Metode Elemen Batas (MEB) untuk Model Konduksi-Konveksi dalam Media Anisotropik ∗ Moh. Ivan Azis† September 13, 2011
Daftar Isi 1 Pendahuluan
1
2 Masalah nilai batas
1
3 Persamaan integral batas
2
4 Hasil numerik
3
5 Konklusi
5
1
Pendahuluan
Berbagai proses fisik dapat dimodel dengan persamaan konduksi-konveksi, seperti transpor dan dispersi polutan dalam tanah dan air permukaan, difusi dalam alat bermaterial semikonduktor, medan magnet bergerak, magnetohydrodynamics, dan lain-lainnya. Dalam bidang perpindahan panas, kita dapat menyebutkan masalah pertumbuhan kristal, pengerasan permukaan benda menggunakan laser, pemotongan logam, dan lain-lainnya (Wrobel dan DeFigueiredo [2]). Model konduksi-konveksi untuk media anisotropik sangat relevan dengan proses fisik transpor dan dispersi polutan dalam tanah, karena medium tanah dapat dianggap sebagai medium anisotropik.
2
Masalah nilai batas
Dengan merujuk pada sistim kordinat Cartesian Ox1 x2 , secara umum persamaan pembangun untuk sistem konduksi-konveksi steady dua dimensi dalam suatu media anisotropik yang homogen, dengan asumsi bahwa tidak terdapat sumber pembangkit dalam media, adalah λij
∂2ϕ ∂ϕ − vi =0 ∂xi ∂xj ∂xi
(1)
dimana ϕ dapat berupa temperatur, λij adalah diffusivitas panas konstan, dan vi adalah komponen vektor kecepatan konstan v. Matriks kofisien [λij ] merupakan matriks bilangan real definit positif dan simetris. Juga, pada (1) penjumlahan untuk index yang berulang (jumlahan dari 1 sampai 2) diberlakukan. Suku pertama di ruas kiri dari (1) merepresentasikan proses konduksi, sedangkan suku keduanya menggambarkan proses konveksi dari sistem. ∗ Disampaikan pada Pelatihan Pemodelan Matematika dan Simulasi beberapa Sistem Fisik untuk Media Anisotropik Menggunakan Metode Elemen Batas, Fak. MIPA Unhas, Pebruari 2002 † Dosen Jurusan Matematika Fak. MIPA Unhas, Makassar Indonesia. Email:
[email protected]
1
Solusi dari (1) dicari dimana solusi tersebut valid dalam daerah Ω di R2 dengan batas Γ yang terdiri dari sejumlah berhingga kurva halus bagian demi bagian. Pada Γ salah satu dari ϕ(x) atau P (x) = λij
∂ϕ(x) nj (x) ∂xi
diberikan, dimana x = (x1 , x2 ), n = (n1 , n2 ) melambangkan vektor normal satuan mengarah ke luar di batas Γ. Metode solusi yang dipakai akan bekerja dengan cara menurunkan suatu persamaan integral batas, yang relevan untuk persamaan differensial (1), darimana nilai numerik ϕ dan P dapat ditentukan untuk semua titik dalam daerah Ω.
3
Persamaan integral batas
Bila (1) diperkalikan dengan ϕ∗ lalu diintegralkan, maka ) ∫ ( ∂ϕ ∂2ϕ − vi λij ϕ∗ dΩ = 0 ∂x ∂x ∂x i j i Ω Dengan menggunakan Teorema Divergensi Gauss pada (2) kita akan memperoleh ) ) ∫ ( ∫ ( ∂ϕ ∂ϕ ∂ϕ∗ ∂ϕ∗ λij nj − ϕvi ni ϕ∗ dΓ − λij − ϕ vi dΩ = 0 ∂xi ∂xi ∂xj ∂xi Γ Ω
(2)
(3)
Penggunaan Teorema Divergensi Gauss sekali lagi pada integral domain di persamaan (3) untuk integran pertamanya akan menghasilkan ) ∫ ( ∫ ∂ϕ ∂ϕ∗ λij nj − ϕvi ni ϕ∗ dΓ − ϕλij nj dΓ ∂xi ∂xi Γ Γ ) ∫ ( ∂ϕ∗ ∂ 2 ϕ∗ + ϕ vi dΩ = 0 (4) + ϕλij ∂xi ∂xj ∂xi Ω Atau
) ∫ ( ∫ ∂ 2 ϕ∗ ∂ϕ∗ λij + vi ϕ dΩ = − [P ϕ∗ − (Pv ϕ∗ + P ∗ )ϕ] dΓ ∂xi ∂xj ∂xi Ω Γ
dimana Pv (x) = vi ni (x)
dan
P ∗ (x, ξ) ≡ λij
(5)
∂ϕ∗ (x, ξ) nj (x) ∂xi
Fungsi ϕ∗ , yang biasa disebut sebagai solusi fundamental, kita definisikan sebagai λij
∂ 2 ϕ∗ ∂ϕ∗ + vi = −δ(x − ξ) ∂xi ∂xj ∂xi
(6)
dimana ξ = (ξ1 , ξ2 ), dan δ adalah fungsi delta Dirac. Solusi fundamental ϕ∗ ini dapat dituliskan sebagai berikut (lihat Azis [1] untuk penurunan ϕ∗ ) ( ) ( ) ˙ ˙ R v. v˙ R˙ K ∗ exp − K0 (7) ϕ (x, ξ) = 2π 2D 2D ˙ x˙ = (x1 + ρx ˙ ˙ = x− dimana D = [λ11 +λ12 (ρ+ρ)+λ ¨/D, R ˙ ξ, ˙ 2 , ρ¨ξ2 ), 22 ρρ]/2, K = ρ √ √˙ 2 , ρ¨x2 ), ξ = (ξ1 + ρξ 2 2 ˙ v˙ = (v1 + ρv ˙ 2 , ρ¨v2 ), R = (x1 + ρx ˙ 2 − ξ1 − ρξ ˙ 2 ) + (¨ ρx2 − ρ¨ξ2 ) , v˙ = (v1 + ρv ˙ 2 )2 + (¨ ρv2 )2 , ρ˙ dan ρ¨ berturut-turut merupakan bagian real dan imajiner positif dari akar kompleks ρ dari persamaan kuadrat λ11 + 2λ12 ρ + λ22 ρ2 = 0 dan tanda bar (¯) melambangkan operasi konjugat untuk bilangan kompleks, serta K0 adalah fungsi Bessel termodifikasi berorde nol. 2
Sebagai salah satu sifat dari fungsi delta Dirac, persamaan berikut berlaku ∫ ϕ(x) δ(x − ξ) dΩ(x) = η(ξ) ϕ(ξ)
(8)
Ω
dimana η = 12 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 (6) ke dalam ruas kiri dari (5) dan penggunaan persamaan (8) memberikan ∫ η(ξ) ϕ(ξ) = {P (x) ϕ∗ (x, ξ) Γ
− [Pv (x) ϕ∗ (x, ξ) + P ∗ (x, ξ)] ϕ(x)} dΓ(x)
(9)
Persamaan (9) dapat digunakan untuk menentukan solusi ϕ 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 (9). Tetapi secara umum integral batas ini tidak mudah dikalkulasi secara analitik, karena bentuk geometri dari Γ tidak beraturan atau kelakuan dari fungsi ϕ dan P sangat bervariasi. Untuk itu, nilai integral batas ini lalu diapproksimasi dengan cara memenggalmenggal batas domain Γ menjadi segmen-segmen kecil berupa garis lurus dan kelakuan dari fungsi ϕ dan P juga didekati dengan mengasumsikan bahwa fungsi-fungsi ini konstan, atau bervariasi secara linear, kuadratik dan seterusnya pada setiap segmen. 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.
4
Hasil numerik
Pada pasal ini beberapa contoh masalah konduksi-konveksi akan diselesaikan secara numerik dengan menggunakan persamaan integral batas (9). Contoh 1 : Masalah Uji Perhatikan solusi analitik untuk (1) berikut ϕ = exp(α1 x1 + α2 x2 )
(10)
dimana α1 dan α2 adalah bilangan ril yang memenuhi λ22 α22 + (2λ12 α1 − v2 )α2 + (λ11 α12 − v1 α1 ) = 0.
(11)
Geometri medium dan syarat batas dari masalahnya adalah (lihat Gambar 1) P, yang dapat dihitung dari (10), diketahui pada AB, BC dan CD, ϕ, seperti diberikan oleh (10), diketahui pada AD. Kofisien λij dan vi adalah λ11 = 1, λ12 = 1, v1 = 1, v2 = 1.
λ22 = 2,
Juga diambil α1 = 1 dan nilai α2 dihitung dari persamaan kuadrat (11). 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. Contoh 2 Perhatikan masalah transpor dan dispersi polutan dalam tanah yang dibangun oleh persamaan (1) untuk suatu medium tanah yang diasumsikan isotropik homogen dan memiliki geometri dan syarat batas seperti diperlihatkan dalam Gambar 2 dengan kofisien konduktivitas λ11 = 1, λ12 = 0, λ22 = 1, dan kecepatan aliran polutan v1 = 0, v2 = 1 untuk kasus pertama, dan v1 = 2, v2 = 1 untuk kasus 3
x2 6 D(0, 1)
C(1, 1)
- x1 A(0, 0)
B(1, 0)
Gambar 1: Geometri dari masalah uji
Tabel 1: Solusi MEB dan analitik untuk Contoh 1 Posisi (x1 , x2 )
ϕ
MEB ∂ϕ/∂x1
(.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)
1.1039 1.3473 1.6453 2.0097 2.4546
1.1004 1.3436 1.6456 2.0096 2.4552
(.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)
1.1045 1.3484 1.6468 2.0116 2.4569
1.1023 1.3462 1.6473 2.0117 2.4567
(.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5)
1.1048 1.3490 1.6476 2.0125 2.4580
1.1034 1.3475 1.6479 2.0126 2.4576
∂ϕ/∂x2 80 segmen .0011 .0040 .0024 .0020 .0020 160 segmen .0007 .0025 .0013 .0014 .0017 320 segmen .0005 .0016 .0008 .0009 .0013
4
ϕ
Analitik ∂ϕ/∂x1 ∂ϕ/∂x2
1.1052 1.3499 1.6487 2.0138 2.4596
1.1052 1.3499 1.6487 2.0138 2.4596
.0000 .0000 .0000 .0000 .0000
1.1052 1.3499 1.6487 2.0138 2.4596
1.1052 1.3499 1.6487 2.0138 2.4596
.0000 .0000 .0000 .0000 .0000
1.1052 1.3499 1.6487 2.0138 2.4596
1.1052 1.3499 1.6487 2.0138 2.4596
.0000 .0000 .0000 .0000 .0000
y 6 D(0, 1)
P =1
C(1, 1)
P =0
P =0
- x A(0, 0)
ϕ=0
B(1, 0)
Gambar 2: Geometri untuk Contoh 2
Tabel 2: Solusi MEB untuk Contoh 2 Posisi (x1 , x2 ) (.1,.5) (.3,.5) (.5,.5) (.7,.5) (.9,.5) (.5,.1) (.5,.3) (.5,.5) (.5,.7) (.5,.9)
ϕ .2384 .2384 .2384 .2384 .2384 .0385 .1285 .2384 .3726 .5366
v1 = 0, v2 = 1 ∂ϕ/∂x1 ∂ϕ/∂x2 .0000 .6064 .0000 .6063 .0000 .6063 .0000 .6063 .0000 .6064 .0000 .4064 .0000 .4964 .0000 .6063 .0000 .7406 .0000 .9046
ϕ .2380 .2380 .2381 .2382 .2382 .0385 .1283 .2381 .3722 .5361
v1 = 2, v2 = 1 ∂ϕ/∂x1 ∂ϕ/∂x2 .0003 .6054 .0003 .6056 .0003 .6057 .0002 .6059 .0001 .6061 .0000 .4059 .0002 .4958 .0003 .6057 .0005 .7401 .0006 .9044
kedua. Tidak tersedia solusi analitik eksak sederhana untuk contoh masalah ini. Kita tertarik untuk melihat pengaruh perubahan komponen kecepatan v1 , dari v1 = 0 menjadi v1 = 2. Tabel 2 memperlihatkan solusi MEB dengan menggunakan 320 segmen. Dari tabel ini dapat diamati pengaruh perubahan komponen kecepatan v1 (kecepatan ke arah sumbu-x1 ) dari v1 = 0 menjadi v1 = 2, dengan komponen kecepatan v2 tetap (v2 = 1). Dapat dikatakan bahwa perubahan ini menurunkan besarnya nilai (magnitude) dari solusi ϕ dan ∂ϕ/∂x2 dan menaikkan besarnya nilai ∂ϕ/∂x1 . Secara intuitif hasil ini diharapkan (expected) karena keberadaan aliran ke arah sumbu-x1 (v1 ̸= 0) akan mempengaruhi (mengurangi) laju aliran ke arah sumbu-x2 (∂ϕ/∂x2 ) dan konsentrasi polutan ϕ itu sendiri. Dapat pula diamati bahwa perubahan nilai v1 dari v1 = 0 menjadi v1 = 2 mempengaruhi kesimetrian solusi di sepanjang ordinat x2 = 0.5. Ketika v1 = 0 solusi di titik-titik dengan ordinat x2 = 0.5 simetris, tapi untuk v1 = 2 kesimetrian ini tidak terjadi lagi. Secara intuitif hal ini juga diharapkan.
5
Konklusi
Suatu MEB untuk solusi masalah nilai batas untuk model konduksi-konveksi 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
5
akurat. Evaluasi integral secara analitik, penerapan proses refinement untuk penyelesaian sistim persamaan aljabar linear, dan strategi peletakan titik sumber di luar domain akan memberikan hasil yang lebih akurat.
References [1] Azis, M. I. (2001). On the boundary integral equation method for the solution of some problems for inhomogeneous media (PhD Thesis), Department of Applied Mathematics, University of Adelaide. [2] Wrobel, L. C. and DeFigueiredo, D. B. Coupled Conduction-Convection Problems. in BEMs in Heat Transfer, Wrobel, L. C. and Brebbia, C. A. (eds.), Computational Mechanics Publications, Elsevier Applied Science
6
