FOTON, Jurnal Fisika dan Pembelajarannya
Volume 9, Nomor 2, Agustus 2005
Pengkajian Konduksi Panas Tak Tunak 2D Berdasarkan Hasil Tinjauan Komputasi Numerik Hari Wisodo Jurusan Fisika UM, Jl. Surabaya 6, Malang 65145 Tlp. (0341) 552125 E-mail: wisodo
[email protected]
Intisari : Pengkajian masalah konduksi panas tak tunak 2D telah berhasil dilakukan berdasarkan hasil tinjauan komputasi numerik. Dilakukan normalisasi terhadap persamaan konduksi panas yang dapat memberikan beberapa keuntungan komputasi. Diskretisasi dilakukan dengan menggunakan metode beda hingga (finite difference) skema FTCS. Pemilihan nilai step waktu ∆t harus memenuhi syarat stabilitas yang terkait dengan penggunaan tersebut.
Kata kunci : konduksi, FTCS
1
PENDAHULUAN
saja belum mencukupi karena hasil nilai-nilai numerik yang diharapkan baru akan muncul setelah diproses di komputer. Karena itu diperlukan bahasa komputer atau piranti tertentu yang biasanya berupa perangkat lunak berbentuk paket siap pakai sehingga pemakai dapat berkomunikasi dengan komputer dan memerintahkannya untuk melakukan proses komputasi seperti yang diharapkan. Pemahaman tentang watak-watak alur langkah komputasi (algoritma) serta watak komputer itu sendiri tentunya sangat membantu dalam melakukan trik dan manipulasi untuk optimasi proses komputasi.
Fisika Komputasi (Computational Physics) merupakan bidang yang mengkaji masalah fisika berdasarkan hasil tinjauan komputasi numerik [1]. Fisika komputasi sekarang telah diterima secara luas sebagai bentuk pendekatan ketiga fisika, selain fisika teori dan fisika eksperimen, dalam mempelajari fisika. Penghargaan Nobel Kimia untuk Pope dalam pengembangan teori fungsional kerapatan (density functional theory) untuk molekul menunjukkan peranan fisika komputasi yang semakin mantap. Bidang fisika komputasi tidak hanya memerlukan pemahaman fisika tetapi juga metode numerik dan bahasa pemrograman. Fisikawan komputasi perlu memiliki keterampilan untuk mengubah persamaan matematis atau hukum fisika ke dalam bentuk diskret yang sesuai sedemikian sehingga informasi fisis mengenai perilaku sistem dan penyelesaiannya dapat diwakili oleh angka-angka atau nilai-nilai numerik. Metode pengubahan persamaan matematis ke bentuk diskret beserta berbagai metode penyelesaiannya ini biasa disebut metode numerik. Pemahaman mengenai metode numerik itu ISSN 1410-3273
Persamaan difusi dipilih mengingat persamaan ini sering muncul di dalam persoalan fisika. Persamaan difusi digunakan untuk mendeskripsikan beberapa fenomena fisika yaitu difusi molekuler, konduksi termal, viskositas [2], dan superkonduktivitas [3]. Tulisan ini menyajikan penentuan penyelesaian persamaan difusi pada masalah konduksi termal dua dimensi dengan metode beda hingga (finite difference) untuk sekema Forward Time Centere Space (FTCS). 67
c 2005 Jurusan Fisika UM
Pengkajian Konduksi Panas Tak Tunak 2D . . . 2
SISTEM FISIS
Ditinjau penghantar panas persegi dengan luas penampang L × L yang terletak di bidang xy dengan salah satu sudutnya terletak di pusat koordinat seperti disajikan pada Gambar 1. Penghantar ini dipandang sebagai tampang lintang (cross section) penghantar tiga dimensi yang panjangnya tak hingga dan homogen dalam arah panjangnya. Mula-mula penghantar bertemperatur T0 . Pada saat t = 0, sisi penghantar yang terletak di x = L dan di y = L secara bersamaan ditempelkan pada reservoir bertemperatur Tb dimana Tb > T0 . Sementara pada sisi penghantar yang terletak di x = 0 dan di y = 0 dijaga sedemikian rupa sehingga tidak ada aliran kalor. Selanjutnya akan dicari keadaan sistem pada saat t.
Gambar 1: Sistem fisis yang ditinjau.
Persamaan yang sesuai dengan keadaan sistem fisis tersebut adalah persamaan difusi yang berbentuk ∂T =D ∂t
∂2T ∂2T + ∂x2 ∂y 2
(1)
dengan D = κ/ρc adalah tetapan difusi, κ adalah konduktifitas termal penghantar, ρ adalah massa jenis penghantar, dan c adalah panas jenis penghantar. Persamaan diferensial parsial (1) termasuk bertipe persamaan diferensial parsial parabolik. Syarat batas yang harus dipenuhi agar didapatkan penyelesaian yang unik adalah berjenis Derichlet untuk T (L, y, t) = Tb , T (x, L, t) = Tb ,
(2) (3)
68
dan berjenis Neumann untuk ∂T (x, y, t) = 0, ∂x (0, y) ∂T (x, y, t) = 0. ∂y (x, 0) 3
(4) (5)
BENTUK TAK BERDIMENSI DAN DISKRETISASI
Persamaan (1) diselesaikan melalui pendekatan numerik dengan menggunakan pendekatan beda hingga untuk derivatifnya. Terlebih dahulu persamaan tersebut dinormalisir. Cara ini setidaknya dapat memberikan tiga keuntungan. Pertama, nilai yang terlibat dalam komputasi dapat dijamin tidak terlalu besar atau terlalu kecil. Kedua, persamaan yang terlibat menjadi berbentuk sederhana. Ketiga, dimungkinkannya diperoleh ketelitian proses komputasi yang tinggi mengingat angka numerik yang terlibat berorde besar sesuai batas ketelitian komputer. Dalam upaya untuk menormalisir persamaan (1), didefinisikan satuan universal bagi semua variabel yang terlibat, yaitu x y , y0 = , L L t = , L2 /D T − T0 = . Tb − T0
x0 =
(6)
t0
(7)
T0
(8)
Pada persamaan (6) tampak bahwa x dan y dinormalisir terhadap L yang berarti bahwa kedua variabel ini diukur dalam satuan universal L. Pemilihan ini akan menjamin x0 dan y 0 bernilai [0, 1]. Mengingat sistem yang ditinjau, temperatur penghantar T (x, y, t) hanya akan bernilai T0 ≤ T ≤ Tb . Penormalisiran temperatur T seperti yang diberikan pada persamaan (8) menjamin T 0 bernilai [0,1]. Penormalisiran x, y, dan T semacam ini mengharuskan variabel waktu t dinormalisir terhadap L2 /D yang juga menjadi satuan universal bagi variabel t. Selanjutnya jika persamaan (6), (7), dan (8) disubstitusikan ke persamaan (1) untuk setiap variabel yang sesuai, diperoleh
FOTON/Vol. 9 No. 2/Agustus 2005
69
Hari Wisodo
persamaan (1) yang ternormalisir, yaitu ∂T (x, y, t) ∂ 2 T (x, y, t) ∂ 2 T (x, y, t) = + (9) ∂t ∂x2 ∂y 2 dengan tanda aksen telah dihilangkan. Sedangkan syarat batas yang ternormalisir diperoleh dengan cara mensubstitusikan T yang diperoleh dari persamaan (8) ke persamaan (2) dan (3) yang menghasilkan T (1, y, t) = 1, T (x, 1, t) = 1,
(10) (11) Gambar 2: Sistem grid yang ditinjau.
dan mensubstitusikan x, y, T yang diperoleh dari persamaan (6) dan (8) ke persamaan (4) dan (5) yang menghasilkan ∂T (x, y, t) = 0, (12) ∂x (0, y) ∂T (x, y, t) = 0. (13) ∂y (x, 0) dengan tanda aksen juga telah dihilangkan. Dalam upaya memperoleh bentuk diskret persamaan (9), terlebih dahulu dilakukan diskretisasi terhadap variabel-variabel yang terlibat, yaitu variabel bebas x, y, t dan variabel terikat T (x, y, t). Grid komputasi yang ditinjau diperoleh dengan cara membagi sistem fisis pada Gambar 1 ke dalam Nx ×Ny kisi yang masing-masing kisi memiliki luas hx ×hy . Pembagian semacam ini menghasilkan grid komputasi yang uniform seperti ditunjukkan pada Gambar 2. Karena itu bentuk diskret variabel bebas x dan y adalah
suku derivatif terhadap waktu, digunakan pendekatan beda maju Euler yaitu [1] n+1 n Ti,j − Ti,j ∂T (x, y, t) = ∂t ∆t (xi , yj , tn ) + O(∆t) (17) yang hanya teliti sampai orde pertama dalam ∆t. Sedangkan derivatif terhadap ruang digunakan pendekatan beda terpusat yang hanya menggunakan kuantitas T pada langkah waktu n, yaitu [1] ∂ 2 T (x, y, t) = ∂x2 (xi , yj , tn ) n n n Ti−1,j − 2Ti,j + Ti+1,j + O(h2x ), (18) h2x ∂ 2 T (x, y, t) = ∂y 2 (xi , yj , tn ) n n n Ti,j−1 − 2Ti,j + Ti,j+1 + O(h2y ). (19) h2y Pengimplementasian ini menghasilkan
xi = (i − 1)hx , i = 1, 2, · · · , Nx + 1 yj = (j − 1)hy , j = 1, 2, · · · , Ny + 1
(14) (15)
n+1 n Ti,j − Ti,j ∆t
(16)
n n n Ti−1,j − 2Ti,j + Ti+1,j h2x n n n Ti,j−1 − 2Ti,j + Ti,j+1 + . h2y (20)
Sedangkan variabel terikat T (x, y, t) di titik (xi , yj ) pada saat tn dinyatakan oleh n T (xi , yj , tn ) ≡ Ti,j . Bentuk diskret persamaan (9) diperoleh dengan cara mengimplementasikan pendekatan beda hingga pada setiap derivatifnya. Untuk
Persamaan (20) disebut representasi FTCS. Persamaan ini dapat diungkapkan dalam bentuk ∆t n n+1 n n Ti,j = Ti,j + 2 Ti−1,j + Ti+1,j hx n n n + Ti,j+1 − 4Ti,j , (21) + Ti,j−1
dan bentuk diskret variabel bebas t adalah tn = n∆t, n = 0, . . . , tmax .
=
FOTON/Vol. 9 No. 2/Agustus 2005
Pengkajian Konduksi Panas Tak Tunak 2D . . . dengan telah dipilih hx = hy . Persamaan (21) digunakan hanya untuk menghitung nilain+1 nilai Ti,j yang terletak di interior, yaitu n+1 Ti,j : i = 2, . . . , Nx ; j = 2, . . . , Ny . n+1 Nilai-nilai Ti,j yang terletak di batas dihitung dengan rumus yang sesuai dengan syarat batas yang telah ditetapkan. Temperatur pada sisi penghantar yang terletak di x = 1 dan y = 1 dihitung berdasarkan rumus yang diperoleh dari syarat batas persamaan (10) dan (11), yaitu TNn+1 = 1, j = 1, . . . , Ny (22) x +1,j n+1 = 1, i = 1, . . . , Nx (23) Ti,N y +1
= 1, TNn+1 x +1,Ny +1
(24)
Temperatur pada sisi penghantar yang terletak di x = 0 dihitung berdasarkan rumus yang diperoleh dari persamaan (9) dengan memasukkan syarat batas persamaan (12). Berdasarkan Gambar 2, ditinjau T2,j pada saat n t = tn , T2,j . Selanjutnya dengan menggunan kan ekspansi Taylor, T2,j diekspansikan di sekitar x = x1 , diperoleh ∂T (x, y, t) n n T2,j = T1,j + hx ∂x (x1 , yj , tn ) 2 2 h ∂ T (x, y, t) + x 2! ∂x2 (x1 , yj , tn ) 3 + O(hx ) (25) dimana hx = x2 − x1 . Berdasarkan persamaan (12) diketahui suku kedua ruas kanan persamaan (25) bernilai nol. Sehingga persamaan ini dapat dituliskan menjadi ∂ 2 T (x, y, t) 2 n n = T − T 2,j 1,j ∂x2 h2x (x1 , yj , tn ) + O(hx ). (26) Jika di dalam persamaan (21) i = 1, suku ke dua ruas kanan diganti dengan persamaan (26), dan hx = hy , maka dapat diperoleh n+1 n T1,j = T1,j +
+
2∆t n n T2,j − T1,j 2 hx
∆t n n n T − 2T1,j + T1,j+1 (27) . h2x 1,j−1
70
Persamaan inilah yang selanjutnya digunakan untuk menghitung nilai-nilai temperatur n+1 yang terletak di x = 0, yaitu T1,j : j = 2, . . . , Ny . Temperatur pada sisi penghantar yang terletak di y = 0 dihitung berdasarkan rumus yang diperoleh dengan cara yang sama seperti sebelumnya yaitu ∆t n n+1 n n n − 2Ti,1 + Ti+1,j Ti,1 = Ti,1 + 2 Ti−1,1 hx 2∆t n n Ti,2 − Ti,1 , (28) + 2 hx untuk i = 2, . . . , Nx . Temperatur yang terletak di (x1 , y1 ) diperoleh berdasarkan nilai rerata temperatur dari dua tetangga terdekatnya yaitu, n+1 T1,1
4
n+1 n+1 T2,1 + T1,2 . = 2
(29)
SYARAT STABILITAS
Penggunaan persamaan (21) harus memenuhi syarat stabilitas yang ungkapannya diperoleh dengan cara sebagai berikut. Pertama mensubstitusikan [4] n Ti,j = ξ n eikx ihx eiky jhy
(30)
untuk hx = hy ke dalam persamaan (21). √ Dalam persamaan tersebut i = −1, kx dan ky adalah bilangan gelombang spasial real, dan ξ = ξ(kx , ky ) adalah bilangan kompleks yang bergantung pada kx dan ky . ξ disebut faktor amplifikasi (amplification factor ). Berikutnya membagi kedua ruasnya dengan ξ n dan menggunakan kaitan e−ikx hx + eikx hx − 2 = −4 sin2 (kx hx /2) (31) yang menghasilkan 4∆t ξ = 1 − 2 sin2 (kx hx /2) + sin2 (ky hx /2) hx (32) Persamaan (21) akan stabil jika dipenuhi |ξ| ≤ 1. Suku dalam tanda kurung {. . .} dalam persamaan (32) maksimum bernilai 2. Karena itu syarat stabilitas terpenuhi jika 4∆t ≤ 1. h2x
(33)
Untuk pilihan hx = 0, 1 maka ∆t ≤ 0, 0025.
FOTON/Vol. 9 No. 2/Agustus 2005
71 5
Hari Wisodo PROGRAM DAN HASIL
Program untuk mengimplementasikan bentuk diskret persamaan konduksi panas tak tunak 2D ditulis dalam MATLAB 6.5. Kode programnya ditunjukkan dalam Lampiran A. Dalam simulasi digunakan hx = 0, 1 dan ∆t = 0, 0025. Hasilnya ditunjukkan pada Gambar 3. Warna bergerak mulai dari putih untuk T 0 = 0 atau T = T0 sampai hitam untuk T 0 = 1 atau T = Tb (diperoleh berdasarkan persamaan (8)). Berdasarkan Gambar 3 dapat diamati pola perambatan kalor. Kalor merambat mulai dari sisi-sisi yang temperaturnya Tb , Gambar 3 (b). Pada saat t = 3L2 /D, hampir dicapai kesetimbangan termal di seluruh bahan.
Gambar 3: Kontur konduksi panas pada plat homogen pada saat (a) t = 2 (b) t = 10 (c) t = 15 dan (d) t = 30 dalam satuan L2 /D.
6
KESIMPULAN
Konduksi panas tak tunak 2D pada plat homogen dapat disimulasikan berdasarkan hasil tinjauan komputasi numerik. Terdapat langkah-langkah yang harus dilakukan agar hasil komputasi numerik dapat memberikan informasi fisis mengenai perilaku sistem yang ditinjau. Pemilihan lebar kisi dan lebar langkah waktu harus memenuhi syarat stabilitas sebagai konsekuensi penggunaan metode beda hingga skema FTCS.
LAMPIRAN A :
KODE PROGRAM
function difusi L = 1; dt = 0.0025; hx = 0.1; perhx = 1/hx; perhx2 = perhx*perhx; Nx = L/hx; Ny = Nx; tmax = 20; for i = 1:Nx for j = 1:Ny T(i,j)=0; end end for j = 1:Ny T(Nx+1,j)=1; end for i = 1:Nx T(i,Ny+1)=1; end T(Nx+1,Ny+1)=1; for t=1:tmax for i = 2:Nx for j = 2:Ny T(i,j) = dt*perhx2*(T(i-1,j) +T(i+1,j)+T(i,j-1) +T(i,j+1)-4*T(i,j))+T(i,j); end end for j = 2:Ny T(1,j) = T(1,j) + 2*dt*perhx2*(T(2,j)-T(1,j)) + dt*perhx2*(T(1,j-1)-2*T(1,j)+T(1,j+1)); end for i = 2:Nx T(i,1) = T(i,1) + dt*perhx2*(T(i-1,1)-2*T(i,1)+T(i+1,1)) + 2*dt*perhx2*(T(i,2)-T(i,1)); end T(1,1)=0.5*(T(2,1)+T(1,2)); save suhu20.dat T -ascii clf pcolor(T) colorbar vert cxs=max(max(abs(T))); caxis([0 1]); shading interp drawnow end
PUSTAKA [1] Vesely, F. J. 1994. Computational Physics: An Introduction. New York: Plenum Press. hal. 14, 17.
FOTON/Vol. 9 No. 2/Agustus 2005
Pengkajian Konduksi Panas Tak Tunak 2D . . . [2] Alonso, M. dan Finn, E. J. 1992. DasarDasar Fisika Universitas Edisi Kedua Jilid 1 Mekanika Dan Termodinamika. Penerbit Erlangga. hal. 356-369. [3] Winiecki T. dan Adams C. S. 2002. A Fast Semi-Implicit Finite-Difference Method for the TDGL Equations. Journal of Computational Physics. Vol. 179 hal. 127 - 139. [4] Press, W.H., Teukolsky S.A., Vetterling, W.T., dan Flannery, B.P. 1992. Numerical Recipes in FORTRAN: The Art of Scientific Computing 2ed. New York: Cambridge University Press. hal. 845.
FOTON/Vol. 9 No. 2/Agustus 2005
72
