Malang, September Penulis

KATA PENGANTAR Dengan mengucap syukur alhamdulillah kehadirat Allah SWT yang telah

melimpahkan

rahmatNya,

sehingga

dapat

terselesaikan

pembuatan diktat kuliah Metode Elemen Hingga ini. Diktat ini disusun dimaksudkan untuk membantu serta menunjang matakuliah Metoda Elemen Hingga sebagai pegangan dasar. Buku ini disusun berdasarkan beberapa buku acuan serta pengalaman penulis selama mengajar

matakuliah ini. Dalam kesempatan ini penulis

mengucapkan pada semua fihak yang telah membantu hingga tersusunnya diktat kuliah ini. Akhirnya penulis menyadari bahwa diktat ini masih banyak kekurangan, untuk itu adanya kritik dan saran yang membangun sangat diharapkan agar karya-karya selanjutnya lebih sempurna lagi.

Malang, September 2003

Penulis

DAFTAR ISI PENDAHULUAN DAFTAR ISI

I

II

BAB I : DASAR-DASAR METODE ELEMEN HINGGA 1.1 Pendahuluan

1

1.2 Sistem Koordinat

2

1.1 Sistem koordinat 2-D/Sistem Koordinat Luasan 1.2 Sistem Koordinat 3-D (Elemen Tetrahedral) 1.3 Transformasi Koordinat

1.5 Konsep Dasar Analisis MEH

6

7

1.6 Metoda Untuk Formulasi Integral

8

1.7 Analisis Prinsip Energi Potensial Minimum 1.8 Konsep Elemen Hingga 2-Dimensi 1.9 Elemen Segitiga Isoparametrik

31

2.2 B E A M

41

2.3 F R A M E

10

18

26

29

BAB II : APLIKASI PADA STRUKTUR 2.1 T R U S S

4

4

1.4 Hubungan Tegangan-Regangan

1.10 Elemen Segiempat

3

31

47

BAB III : INTERPOLASI DAN INTEGRASI NUMERIK

51

BAB IV : APLIKASI PADA PERPINDAHAN PANAS

54

4.1 Steady State Uniaxial Heat Flow

54

4.2 Model Elemen Hingga Aliran Panas 1-Dimensi

56

4.3 One Dimensional Heat Flow With Convection

58

4.4 Perpindahan Panas dan Aliran Fluida 2-Dimensi

62

1

BAB V : ANALISA TEGANGAN AXISYMMETRIC 5.1 Persamaan Dasar untuk Elemen

66

5.2 Persamaan Elastisitas Axisymmetric

DAFTAR PUSTAKA

71

67

64

DIKTAT METODE ELEMEN HINGGA Oleh : Ir. A. As’ad Sonief, MT.

BAB I

DASAR-DASAR METODE ELEMEN HINGGA Pada bab ini dibahas mengenai dasar-dasar analisa elemen hingga, yang didalamnya meliputi sistem koordinat, transformasi koordinat, hubungan tegangan-regangan, prinsip energi potensial minimum, dan juga konsep model untuk elemen 2 dimensi. 1.1 Pendahuluan Perkembangan

dunia

komputer

telah

begitu

cepatnya

mempengaruhi bidang-bidang penelitian dan industri, sehingga impian para ahli dalam mengembangkan ilmu pengetahuan dan industri telah menjadi kenyataan. Pada trend sekarang ini, metoda dan analisa desain telah banyak menggunakan perhitungan metematis yang rumit dalam penggunaan sehari-hari. Metode elemen hingga (finite element method) banyak memberikan andil dalam melahirkan penemuan-penemuan bidang riset dan industri, hal ini dikarenakan dapat berperan sebagai research tool pada eksperimen numerik. Aplikasi banyak dilakukan pada problem kompleks diselesaikan dengan metode elemen hingga seperti rekayasa struktur, steady state dan time dependent heat transfer, fluid flow,

dan

electrical

potential

problem,

aplikasi

bidang

medikal.

Gambaran dasar sebagai berikut.

Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

1


Proses Analisa M E H

Change of physical problem

Physical problem

Mathematic model Governed by differential equations Assumptions on • Geometry • Kinematics • Material law • Loading • Boundary conditions, etc.

Finite element solution of mathematical model

Finite element solution Choice of • Finite elements • Mesh density • Solution parameters Representation of • Loading • Boundary conditions, etc.

Improve mathematical model

Refine mesh, solution parameter etc.

Assessment of accuracy of finite element solution of mathematical model Interpretation result

Refine analysis

Design improvements Structural optimization

1.2 SISTEM KOORDINAT -

Sistem koordinat global → koordinat struktur untuk sebuah titik pada continum

-

-

Ref untuk seluruh continum

-

Ref untuk seluruh struktur

Sistem koordinat lokal → Sistem koordinat yang dipasang pada elemen (acuan pada elemen yang bersangkutan)


2


-

-

dipasang elemen

-

Ref untuk titik-titik yang ada di elemen

Sistem koordinat natural → Terdiri atas koordinat tanpa dimensi untuk identifikasi posisi, dengan tanpa terpengaruh oleh keluaran elemen.

→ Merupakan nisbah koordinat tersebut terhadap ukuran elemen Sistem koordinat Natural 1-D (elemen garis)

Koordinat global P(xp) Koordinat lokal P (xs) Koordinat natural P(L1,L2)

L1 = 1 − 1.2.1

S S ; L2 = L L

Sistem Koordinat 2-D / Sistem Koordinat Luasan

(elemen segitiga) P (L1, L2, L3)

Dimana

L1 + L2 + L3 = 1 L1 =

Luas ∆P − 2 − 3 S1 = Luas ∆1 − 2 − 3 t1

Luas ∆P − 1 − 3 S 2 = Luas ∆1 − 2 − 3 t 2 Luas ∆P − 1 − 3 S 3 L3 = = Luas ∆1 − 2 − 3 t 3

L2 =


3


1.2.2

Sistem koordinat 3-D (elemen tetrahedral)

P (L1, L2, L3, L4)

Dimana

L1 =

Vol ∆ P − 2 − 3 − 4 Vol ∆ 1 − 2 − 3 − 4

L2 =

Vol ∆ P − 1 − 3 − 4 Vol ∆ 1 − 2 − 3 − 4

L3 =

Vol ∆ P − 1 − 2 − 4 Vol ∆ 1 − 2 − 3 − 4

L1 =

Vol ∆ P − 1 − 2 − 3 Vol ∆ 1 − 2 − 3 − 4

L1 + L2 + L3 + L4 = 1 1.3

TRANSFORMASI KOORDINAT Koordinat yang banyak digunakan dalam metode elemen hingga

adalah koordinat kartesian, dan koordinat

sering dinyatakan dalam

bentuk vektor yang dijabarkan sebagai berikut :


4


X  p =  p  Yp 

X p  P=   Yp 

X p = X p Cosθ + Y p Sinθ Y p = − X p Sinθ + Y p Cosθ  X  − Cosθ  = Y   − Sinθ

Sinθ   X  .  Cosθ   Y  − Cosθ

Matrik transformasi [T] =   − Sinθ

 X  Cosθ  =  Y   Sinθ

Sinθ  Cosθ

− Sinθ   X  .  Cosθ  Y 

[T]-1 = [T]T → orthogonality Koordinat dinyatakan dalam 3 Dimensi Orientasi X (l1, m1, n1) Orientasi Y (l2, m2, n2) Orientasi Y (l3, m3, n3)

 X   l1    Y  = l 2  Z  l   3

m1 m2 m3

n1   X    n2   Y  n3   Z 

[T]


5


1.4

HUBUNGAN TEGANGAN – REGANGAN

εx =

σ x − v.σ y − v.σ z E

εy =

σ y − v.σ x − v.σ z E

εz =

σ z − v.σ x − v.σ y E

γ xy =

τ xy τ yz τ ; γ yz = ; γ zx = zx G G G

dimana : G =

E 2(1 + v)

E = Modulus Elastisitas ν = poisson ratio

{ε } = [C ].{σ }  1 − v  1 − v c = . E 0   0  0 

−v 1 −v 0 0 0

{ε }T = {ε x −v 0 1 0 0 0

εy

εz

ε xy

ε yz

ε zx }

0 0 0  0 0 0   0 0 0   2.(1 + v ) 0 0  0 2.(1 + v ) 0  0 0 2.(1 + v )

Selanjutnya :

{σ } = [E ].{ε } Dimana ;

a b  E b [E ] =  1 + V 0 0  0

b a b 0 0 0 a=

b b a 0 0 0

0 0 0 c 0 0

0 0 0 0 c 0

0 0 0  0 0  c 

1− V V 1 ; b= ; c= 1 − 2V 2 1 − 2V

[ E ] = [C ] −1 Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

6


1.5

KONSEP DASAR ANALISIS MEH.

Dua kategori model matematik : -

lumped-parameter models (“discrete-system”)

-

continuum-mechanics-based models (“continuous-ystem”).

Kondisi Problem : 1. Steady -State Problems. K.U =R 2. Propagation Problems/Dynamic Problem. M . Ü + K . U = R(t) 3. Eigenvalue Problems. Konsep Dasar Metode Elemen Hingga 1.

Menjadikan elemen-elemen diskrit untuk memperoleh simpangansimpangan dan gaya-gaya anggota dari suatu struktur.

2.

Menggunakan elemen-elemen kontinum untuk memperoleh solusi pendekatan

terhadap

permasalahan-permasalahan

perpindahan panas, mekanika fluida dan mekanika solid. Dua karakteristik yang membedakan metoda elemen hingga dengan metoda numeric yang lain yaitu : -. Metoda ini menggunakan formulasi integral untuk menghasilkan sistem persamaan aljabar. -. Metoda ini menggunakan fungs-fungsi kontinyu untuk pendekatan parameter-parameter yang belum diketahui. Lima langkah untuk menyelesaikan permasalahan fisik dengan metoda elemen hingga yaitu : 1. Permasalahan fisik dibuat elemen-elemen kecil. Elemen-elemen tersebut ditandai dengan nomor elemen dan nomor titik nodal, termasuk juga harga-harga koordinat.


7


2. Tentukan persamaan pendekatannya, linear atau kuadratik. Persamaan-permsamaan tersebut harus ditulis dalam bentuk harga-harga nodal yang belum diketahui. Ini berlaku untuk setiap elemen, artinya setiap elemen harus didefinisikan sifatnya dalam bentuk persamaan diatas. 3. Bentuklah sistem persamaan diatas dengan metoda Galerkin, Varisional, Formulasi energi potensial, Collocation, Subdomain, dll. Khusus untuk formulasi energi potensial, energi potensial dari sistem ditulis

dalam

bentuk

simpangan

nodal

dan

kemudian

diminimalkan. Dimana akan diberikan satu persamaan setiap simpangan yang belum diketahui. 4. Selesaikan sistem persamaan diatas. 5. Hitung besaran yang dicari. Besaran bisa berupa komponenkomponen tegangan, aliran panas atau kecepatan fluida. 1.6 METODA UNTUK FORMULASI INTEGRAL Metoda Varisional

Π=

H

∫0

  D  dy  2 − Qy   dx  dx 2    

(1)

Harga numeric Π dapat dikalkulasi dengan memberikan persamaan coba-coba y=f(x). Misal persamaan coba-coba yang memberikan harga terkecil Π adalah y=g(x), maka persamaan ini merupakan jawab dari persamaan diferensial berikut :

D

d 2y dx 2

+Q = 0

dengan kondisi batas y(0)=y0 dan y(H)=yH

(2) harga Π minimum adalah

merupakan jawab pendekatan.


8


Weighted Residual Method; Ritz Method

Andaikan bahwa y=h(x) adalah merupakan jawab pendekatan terhadap persamaan (2), dengan subsitusi akan memberikan :

D

d 2 h( x ) dx 2

+ Q = R( x ) ≠ 0

karena y=h(x) tidak memenuhi persyaratan persamaan, WRM mengharuskan : H

∫0 W i ( x )R( x )dx = 0

fungsi residual R(x) ;fungsi pemberat (weighting) Wi(x), Beberapa pilihan fungsi pemberat dengan beberapa metoda yang popular : 1. Metoda Collocation 2. Metoda Subdomain 3. Metoda Galerkin 4. Metoda Least Squares

Formulasi Energi Potensial Integral volume dengan hasil kali komponen tegangan & regangan.

Λ=

∫v

σ xx .ε xx dV . 2

Prinsip energi potensial minimum dan energi regangan banyak digunakan untuk menganalisis masalah-masalah struktur dan mekanika solid.


9


1.7 ANALISIS PRINSIP ENERGI POTENSIAL MINIMUM

Variabel tak bebas → dof Variabel bebas → koordinat Ada syarat kontinuitas → bentuk persamaan tidak ada gabungan Kompatibilitas → berkaitan dengan dof Elemen linear → node diujung, sebagai contoh seperti pada elemen ∆ linear sederhana Dalam domain mekanika solid → harus ada boundary condition (BC’s) yaitu dof yang direstrin/ diberikan kendala. Domain yang terbagi sumbu domain merupakan : -

Kasus per elemen dengan f interpolasi

-

Keseimbangan statis pada elemen dengan kaidah struktur yang dikenai beban akan terdeferensi (prinsip energi potensial minimum)

Keseimbangan terjadi kalau energi potensial minimum dalam suatu sistem. Dalam MEH merupakan suatu teknik numerik dari model matematis suatu sistem yang digambarkan dari suatu fenomena problem. Sebagai gambaran dapat diterapkan pada elemen garis, dan dengan konsep energi potensial minimum (pada solid mekanik) kemudian dilakukan dengan teknik numerik murni sehingga membentuk persamaan diskrit sebagai berikut:

[

]

{φ} = {f}, yaitu suatu matrik dikalikan dengan

vektor dof sama dengan vektor beban. Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

10


Energi potential total = Kerja gaya luar + Energi regangan -

Beban terpusat

-

Beban traksi (bekerja pada permukaan)

-

Body force (centrifugal, gaya magnit gravitasi, gaya elektromaknetik) (Beban/Variabel)

Prinsip Energi Potensial Minimum Analisa tegangan (prob elastisitas benda padat) dengan FEM didasarkan

pada

prinsip

Energi

potensial

minimum

yang

menyatakan : Dari

sekian

kompatibilitas persamaan

persamaan interval

dan

perpindahan

perpindahan memenuhi yang

juga

yang

syarat

memenuhi

batas,

memenuhi

maka kondisi

keseimbangan stabil adalah persamaan perpindahan yang memberikan / menghasilkan energi potensial yang terkecil (minimum). Prinsip tersebut mengimplikasikan hal-hal sebagai berikut : -

Perlunya pendefinisian persamaan perpindahan untuk setiap elemen yang memenuhi syarat kompabilitas antar elemen.

-

Persamaan perpindahan tersebut diatas harus memenuhi semua syarat batas

-

Penjabaran persamaan energi potensial yang dianalisa. Persamaan diumpamakan sebagai fungsi persamaan (dalam hal ini persamaan node) yang akan dicari nilainya (yang tidak diketahui)

-

Minimalisasi energi potensial terhadap persamaan yang tidak diketahui tersebut.

Energi Potensial Energi regangan – kerja yang dilakukan oleh gaya-gaya eksternal yang bekerja pada sistem. Energi Potensial Energi regangan – kerja yang dilakukan oleh gaya-gaya eksternal yang bekerja pada sistem. Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

11


Energi regangan

U (e) = 1

2∫

(σ x .ε x + σ y .ε y + σ z .ε z + τ xy γ xy + τ yz .γ yz + τ zx .γ zx )dv

V

=1

2∫

{σ}T {ε}dv = 1

V

2∫

{d}T [B]T [E][B]dv

V

Kerja yang dilakukan body force

Wbf = ∫ (u.bx + v.b y + w.bz )dV V

Kerja yang dilakukan oleh beban traksi (beban terdistribusi)

Wt = ∫ (u. p x + v. p y + w. p z )dA V

Kerja yang dilakukan oleh beban terpusat

W f = d x .Px + d y .Py + d z .Pz Energi potensial total : n

π = ∑ π e − {d } .{P} T

e =1

Dimana : π e = u e − Wbf − Wt Minimalisasi energi potensial,

∂x = 0 , maka ∂d

∑ [K ]e .{d } = ∑ {f e }+ {P} n

n

e =1

e −1

Merupakan rumus umum. Dimana :

{f }= {f }+ {f } e

e bf

e

t

Contoh penyelesaian MEH dari persamaan diferensial : Persamaan deferensial :

d 2u +u =1 dx 2

Kondisi batas :

u(0)= 0


; u(2π)=0

12


Solusi eksak : u = 1 – cos x. Prosedur Penyelesaian : 1. Diskrititasi region. Dalam region dibagi dalam 4 elemen dan elemen dan nodal diberi nomor.

u~

1

1

2

2 π/2

0

3

3

π

4 4 3π/2

5 2π

2. Buat trial function.

u~ e

i

j x

xi

L

Fungsi asumsi :

u~ = a1 + a2 x u~ = u i = a1 + a2 x i

a1 =

u~ = u j = a1 + a 2 x j

a2 =


x j ui − xi u j x j − xi − ui + u j x j − xi

13


 xj − x u~ =   x j − x i

x − xi x j − xi

 u i     = [N1 N 2 ]{q}  u j 

3. Substitusi trial functions kedalam governing equation.

W

du dx

X2 X1

4

4 4 dW du~ dx + ∑ ∫ W .u~.dx −∑ ∫ W .dx = 0 X e dx dx X X e =1 e e =1 e

− ∑∫ e =1

Weighting function untuk metode Galerkin :

Wi =

∂u~ ∂ai

untuk masing-masing konstanta a1 dan a2 :

W1 =

∂u~ = N1 ∂ai

W1 = N1 = dan :

xj − x x j − xi

dW1 dN1 −1 = = dx dx x j − xi

W2 =

∂u~ = N2 ∂a j

W2 = N 2 =

x − xi x j − xi

dW2 dN 2 1 = = dx dx x j − xi

governing equation dalam bentuk matrik :

 dN1   dN dN 2  u i   − ∫  dx  1  u dx xi dN dx dx 2   j  xi    dx  u i  x j N  x j N  + ∫  1 [N1 N 2 ] dx − ∫  1 dx = 0 xi N xi N  2  2 u j 

 N1  du  .. N 2  dx

xj

xj


14


Pengembangan suku 1 :

  du   du    N1  0   du  − dx   dx xi   du     xj   xi  =  −  dx xi  =  du −   N 2 du   dx x j   0      dx  dx xi   xj  x ji    

 du  N ` dx   N 2 du  dx  Suku 2 :

 dN1 xj  dx ∫xi  dN  2  dx

dN1 dx dN1 dx

dN1 dx dN 2 dx

dN 2  dx  u i dx = i  1 − 1 u i    dN 2  u j  Le − 1 1  u j   dx 

dimana : Le = xj - xi Suku 3 : xj

∫x

i

u i  N1N 2  u i  x j  N1N1  N1  [ ] = dx N N    1 2   u dx ∫ x i N N u N N N j  2 1  2 2 2  j    x 3j − x i3 − 3 x i x j Le 2 1 u i  =  1 2 u  6.L2e  j  

Suku 4 : xj

∫x

i

x 2j + x i2 − 2.x j .x i  N1   dx = 2.Le N 2 

1  1

Secara keseluruhan :

 − du   dx ( x i ) 1  du −  ( x j )  Le  dx 

3 3  1 − 1 u i  x j − x i − 3.x i .x j .Le − 1 1  u  + 6.L2e   j 

2 1 u i  1 2 u    j 

x 2j + x i2 − 2 x j .x i 1 0 −  =  2.Le 1 0 Aplikasi untuk setiap elemen, dengan asumsi Le = L Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

15


Elemen 1 : i = 1, j = 2, x1 = 0 , dan x2 = L

 du  − dx  1  1 − 1 u  L 2 1 u  L 1 0 1 1 x =0   −   =    +   −     du  L − 1 1  u 2  6 1 2 u 2  2 1 0  dx x =L  Elemen 1 : i = 2, j = 3, x2 = L , dan x3 = 2L

 du  − dx   1 − 1 u 2  L 2 1 u 2  L 1 0 x =L  − 1   − 1 1  u  + 6 1 2 u  − 2 1 = 0 du L   3    3        dx x =2L 

dst.

Diasumsikan du/dxIx=L pada elemen 1 sama dengan du/dxI x=L pada elemen 2 maka : Asembly persamaan :

 du  − dx   1 − 1 0 0 0   u1  2 x =0        0  1 − 1 2 − 1 0 0  u 2  L 1      0  −  0 − 1 2 − 1 0  u 3  + 0   L  0 0 − 1 2 − 1 u  6 0 0   4     du  0 0 0 − 1 1  u 5  0    dx x = 4L 

1 4 1 0 0

0 1 4 1 0

0 0 1 4 1

0   u1  1 0 2 0 0 u 2    L     0 u 3  − 2 = 0 2     0 2 1 u 4        1 0 2 u 5 

dengan kondisi batas essential : u1 = 0 ; u5 = 0 maka :

 2 − 1 0  u 2  4 ` 0  u 2  2 0 − 1 L L      u +  1 4 1 u − 2 = 0 − − 1 2 1  3  6   3  2     L    2 0  0 − 1 2  u 4  0 1 4 u 4      disederhanakan dan dengan L = π/2 didapat :


16


0 − 2.1304 8.4674  u 2  1     − 2.1304 8.4674  u 3  = 14.804 1  1  − 2.1304  u 4   u 2  1.130      u 3  = 2.033 u  1.130   4   X π/4 2π/4 3π/4 4π/4 5π/4 6π/4 7π/4

Exact 0.293 1.000 1.707 2.000 1.707 1.000 0.293

4 Elemen 1.130 2.033 1.130 -

8 Elemen 0.332 1.038 1.722 2.003 1.722 1.038 0.332

Gambar hasil yang yang dibandingkan dengan solusi eksak dan MEH dengan beda jumlah elemen sebagai berikut :

u 2

X(rad) 0


π

17


1.8 KONSEP MODEL ELEMEN HINGGA 2 – DIMENSI ELEMEN LUASAN (SEGITIGA , SEGIEMPAT).

Y

1

X •

Sistem koordinat.

! Global Coordinate Fungsi asumsi : U(X,Y) = α1 + α2 X + α3 Y V(X,Y)= β1 + β2X + β3 Y

Y

3 V1

2 U1

1 X

u1  1 X 1 Y1  α 1       u 2  = 1 X 2 Y2  α 2  u  1 X Y  α  3 3  3   3  {q1}= [A1] . {α} Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

18


v 1  1 X 1 Y1   β 1       v 2  = 1 X 2 Y2   β 2  v  1 X Y   β  3 3  3   3  {q2}= [A1] . {β}

[A1 ]

−1

1 X 1 Y1  = 1 X 2 Y2  1 X 3 Y3 

−1

a adjo int of [ A1 ] 1  1 = = b1 det er min ant of [ A1 ] det   c1

a2 b2 c2

a3  b3  c 3 

{α} = [A1] -1 . {q1} {β} = [A1] -1. {q2} {u} = [1

X

Y}.[A1] -1 . {q1}

{v} = [1

X

Y}. [A1] -1. {q2}

u1  {q1} = u 2  u   3 Ekspansi : [1

=

X

v 1  {q 2 } = v 2  v   3 [1

X

Y}.[A1] -1 .

Y}.[A1] -1 .

1 [[a1 + Xb1 + Yc1 ] det

= [N1

N2

[a2 + Xb2 + Yc 2 ] [a3 + Xb3 + Yc 3 ]]

N3]

sehingga :

u = [N1

N2

N3] .

 u1    u 2  u   3


19


v = [N1

N2

N3] .

v 1    v 2  v   3

dalam bentuk matrik

u  N1 0 N 2  = v   0 N1 0

O N2

 u1  v   1 0  u 2    N 3  v 2  u 3    v 3 

N3 0

atau bentuk symbol : {u} = [N] . {q} Koordinat local : u(X,Y) = α1 + α2 x + α3 y v(X,Y)= β 1 + β2x + β3 y

Y

3 y 2

x

1 X

u1  1 0    u 2  = 1 x 2 u  1 x 3  3 

0  α 1    0  α 2  y 3  α 3 

{q1}= [A1] . {α}


20


[A1 ]−1

a adjo int of [ A1 ] 1  1 = = b1 det er min ant of [ A1 ] det   c1

a2 b2 c2

a3  b3  c 3 

{α} = [A1] -1 . {q1} {β} = [A1] -1. {q2} {u} = [1

x

y}.[A1] -1 . {q1}

{v} = [1

x

y}. [A1] -1. {q2}

u  N1 0 N 2  = v   0 N1 0

O N2

 u1  v   1 0  u 2    N 3  v 2  u 3    v 3 

N3 0

atau bentuk symbol : {u} = [N] . {q} dimana :

y 3 (x 2 − x ) + y ( x3 − x2 ) ; x 2 .y 3 yx 2 N3 = x 2 .y 3 N1 =

N2 =

x.y 3 − yx 3 x 2 .y 3

Koordinat Natural

L2

Y

3 L3

2

1

x

L1 X


21


Fungsi asumsi : u = L 1 u1 + L 2 u2 + L 3 u3 Hubungan koordinat natural :

L1 + L 2 + L 3 = 1

u = L1 u1 + L2 u2 + (1 – L1 – L2) u3 v = L1 v1 + L2 v2 + (1 – L1 – L2) v3 Untuk elemen isoparametrik : X = L1 X1 + L2 X2 + (1 – L1 – L2) X3 Y = L1 Y1 + L2 Y2 + (1 – L1 – L2) Y3

Aplikasi solid (mekanik) : -plane stress - plane strain " Elemen segitiga linear (elemen regangan konstan)

Ciri : - 3 node per elemen - 2 dof per node u : displacement arah x v : displacement arah y Q variasinya diasumsikan fungsi linear (pada sub domain bervariasi linear) Pada solid mekanik, konsekuensi linear → regangan konstan di titik manapun di elemen sehingga tegangan juga konstan.


22


Step 1 * membuat fungsi linear Fungsi interpolasi (asumsi) displacement

u ( x , y ) = α 1 + α 2 .x + α 3 . y v ( x, y ) = β 1 + β 2 . x + β 3 . y

ε x ( x, y ) =

∂u = α2 ∂x

ε y ( x, y ) =

∂v = β3 ∂y

γ xy ( x, y ) =

∂u ∂v + = α3 + β2 ∂y ∂x

u dan v → titik sebarang pada elemen (boleh node/tidak) Shape function ; Step 2 Menyatakan hubungan ∈ dengan displacement node {∈} = [B] {d} Step 3

[ K ( e ) ] = ∫ [ B ]T .[ E ].[ B ]dv V

Untuk tebal elemen konstan = h

[ K ( e ) ] = ∫ [ B]T .[ E ].[ B].h.dA A

[ K ( e ) ] = [ B]T .[ E ].[ B].h. A → Untuk : plane stress Plane strain → untuk h = 1 unit yang membedakan [E] " Beban node ekuivalen akibat body force

{f }= ∫ [N ] b  T

bf

V

bx  dV y


23


-

body force → jadi 2 komponen dalam fungsi x dan y

-

batas integral untuk elemen

" Beban node ekuivalen akibat traksi

{f }= ∫ [N ]  p  T

bf

A

px  dA y

" Beban node ekuivalen akibat beban thermal (beban mula)

{∈th }T = [α .∆T

α .∆T

0]

[ f th ] = ∫ [B]T .[E ]{∈th }dV V

Untuk setiap elemen perlu dianalisa

[K ] {f }= {f } (e )

(e)

e bf

Untuk struktur n

[K ] = ∑ [ K ( e) ] e =1 n

{ }+ {P} ← Beban terpusat

{F } = ∑ f e =1

(e)

[K ]{D} = {F } Solusi kasus 2-D

Fungsi interpolasi

u ( x , y ) = α 1 + α 2 .x + α 3 . y


24


U 1 = α 1 + α 2 .x1 + α 3 . y1 U 2 = α 1 + α 2 .x 2 + α 3 . y 2 U 3 = α 1 + α 2 .x 3 + α 3 . y 3

U 1  1 X 1    U 2  = 1 X 2 U  1 X 3  3 

Y1  α 1    Y2 .α 2  Y3  α 3 

α 1  1 X 1    α 2  = 1 X 2 α  1 X 3  3 

Y1  Y2  Y3 

α 1  a1   1 α 2  =  b1 α  J  c  3  1

a2 b2 c2

−1

U 1    .U 2  U   3

a 3  U 1    b3 .U 2  c3  U 3 

a1 = ( x 2 . y 3 − y 2 .x3 ) ; a 2 = ( y1 .x3 − x1 . y 3 ) ; a 3 = ( x1 . y 2 − y1 .x 2 ) b1 = y 2 − y 3 ; b2 = y 3 − y1 ; b3 = y1 − y 2 c1 = x3 − x 2 ; c 2 = x1 − x3 ; c3 = x 2 − x1 J = ( x 2 . y 3 − y 2 .x3 ) + x1 ( y 2 − y 3 ) + y1 ( x3 − x 2 )

{U } = [1 {U } = [1

{U } = [1

x

y ].{α }

x

α 1    y ].α 2  α   3

x

a1 y ]. b1  c1

a2 b2 c2

{U } = 1 [(a1 + b1 x + c1 y) J

a3  U 1    b3 .U 2  c3  U 3  ( a 2 + b2 x + c 2 y )


U 1    (a 3 + b3 x + c3 y )].U 2  U   3

25


{U } = [N1

N2

U 1    N 3 ].U 2  U   3

N1 =

1 (a1 + b1 x + c1 y ) J

N2 =

1 (a 2 + b2 x + c 2 y ) J

N3 =

1 (a 3 + b3 x + c3 y ) J

1.9 ELEMEN SEGITIGA ISOPARAMETRIK Elemen isoparametrik yaitu

fungsi interpolasi untuk koordinat

geometri-identik dengan fungsi interpolasi untuk perpindahan. Pada Elemen segitiga digambarkan sebagai berikut

 X1  Y   1 X   X 2    = [N ].  Y   Y2  X3     Y3  Misal


26


Sehingga yang dibicarakan adalah koodinat natural, tidak hanya :

U = L1 .U 1 + L2 .U 2 + L3 .U 3 V = L1 .V1 + L2 .V2 + L3 .V3 Tetapi

X = L1 . X 1 + L2 . X 2 + L3 . X 3

X1, Y1 → koordinat node

Y = L1 .Y1 + L2 .Y2 + L3 .Y3 L1, L2, L3 = koordinat natural (luasan) L1, L2, L3 = 1 Interpolasi Formula

X ( s, t ) = N 1 . X 1 + N 2 . X 2 + N 3 . X 3 + N 4 . X 4 Y ( s, t ) = N 1 .Y1 + N 2 .Y2 + N 3 .Y3 + N 4 .Y4 N i ( s, t ) Dengan formula interpolasi lagrange Dalam arah x

dalam arah y

Untuk n = 2

l1 ( x) = 1

x − x2 x1 − x 2

l1 ( y ) = 1

y − y4 y1 − y 4

Elemen shape function N1e

N 1 ( x, y ) = l11 ( x).l1 ( y ) = 1

e

x − x2 y − y 4 . x1 − x 2 y1 − y 4

Untuk :

N 1 ( s, t ) =

s − s2 t − t 4 . s1 − s 2 t1 − t 4


27


Node 1 : s1 = -1 ; t1 = -1

Node 3 : s3 = 1

; t3 = 1

Node 2 : s2 = 1 ; t2 = -1

Node 4 : s4 = -1 ; t4 = 1

N 1 ( s, t ) =

( s − 1) (t − 1) (1 − s ) (1 − t ) = . . −1−1 −1−1 2 2

N 1 ( s, t ) =

(1 − s ).(1 − t ) 4

N 2 ( s, t ) =

( s + 1) (t − 1) (1 + s ) (1 − t ) s − s1 t − t3 . = . = . s2 − s1 t 2 − t3 1+1 −1−1 2 2

N 2 ( s, t ) =

(1 + s ).(1 − t ) 4

N 3 ( s, t ) =

s − s4 t − t 2 ( s + 1) (t + 1) (1 + s ).(1 + t ) . = . = s3 − s 4 t 3 − t 2 1+1 1+1 4

N 4 ( s, t ) =

s − s 3 t − t1 ( s − 1) (t + 1) (1 − s ).(1 + t ) = . = . 4 s 4 − s 3 t 4 − t1 − 1 − 1 1 + 1

Kelemahan elemen linear -

Berawal dari asumsi U = a1 + a 2 x + a3 . y

→ regangan konstan maka kalau membahas defleksi tegangan baik hanya ditengah perbaikan dengan membentuk elemen nonlinear untuk

U = N 1 .U 1 + N 2 .U 2 + N 3 .U 3

N1=Li

i = 1, 2, 3

V = N1.V1 + N 2 .V2 + N 3.V34 dengan asumsi :

U = a1 + a 2 x + a3 . y V = b1 + b2 x + b3 . y


28


1.10

ELEMEN SEGI EMPAT

Keuntungan : pada FEM yang didapat → distribusi Pada konvensional → yang didapat pada titik tertentu " Elemen Isoparametrik n

U = ∑ N i .U i i =1 n

V = ∑ N i .Vi i =1

n

X = ∑ N i .X i 1

i =1

n

Y = ∑ N i .Yi 1

i =1

Ni = Ni1 → isoparametrik " Elemen Isoparametrik

Linear → hanya mempunyai node diujung-ujungnya Penomoran : sebarang, tapi analisanya dimulai dengan CCW Dimapping ke koordinat s. t → ke koordinat natural Isoparametrik

U ( s, t ) = N 1 .U 1 + N 2 .U 2 + N 3 .U 3 + N 4 .U 4 V ( s, t ) = N 1 .V1 + N 2 .V2 + N 3 .V3 + N 4 .V4


29


X ( s, t ) = N 1 . X 1 + N 2 . X 2 + N 3 . X 3 + N 4 . X 4 Y ( s, t ) = N 1 .Y1 + N 2 .Y2 + N 3 .Y3 + N 4 .Y4 N1 =

(1 − s )(1 − t ) 4

N3 =

(1 + s )(1 + t ) 4

N2 =

(1 + s )(1 − t ) 4

N4 =

(1 − s )(1 + t ) 4

Asumsi fungsi Interpolasi untuk perpindahan

U = L1 .U 1 + L2 .U 2 + L3 .U 3 V = L1 .V1 + L2 .V2 + L3 .V3

U   L1  = V   0

0 L1

L2 0

0 L2


L3 0

U 1  V   1 0  U 2  .  L3  V2  U 3     V3 

30


BAB II

APLIKASI PADA STRUKTUR Aplikasi elemen hingga untuk analisa struktur, yaitu untuk struktur truss, beam, dan frame. Juga dijelaskan mengenai ciri-ciri masingmasing stuktur tersebut, kelebihan dan kekurangannya masingmasing 2.1 T R U S S Adalah struktur yang istimewa, dimana joint yang dirancang tidak untuk mendukung momen, dan dapat dikatakan merupakan elemen 2 – Force member yang seolah-olah merupakan sambungan pin.

Konsekuensi Karena tidak mendukung momen dalam keseimbangannya → batang sebagai 2- force member sehingga beban selalu dikerjakan di joint. Sehingga gaya-gaya berimpit dengan sumbu aksial batang.

Dalam MEH → diskritisasi dengan setiap batang sebagai elemen dengan membuat node-node, dengan berat sendiri diabaikan. Struktur yang dilas bisa didekati dengan truss asal fabrikasinya baik yaitu sumbu aksial bertemu di satu titik. Elemen garis dapat berupa truss, beam, frame. Metoda langsung → Hubungan displacement dan kekakuan

P


31


∆=

PL P = AE AE

= L

P K

Derajat kebebasan (dof) → displacement (dalam struktur) → variable analisa Per node → memiliki 1 dof Elemen truss yang terletak pada sumbu x

Hubungan → gaya, displacement, stifness Bagaimana dengan display yang ditengah → Fungsi interpolasi (pendekatan) untuk displacement : dipilih polynomial (karena mudah didefferensialkan / diintegrasikan) Syarat : - Kontinuitas Kompabilitas

-

U ( x) = a1 + a 2 .x (asumsi) E ( x) =

du ( x) = a 2 (konstanta) dx

T ( x) = Eε ( x) = E.a 2 (konstanta) pada x = 0 U1 = a1

a1 = ui

pada x = L

a2 =

U2 = a1 + a2 L

U ( x) = U 1 +

u 2 − u1 L

U 2 − U1 x X  X = 1 − U 1 + U 2 L L  L

U ( x) = f 1 ( x)U 1 + f 2 ( x)U 2 E ( x) = f 1 ( x)U 1 + f 2 ( x)U 2 1

1


32


N 1 ( x) = f 1 ( x ) = 1 − N 2 ( x) = f 2 ( x ) =

x L

x L

Shape Function

(Sebagai pola umum perpindahan sebagai fungsi dari Shape function dengan dof)

 L 1 1   L 1 1    X 1 =  EA∫ f1 . f 2 .dx u1 +  EA∫ f1 . f 2 .dx u 2  0   0  ditulis dalam bentuk vektor [k]

{d}

↓

↓

Stiffness matrix

=

{f} ↓

vektor

vektor

disp. node

load node

[k] = matrik kekakuan elemen L

k ij = EA∫ f i ( x). f j ( x)dx 1

1

0

1 − 1 L − 1 1 

[k ] = EA 

1 − 1  − 1 1 

[k ] = k 

Persamaan kekakuan dengan Metode Energi : axial force :

S = T . A( x) = Eε . A( x) = EA

∂u ∂x

= E. A( x)[ f 1 ( x)U 1 + f 2 ( x)U 2 ] 1

1

{d }= [T ]{d } dengan cara sama :

{f }= [T ]{f } Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

33


{}

[ K ] d = [T ]{f } [ K ].[T ].{d } = [T ]{f } [T ]T .[ K ].[T ].{d } = {f } [ K ].{d } = { f } dimana

 Cos 2θ  AE  Sinθ .Cosθ [K ] = L  − Cos 2θ  − Sinθ .Cosθ

Sinθ .Cosθ Sin 2θ − Sinθ .Cosθ − Sin 2θ

− Cos 2θ − Sinθ .Cosθ Cos 2θ Sinθ .Cosθ

− Sin 2θ .Cosθ   − Sin 2θ  Sinθ .Cosθ   Sin 2θ 

model matematis

 u1   x1  v   y   1  1 [K ]   =   u 2   x 2  v 2   y 2  Elemen truss dengan orientasi sembarang

Model matematis (Persamaan keseimbangan node)

AE  1 − 1 u1   X 1   =  L − 1 1  u 2   X 2  [K]

{d} = {f}

Spesifikasi elemen : Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

34


-

2 node pe elemen

-

2 dof per node (u dan v)

Data teknis yang diperlukan : E, A, L, θ 2 node per elemen dengan asumsi perpindahan yang terjadi sepanjang → merupakan variasi linear

X , U , Y , V → Koordinat lokal Dalam sistem sumbu lokal

AE  1 − 1 U 1   X 1  .  =   L − 1 1  U 2   X 2  Dikembangkan dengan 2 persamaan : nol = nol

1  AE  0 L − 1  0

0 −1 0 0 0 1 0 0

0 u1   X 1  0  v1   Y1    =   Atau 0 u 2   X 2   0 v 2   Y2 

[K ].{d } = { f } Dimana

u1   Cosθ  v  − Sinθ  1   = u 2   0 v 2   0

Sinθ Cosθ 0 0

0 0 Cosθ − Sinθ

0   u1  0   v1  .  Sinθ  u 2   Cosθ  v 2 

Resume Truss → digunakan tidak untuk mendukung momen * Steps : 1. Diskritisasi dengan setiap batang sebagai elemen dengan membuat node-node dan diberi nomor.


35


2. Membuat tabel, data yang diketahui dan Cos dan Sin arah setiap elemen 3. Buat model matematis elemen / K elemen 4. Beri notasi pada K elemen sesuai dengan dof 5. Susun nomor notasi dari K elemen pada susunan K total / assembly 6. Identifikasi B . C 7. Temukan dof aktifnya 8. Temukan problem yang ditanyakan (reaksi pada tumpuan, tegangan pada batang, dsb) * Ciri [K] struktur / assemble -

Elemen matriknya : 2 x joint

-

Simetris matrik

-

Singular matrik

-

Tidak semua persamaan independent (hanya 2 persamaan independent)

* Konsep K Struktur / Assemble Gaya node di tiap-tiap node pada struktur merupakan sigma gaya node elemen yang dikontribusikan masing-masing nodenya. * Konsep keseimbangan truss Gaya node pada setiap node sama dengan gaya luar (beban / reaksi tumpuan) dalam arah yang sama. Contoh


36


Tabel i

j

E

A

L

θ

1

2

E

A

L

0o

1

3

E

A

L

60o

2

3

E

A

L

120o

Cos θ

Sin θ

* K elemen / model matematis elemen

 Cos 2θ  AE  Sinθ .Cosθ [K ] = L  − Cos 2θ  − Sinθ .Cosθ

Sinθ .Cosθ Sin 2θ − Sinθ .Cosθ − Sin 2θ

− Cos 2θ − Sinθ .Cosθ Cos 2θ Sinθ .Cosθ

− Sin 2θ .Cosθ   − Sin 2θ  Sinθ .Cosθ   Sin 2θ 

(1–2):

0 −1 0 0 0 1 0 0

1  AE  0 [K ] = L − 1  0 1  AE  0 L − 1  0

0 −1 0 0 0 1 0 0

0 0 0  0

0  u1   x1  0  v1   y1  .  =   0 u 2   x 2   0 v 2   y 2 

(1–3):

 14  AE  3 4 L  − 1/ 4  − 3 4

34 3/ 4 − 34 − 3/ 4

−1 4 − 34 1/ 4 34

− 3 4  u1   x1      − 3 / 4   v1   y1  .  =   3 4  u 3   x3   3 / 4  v3   y 3 

 14  AE − 3 4 L  − 1/ 4   3 4

− 34 3/ 4 34 − 3/ 4

−1 4 34 1/ 4 − 34

3 4  u 2   x 2      − 3 / 4  v 2   y 2  .  =   − 3 4  u 3   x3   3 / 4   v3   y 3 

(2–3):


37


* K struktur

X 1 (1 − 2) =

AE [1.u1 + 0.v1 − 1.u 2 + 0.v 2 ] L

X 1 (1 − 3) =

1  AE  1 1 1 3.v1 − .u 3 − .v3   u1 + L 4 4 4 4  +

1 3  AE  5 3 .v1 − 1.u 2 + 0.v 2 − u 3 − .v3   u1 + L 4 4 4 4 

X1 =

Y1 (1 − 2) =

AE [0.u1 + 0.v1 + 0.u 2 + 0.v 2 ] L

Y1 (1 − 3) =

AE  3 3 3 3  .u 3 − .v3   .u1 + .v1 − L  4 4 4 4  +

Y1 =

AE  3 3 3 3  .u 3 − v3   .u1 + .v1 + 0.u 2 + 0.v 2 − L  4 4 4 4 

X 2 (1 − 2) =

AE [−1.u1 + 0.v1 + 1.u 2 + 0.v 2 ] L

X 2 (2 − 3) =

AE  1 3 1 3  .v 2 − .u 3 + .v3   u2 − L 4 4 4 4  +

X2 =

AE  5 3 1 3  .v 2 − u 3 + .v3  − 1.u1 + 0.v1 + .u 2 − L  4 4 4 4 

Y2 (1 − 2) =

AE [0.u1 + 0.v1 + 0.u 2 + 0.v 2 ] L

Y2 (2 − 3) =

AE  3 3 3 3  .u 2 + .v 2 + .u 3 − .v3  − L  4 4 4 4  +


38


Y2 =

AE  3 3 3 3  ..u 2 + .v 2 + .u 3 − v3  0.u1 + 0.v1 − L  4 4 4 4 

X 3 (1 − 3) =

3  AE  1 3 1 .v1 + .u 3 + .v3  − .u1 − L  4 4 4 4 

X 3 (2 − 3) =

3  AE  1 3 1 .v 2 + .u 3 − .v3  − .u 2 + L  4 4 4 4  +

 1 AE  1 3 1 3 .v1 − .u 2 + .v 2 + u 3 + 0.v3  − .u1 − L  4 4 4 4 2 

X3 =

Y3 (1 − 3) =

AE 3 3 3 3 [− .u1 − .v1 + .u 3 + .v3 ] L 4 4 4 4

Y3 (2 − 3) =

AE  3 3 3 3  .u 3 + .v3   .u 2 − .v 2 − L  4 4 4 4  +

Y2 =

3 3 3 6  AE  3 .u1 − .v1 + ..u 2 − .v 2 + 0.u 3 + v3  − L  4 4 4 4 4 

* Model matematis struktur

 5/ 4   3  4   −1   0   −1/ 4  − 3 4

3

4 3/ 4

−1

0

− 1/ 4

0

0

− 3

0

5/ 4

0

− 3

− 3

4 − 3/ 4

4 − 1/ 4 3

4

− 3

4 3/ 4 3

4 − 3/ 4

[K]


− 3 4  U 1 = 0   X 1 = R x1  − 3 / 4     Y =0    V1 ?   1  3   X = P 4 . U 2 ?  =  2      − 3 / 4  V2 = 0   Y2 = R y 2   U 3 = 0  X 3 = R x 3     0  V3 = 0   Y3 = R y 3   6/4  

4 − 1/ 4 3

4 1/ 2 0 .

{D} =

{f}

39


* Identifikasi B.C U1 = V2 = U3 = V3 = 0 (kondisi tumpuan pada joint) * Dof aktif

AE 3 / 4 0   V1   0  .  =   L  0 5 / 4 U 2   P  V1 = 0 ; U 2 =

 V1  4 P.L 0 4 PL →  =   5 EA U 2  5 EA 1

* Gaya reaksi

 3 4   R x1   0 R   y 2  EA   =   R x 3  L − 3 4  R y 3    3 / 4

−1    3 − 4 . V1    − 1 / 4  U 2  3  4 

 3 4   R x1   R   y 2  EA  0  =   R x 3  L − 3 4  R y 3   3/ 4 

−1    −1     3 − 4 . 4 P.L 0 = 4 P − 3 4       − 1 / 4  5 EA 1 5  − 1 / 4   3 4  3  4 

* Gaya Aksial

S = X 2 .Cosθ + Y2 .Sinθ S = EA[Cosθ X2 = = Y2 =

U 2 − U 1  Sinθ ].   V2 − V1 

(

EA − C 2U 1 − SCV1 + C 2U 2 + SCV2 L

)

EA [(C 2 (U 2 − U 1 ) + SC (V2 − V1 )] L

(

EA − SCU 1 − S 2V1 + S (U 2 + S 2V2 L


) 40


= S=

S=

EA [ SC (U 2 − U 1 ) + S 2 (V2 − V1 )] L EA [C{C 2 (U 2 − U 1 ) + SC (V2 − V1 )} + S{SC (U 2 − U 1 ) + S 2 (V2 − V1 )}] L =

EA 2 [C {C (U 2 − U 1 ) + S (V2 − V1 )} + S 2 {C (U 2 − U 1 ) + S (V2 − V1 )}] L

=

EA 2 (C + S 2 )[C (U 2 − U 1 ) + S (V2 − V1 )] L

EA [C (U 2 − U 1 ) + S (V2 − V1 )] L

S1− 2 =

EA [Cθ L

U 2 − U 1  Sθ ].   V2 − V1 

* Gaya batang / axial

S (1− 2) =

=

S (1−3) = =

S ( 2−3)

EA [Cosθ 1−2 L

U 2 − U 1  Sinθ 1− 2 ].   V2 − V1 

 4 P.L  EA [1 0]. 5 . EA  = 4 .P (tension) L  0  5 EA [Cosθ 1−3 L

[

EA 1/ 2 L

U 3 − U 1  Sinθ 1−3 ].   V3 − V1 

]

0 3 / 2 .  = 0 0

EA [Cosθ 2−3 = L

[

U 3 − U 2  EA Sinθ 2−3 ]. 1/ 2 = L  V3 − V2 

 4 P.L  2 − .  3 / 2 . 5 EA  = .P  0  5

]

(tension) 2.2 B E A M Struktur yang dirancang untuk mendukung beban lateral. Sehinngga utamanya dapat meneruskan bending, meskipun ada shear (sebagai konsekuensi logis) Tegangan Bending → Tegangan normal Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

41


Data teknis : E, I, L Pola model matematis → titik diluar node bagaimana defleksi (asumsi dengan interpolasi) Pada elemen ada 4 yang tidak diketahui → 4 suku Fisik Justifikasi : truss → dapat menurunkan ∈ yang konstan → sehingga T yang konstan. Beam Fungsi interpolasi (asumsi) : → Upaya untuk mendukung yang sebenarnya (yang didekati bukan fungsinya tetapi nilai numeriknya)

V ( x) = a1 + a 2 x + a3 x 2 + a 4 x 3 Justifikasi : di Beam

d 2 v M ( x) d 2v = → EI 2 = M ( x) EI dx 2 dx


42


Keseimbangan

dV = W .dx

W=

Keseimbangan

( M + dM ) − M = V .dx

dV dx

W = EI

dM = V .dx

d 4V dx 4

V=

dM dx

V = EI

d 3V dx 3

Pemisalan harus bisa memodelkan daerah beam tidak ada beban merata sehingga fungsi interpolasi turunan ke IV nya = nol Model umum ; Displacement =

∑ fi ( x)..di

Dimana fi(x) merupakan fungsi bentuk dan di merupakan Displacement dari node. Fungsi Interpolasi (asumsi)

V ( x) = a1 + a 2 x + a3 x 2 + a 4 x 3 ↓

V ( x) = f 1 ( x)V1 + f 2 ( x)θ 1 + f 3 ( x)V2 + f 4 ( x)θ 2

θ ( x) =

dV ( x) 1 1 1 1 = f 1 ( x)V1 + f 2 ( x)θ 1 + f 3 ( x)V2 + f 4 ( x)θ 2 dx

Gambaran penyelesaian pada aplikasi Beam digambarkan sebagai berikut : Suatu struktur Beam dengan berbagai beban .


43


W2(x) W1(x)

F

M

Langkah yang dilakukan sebagai berikut : 1. Diskrititasi (minimal) dengan cara sebagai berikut : -

Pada ujung-ujung beam diberi nodal

-

Pada setiap tumpuan diberi nodal

-

Pada diskontinuitas geometri diberi nodal

-

Pada beban terpusat diberi nodal

-

Pada diskontinuitas beban merata diberi nodal

2. Memberikan nomor nodal dan elemen dilakukan dari kiri ke kanan 3. Membuat tabel spesifikasi dari model yang dianalisa 4. Membuat model matematik atau persamaan kekakuan per elemen Dengan memberikan penomoran dof : Elemen K :

elemen

nomor dof

1

2

1 2

3 4

2

3

3 4

5 6

dan seterusnya. 5. Membuat matrik kekakuan total dengan mengasembly masing elemen


44


6. Dengan adanya beban merata, maka harus dibuat dulu beban ekivalensinya dengan cara sebagai berikut :

Y,V

P(x) M2 X,U

M1 Y1

1

2 L

Y2

Bentuk beban ekivalen :

 Y1  M  L {Fi } =  1  = ∫ P( x ).fi ( x ).dx  Y2  0 M 2  7. Indentifikasi kondisi batas menjadi dof aktif dan dof non aktif


45


8. Dengan persamaan kesimbangan total , tentukan dof aktif dengan metoda gauss eliminasi. 9. Menjawab pertanyaan dari problem. Prosedur yang dilakukan dalam struktur beam sebagai berikut : " Elemen Beam Spesifikasi -

2 node/elemen

-

2 dof / node

" Fungsi Interpolasi

V ( x) = f 1 ( x)V1 + f 2 ( x)θ 1 + f 3 ( x)V2 + f 4 ( x)θ 2

θ ( x) =

dV ( x) 1 1 1 1 = f 1 ( x)V1 + f 2 ( x)θ 1 + f 3 ( x)V2 + f 4 ( x)θ 2 dx

Shape Function 2

X X f1( x ) = 1 − 3  + 2  L L

3

 X 2  X3 + f 2 ( x ) = X − 2  L  L2   2

3

X X f3 ( x ) = 3  − 2  L L f4 (x) = −

X 2 X3 + L L2

Persamaan keseimbangan struktur : {f} = [K] {d} L

∫

dengan Elemen stiffness : k ij = EI fi" ( x ).f "j ( x ).dx 0


46


2.3 F R A M E → masing-masing elemen bisa menerima gaya kearah x dan y dan mampu mendukung momen sehingga dof = 3

mampu menerima : -

Beban lateral (bending)

-

Beban aksial

-

Beban terpusat/merata

-

Beban momen

Data teknis E, A, I, L, φ 2 node per elemen 3 dof pernode (u, v, θ) Konsep Seperti beam yang berorientasi φ terhadap x Dalam pemodelan matematis → kombinasi elemen truss dan beam I. Analisa → elemen tersebut terletak pada sumbu x (tapi bukan beam) (merupakan ide frame = truss + beam)

θ lokal = θ global


47


(karena diputar pada sumbu yang sama)

-

Diskritisasi

-

K (6 x 6) elemen

-

Assemble

-

Beban node ekivalen (karena ada beban merata)

-

B.C

-

Dof aktif

-

Jawab pertanyaan

•

Tidak ada tumpuan (dari soal terlihat kesetimbangan statis)

•

Tidak ada rigid body motion

•

Tumpuan → jadi B.C

•

Simetri

•

Sumbu simetri

BC dengan kesimetriannya (dari bentuk defleksi) V1 = θ1 = U3 = θ3 = 0 U2 = 0 V2 = ? (tidak nol/hampir nol) Penentuan BC -

BC yang lebih / kelewatan bisa membuat K tetap singular

-

Atau kalau tidak singular → maka proses kalkulasinya lebih panjang

Bidang simetri tengah Dua buah titik yang berjarak sama terhadap bidang simetri Pada bidang simetri → syarat : -

Struktur simetri


48


-

Beban simetri BC’s u=0 θy = 0 θz = 0

contoh soal

Analisa

Diskritisasi node 1 anggota frame aslinya (v, u, θ) sebagai truss hanya punya (u, v) dof aktif

[K

frame

]

12 / L2 x  x x EA   x = x L  − 6 / L x  0 −x 

x − 6/ L x x x x x x

8 6/ L

  u1     u 2  .v 2   − 6 / L  θ 2  12 / L2   v3  0 x x

Dof aktif

1 / 2 − 1 / 2 2L − 1 / 2 1 / 2 

[K truss ] = EA 


49


12 EI + EA x  L3 4L   x x  x x [ K struktur ] =   −6 x  L   0 − EA x  4L 

x −6 L x x x x x 8

0 − EA   4L   x  x   −6  L  x − 6 12 EI + EA  L3 4L 


50


BAB III

INTERPOLASI DAN INTEGRASI NUMERIK Interpolasi Lagrange merupakan pendekatan fungsi polynomial. Sedangkan Integrasi Gauss Quadrature merupakan suatu proses integrasi numerik dimana batas integral harus sudah dilihat melalui analisa numerik. Shape function → hubungan matematik dari fungsi interpolasi

φ = C 0 + C1θ + C 2θ 2 C 0    φ = 1 θ θ .C1  C   2

[

2

]

Tiga titik di

θ = θ 1 → φ = φ1 θ = θ 2 → φ = φ2 θ = θ 3 → φ = φ3 φ1 = 1 θ 1

[

C 0    θ 1 .C1  C   2

[

C 0   2  θ 2 . C1  C   2

φ2 = 1 θ 2

[

2

φ3 = 1 θ 3 θ 3

]

]

2

2 φ1  1 θ 1 θ 1  C 0       → φ 2  = 1 θ 2 θ 22 . C1  φ  1 θ θ 2  C  3 3  2  3 

C 0    .C1  C   2

]

−1

2 C 0  1 θ 1 θ 1  φ1       2 C1  = 1 θ 2 θ 2  .φ 2  C  1 θ θ 2  φ  3 3   2   3


51


φ = N 1 (θ )φ1 + N 2 (θ )φ 2 + N 3 (θ )φ 3 Curve fitting → suatu pendekatan Lagrange’s interpolation → pendekatan f polynomial FEM→ yang didekati bukan fungsinya karena kompleksnya tapi nilainya " 2 independent variables φ1, φ2 . . . . . . φ9 → diketahui

φ I ( x, y1 ) = N 1 ( x).φ1 + N 2 ( x).φ 2 + N 3 ( x).φ 3 N1 =

( x − x 2 ).( x − x3 ) ( x1 − x 2 ).( x1 − x3 )

N2 =

( x − x1 ).( x − x3 ) ( x 2 − x1 ).( x 2 − x3 )

N3 =

( x − x1 ).( x − x 2 ) ( x3 − x1 ).( x3 − x 2 )

φ II ( x, y 2 ) = N 4 ( x).φ 4 + N 5 ( x).φ 5 + N 6 ( x).φ 6 N1 =

( x − x5 ).( x − x 6 ) ( x 4 − x5 ).( x 4 − x6 )

φ III ( x, y 3 ) = N 7 ( x).φ 7 + N 8 ( x).φ 8 + N 9 ( x).φ 9 N 7 ( x) = Shape kurva :

φ ( x, y ) = N 1 ( y )φ I ( x, y1 ) + N 2 ( y )φ II ( x, y 2 ) + N 3 ( y )φ III ( x, y 3 ) N1 ( y ) =

( y − y 2 ).( y − y 3 ) ( y − y1 ).( y − y 3 ) ; N 2 ( y) = ( y1 − y 2 ).( y1 − y 3 ) ( y 2 − y1 ).( y 2 − y 3 )

N 3 ( y) = " Integrasi numerik Pada software → yang dipakai integrasi Gauss * GAUSS QUADRATURE Batas integrasi :harus sudah lihat : Analisa Numerik Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

52


" Mapping → merubah batas integral dengan menggunakan determinan Jacobi

4 titik

3 titik

titik gauss → dinyatakan dengan koordinat natural koordinat natural

faktor bobot

(½, ½, 0)

1/ 3

A

(0, ½, ½)

1/ 3

A

(½, 0, ½)

1/3

A

Koordinat natural

faktor bobot

( 1/3, 1/3, 1/3 )

-27/48

A

( 3/5, 1/5, 1/5 ) ( 1/5, 3/5, 1/5 )

25/48

A

( 1/5, 1/5, 3/5 ) Hubungan antara x dan interpolasi dalam natural : X = L1 X1 + L2 X2 + L3 X3 Kalau ada y Y = L1 Y1 + L2 Y2 + L3 Y3 Shape function pada elemen segitiga = koordinat natural Ni = Li Dalam pengertian koordinat natural sebagai interpolasi.


53


BAB IV

APLIKASI PADA PERPINDAHAN PANAS Disamping aplikasi untuk struktur, metode elemen hingga dapat juga diterapkan untuk perpindahan panas. Disini akan dibahas mengenai perpindahan aliran panas untuk 1-Dimensi dan juga untuk 2-Dimensi.

4.1 Steady State Uniaxial Heat Flow.

H(x)

Q1

Q2 x H(x)

X1

X2

q+(dq/dx).dx

q A

A+(dA/dx).dx dx Suatu daerah dengan luas penampang variable A(x) dengan aliran panas Q (energy/time) pada ujung dan sumber fluks panas, H(x) (energy/time-length), didistribusikan sepanjang arah x. Kesetimbangan energi dari differential element :

d A.q − H ( x ) = 0 dx


54


Fourier’s Law :

q = −k .

dT dx

k: thermal conductivity. ; T : Temperature

Substitusi Fourier Law ke differential equation :

d dT + H(x) = 0 A.k dx dx Bentuk varisional ekivalen dari persamaan diferensial : x2  d dT  δΠ = 0 = ∫  A.k + H ( x ).δT .dx x1 dx dx  

=

x2

∫x

1

x2 dT  d A.k  δT .dx + ∫x H ( x ).δT .dx 1 dx   dx

Integrasi suku pertama dan dikalikan dengan –1 didapat :

dT δΠ = − Ak δT dx

x2 x1

+∫

x2

A.k

x1

x2 dT d δT .dx − ∫ H ( x ).δT .dx x1 dx dx

dengan T : essential boundary condition (Dirichlet Boundary Condition) dT/dx: natural boundary value (Neumann Boundary Condition) untuk : Q =-A.k.dT/dx, maka

δΠ = δ (Q2 .T2 − Q1.T1 ) + ∫

x2

x1

A.k

x2 dT d δT .dx − ∫ H ( x ).δT .dx x1 dx dx

Functional untuk 1 dimensi problem perpindahan panas adalah :


55


Π = (Q2 .T 2 − Q1.T1 ) + ∫

x2

x1

2

x2  dT  A.k   .dx − ∫x H ( x ).T ( x ).dx 1  dx 

Newton’s Law of cooling, aliran panas konveksi pada batas 1 dan 2: Q1 = Q1c = h.A.(T∝ - T1)

dan

Q2 = Q2c = h.A.(T2 - T∝ )

T∝ : temperatur ambient ; h : koefisien perpindahan panas konveksi. Energi yang ditambahkan dengan konveksi pada daerah panjang dx : H(x).dx = h.(P.dx)(T∝ - T(x)) H(x) = h.P.(T∝ - T(x)) 4.2 MODEL ELEMEN HINGGA UNTUK ALIRAN PANAS 1-DIMENSI. Functional :

Π = (Q2 .T 2 − Q1.T1 ) + ∫

x2

x1

2

x2  dT  A.k   .dx − ∫x H ( x ).T ( x ).dx 1  dx 

model elemen : dua nodal heat flow element.

L

Q1

i

j

Q2 X

1

2

1. Asumsi fungsi yang menyatakan variable dependen melalui elemen. Variasi linear temperatur : T = [N] . {qt}

[N] = [N1

N2] =

 Xj − X .   X j − X i


X − Xi X j − Xi

 .  56


Xj – Xi = L

;

{qt} = [Ti

dT 1 = [− 1 1].{q t } dx L

atau

Tj]T

dT = [B ].{q t } dx

Substitusi :

Π = −[Qi

− Q j ].{q t } +

A.k 2

− Q j ].{q t } +

1 − 1 xj A.k {q t }T  .{q t } − ∫ H ( x ).[N ]{q t }.dx  xi 2.L − 1 1 

atau :

Π = −[Qi

Dengan Ritz procedure

T T ∫x {q t } [B ] [B].{q t }.dx − ∫x xj

xj

i

i

H ( x ).[N ]{q t }.dx

dΠ/d{qt} = 0, maka governing equation for the

single element :

A.k  1 − 1 .{q t } = .L − 1 1 

 Qi  x j  + ∫x H ( x ).[N ].dx  − i Q j  

atau : [ kcd ] . {qt} = {Qt }N + {Qt}H. dimana : [ kcd ] = element conduction matrix ; { qt }

= nodal temperature vector

{ Qt }N = nodal heat flow vector; { Qt }H = nodal heat flow vector equivalent to the distributed flux. Assembly elemen, dgn Rayleigh-Ritz Procedure thd functional seluruh region : [ Kcd ] . {rt} = {Rt }N + {Rt}H. dimana :

[ Kcd ] = assembled conduction matrix; { rt }

= assembled nodal temperature vector

{ Rt }N = nodal heat flow at boundary and node sources Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

57


{ Rt }H = distributed heat flux vector. 4.3 ONE-DIMENSIONAL HEAT FLOW WITH CONVECTION

T∝H , hH T∝L

T∝R

x

hL

hR L

Persamaan kesetimbangan : [ kcd ] . {qt} = {Qt }N + {Qt}H. asumsi konveksi terjadi hanya pada nodal local 1.

{Qt }N

h.A(T∞ L − T1 ) = = − Q2  

h.A.T∞ L  h.A.T1   −   − Q2   0 

h.A.T∞ L  1 0T1  =  − h.A.   0 0T 2   − Q2  atau :

{Qt }N = {Qcv}L.- [ kcv ]L .{qt}

Jika ujung kanan mempunyai konveksi., kemudian dengan subtitusi Q2 = h.A. (T2 - T∝R) didapat :

{Qt }N atau :

 Q1  0 0T1  =  − h.A.   0 1T2  h.A.T∞ R 

{Qt }N = {Qcv}R.- [ kcv ]R .{qt} dimana :

{Qcv} : Vektor aliran panas konveksi; [Kcv] : Matrik konveksi


58


Fluks panas terdistribusi :

{Qt }H = ∫0 H ( x )[N ]T .dx = ∫0 h.P.(T∞H L

{Qt }H

L

L

L

0

0

− T ( x )).[N ]T .dx

atau :

= h.P.T∞H ∫ [N ]T .dx − h.P ∫ [N ]T .[N ].dx.{q t }

Matrik fungsi bentuk dalam koordinat local :

 x [N ] = 1 −  L

x L 

Fluks terdistribusi :

{Qt }H

=

atau :

h.P.L.T∞ 2 H

1 h.P.L 2 1 T1   −   6  1 2 T2  1

{Qt}H = {Qcv}H – [kcv]H .{qt}.

Asumsi single elemen dengan konveksi pada sisi batas kiri dan sepanjang elemen dan aliran panas Q2 pada batas kanan. [ kcd ] . {qt} = {Qt }N + {Qt}H. = {Qcv}L.- [ kcv ]L .{qt} + {Qcv}H – [kcv]H .{qt}.

[k cd ]

T1  h.A.T∞L  1 0 T1  h.P.L.T∞H =  − h.A   T  + − T Q 0 0 2   2   2  2 

1 h.P.L 2 1 T1   −   6  1 2 T2  1

direorganisir : (konveksi pada sisi kiri)

[[k cd ] + [k cv ]L + [k cv ]H ]{q t } = {Qcv }L + {Qcv }H (Konveksi pada sisi kanan ) :

[[k cd ] + [k cv ]R + [k cv ]H ]{q t } = {Qcv }R + {Qcv }H Contoh : Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

59


Aliran panas dalam sirip segiempat seperti pada gambar dimodelkan sebagai problem 1 dimensi. Sisi kiri sirip dipertahankan pada temperatur 2000C dan semua permukaan diekspos pada temperatur ambien 500C. Koefisien konveksi untuk semua permukaan 0.02 W/cm2.0C. konduktifitas termal bahan 4 W/cm.0C. Pertama menggunakan model elemen tunggal dan kemudian model dua-elemen , estimasikan temperatur pada ujung sirip dan panas yang hilang. Penyelesaian :

20 cm 5 cm

20 cm

T∞=50 0C 2000C

Q1

X

1

Model satu-elemen Matrik konduksi :

1 − 1 100.4  1 − 1  20 − 20 = = L − 1 1  20 − 1 1  − 20 20 

[k cd ] = Ak 

Elemen dengan konveksi pada sisi kanan. Matrik konveksi untuk aliran panas dari sisi kanan adalah : Program Semi-Que IV Fakultas Teknik—Jurusan mesin Universitas Brawijaya

60


[k cv ]R

0 0  0 0  0 0  = h.A  = 0.02(100)   =  0 1  0 1  0 2 

Matrik konveksi untuk aliran panas dari semua sisi :

[k cv ]H = hPL 

2 1 0.02(50)(20) 2 1 6.7 3.3 = 1 2 = 3.3 6.7 6 1 2 6    

Vektor konveksi untuk konveksi sisi kanan :

{Qcv }R

Q1  Q1     Q1  = = =  h.A.T∞R  0.02(100)(50) 100

Matrik konveksi untuk sisa sisi bebas :

{Qcv }H = h.P.L.T∞H   = 0.02(50)(20)(50)   =  2

1 1

1 1

2

500  500

Asembly persamaan matrik aliran panas komplit :

[[k cd ]+ [k cv ]R + [k cv ]H ]{q t } = {Qcv }R + {Qcv }H  26.7 − 16.7   T1  Q1 + 500 − 16.7 28.7  T  =  600    2   

kondisi batas esensial, T1 =200

0   T1  200  1 − 16.7 28.7  T  = 600   2    solusi untuk T2 : -16.7 (200) + 28.7 T2 = 600

T2 = 137.3 0C.

Aliran panas Q1 dalam sisi kiri didapatkan : 26.7T1 - - 16.7T2 = Q1 + 500 26.7(200) – 16.7(137.3) = Q1 + 500

Q1 = 2547 W.

Aliran panas rata-rata dalam elemen : Q1 = -k dT/dx = - k[B] {qt} = −

=−

T k [− 1 1] 1  L T2 

200  4 2 [− 1 1]  = 26.3W / cm 137 . 3 20  


61


4.4 PERPINDAHAN PANAS DAN ALIRAN FLUIDA 2-DIMENSI •

Governing equation.

Y Laplace eq. :

∇ 2T = 0

Fourier eq. :

q cd X = −k

nˆ

∂T ∂x

T∞

T(x,y)

∂T ∂y ∂T q n = −k ∂n q cdy = −k

atau

:

Newton’s Law

of cooling :

S

X

qCV = h.A (T - T∞ )

Galerkin Approximation :

∫AWi .∇ T .t.dx.dy = 0 2

dalam bentuk lain :

∇.W i ∇ T = ∇ W .∇ T + W i ∇ 2T

sehingga disubsitusi menjadi :

t ∫ ∇.W i ∇ T .dx.dy − t ∫ ∇ W i .∇ T .dx.dy = 0 A

A

Dengan Gauss Theorem :

! t ∫ ∇.W i ∇ T .dx.dy = t ∫ W i .∇ T .nds A

S

, maka

! t ∫ W i .∇ T .nds − t ∫ ∇ W i .∇ T .dx.dy = 0 S

atau :

t ∫ Wi . S

A

 ∂W i ∂T ∂W i ∂T  ∂T + ds − t ∫  .dx.dy = 0 A ∂y ∂y  ∂n  ∂x ∂x

Interpolation formula : T = [N] . {qi}

dan

Wi = Ni

Sehingga :

 ∂[N ]T ∂[N ] ∂[N ]T ∂[N ] ∂T ∫Se [N ] . ∂n ds − t ∫Ae  ∂x ∂x + ∂y ∂y .dx.dy .{q i } = 0   T


62


Persamaan Elemen :

Y

Keseimbangan energi : qcdn = qcv

(q

ˆ+q

cd x .i

cd y

T∞

qb

qcv

)

Se’’

. ˆj .nˆ = q cv

qcdn

∂T .ds = h.t.ds.(T − T∞ ) ∂n ∂T h − .ds = (T − T∞ ).ds ∂n k

Se’

− t .k.

dan

∂T/∂n=0

untuk : T = [N] . {qt}

− subsitusi :

X

∂T h h .ds = [N ].{q t }.ds − T∞ .ds ∂n k k

T [ N ] .[N ].ds.{q t }+ h.∫ [N ] Se ' Se '

T

− h∫

[N ] Se ''

.T∞ .ds + k.∫

 ∂[N ]T ∂[N ] ∂[N ]T ∂[N ] − k.∫  + . . .dx.dy .{q t } = 0 Ae ∂x ∂y ∂y   ∂x

T

.

∂T .ds ∂n

dalam bentuk persamaan elemen : [KT] . {qt} = {Qcv} + {Qb} “Thermal stiffness” matrix : [KT] = [kcdx] + [kcdy] + [kcv].

 ∂[N ]T ∂[N ] [k cdx ] = k.∫Ae  . .dx.dy ∂x   ∂x  ∂[N ]T ∂[N ] [k cdy ] = k.∫Ae  ∂y . ∂y .dx.dy   [k cv ] = h.∫ [N ]T .[N ].ds Se '

Convection boundary vector Se’ :

[Qcv ] = h.∫Se' [N ]T .T∞ .ds

Applied heat boundary vector, Se’’ :


[Qb ] = k.∫Se'' [N ]T . ∂T .ds ∂n

63


BAB V

ANALISA TEGANGAN AXISYMMETRIC Sekelompok problem yang ada pada kenyataannya meliptui gaya dan domainnya dalam tiga dimensi, tetapi akan diupayakan mereduksi secara matematik menjadi dua dimensi. Problem-problem tersebut disebut dengan axisymmetric problems, dan dikarakteristikan dengan putran solid dan sifat-sifat material dan beban yang tak berubah sepanjang sekeliling putaran.Gambar berikut adalah putaran solid, dengan elemen yang akan digunakan pada diskrititasi dari kontinum yaitu toroid dengan penampang segitiga. Suatu hal yang penting untuk merealisasikan pada axisymmetric

problems, perpindahan dalam kontinum dapat terjadi hanya dalam arah radial dan aksial; perpindahan tidak dapat terjadi

dalam

arah

sirkumferensial, sebagai akibat hal tersebut, menjadi biasa menggunakan sistem koordinat silinder dalam mengembangkan persamaan elemen umum, seperti pada gambar berikut.

Sumbu putaran

Program Semi-Que IV Fakultas Teknik Jurusan Mesin Universitas Brawijaya

64


Putaran benda dari elemen toroidal. Sistem Koordinat

w1 = w1.kˆ

z

2

u1 = u1.eˆ r

1

θ

kˆ

eˆ θ

r

3

eˆ r

Komponen tegangan koordinat silinder untuk keadaan axisymmetric.

z

σz

θ σθ

dθ

τrz

dz

σr

dr

r


65


5.1 Persamaan dasar untuk elemen Persamaan elemen secara umum untuk analisa tegangan kontinum tiga dimensi identik dengan bentuk : T T ∫ [B] .[C][. B].dΩ.{q} = ∫Ω [B] .[C].{εT }.dΩ + {Q}NF + {Q}T + {Q}BF

Ω

walaupun aplikasi persamaan ini untuk elemen tiga dimensi adalah identik dengan konsep elemen dua dimensi, upaya lebih besar karena perpindahan tambahan pada setiap nodal dan dimensi dalam tiga variabel. Integral garis dan luasan dari elemen problem bidang sekarang menjadi integral permukaan dan volume. Dalam persamaan diatas, jika diaplikasikan ke kontinum tiga dimensi didefinisikan kembali sebagai berikut : Matrik kekakuan

[k ] = ∫Ω [B]T .[C][. B].dΩ Vektor beban nodal temperatur :

[Q]temp = ∫Ω [B]T .[C].{εT }.dΩ Vektor gaya nodal {Q}NF = gaya-gaya aplikasi pada nodal Vektor traksi permukaan

[Q]T = ∫A [N]T .{T}.dΩ Vektor Gaya bodi

[Q]BF = ∫Ω [N]T ..{Bf }.dΩ


66


5.2 Persamaan Elastisitas Axisymmetric Pada Axisymmetric, semua persamaan harus menjadi bebas dari θ dan semua perpindahan harus berada dalam bidang rz. Hubungan perpindahan regangan dalam koordinat silinder pada problem khusus sebagai berikut.

∂u εr = ∂r

∂w ; εz = ∂z

u ; εθ = r

;

γ rz =

∂w ∂w + ∂r ∂r

dalam bentuk matrik :

∂  ∂r  εr   ε   0 {ε} =  z  =  1  εθ   γ rz   r ∂   ∂z

 0 ∂  ∂z . u  = [℘]. u       w  w  0 ∂  ∂r 

Hubungan untuk material isotropik :

ν ν 0   ε  1 − ν  σr  1     r  σ  1  ν − ν ν 1 0 E  z       εz  − α ∆ x . T  =       ν ν 1− ν 0  εθ   σθ  (1 + ν)(1 − 2ν)   1   1 − 2ν   0 0 τrz  0   0  γ    2    rz  atau {σ} = [C] . ({ε} - {ε}T) Vektor regangan termal didefinisikan sebagai

 εr  1 ε  1  z {ε}T =   = α∆T    εθ  1 γ rz  0 Program Semi-Que IV Fakultas Teknik Jurusan Mesin Universitas Brawijaya

67


Fungsi perpindahan elemen Nodal dari elemen toroidal sebenarnya adalah lingkaran konsentrik yang lewat melalui puncak penampang segitiga. Koordinatnya adalah r dan z. Spesifikasi perpindahan radial , u, perpindahan aksial, w, posisi radial, r, dan posisi aksial, z dari suatu toroidal yang akan didefinisikan dengan formulasi interpolasi linear dalam

koordinat natural dan sifat-sifat nodal.

u = L1u1 + L2u2 + L3u3 w = L1w1 + L2w2 + L3w3 r = L1r1 + L2r2 + L3r3 z = L1z1 + L2z2 + L3z3 dimana :

L1+ L2 + L3 = 1

dalam bentuk matrik

1   r  = z   

1 1 r r 1 2 z1 z 2

1   L1    r3  L 2  z3   L3 

invers matrik :

 L1   a1 1    L 2  = a 2 det L  a 3  3

b1 b2 b3

c1  1   c 2   r  c3  z 

dimana : a1 = r2z3 – r3z2 ;

a2 = r3z1 – r1z3 ;

a1 = r1z2 – r1z2 ;

b1 = z2 – z3

;

b2 = z3 – z1

;

b3 = z1 – z2

;

c 1 = r3 – r2

;

c 2 = r1 – r3

;

c 3 = r2 – r1

;

dan det = (r1 – r3)( z2 – z3) - (r2 – r3)( z1 – z3) = 2 x luas segitiga.


68


Vektor fungsi perpindahan :

 u  L1 0 L 2  = w   0 L1 0

0 L2

 u1  w   1 0   u 2  .  = [N ].{q} L3  w 2   u3     w 3 

L3 0

Hubungan regangan dengan vektor dof :

{ε} = [℘][N].{q} = [B].{q} derivatif koordinat natural :

∂L1 ∂  a1 + r.b1 + z.c1  b1 =  = ∂r ∂r  det  det

dan seterusnya.

Selanjutnya matrik [B] menjadi :

 b1 0 1  [B] =  L*1 det   cr  1

0 c1

b2 0

0 c2

b3 0

0

L*2 r c2

0

L*3 r c3

b1

b2

0 c3   0 b3 

dimana : L1* = a1 + r.b1 + z.c1 ;

L2* = a2 + r.b2 + z.c2 ;

L3* = a3 + r.b3 + z.c3

Matrik kekakuan

[k ] = ∫Ω [B]T .[C][. B].dΩ Metode pendekatan yang sederhana [Zienkiewics] dinyatakan sebagai berrikut :


69


[B] = [B(r, z )] r +r +r dimana ; r = 1 2 3

z + z 2 + z3 ; z= 1

3

3

volume : V = 2.π.r.A Matrik kekakuan elemen :

[k ] = 2.π.r.A[B]T .[C].[B] Vektor beban nodal temperatur :

1 1 T [Q]temp = ∫Ω B .[C].{εT }.dΩ = 2.π.r.A.α.∆T. B .[C]  1 0

[]

[]

Vektor gaya nodal {Q}NF = gaya-gaya aplikasi pada nodal {Q}NF = [F1r

F1z

F2r

F2z

F3r

F3z]T

Vektor traksi permukaan

[Q]T = ∫A [N]T .{T}.dΩ = 2.π.∫S r.[N]T .

Tr  .ds Tz 

 rL1  0  rL [Q]T = 2.π.∫S  2  0  rL3   0

0  rL1  0   Tr  . .ds rL 2  Tz  0   rL3 

r dalam istilah koordinat natural : r = L1r1+ L2r2+ L3r3 Vektor Gaya bodi

[Q]BF = ∫Ω [N]T ..

Br  T  Br  .dΩ = 2.π.∫A r.[N ] . .dA B z   Bz 


70


REFERENSI 1. Grandin Hartley, Jr.,1986, “ Fundamentals of the Finite Element Method”, Macmillan Publishing Company, New York. 2. Yang, T.Y., 1986,”Finite Element Structural Analysis”, PrenticeHall,Inc,Englewood Cliffs. 3. Buchanan, George R.,1995, “Finite Element Analysis, Schaum’sOutline Series, McGraw-Hill International Editions 4. Bathe Klaus-Jurgen, 1996, “Finite Element Procedures”, Prentice Hall International Editions, Inc, USA. 5. Hughes Thomas J.R.,1987, “ The Finite Element Method”, Prentice-Hall Inc, New Jersey 6. Segerlind L., J., “Applied Finite Element Analysis”, John Willey & Son,Inc.


71

Malang, September Penulis

Recommend Documents