Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
Pengenalan Program Komputer LD-FEM (Large Deformation-Finite Element Method) Sebagai Alat Bantu Pemodelan dan Simulasi Material Hyperelastic Sugeng Waluyo1) Jurusan Teknik, Universitas Jenderal Soedirman Jl. Mayjen Sungkono Km.5 Blater, Purbalingga, Indonesia 53371
LD-FEM is an open source computer program working on the basis of finite element method (FEM) which is aimed to model and simulate hyperelastic materials. It is mainly designed as a solver running on Fortran without embedded pre- or/and post-processing units. Recently, the solver is compatible only for reading and writing in the Gmsh pre- and post-processing software environment (www.guez.org). The FEM kinematics formulation applied here is relied on the 4-nodes tetrahedral element using linear shape functions. Meanwhile, the Newton-Raphson iterative scheme is performed to solve the nonlinear equation arises from the use of Total Lagrange framework on the basis of large deformation theory. To obtain the material tangent stiffness directly from the strain energy density function, the Gill-Murray theory of numerical 2nd derivative is used here. Finally, in order to validate the LD-FEM performance, an example of simple uniaxial compressive test is simulated using the Saint Venant-Kirchhoff strain energy density function. The results confirm the capability of LD-FEM to capture nonlinear behavior in the large deformation either with analytical or numerical 2-nd derivative of the Saint Venant-Kirchhoff model . Keywords: LD-FEM, Gmsh, finite element, linear tetrahedral element, hyperelastic material, Total Lagrange formulation
1. PENDAHULUAN Proses penyusunan kode program LD-FEM dimulai pada pertengahan Juni 2009. Tujuan utamanya adalah menghasilkan produk berupa program komputer sumber terbuka (open source) yang mampu mensimulasikan material hyperelastic, sebagai contoh karet dan jaringan tubuh lunak (soft tissue), secara 3dimensi. Sejauh pengetahuan penulis, hingga saat ini belum ada program FEM sejenis yang bersifat sumber terbuka di Indonesia. Saat ini program LD-FEM masih diperuntukan untuk kalangan akademisi sebagai sarana pembelajaran dan penelitian yang behubungan dengan pemodelan dan simulasi material hyperelastic seperti karet dan soft tissue. Struktur program LD-FEM didesain mengikuti teori dasar large deformation FEM (Taylor, 1999), tanpa mengalami optimasi pemrograman. Dengan struktur tersebut diharapkan pengguna dapat dengan mudah melihat kesesuaian antara teori dan kode program. Sesuatu yang saat ini sulit diperoleh dari open source lain, sebagai contoh FEAP (Taylor, 2000), akibat proses evolusi dan optimalisasi kode program. Harapannya ada dorongan bagi pengguna untuk memodifikasi program ini sesuai kebutuhan mereka. Pada artikel ini akan disajikan secara komprehensif struktur utama program LD-FEM. Porsi terbesar ada pada pembahasan teknis formulasi kinematika deformasi besar (large deformation kinematics) yang digunakan. Struktur program sendiri terdiri dari 3 Fortran subrutin (subroutine) utama dan didukung oleh
setidaknya 5 subrutin pendukung. yang akan dijelaskan secara singkat kegunaan dan parameter input-outputnya. Bagian terbesar makalah ini akan didedikasikan guna menurunkan formulasi FEM pada kasus deformasi besar dengan menggunakan kerangka deformasi Total Lagrange. Seluruh formulasi akan disajikan dalam bentuk notasi matriks-vektor meskipun di beberapa bagian tidak dapat dihindari penggunaan indeks tensor. Namun demikian, perlu diperhatikan bahwa sejatinya indeks atau tensor notasi lebih umum digunakan dalam kasus ini. Pembaca yang tertarik akan hal ini dapat merujuk pada Taylor (1999) dan Sussman dan Bathe (1987). Selanjutnya, untuk memperoleh solusi persamaan nolinear digunakan Metoda Newton-Raphson dengan kendali beban (load-control). Nantinya formulasi FEM yang diperoleh akan diterapkan pada elemen tetrahedral 4 nodal dengan fungsi bentuk (shape function) linear. Fungsi energi densitas regangan (strain energy density) diberikan sebagai subrutin tersendiri untuk memudahkan pengguna mengkaji model material hyperelastic yang mereka miliki. Selain evaluasi secara analitis pada tangent stiffness atau turunan kedua fungsi strain energy, pada makalah ini juga diberikan teknik evaluasi turunan kedua secara numerik terhadap fungsi tersebut. Sebagai ilustrasi kinerja akan digunakan model hyperelastic paling sederhana yaitu Saint Venant-Kirchoff hyperelastic model. Pada penutup, pengujian tekan pada spesimen silinder akan dimodelkan dan disimulasikan dengan
SNTTM X | 825
Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
LD-FEM. Hasil pemodelan dan simulasi akan dibahas secara lengkap untuk membuktikan kinerja secara kuantitatis kode program LD-FEM. Dua hal akan dijadikan referensi disini yaitu kemampuan mewakili hubungan nolinear beban-perpindahan dan kecepatan konvergensi pada saat mengevaluasi turunan kedua dengan metoda numeris. 2. STRUKTUR UTAMA PROGRAM LD-FEM Struktur program LD-FEM tersusun oleh 3 (tiga) subrutin atau sub program utama yaitu PRE_SOLID_FSFEM(), NONLIN_SOLVING(), dan POST_SOLID_FSFEM(). Disamping ketiga subrutin tersebut terdapat pula sejumlah subrutin lain yang melakukan tugas operasii matematika seperti aljabar matriks dan LU Decomposition. 2.1
Subroutine PRE_SOLID_FSFEM() Tugas utama dari subroutin ini adalah mendapatkan data dari proses pembentukan mesh yang dilakukan Gmsh terhadap geometri benda pejal pengguna. Dalam melakukan tugasnya, subrutin ini akan membaca file hasil proses mesh tersebut yang berekstensi *.msh. Empat data utama yang diperoleh berupa jumlah nodal, koordinat nodal, jumlah elemen dan nodal penyusun elemen. Subrutin ini dirancang untuk tidak mengalami perubahan apabila pengguna berhendak menambahkan fungsi perpindahan internal (bubble function) pada elemen tetrahedral. Referensi akan jenis fungsi perpindahan ini dan aplikasinya pada elemen tetrahedral dapat dilihat pada Taylor (1999). 2.2
Subroutine NONLIN_SOLVING() Subrutin ini tersusun dari subroutin-subroutin pendukung seperti TET4nodes(), SOLVER(), LOAD_SOLID_FSFEM(), dan CONSTRAINT_SOLID_FSFEM(). Subroutin ini bertugas melakukan tahapan yang diperlukan bagi iterasi menggunakan Metoda Newton-Raphson dengan beban sebagai parameter yang mengalami penambahan (increments). Secara umum tahapan tersebut berlangsung sebagai berikut secara berurutan: 1. Memanggil LOAD_SOLID_FSFEM() untuk mendapatkan informasi beban yang diaplikasikan pada setiap nodal. Beban tersebut kemudian dibagi menjadi beberapa tahapan beban yang lebih kecil sesuai keinginan pengguna. 2. Memanggil TET4nodes() untuk setiap tahapan pembebanan dan iterasi. Data dari TET4nodes() digunakan untuk menyusun matriks kekakuan elemen berdasarkan formulasi elemen tetrahedral 4-nodal menggunakan fungsi perpindahan linear. Matriks tersebut selanjutnya akan disusun menjadi matriks kekakuan global 3. Memanggil data spesifikasi tumpuan (constraint) pada masing-masing nodal untuk setiap tahapan beban dan iterasi melalui
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
CONSTRAINT_SOLID_FSFEM(). Data tersebut digunakan sebagai dasar modifikasi matriks kekakuan global yang telah disusun pada Tahap 2. 4. Hasil modifikasi matriks kekakuan global pada Tahap 3 selanjutnya diselesaikan oleh SOLVER() menggunakan algoritma LU Decomposition. Hasil ini selanjutnya dicek norm-nya untuk menghitung kriteria konvergensinya. 5. Jika norm pada Tahap 4 konvergen maka dilanjutkan ke tahapan beban selanjutnya. Sebaliknya jika tidak konvergen maka dilanjutkan ke iterasi berikutnya. 2.3
Subroutine POST_SOLID_FSFEM() Bagian ini mempunyai tugas untuk menuliskan pada sebuah file data hasil simulasi dari LD-FEM sehingga dapat dibaca oleh Gmsh. File yang digunakan mengunakan ekstensi *.msh. File tersebut nantinya akan berisi harga deformasi untuk setiap nodal pada setiap tahapan beban. Dalam hal ini vektor deformasi didefinisikan melalui koordinat kartesian global x = x ( x, y, z ) .
3. PERSAMAAN ENERGI MEDIA KONTINUUM UNTUK KASUS DEFORMASI BESAR (LARGE DEFORMATION) Pada kerangka koordinat acuan awal (references coordinat) tetap x = x ( x, y, z ) , formulasi kinematika finite element (FE) untuk kasus deformasi besar dapat diturunkan dari persamaan energi pada sebuah volume dV berikut
∂ν T FM T dV = ν T P fe x t . ∂ x V
∫
(1)
T
Disini ν = ν x ν y ν z dan fe x t berturut-turut adalah vektor perpindahan maya dan gaya ekternal. Sedangkan definisi masing-masing matriks dijabarkan sebagai
∂ ∂ ∂x ∂y ∂ = 0 0 ∂x 0 0
0 ∂ ∂ ∂ 0 ∂x ∂y ∂z 0 0 0 ∂ ∂ ∂ 0 0 0 0 ∂x ∂y ∂z
∂ ∂z
0
0
0
0
0
(2)
SNTTM X | 826
Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
1 + F11 0 0 F12 FM = 0 0 F 21 0 0
0
0
0
0
F21
0
F21
0
0
F31
F21
F21
0
0
0
F21
F21
0
F21
0
0
F21
F21
F21
0
0
0
F21
F21
0
F21
0
0
F21
F21
F13
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
0
Ct =
0 F21 F21 0 F21 F13 0 F21
T = T11 T22 T33 T23 T13 T12
(4)
dengan komponen FM berkorelasi dengan komponen
F berikut ∂u x 1 + ∂x ∂u y F= ∂x ∂u z ∂x
∂u x ∂u x ∂y ∂z ∂u ∂u y 1+ y ∂y ∂z . ∂u z ∂u 1+ z ∂y ∂z
(8)
Perlu dicatat bahwa perhitungan C t dapat dilakukan secara analitik maupun numerik untuk menghasilkan matriks dimensi (9x9). Pembaca yang tertarik pada proses penyederhanakan menjadi bentuk (6x6) pada Persamaan (6) dapat merujuk pada berbagai referensi, salah satunya Bower (2008).
(3) T
∂ W ( E) . ∂E ∂E
3. IMPLEMENTASI METODA NEWTON RAPHSON Metoda Newton-Raphson yang digunakan pada LD-FEM menggunakan pembebanan bertahap (load increment). Berbeda dengan perpindahan bertahap (displacement increment), pembebanan bertahap membagi beban total yang didefinisikan pada obyek dengan interval beban tertentu. Sebagai ilustrasi dapat dilihat pada grafik proyeksi 1-D dari residual vektor r pada Gambar 1. Linearisasi persamaan (1) dapat dilakukan dengan menggunakan persamaan residual
∂ν T FmCt EvdV − ν T Pfext . ∂ x V
r=∫ (5)
(9)
Berdasarkan (9), solusi Newton-Raphson dapat didefinisikan pada koordinat global sebagai berikut ∂r n rn + ∆ ( ∂u ∂x ) = 0 (10) n ∂ ( ∂u ∂x )
Perpindahan titik material (material point) dalam masing-masing arah didalam kerangka acuan
dengan
koordinat
∆ ( ∂u ∂x ) = ( ∂u ∂x )n +1 − ( ∂u ∂x )n
adalah
u = u(ux, uy , uz )
.
Matriks
P adalah matriks berukuran (3x9) yang digunakan sebagai penghubung agar traksi yang ada dipermukaan domain bersesuaian dengan vektor perpindahan maya υ . Sementara itu, analogi dengan permasalahan deformasi kecil (small deformation), tegangan 2nd Piola-Kirchhoff T dapat pula didekati sebagai perkalian antara sebuah modulus kekakuan dan regangan sebagai berikut
T ≈ C t Ev
(6)
(11)
Dalam konteks waktu t , perubahan beban pada (10) diasumsikan berjalan dari kondisi awal t = n menuju kondisi terkini pada t = n + 1 . Secara prinsip turunan terhadap residual dapat diberikan dengan aturan berikut
∂ r ( FM ) ∂ r ( E v ) ∂r = + ∂ u ∂u ∂u ∂ ∂ ∂ ∂x ∂x ∂x
(12)
t
dimana C adalah matrik modulus kekakuan (tangent stiffness) berukuran (6x6) dan
Ev = [ E11 E22 E33 E23 E13 E12 ]
T
(7) Disini
komponen E v bersesuaian dengan komponen dari matriks regangan Green-Lagrange E = 0.5 ( F T F − I ) dengan I adalah matriks satuan. Matriks modulus kekakuan sendiri dapat diperoleh dari turunan kedua fungsi densitas energi regangan (strain energy density function) W ( E ) yakni
SNTTM X | 827
Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
4. FORMULASI FINITE ELEMENT Persamaan (14) selanjutnya dapat digunakan untuk memfasilitasi elemen tetrahedral 4-nodal dengan fungsi bentuk (shape function) linear.
fext solution for time tn+1
(x 4 , y4 , z4 ) 4
( fext )n +1
( x 3 , y3, z3 )
∂r n ∂ ( ∂u ∂x )
(fext )n
3
( x1, y1 ,z1 )
n
1
Z Y
∂u
∂u ∂u ∂x ∂x ∂x n +1 n Gambar 1. Ilustrasi 1-D metoda NewtonRaphson
dimana masing-masing suku di kanan tanda sama dengan diturunkan berdasarkan variabel fungsi. Jadi, pada suku pertama hanya kontribusi dari komponen FM yang dianggap variabel dengan komponen EV dianggap konstanta dan sebaliknya. Perlu dicatat bahwa perhitungan turunan harus dilakukan berdasarkan sistem persamaan yang terbentuk pada (9). Hal ini dilakukan karena aturan turunan pada tensor hanya berlaku pada matriks tensor yang bersangkutan dan tidak pada matriks persegi panjang seperti FM . Lebih jauh lagi, Persamaan (12) masih membutuhkan modifikasi guna memperoleh vektor ∆u . Berdasarkan prinsip dasar kerangka acuan awal tetap, ∂x tidak berubah sepanjang proses iterasi dalam n . Konsekuensi dari hal ini, Persamaan (12) dapat ditulis sebagai
∂u n +1 ∂un ∂∆u ∆ ( ∂u ∂x ) = − = . ∂x ∂x ∂x
(13)
Sekarang, menggunakan aturan (12), dan definisi pada (13), hasil akhir perhitungan (10) adalah
∂ν T * ∂∆u ∂ν T dV T FM C t FM dV + = ∫ M ∫ V ∂x V ∂x ∂x ∂ν T FM TdV + ν T Pf ext −∫ x ∂ V (14) dengan T * adalah matriks (9x9) berisi komponen dari T yang mewakili hubungan perkalian matriks
Ct E v ( ∂∆u ∂x ) .
(x 2 , y2 , z2 ) 2
X Gambar 2. Elemen tetrahedral 4-nodal dengan formulasi linear
Untuk elemen tetrahedral 4-nodal pada Gambar 2, koordinat pada setiap domain x = x(x, y, z) diwakili oleh koordinat keempat nodalnya berikut
x = N1x1 + N 2 x 2 + N 3x 3 + N 4x 4
y = N1y1 + N 2 y 2 + N3 y3 + N 4 y 4 . z = N1z1 + N 2z 2 + N 3z 3 + N 4 z4
(15)
Dalam hal ini fungsi bentuk didefinisikan oleh N1 = ξ ,
N 2 = η , N 3 = ϕ , dan N 4 = 1 − ϕ − η − ξ . Koordinat parametrik ini mendefinisikan koordinat lokal pada seluruh tetrahedral x = x (ξ, η, ϕ) . Menggunakan prinsip isoparametrik, perpindahan u didefinisikan dalam orde yang sama dengan Persamaan (15) sebagai berikut
u x = N1u x1 + N 2 u x 2 + N3u x3 + N 4u x 4
u y = N1u y1 + N2u y2 + N3u y3 + N4u y4
(16)
u z = N1u z1 + N2 u z2 + N3u z3 + N4 u z4 atau
u = Nu nodal
(17)
u = u x u y u z
(18)
dengan
unodal = ux1 uy1 uz1 ux2 uy2 uz2 T
(19)
ux3 uy3 uz3 ux4 uy4 uz4 dan N merupakan matriks berukuran (3x9). Dengan menggunakan matriks fungsi bentuk yang sama pada (15), ν dan ∆u didefinisikan sebagai
SNTTM X | 828
Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
ν = Nν nodal
(20)
∂W ( E n )
∆u = N∆u nodal
(21)
∂E
Adanya definisi koordinat lokal dan global membutuhkan persamaan matriks transformasi dalam kedua arah x(x, y,z) = J x(ξ, η, ϕ) dengan J adalah matriks transformasi. Pada umumnya matriks tersebut didefinisikan sebagai Jacobian berikut
x1 − x 4 J = y1 − y 4 z1 − z4
x 2 − x4 y2 − y4 z2 − z4
x3 − x4 y 3 − y4 z 3 − z4
W ( En + ε ) − W ( E n )
(30)
ε
dengan ε adalah parameter scalar perturbation. Selanjutnya turunan dari Persamaan (30) dapat diberikan sebagai
∂ 2W (En ) ∂ E ∂E
=
1 ∂W ( E n + ε , E n + ε ) ε ∂E
−
(22)
∂W ( E n , E n ) ∂E (31)
Berdasarkan fungsi bentuk pada (15), (20), dan (21) dapat diberikan definisi-definisi berikut
∂u ∂N u nodal = ∂x ∂x ∂υ ∂N = υnodal ∂x ∂x ∂∆u ∂N = ∆u nodal ∂x ∂x
=
(23) (24) (25)
dengan
Untuk menuliskan Persamaan (31) dalam bentuk algoritma pemrograman diperlukan indeks. Salah satu alternatif pemberian indeks adalah sebagai berikut
∂ ∂E kl
∂W ( E ) 1 ∂W ( E ij + ε, E kl + ε ) = ∂E ij ∂E ij ε ∂W ( E ij , E kl ) − ∂E ij (32)
∂N ∂N ∂ξ ∂N ∂η ∂N ∂ϕ = + + ∂x ∂ξ ∂x ∂η ∂x ∂ϕ ∂x
(26)
Pada akhirnya, menggunakan persamaan (23) hingga (25), persamaan (14) dapat di tulis kembali menjadi
K ∆u nodal = f
(27)
(29)
5. EVALUASI MODULUS TANGEN SECARA NUMERIK Pada umumnya energi deformasi pada Persamaan (8) adalah fungsi nonlinear dari regangan E . Fakta ini seringkali menyulitkan ketika harus mengevaluasi turunan kedua untuk memperoleh matriks modulus t
t
mendapatkan komponen C11 . Dari proses simetrisasi dan pengurangan komponen matriks pada Persamaan t
dengan ∂N T * ∂N ∂N T ∂N dV + ∫ dV K =∫ TM FMCt FM* ∂ ∂ ∂ ∂x x x x V V (28)
∂N T f = −∫ FM TdV + ν T Pfext ∂ x V
dimana indeks i, j, k, dan l menunjukan lokasi baris dan kolom dari matriks Green-Lagrange E . Sebagai ilustrasi, akan diberikan contoh untuk
kekakuan C . Untuk mengatasi kesulitan tersebut, pada LD-FEM diaplikasikan metoda Gill-Murray (Gill et al, 1980) untuk menghitung turunan numerik kedua dari W ( E ) terhadap E . Pada satu dimensi, turunan numerik dengan forward difference didefinisikan sebagai
t
(8), C11 adalah komponen matriks berindeks C1111 . Untuk memperoleh harga komponen tersebut Persamaan (32) akan menjadi
∂ ∂W ( E ) 1 ∂W ( E11 + 2ε ) ∂W ( E11 ) − = . ∂ E11 ∂ E11 ε ∂E11 ∂ E11 (33) Ilustrasi lain juga dapat diberikan untuk komponen
C t45 yang bersesuaian dengan C t2313 yaitu ∂ ∂W ( E ) 1 ∂W ( E 23 + ε, E13 + ε ) = ∂E13 ∂E 23 ε ∂E 23 −
∂W ( E 23 , E13 ) ∂E 23 (34)
6. CONTOH SIMULASI Sebagi gambaran kinerja program LD-FEM, akan disajikan contoh simulasi pengujian tekan karet. Model specimen yang digunakan adalah silinder karet dengan kedua sisi direkatkan kepada dua buah pelat.
SNTTM X | 829
Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
Pembebanan dilakukan pada salah satu sisi dengan sisi yang lain tetap, tidak berotasi dan bertranslasi. Lebih jelas mengenai ilustrasi tumpuan dan modal dapat dilihat pada Gambar 3. Untuk memodelkan karet, LD-FEM menggunakan fungsi energi regangan Saint Venant-Kirchhoff. Fungsi ini diberikan sebagai
W (E) =
2 λ tr ( E ) + 2 µ tr ( E2 ) 2
∆ u T ∆ u ≤ 10 − 12 . Dalam hal ini penggunaan besaran 10-12 adalah kompromi antara ketelitian perhitungan dan keterbatasan presisi komputer.
(35)
F
13 mm
φ = 28 mm
Gambar 3. Model pembebanan dan tumpuan dengan F = 97.5 N
Gambar 5. Hubungan nonlinear antara beban dan perpindahan pada lokasi nodal nomor 87 pada Gambar 4. Sebagai pembanding adalah solusi 1D secara analitik dan MSC.Nastran.
Gambar 6. Ilustrasi tampilan Gmsh untuk deformasi silinder hasil perhitungan LD-FEM pada tahapan beban ke-8 Gambar 4. Solid mesh yang diperoleh dari Gmsh. Gaya bekerja pada permukaan atas dalam arah sumbu Z negatif.
dengan λ dan µ adalah konstanta Lame. Masingmasing konstanta diberikan sebagai fungsi poisson ratio ν dan modulus elastisitas E berikut Eν λ= (36) 1+ ν ( )(1 − 2ν )
µ=
E 2 (1+ν ) .
(37)
Pada contoh simulasi ini material karet diasumsikan memiliki harga E = 1 MPa dan ν = 0.3 . Sementara itu, kriteria konvergensi pada proses iterasi Newton-Raphson adalah norm
Gambar 5 menunjukan hubungan yang tidak linear antara beban dan perpindahan. Perlu ditekankan disini bahwa ketidaklinearan tersebut disumbang oleh perilaku regangan Green-Lagrange. Perilaku tersebut dapat sesuai dengan persamaan konstitutif yang merupakan turunan pertama fungsi kerapatan energi
T ( E ) = λ tr ( E ) I + 2 µE .
(38) Untuk memvalidasi harga deformasi keluaran LDFEM, solusi analitik 1-D dan MSC.Nastran disajikan pula di Gambar 5 sebagai perbandingan. Solusi dari MSC.Nastran menggambarkan spesimen lebih fleksibel dibandingkan dari LD-FEM. Hal ini mudah dimengerti karena formulasi elemen tetrahedral pada MSC.Nastran menggunakan fungsi tambahan internal yang dikenal sebagai bubble function. Sementara itu solusi analtik 1-D mengasumsikan kondisi tegangan 1D tanpa melibatkan pengaruh poisson ratio sehingga flesksibilitas masih diatas hasil MSC.Nastran. Perlu
SNTTM X | 830
Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
diperhatikan disini bahwa MSC.Nastran memberikan grafik linear karena solusi yang digunakan adalah formula defomasi kecil. Ini dikarenakan modul simulasi deformasi besar pada MSC.Nastran tidak menggunakan fungsi regangan Saint Venant-Kirchoff. Dilain pihak kinerja algortima Newton-Raphson pada LD-FEM terbukti cukup meyakinkan. Hal ini ditandai dengan trend penurunan nilai norm
∆u T ∆u secara kuadratik (lihat Gambar 7). Hanya pada iterasi terakhir tidak tercapai konvergensi akibat grafik mulai melandai yang menandai modulus tangen residual mendekati 0 (nol). Gambar 8. Perbandingan jumlah iterasi pada tahapan beban ke-9 antara solusi turunan ke-2 secara analitik fungsi Saint Venant-Kirchhoff dengan solusi numerik untuk harga ε yang berbeda. Target perpindahan pada nodal 87 sebesar 1.84 mm tercapai untuk kedua harga ε .
Gambar
7.
Penurunan
nilai
norm
∆u T ∆u berlangsung secara kudratik hingga tahapan beban ke-7. Mulai tahapan beban ke-9 modulus tangen mulai melandai dan harganya mendekati nol. Akibatnya laju konvergensi menjadi lambat bahkan pada thapan ke-10 tidak konvergen.
Permasalahan pada tahapan ke-10 merupakan kelemahan metoda Newton-Raphson dengan kendali beban (load control). Metoda ini mengalami kesulitan menemukan solusi karena melibatkan pembagian dengan harga ∂r ∂ ( ∂u ∂x ) → 0 (lihat Pers.10). Salah satu alternatif untuk mengatasi hal ini adalah modifikasi Newton-Raphson dengan kendali deformasi . Sementar itu, kinerja turunan numeris kedua pada LD-FEM cukup menjanjikan. Hal ini terlihat pada Gambar 8 yang menunjukan perbandingan laju norm terhadap jumlah iterasi pada tahapan beban ke-9. Pengujian dengan tahapan beban ke-9 dilakukan karena tahapan ini merupakan batas sensitif konvergensi (lihat n
n
Gambar 7). Dalam rentang harga ε antara 10000 hingga 0.001 laju norm berimpit dengan turunan analitik dan menghasilkan perpindahan pada nodal 87 sebesar 1.84 mm. Sebagai sandaran teori, dapat dilihat pada GillMurray (Gill et al, 1980). Pada artikel tersebut dinyatakan bahwa pilihan ε yang terbaik dapat ditentukan berdasarkan formula ε≈
3
Ef W ( E)
(39)
K
dengan E f adalah kesalahan dari floating point komputer dan
K =
∂ 3 W (E ) ∂ E ij 2 ∂ E k l
+
∂ 3 W (E ) ∂ E ij ∂ E k l 2
.
(40)
Menggunakan Persamaan (40) dapat dihitung bahwa harga K = 0 pada kasus fungsi Saint Venant-Kirchhoff pada (35). Untuk itu, berdasarkan formula (39) pilihan ε yang terbaik adalah ε → ∞ . 7. KESIMPULAN Berdasarkan hasil simulasi dapat disimpulkan bahwa kinerja LD-FEM cukup memuaskan dilihat dari dua aspek. Aspek pertama adalah LD-FEM mampu mewakili respon nolinear pada kasus deformasi besar (large deformation), dengan prediksi harga deformasi tidak terlampau jauh dengan prediksi perhitungan analitik 1-D dan MSC.Nastran. Kedua adalah laju konvergensi LD-FEM terbukti memiliki karakteristik kuadratik yang merupakan ciri khas Metoda Newton-Raphson. Selain itu juga dapat disimpulkan bahwa implementasi metoda Gill-Murray pada LD-FEM untuk kasus fungsi energi regangan
SNTTM X | 831
Seminar Nasional Teknik Mesin X Jurusan Mesin Fakultas Teknik UB
2-3 November 2011 ISBN 978 – 602 – 19028 – 0 – 6
Saint Venant –Kirchhoff memperoleh hasil yang sesuai dengan prediksi analitik. 8. DAFTAR PUSTAKA [1] Bower, A. F. (2008). Applied Mechanics of Solids at http://solidmechanics.org diakses Janurai 2010 [2] Gill, P. E., Murray, W., Saunders, M. A. and Wright, M. H. (1980). Computing the Finite Difference Approximations to Derivatives for Numerical Optimization. Technical Report, Department of Operations Research. Stanford University [3] Geuzaine, C., Remacle, JF. (2009). Gmsh: a three dimensional finite element generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 0:1-24 [4] MSC. Software Corporation Inc. (1999). MSC. Nastran Documentation. MSC. Nastran for Windows 4.5 [5] Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (1986). Numerical Recipies in Fortran 77: The Art of Scientific Computing. Cambridge University Press, 2nd Edition [6] Strauss, W. (2008). Partial Differential Equations: An Introduction. Wiley [7] Sussman T., Bathe, K. J. (1987). A Finite Element Formulation for Nonlinear Incompressible Elastic and Inelastic Analysis. Journal Computers & Structures, Vol. 26, pp.357-409 [7] Taylor, R. L. (1999). A Mixed-Enhanced Formulation for Tetrahedral Finite Elements. Report No. UCB/SEMM-99/02. Berkeley: University of California [8] Taylor, R. L. (2000). FEAP – A Finite Element Analysis Program – Version 7.3. Berkeley: University of Califo
SNTTM X | 832
