Metode Elemen Batas (MEB) untuk Model Konduksi Panas

Metode Elemen Batas (MEB) untuk Model Konduksi Panas ∗ Moh. Ivan Azis† October 14, 2011

Abstrak Metode Elemen Batas untuk masalah konduksi panas pada media ortotropik berhasil ditemukan pada tulisan ini. Solusi numerik yang diperoleh memperlihatkan kecocokan dan keakuratan pemakaian metode elemen batas untuk masalah ini.

Daftar Isi 1 Persamaan pembangun

1

2 Masalah nilai batas

3

3 Solusi fundamental

4

4 Persamaan integral batas

5

5 Diskritisasi

5

6 Hasil numerik

8

7 Konklusi

9

1

Persamaan pembangun

Di pasal ini akan ditunjukkan penurunan persamaan pembangun untuk model konduksi panas dalam media anisotropik. Lihat Holman [4] untuk penjelasan lebih detail. Ketika gradien suhu terjadi dalam suatu medium, maka pengalaman menunjukkan bahwa akan terdapat perpindahan energi dari daerah dengan suhu tinggi ke daerah dengan suhu rendah. Kita katakan bahwa energi dipindahkan secara konduksi, dan ∗

Dipublikasikan pada Jurnal Fisika FUSI Jurusan Fisika FMIPA Unhas, Vol. 6, Nomor 8, 2002, Halaman 35 - 43 † Jurusan Matematika Universitas Hasanuddin, Indonesia. mailto:[email protected]

1

bahwa laju perpindahan panas (heat-transfer rate) per satuan area sebanding dengan gradien suhu dalam arah normal: q ∂T ∼ A ∂x Bila konstanta kesebandingan dimasukkan, q = −kA

∂T ∂x

(1)

dimana q adalah laju perpindahan panas dan ∂T /∂x gradien suhu dalam arah aliran panas. Konstanta positif k disebut konduktivitas panas dari material, dan tanda minus “-” disisipkan sehingga prinsip kedua dari thermodynamics terpenuhi; yakni bahwa panas pasti berpindah dalam arah menurun dalam skala ukuran suhu. Persamaan (1) disebut hukum Fourier untuk konduksi panas. Dan dalam (1) k memiliki satuan watt per meter per derajat Celcius untuk suatu sistem satuan tertentu dimana q bersatuan watt. Bila suhu dapat berubah dengan waktu dan terdapat suatu sumber panas dalam medium, maka untuk elemen dengan ketebalan dx kesetimbangan energi berikut dapat diperlakukan: “Energi masuk + panas yang dibangkitkan dalam elemen = perubahan energi dalam + energy keluar” Kwantitas energi-energi ini diberikan sebagai berikut: Energi masuk

=

qx = −kA ∂T ∂x

Energi dibangkitkan dalam elemen

=

qA ˙ dx

Perubahan energi dalam

=

ρcA ∂T dx ∂τ

Energi keluar

= qx+dx = −kA  + = −A k ∂T ∂x

∂ ∂x

∂T ∂x x+dx




 k ∂T dx ∂x

dimana q˙ = energi dibangkitkan per satuan volume (W/m3 ), c = panas spesifik dari material (J/kg.◦ C), ρ = kepadatan (kg/m3 ). Penggabungan relasi-relasi di atas memberikan     ∂T ∂T ∂ ∂T ∂T + qAdx ˙ = ρcA dx − A k + k dx −kA ∂x ∂τ ∂x ∂x ∂x atau

∂ ∂x

  ∂T ∂T k + q˙ = ρc ∂x ∂τ

(2)

Persamaan ini adalah persamaan konduksi panas satu-dimensi. Untuk kasus tiga-dimensi, persamaan kesetimbangan energi dapat dituliskan sebagai ∂E qx + qy + qz + qgen = qx+dx + qy+dy + qz+dz + ∂τ 2

dan kuantitas-kuantitas energi diberikan sebagai ∂T qx = −k dy dz ∂x     ∂T ∂ ∂T qx+dx = − k + k dx dy dz ∂x ∂x ∂x ∂T qy = −k dx dz ∂y     ∂ ∂T ∂T + k dy dx dz qy+dy = − k ∂y ∂y ∂y ∂T qz = −k dx dy ∂z     ∂T ∂ ∂T qz+dz = − k + k dz dx dy ∂z ∂z ∂z qgen = qdx ˙ dy dz ∂E ∂T = ρc dx dy dz ∂τ ∂τ sehingga bentuk umum dari persamaan konduksi panas tiga-dimensi adalah       ∂ ∂T ∂ ∂T ∂ ∂T ∂T k + k + k + q˙ = ρc ∂x ∂x ∂y ∂y ∂z ∂z ∂τ

(3)

Perhatikan bahwa pada penurunan dari persamaan di atas, ekspresi derivatif pada x + dx telah ditulis dalam bentuk ekspansi deret Taylor dengan pengambilan dua suku pertama saja. Untuk media yang anisotropik nilai konduktivitas ke arah x, y dan z tidak semuanya sama. Namakan konduktivitas ke arah x, y dan z masing-masing sebagai kx , ky dan kz . Sehingga (3) dapat ditulis sebagai       ∂T ∂ ∂T ∂ ∂T ∂T ∂ (4) kx + ky + kz + q˙ = ρc ∂x ∂x ∂y ∂y ∂z ∂z ∂τ Pembicaraan selanjutnya dibatasi hanya untuk kasus steady (suhu tidak berubah dengan waktu τ , sehingga ∂T /∂τ = 0) dan dengan asumsi bahwa di dalam medium tak terdapat suatu sumber atau pembangkit panas, (q˙ = 0). Lebih jauh diasumsikan bahwa medianya homogen (konduktivitas panas konstan), serta kita hanya akan memandang kasus untuk dua-dimensi saja. Sehingga persamaan yang akan kita hadapi adalah ∂ 2T ∂ 2T kx 2 + ky 2 = 0 ∂x ∂y

2

(5)

Masalah nilai batas

Dengan merujuk pada sistim kordinat Cartesian Oxy solusi dari (5) dicari dimana solusi tersebut valid dalam daerah Ω di R2 dengan batas Γ yang terdiri dari sejumlah 3

berhingga kurva muls bagian demi bagian. Pada Γ salah satu dari suhu T (x, y) atau flux ∂T ∂T ∂T P (x, y) = = kx nx + ky ny ∂n ∂x ∂y diberikan, dimana n = (nx , ny ) melambangkan vektor normal satuan mengarah ke luar di batas Γ. Metode solusi yang dipakai akan bekerja dengan cara menurunkan suatu persamaan integral batas yang relevan dengan (5), darimana nilai numerik T dan P dapat ditentukan untuk semua titik dalam daerah Ω.

3

Solusi fundamental

Persamaan integral batas yang disebutkan pada Pasal 2 melibatkan suatu fungsi solusi fundamental T ∗ yang didefinisikan sebagai kx

∂ 2T ∗ ∂ 2T ∗ + k = δ(x − ξ) y ∂x2 ∂y 2

(6)

dimana x = (x, y), ξ = (a, b) dan δ adalah fungsi delta Dirac. Solusi fundamental T ∗ ini dapat dituliskan sebagai berikut (lihat Azis [2] untuk penurunan T ∗ ) T ∗ (x, ξ) =

K ln R 2π

(7)

dimana D = [kx + ky ρρ]/2 K = ρ¨/D p R˙ = (x + ρy ˙ − a − ρb) ˙ 2 + (¨ ρy − ρ¨b)2 ρ˙ dan ρ¨ berturut-turut merupakan bagian real dan imajiner positif dari akar kompleks ρ dari persamaan kuadrat kx + ky ρ2 = 0 dan tanda bar (¯.) melambangkan operasi konjugat untuk bilangan kompleks. ∗ ≡ Selain T ∗ kita juga memerlukan fungsi P ∗ , yang didefinisikan sebagai P ∗ = ∂T ∂n ∗ ∗ kx (∂T /∂x)nx +ky (∂T ∂y)ny , untuk evaluasi persamaan integral batas tersebut di atas. Fungsi P ∗ ini adalah   K 1 ∂R ∂R ∗ P (x, ξ) = kx nx + ky ny (8) 2π R ∂x ∂y Turunan ∂R/∂x dan ∂R/∂y yang diperlukan untuk evaluasi nilai P ∗ diberikan oleh ∂R 1 = (x + ρy ˙ − a − ρb) ˙ ∂x R ∂R 1 = [ρ(x ˙ + ρy ˙ − a − ρb) ˙ + ρ¨(¨ ρy − ρ¨b)] ∂y R Perlu dicatat bahwa P ∗ memiliki titik singular pada x = ξ. 4

4

Persamaan integral batas

Identitas Green kedua memberikan     Z  Z   ∂ 2T ∂T ∗ ∂ 2T ∗ ∂ 2T ∗ ∂ 2T ∗ ∂T ∗ T −T dΓ = T kx − T kx 2 + ky 2 dΩ + ky ∂n ∂n ∂x2 ∂y 2 ∂x ∂y Γ Ω (9) Sebagai salah satu sifat dari fungsi delta Dirac, persamaan berikut berlaku Z T (x) δ(x − ξ) dΩ(x) = η(ξ) T (ξ) (10) Ω 1 2

dimana η = bila ξ berada pada batas domain Γ dan Γ mempunyai kemiringan yang berubah secara kontinyu pada ξ, η = 1 bila ξ berada di dalam domain Ω, η = 0 bila ξ berada di luar domain Ω. Substitusi persamaan (5), (6), (10) ke dalam persamaan (9) akan diperoleh Z η(ξ) T (ξ) = [P ∗ (x, ξ) T (x) − P (x) T ∗ (x, ξ)] dΓ(x) (11) Γ

Persamaan (11) dapat digunakan untuk menentukan solusi T dan P di setiap titik x di batas Γ dan di dalam domain Ω. Dan kalkulasi solusi ini sepenuhnya hanya memerlukan kalkulasi integral batas pada ruas kanan persamaan (11). Tetapi secara umum integral batas ini tidak mudah dikalkulasi secara analitik, karena bentuk geometri dari Γ tidak beraturan atau kelakuan dari fungsi T dan P sangat bervariasi. Untuk itu, nilai integral batas ini lalu diapproksimasi dengan cara memenggal-menggal batas domain Γ menjadi segmen-segmen kecil berupa garis lurus dan kelakuan dari fungsi T dan P pada setiap segmen juga didekati dengan mengasumsikan bahwa fungsi-fungsi ini konstan, atau bervariasi secara linear, kuadratik dan seterusnya. Lalu integral dihitung untuk setiap segmen dan kemudian menjumlahkannya. Dengan kata lain batas domain Γ diapproksimasi oleh suatu poligon yang jumlah sisinya diambil sebanyak mungkin sehingga nilai pendekatan akurat dapat diperoleh.

5

Diskritisasi

Misalkan batas domain Γ didekati oleh suatu poligon dengan sejumlah J sisi, sehingga Γ terdiri atas segmen-segmen garis lurus Γj ≡ [qj−1 , qj ], j = 1, 2, . . . , J dimana qj−1 dan qj adalah titik-titik ujung awal dan akhir dari segmen Γj , maka persamaan (11) dapat ditulis sebagai J Z X η(ξ) T (ξ) = [P ∗ (x, ξ) T (x) − P (x) T ∗ (x, ξ)] dΓ(x) (12) j=1

Γj

Selanjutnya, bila kita mengasumsikan bahwa pada setiap segmen Γj nilai T dan P konstan, dan masing-masing diwakili oleh nilainya pada titik-tengah qj = (qj−1 +qj )/2 dari segmen tertentu Γj , maka persamaan (12) dapat ditulis sebagai " # Z qj Z qj J X η(ξ) T (ξ) = T (qj ) P ∗ (x, ξ) dΓ(x) − P (qj ) T ∗ (x, ξ) dΓ(x) (13) j=1

qj−1

qj−1

5

Hasil penelitian telah menunjukkan bahwa pengambilan nilai T dan P pada titiktengah qj untuk setiap segmen Γj menghasilkan keakuratan maksimal. Sebagaimana disebutkan pada Pasal 2, pada suatu segmen Γj hanya salah satu dari T dan P diketahui. Bila nilai T (qj ) diberikan maka nilai P (qj ) menjadi unknown di Γj . Sebaliknya, bila pada segmen Γj nilai P (qj ) diberikan maka nilai T (qj ) menjadi unknown. Untuk penentuan nilai unknown di batas domain Γ, hanya ada dua kemungkinan pengambilan posisi titik ξ, yakni diletakkan di batas domain Γ (yang mengimplikasikan bahwa η = 21 ) atau di luar domain Ω (mengimplikasikan η = 0). Kita tidak dapat meletakkan ξ di dalam domain Ω (untuk mana η = 1) untuk penentuan nilai unknown di batas domain Γ, kecuali bila kita mempunyai informasi tambahan mengenai nilai T di titik dalam ξ ini. Sementara itu, peletakan titik ξ di luar domain Ω akan menghindari titik singular dari P di x = ξ dan hal ini tentu akan memiliki advantage untuk hasil evaluasi integral. Dan telah ada beberapa kajian di dalam beberapa paper yang telah terpublish, yang memakai analisis peletakan titik ξ di luar domain Ω. Umumnya kajian-kajian ini telah berhasil menentukan jarak ideal dari titik ξ ke batas domain Γ untuk tingkat keakuratan yang cukup bagus. Namun penentuan jarak optimal ini masih sebatas cara coba-coba (trial and error), dan belum dilandasi oleh dan belum dibuktikan keabsahannya secara analitik matematik. Untuk itu, pada tulisan ini kita akan memposisikan titik ξ pada batas domain Γ. Sehingga untuk penentuan nilai unknown di batas domain Γ, nilai T (ξ) pada ruas kiri (13) akan mengambil nilai T (ql ) dan η(ξ) = 21 . Persamaan (13) kemudian dapat dituliskan sebagai " # Z qj Z qj J X 1 P ∗ (x, ql ) dΓ(x) − P (qj ) T ∗ (x, ql ) dΓ(x) T (ql ) = T (qj ) (14) 2 q q j−1 j−1 j=1 untuk l = 1, 2, . . . , J. Persamaan ini dapat dituliskan dalam bentuk matriks J

J

X X 1 b lj Tj = − Tl + H Glj Pj 2 j=1 j=1

(15)

dimana Tj = T (qj ), Pj = P (qj ), dan Z

qj

b lj = H q Z j−1 qj

Glj =

P ∗ (x, ql ) dΓ(x)

(16)

T ∗ (x, ql ) dΓ(x)

(17)

qj−1

Evaluasi integral pada persamaan (16) dan (17) dapat dilakukan secara analitik maupun numerik. Tentunya evaluasi analitik (eksak) akan memberikan hasil yang lebih memuaskan (akurat) daripada evaluasi numerik (pendekatan). Namun perlu diperhatikan bahwa untuk j = l selang integral dalam (16) memuat titik singular ql dari integran P ∗ (x, ql ). Untuk itu nilai prinsipal Cauchy (Cauchy principal value) dari integral ini biasanya diambil untuk evaluasi analitik. Di lain hal, dangan evaluasi 6

numerik dari kedua integral ini, strategi pemilihan metode kuadratur (pengintegralan numerik) sangatlah penting, sebab terdapat beberapa metode kuadratur yang melibatkan dan ada pula yang tidak melibatkan (misalnya aturan trapezoidal) kalkulasi nilai fungsi integran pada titik tengah dari selang integral. Dan metode kuadratur yang terakhir inilah yang dikehendaki. Pada tulisan ini, untuk hasil numerik dari setiap contoh masalah yang akan dibicarakan pada Pasal 6, evaluasi integral dilakukan secara numerik dengan menggunakan aturan trapezoidal termodifikasi enam titik (lihat Abramowitz and Stegun [1]). Lebih kompak, persamaan (15) dapat ditulis sebagai J X

Hlj Tj =

j=1

dimana

Glj Pj

(18)

j=1

( Hlj =

J X

b lj H b lj − H

1 2

bila l 6= j bila l = j

Persamaan matriks (18) dapat diurutkan ulang dengan meletakkan unknown di ruas kiri dan known-nya di ruas kanan, dalam bentuk AX = B

(19)

dimana X adalah vektor unknown T dan/atau P . Persamaan ini merupakan suatu sistem persamaan aljabar linear dengan J persamaan dan J unknown. Penyelesian sistem persamaan aljabar linear (19) dapat dilakukan dengan berbagai metode, antara lain dengan metode eliminasi Gauss. Namun, pada tulisan ini, untuk hasil numerik dari setiap contoh masalah yang akan dibicarakan pada Pasal 6, solusi sistem persamaan aljabar linear (19) ditentukan dengan menggunakan metode gradien konjugat (lihat Coleman [3]), yang secara empiris telah diketahui lebih stabil ketimbang metode eliminasi Gauss. Solusi dari persamaan (19) ini dapat ditentukan untuk unknown T dan P di batas domain Γ. Sekali nilai T dan P pada batas domain Γ telah diketahui, maka kita bisa menentukan nilai T dan P pada sebarang titik dalam ξ dengan menggunakan persamaan (13), yakni # " Z qj Z qj J X P ∗ (x, ξ) dΓ(x) − P (qj ) T ∗ (x, ξ) dΓ(x) (20) T (ξ) = T (qj ) j=1

qj−1

qj−1

Selain itu dapat pula ditentukan nilai turunan ∂T /∂a dan ∂T /∂b melalui persamaan berikut " # Z qj Z qj J X ∂T ∂P ∗ (x, ξ) ∂T ∗ (x, ξ) (ξ) = T (qj ) dΓ(x) − P (qj ) dΓ(x) (21) ∂a ∂a ∂a q q j−1 j−1 j=1 " # Z qj Z qj J X ∂T ∂P ∗ (x, ξ) ∂T ∗ (x, ξ) (ξ) = T (qj ) dΓ(x) − P (qj ) dΓ(x) (22) ∂b ∂b ∂b qj−1 qj−1 j=1 7

y 6

D(0, 1)

C(1, 1)

-

A(0, 0)

x

B(1, 0)

Gambar 1: Geometri dari masalah uji

6

Hasil numerik

Pada pasal ini beberapa contoh masalah konduksi panas dalam media anisotropik akan diselesaikan secara numerik dengan menggunakan persamaan integral batas yang telah diturunkan pada Pasal 4. Contoh 1 : Masalah Uji Perhatikan solusi analitik untuk (5) berikut T = kx x + ky y P = kx2 nx + ky2 ny Geometri medium dan syarat batas dari masalahnya adalah (lihat Gambar 1) P diketahui pada sisi AB, BC dan CD, T diketahui pada sisi AD. Kofisien konduktivitas kx = ky = 1. Tabel 1 memperlihatkan perbandingan antara solusi MEB dan solusi analitik. Dapat diamati bahwa solusi MEB konvergen ke solusi analitik sejalan dengan meningkatnya jumlah segmen dari 80, 160 dan 320. Hasil ini sesuai dengan yang diharapkan, dengan alasan bahwa semakin kecil selang integral yang digunakan maka semakin akurat pendekatan integrasi numerik yang akan diperoleh. Contoh 2 Perhatikan masalah konduksi panas yang dibangun oleh persamaan (5) untuk suatu medium isotropik (kx = ky = 1) homogen yang memiliki geometri seperti diperlihatkan dalam Gambar ??.

8

Table 1: Solusi MEB dan analitik untuk Contoh 1

7

Posisi (x, y)

T

(.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)

.5977 .7970 .9963 1.1957 1.3952

(.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)

.5989 .7986 .9983 1.1980 1.3978

(.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)

.5995 .7993 .9992 1.1991 1.3989

MEB Analitik ∂T /∂x ∂T /∂y T ∂T /∂x ∂T /∂y J = 80 segmen .9959 .9994 .6000 1.0000 1.0000 .9965 .9985 .8000 1.0000 1.0000 .9969 .9980 1.0000 1.0000 1.0000 .9972 .9976 1.2000 1.0000 1.0000 .9978 .9976 1.4000 1.0000 1.0000 J = 160 segmen .9982 .9997 .6000 1.0000 1.0000 .9985 .9993 .8000 1.0000 1.0000 .9986 .9991 1.0000 1.0000 1.0000 .9987 .9989 1.2000 1.0000 1.0000 .9989 .9989 1.4000 1.0000 1.0000 J = 320 segmen .9992 .9999 .6000 1.0000 1.0000 .9993 .9997 .8000 1.0000 1.0000 .9994 .9996 1.0000 1.0000 1.0000 .9994 .9995 1.2000 1.0000 1.0000 .9995 .9995 1.4000 1.0000 1.0000

Konklusi

Suatu MEB untuk solusi masalah nilai batas untuk model konduksi panas steady dalam suatu medium anisotropik telah ditemukan. MEB ini secara umum cukup mudah untuk diimplementasikan untuk memperoleh solusi numerik untuk masalah tertentu. Hasil numerik yang diperoleh dengan menggunakan MEB ini mengindikasikan bahwa MEB ini dapat menghasilkan solusi numerik yang cukup akurat. Evaluasi integral secara analitik, penerapan proses refinement untuk penyelesaian sistim persamaan aljabar linear, dan strategi peletakan titik ξ di luar domain Ω akan memberikan hasil yang lebih akurat.

References [1] Abramowitz, M. and Stegun, A. Handbook of Mathematical Functions, Dover, New York, 1970. [2] Azis, M. I. On the boundary integral equation method for the solution of some problems for inhomogeneous media (PhD Thesis). Department of Applied Mathematics, University of Adelaide, 2001. [3] Coleman, C. J. University of Wollonggong, Australia.

9

[4] Holman, J. P. Heat Transfer, 8th edition. Mc. Graw Hill pp. 2–7, 1999.

10