AIP/123-QED
Metode Lattice-Boltzmann, Aplikasi pada Kasus Difusi Kalor Ridho Muhammad Akbar∗ Departemen Fisika, Institut Teknologi Bandung, Indonesia (Dated: December 9, 2015)
Abstract Abstrak - Metode Lattice-Boltzmann adalah metode penyelesaian kasus fenomena transfer menggunakan cara pandang mesoskopik, melihat sistem secara statistik menggunakan fungsi distribusi. Pada tlisan ini, metode Lattice-Boltzmann digunakan untuk menyelesaikan fenomena transfer paling sederhana yaitu fenomena difusi kalor. Simulasi dilakukan dalam 1D dan 2D menggunakan model D1Q2 dan D2Q9. Simulasi menunjukkan bahwa metode Lattice-Boltzmann memberikan hasil yang cocok dengan solusi analitik dan solusi numerik menggunakan metode beda hingga finite difference.
∗
[email protected]
1
I.
PENDAHULUAN
Ada dua jenis metode pendekatan umum yang digunakan dalam simulasi fenomena transfer, baik itu transfer kalor, transfer massa, maupun transfer momentum. Kedua jenis metode tersebut adalah metode disrit dan metode kontinum. Pada pendekatan kontinum, sistem, misalnya fluida, dianggap sebagai benda yang kontinu walaupun sebenarnya jarak antar sati partikel fluida dengan molekul fluida yang lain sangat besar dibandingkan dengan ukuran molekul itu sendiri. Persamaan diferensial dapat diselesaikan dengan mengaplikasikan kekekalan energi, massa, dan momentum pada suatu area kecil (infinitesimal ) yang kita sebut dengan kontrol volume. Karena akan sangat sulit untuk menyelesaikan persamaan diferensial secara langsung, maka digunakanlah skema-skema penyelesaian seperti finite volume, finite difference, atau finite element untuk mengonversi persamaan diferensial menjadi persamaan aljabar yang mudah diselesaikan. Pada metode kontinum, hal pertama yang dilakukan adalah mendefinisikan persamaan diferensial yang menjelaskan keadaan sistem. Kemudian dilakukan diskritisasi daerah domain penyelesaian menjadi grid atau elemen. Besaran-besaran fenomenologi seperti temperatur, massa jenis, dan sebagainya direpresentasikan pada setiap grid dan bervariasi secara linier antara grid satu dengan grid yang lain. Pendekatan kontinum menyelesaikan sistem secara makroskopis. Kita dapat pandang bahwa nilai besaran tersebut sebagai nilai ratarata semua partikel yang ada di dalam grid tersebut. Masalah yang dijumpai pada metode ini adalah sulitnya menemukan solusi untuk sistem-sistem dengan geometri yang kompleks atau dengan syarat-syarat batas yang kompleks. Di sisi lain, pendekatan diskrit melihat sistem sebagai kumpulan partikel yang berinteraksi satu sama lain. Dengan demikian, kita diharuskan mendefinisikan gaya-gaya yang bekerja tiap partikel dan menyelesaikan hukum II Newton untuk mengetahui posisi dan kecepatan partikel. Pada metode ini, kita tidak bisa mendapatkan nilai besaran fenomenologi secara langsung namun dapat dihitung dari posisi dan kecepatan partikel. Misalnya, besaran temperatur berhubungan dengan energi kinetik rata-rata pertikel yang tidak lain adalah massa partikel dikalikan dengan setengah dari kuadrat kecepatan rata-rata pertikel. Begitu juga dengan besaran lain seperti tekanan, massa jenis, konsentrasi, dan sebagainya. Simulasi ini disebut dengan Molecular Dynamics (MD) dan melihat sistem secara mikroskopis kemudian menerjemahkannya menjadi besaran-besaran fenomenologi (makroskopis) melalui 2
mekanika statistik. Metode MD menjawab kekurangan-kekurangan yang dialami oleh metode kontinum. Karena MD berangkat dari persaman yang paling fundamental, hukum II Newton, dan melihat keadaan masing-masing partikel dalam sistem. maka kerumitan geometri dan syarat batas bukanlah masalah. MD bahkan dapat menangani sistem multifasa dengan atau tanpa transisi fasa dengan baik. Namun kelemahan dari MD adalah biaya komputasi yang panjang, baik waktu komputasi maupun memori. Bayangkan saja, dalam satu mol gas saja terdapat 1023 partikel dan harus dihitung interaksi antar satu partikel dengan masing-masing partikel yang lain. Muncul pula pertanyaan pada metode MD, untuk mengetahui keadaan makroskopis suatu sistem apakah perlu ktia mengetahui keadaan mikroskopis seluruh partikel penyusunnya? Muncul sebagai alternatif, metode Lattice-Boltzmann (LBM) memandang sistem di antara skala mikroskopis dan makroskopis. LBM tidak memandang sistem sebagai benda yang kontinum, tidak juga memperhitungkan sifat-sifat setiap partikel, melainkan memperhatikan sifat kolektif beberapa partikel sebagai suatu kesatuan yang direpresentasikan dengan fungsi distribusi. Setiap grid pada domain penyelesaian memiiki fungsi distribusi yang berbedabeda yang merepresentasikan keadaan mikroskopis kolektif pada grid tersebut dan dapat ditransformasi sehingga dapat dihitung nilai besaran makroskopisnya. Cara pandang ini dikenal dengan istilah pendekatan mesoskopis (di antara makroskopis dan mikroskopis).
Gambar 1. Ilustrasi perbandingan cara pandang sistem untuk metode kontinum (kiri), LBM (tengah), dan MD (kanan). Pada metode kontinum sistem dianggap sebagai benda yang kontinu sehingga solusi ada di semua grid tanpa mempertahikan tingkah laku masing-masing partikel. Pada metode LBM, tingkah laku partikel di setiap grid diperhitungkan secara kolektif sebagai fungsi distribusi. Pada metode MD, tingkah laku setiap partikel diperhitungkan.
3
II.
METODE LATTICE-BOLTZMANN
Kebanyakan fenonema transfer dideskripsikan secara matematis melalui persamaan diferensial parsial, misalnya persamaan Navier-Stokes pada kasus aliran fluida. Karena LBM tidak menyelesaikan sistem secara makroskopis, maka kita tidak perlu menyelesaikan persamaan Navier-Stokes (pada kasus aliran fluida) melainkan menyelesaikan persamaan transfer Boltzmann (Boltzmann Transport Equatoin).
A.
Persamaan Transfer Boltzmann
Ludwig Eduard Boltzmann, fisikawan Austria yang terkena dengan inovasinya dalam mekanika statistik menyatakan bahwa setiap sistem dapat dideskripsikan secara statistik melalui fungsi distribusi f (r, c, t), jumlah partikel dengan kecepatan antara c dan c + dc di antara posisi r dan r + dr dalam rentang waktu t sampai t + dt pada suatu asembli. Ketika gaya external F diberikan pada sistem, kita dapat menyatakan fungsi distribusinya F dt, t+dt) dengan m adalah massa partikel. Jika tidak ada tumbukan menjadi; f (r+cdt, c+ m
antar partikel maka tidak akan ada perubahan jumlah partikel pada interval keadaan drdc. Kita dapat tuliskan sebagai; f (r + cdt, c +
F dt, t + dt)drdc − f (r, c, t)drdc = 0 m
(1)
Akan tetapi, tentu saja tumbukan terjadi antar partikel sehingga ruas kanan pada persamaan (1) diganti dengan operator tumbukan Ω(f ) yang menyatakan laju perubahan fungsi distribusi per satuan waktu. Persamaan 1 dapat ditulis kembali menjadi; f (r + cdt, c +
F dt, t + dt)drdc − f (r, c, t)drdc = Ω(f )drdcdt m
(2)
bagi kedua ruas dengan drdcdt didapat; df ∂f ∂f F ∂f = c+ + = Ω(f ) dt ∂r ∂c m ∂t
(3)
atau dapat kita tulis menjadi; ~ ∂f F~ ∂f + ~c · ∇f + · = Ω(f ) ∂t m ∂c
(4)
Persamaan 4 adalah persamaan transfer Boltzmann yang perlu diselesaikan untuk mendapatkan solusi dari kasus fenomena transfer. Keindahan dari metode Lattice-Boltzmann 4
ini adalah pada kesederhanaan persamaan yang digunakan. Hanya dengan menggunakan persamaan transfer Boltzmann tersebut kita dapat menyelesaikan berbagai sistem fenomena transfer yang kompleks.
B.
Aproksimasi BGKW
Dibalik kesederhanaan persamaan transfer Boltzmann, permasalahan timbul pada operator Ω(f ) dimana sulit menentukan ekspresi eskplisit untuk operator tersebut. Akan tetapi pada tahun 1954, empat orang fisikawan, Bhatnagar, Gross, Krook, dan Welander (independen di tempat yang berbeda) memperkenalkan ekspresi untuk menyatakan operator Ω(f ) yang kemudian dikenal dengan nama Approksimasi BGKW; 1 (5) Ω(f ) = ω(f eq − f ) = (f eq − f ) τ dimana τ disebut sebagai waktu relaksasi. Besaran τ inilah yang menjadi jembatan antara dunia makroskopis (fenomenologi) dan dunia mikroskopis. Variabel f eq adalah fungsi distribusi equilibrium yang tidak lain adalah fungsi distribusi Maxwell-Boltzmann untuk kasus klasik. Aproksimasi tersebut terbukti sangat baik saat dicoba dalam kasus-kasus fenomean transfer partikel-partikel klasik. Dengan approksimasi ini, maka persamaan (4) dapat ditulis kembali menjadi; ~ ∂f F~ ∂f 1 + ~c · ∇f + · = (f eq − f ) (6) ∂t m ∂c τ Ide dari LBM adalah menyelesaikan persamaan (6) melalui dua tahap, collision/tumbukan dan streaming/aliran. Tahap tumbukkan adalah proses mencari fungsi distribusi equilibrium baru sedangkan proses aliran adalah proses perambatan fungsi distribusi ke grid berikutnya sesuai dengan arah rambat yang didefinisikan. Kedua tahap ini dilakukan di setiap selang waktu sampai batas waktu atau iterasi tertentu.
C.
Fungsi Distribusi Equilibrium
Untuk suatu sistem partikel dengan densitas ρ yang berada dalam suatu medium dengan kecepatan makroskipis ~v (kecepatan rata-rata tiap partikel), maka distribusi kecepatan tiap partikel (~c) dalam sistem tersebut memenuhi distribusi Maxwell-Boltzmann; " # 2 n −(~ c − ~ v ) f eq (~c) = exp (2πθ)2/3 2θ 5
(7)
dimana θ = kT , T adalah temperatur dan k adalah konstanta Boltzmann. Persamaan tersebut dapat ditulis kembali menjadi; " # h −c2 i 2 ~ c · ~ v − v n exp f eq (~c) = exp (2πθ)3/2 2θ 2θ
(8)
Mengingat expansi Taylor untuk fungsi ex ; ex = 1 + x +
x2 x3 xn + + ··· + 2 6 n!
(9)
kita dapat tulis persamaan (8) menjadi; " # h −c2 i 2 2 2 n ~ c · ~ v − v (~ c · ~ v − v ) f eq (~c) = exp 1+ + + ··· (2πθ)3/2 2θ 2θ 8θ2
(10)
dengan mengabaikan suku v n untuk n > 2, maka persamaan (10) dapat disusun kembali menjadi; f eq (~c) =
h −c2 ih i n 2 2 exp A + B(~ c · ~ v ) + C(~ c · ~ v ) + Dv (2πθ)3/2 2θ
(11)
dimana A, B, C,dan D adalah konstanta yang harus dihitung berdasarkan sistem yang ditinjau. Densitas n ternyata dapat dipandang sebagai suatu fungsi skalar yang mendeskripsikan keadaan makroskopis sistem, misalnya densitas massa (massa jenis), densitas energi (temperatur), densitas populasi (konsentrasi), dan sebagainya tergantung sistem yang ditinjau. h i Kemudian, komponen 1/(2πθ)3/2 exp − c2 /2θ dapat kita pandang sebagai suatu fungsi bobot yang bergantung pada kecepatan, misal w(~c). Kita dapat melakukan diskritisasi pada persamaan (11) untuk setiap arah vektor menjadi; h i 2 2 f (~c) = Φw(~c) A + B(~c · ~v ) + C(~c · ~v ) + Dv eq
(12)
dimana Φ adalah fungsi skalar (besaran makroskopis) dan wi adalah fungsi bobot.
III. A.
ALGORITMA DAN PEMROGRAMAN Diskritisasi Persamaan Transfer Boltzmann
Sebelum mengaplikasikan persamaan (6) dalam program komputer, kita perlu melakukan diskritisasi terlebih dahulu. Dalam kasus difusi dimana tidak ada gaya eksternal yang bek6
erja, maka persamaan (6) dapat didiskritisasi sebagai berikut; i fk (x, t + ∆t) − fk (x, t) fk (x + ∆x, t + ∆t) − fk (x, t + ∆t) 1h + ck · = − fk (x, t) − fkeq (x, t) ∆t ∆x τ (13) Diketahui bahwa ∆x = ck ∆t maka persamaan (13) dapat ditulis kembali menjadi; fk (x + ∆x, t + ∆t) = fk (x, t) 1 − ω + ωfkeq (x, t)
(14)
dimana ω = ∆t/τ disebut dengan frekuensi tumbukan. Dalam LBM, persamaan (14) diselesaikan melalui dua tahap yaitu tumbukkan dan aliran. Tahap tumbukkan dapat ditulis sebagai; fk (x, t + ∆t) = fk (x, t) 1 − ω + ωfkeq (x, t)
(15)
yang dapat dipandang sebagai proses mencari fungsi equilibrium baru. Sementara proses aliran dapat ditulis sebagai; fk (x + ∆x, t + ∆t) = fk (x, t + ∆t)
(16)
Fungsi distribusi equilibrium (12) sendiri dapat ditulis dalam bentuk diskrit menjadi; h i fkeq = Φwk A + B(ck · ~v ) + C(ck · ~v )2 + Dv 2 (17) kita pilih wk sedemikian sehingga
P
k
wk = 1. Dapat ditunjukkan dengan menggunakan inte-
gral Gauss dari persamaan (7) bahwa hubungan antara besaran makroskopis dan mesoskopis memenuhi hubungan; X
fk = Φ
(18)
fk ck = Φ~v
(19)
k
X k
B.
Model Kecepatan
Persamaan 13-17 menyatakan probabilitas partikel mengalir pada arah k dengan kecepatan ck . Pada LBM, beberapa jumlah kemungkinan arah kecepatan ditentukan dan biasa dinyatakan dalam dimbol DmQn (m dimensi dengan n kemungkinan arah kecepatan). Dalam tulisan ini, akan ditinjau model kecepatan D1Q2 untuk menyelesaikan kasus difusi 1D dan D2Q9 untuk menyelesaikan kasus 2D. Banyak model-model kecepatan lain yang dapat digunakan untuk menyelesaikan kasus-kasus fenomena transfer. 7
Gambar 2. Model kecepatan D2Q9 (a) dan model kecepatan D1Q2 (b). Pada model D2Q9, arah no.9 adalah keadaan partikel tidak bergerak, yang tidak ada pada model D1Q2. Dapat pula dibuat model D1Q3 dan mempertimbangkan partikel diam untuk kasus 1D, namum model D1Q2 sudah cukup baik menyelesaikan kasus difusi 1D.
Nilai wk untuk masing-masing arah kecepatan pada model D2Q9 adalah; wk = 4/9; k = 9
(20)
wk = 1/9; k = 1, 2, 3, 4
(21)
wk = 1/36; k = 5, 6, 7, 8
(22)
sementara untuk model D1Q2, nilai wk nya adalah; wk = 1/2; k = 1, 2
(23)
Gambar 3 menunjukkan ilustrasi proses aliran fungsi distribusi (persamaan (16)) dalam 1D sesuai dengan model kecepatan D1Q2.
Gambar 3. Ilustrasi proses aliran 1D. f1 mengalir ke arah kanan dan f2 mengalir ke arah kiri sesuai dengan model kecepatan pada Gambar 2
8
Gambar 4. Diagram Alir metode LBM. Nomor di dalam tanda kurung menyatakan nomor persamaan yang digunakan dalam proses yang bersangkutan. Untuk tahap syarat batas akan dibahas pada bagian berikutnya C.
Diagrm Alir Algoritma
Gambar 4 menunjukkan diagram alir metode LBM untuk menyelesaikan kasus-kasus fenomena transfer. Algoritma LBM terlihat sangat sederhana dan mudah diaplikasikan untuk berbagai kasus, meskipun dalam geometri yang rumit. Pada bagian berikutnya, diagram alir ini akan diaplikasikan dalam penyelesaian kasus difusi termal dalam 1D dan 2D.
IV.
STUDI KASUS PERSAMAAN DIFUSI
Persamaan difusi termal dinyatakan sebagai beriku; ρC
∂T = ∇ · k∇T ∂t 9
(24)
dimana T, ρ, C, k masing-masing adalah temperatur, massa jenis, kalor jenis, dan konduktifitas medium. Jika ρ, C dan k konstan, maka persamaan (24) menjadi; ∂T = α∇2 T ∂t dimana α adalah difusifitas termal medium.
(25)
Kita tidak akan menyelesaikan persamaan difusi (25) secara langsung, melainkan melalui persamaan (13) dimana besaran α terkait dengan besaran τ melalui hubungan; (∆x)2 1 α= τ− (26) ∆tD 2 dimana ∆x adalah ukuran partisi bidang (ukuran grid), ∆t adalah ukuran partisi waktu (selang waktu), dan D adalah dimensi sistem yang ditinjau. Misalka untuk model kecepatan D2Q9 maka D = 2 dan D = 1 untuk model kecepatan D1Q2. Hubungan antara α (besaran makroskopis) dengan τ (besaran mikroskopis) ini dapat ditunjukkan melalui expnasi Chapman-Enskog yang tidak dijelasakan dalam tulisan ini. Pada tulisan ini, semua simulasi dilakukan dengan menetapkan nilai α = 0.25, Pada kasus difusi, tidak ada perpindahan massa dari medium, artinya medium tidak bergerak dan suku-suku yang mengandung kecepatan pada persamaan (17) bernilai nol. Sehingga fungsi distribusi equilibrium untuk kasus difusi adalah; fkeq = Awk Φ
(27)
persamaan (27) harus memenuhi persamaan (18) sehingga kita dapatkan nilai A = 1. Setiap persoalan sistem fisis memiliki syarat batas yang berbeda-beda. Pada sub bagian berikut akan dibahas beberapa permasalahan sistem fisis beserta syarat-syarat batasnya. Perlu dicatat bahwa semua besaran yang digunakan dalam tulisan ini tidak memiliki satuan (dimensionless unit). Hal ini biasa digunakan pada simulasi-simulasi fisika, transformasi sederhana dapat digunakan untuk mengonversi dimensionless unit tersebut menjadi satuan fisis yang berarti. Namun bagaimanapun, tujuan dari tulisan ini hanyalah untuk menunjukkan bahwa LBM dapat digunakan untuk menyelesaikan persamaan difusi, sehingga transformasi demikian tidak dilakukan dalam tulisan ini.
A.
Difusi 1D, Syarat Batas Temperatur Konstan (Syarat batas Dirichlet)
Anggaplah ada suatu batang dengan panjang 100 unit dengan temperatur mula-mula 0. Pada waktu t = 0 sisi kanan batang dihubungkan dengan suatu reservoir panas dengan 10
Gambar 5. Kasus Difusi 1D dengan syarat batas temperatur konstan
temperatur konstan 1. Kasus ini ditunjukkan oleh Gambar 5. Misalkan nilai α = 0.25, maka dengan LBM akan kita tentukan distribusi temperatur sepanjang batang. Pada sisi kiri (x = 0), nilai f1eq (0) = w1 Tw = 0.5Tw dan f2eq (0) = w2 Tw = 0.5Tw . Nilai f2 (0) dapat diketahui dari proses aliran, sehingga nilai f1 (0) harus kita tentukan. Diketahui bahwa T (0) = f1 (0) + f2 (0) = Tw maka kita dapatkan; f1 (0) = Tw − f2 (0). Pada sisi kanan (x = N ), nilai f1eq (N ) = f2eq (N ) = 0.5T (N ) = 0. Nilai f1 (N ) dapat diketahui dari proses streaming, sehingga nilai f2 (N ) perlu kita tentukan. Dengan cara yang sama dengan sebelumnya kita dapatkan bahwa f2 (N ) = −f1 (N ).
B.
Difusi 1D, Syarat Batas Adiabatik (Syarat Batas Neumann)
Persoalan yang sama seperti kasus 1, namun sisi kanan batang ditutup dengan isolator sehingga ujung kanan batang bersifat adiabatik, tidak ada kalor yang keluar atau masuk dari batang, seperti digambarkan pada Gambar 6. Hukum Fourier tentang transfer kalor menyatakan bahwa flux kalor berbanding lurus dengan gradien temperatur; 00
q = −k
Ti − Ti−1 ∂T = −k ∂x ∆x
(28)
pada kondisi adiabatik, tidak ada transfer kalor sehingga didapatkan Ti = Ti−1 atau pada kasus kedua ini, T (N ) = T (N − 1). Karena T (N ) = f1 (N ) + f2 (N ) dan f1 (N ) = f1 (N − 1) (proses aliran) maka dapat kita simpulkan bahwa f2 (N ) = f2 (N − 1).
Gambar 6. Kasus Difusi 1D dengan syarat batas adiabatik
11
C.
Difusi 2D pada Plat Persegi
Pada kasus ini, kita tinjau peristiwa difusi 2D pada pelat persegi. Misalkan sebuah pelat persegi seperti pada Gambar 7. Dua variasi akan ditinjau pada kasus ini. Pertama adalah apabila semua sisi pelat dibuat memiliki temperatur konstan seperti pada Gambar 7(a). Variasi kedua adalah sisi selatan pelat dibuat adiabatik dengan menempelkannya pada isolator seperti gambar 7(b).
Gambar 7. Kasus Difusi 2D dengan syarat batas (a) temperatur konstan di semau sisi, syarat batas Dirichlet, dan (b) salah satu sisi dibuat adiabatik, syarat batas Neumann.
Untuk kasus (a), kita dapat menggunakan prinsip yang sama seperti yang kita lakukan pada kasus difusi 1D dengan syarat batas Dirichlet. Pada sisi barat pelat, f2 , f3 , f4 , f6 , dan f7 dapat diketahui melalui proses aliran, maka kita perlu menghitung f1 , f5 , dan f8 . Berdasarkan persamaan proses tumbukkan (15) dan aliran (16), kita dapat tuliskan hubungan; f1 − f1eq = f3 − f3eq , f5 − f5eq = f7 − f7eq , dan f8 − f8eq = f6 − f6eq . Kita ketahui bahwa fkeq (x, y) = wk T (x, y). Dari sini, nilai f1 , f5 , dan f8 dapat dihitung melalui hubungan; f1 (0, y) = w1 TW + f3 (0, y) − w3 TW
(29)
f5 (0, y) = w5 TW + f7 (0, y) − w7 TW
(30)
f8 (0, y) = w8 TW + f6 (0, y) − w6 TW
(31)
dimana index W melambangkan sisi barat (west). Dengan cara yang sama, kita dapatkan persamaan untuk menghitung fungsi distribusi yang tidak diketahui untuk sisi timur, utara, 12
dan selatan; f3 (Lx , y) = w3 TE + f1 (Lx , y) − w3 TE
(32)
f7 (Lx , y) = w7 TE + f5 (Lx , y) − w7 TE
(33)
f6 (Lx , y) = w6 TE + f8 (Lx , y) − w6 TE
(34)
f4 (x, Ly ) = w4 TN + f2 (x, Ly ) − w2 TN
(35)
f7 (x, Ly ) = w7 TN + f5 (x, Ly ) − w5 TN
(36)
f8 (x, Ly ) = w8 TN + f6 (x, Ly ) − w6 TN
(37)
f2 (x, 0) = w2 TS + f4 (x, 0) − w4 TS
(38)
f5 (x, 0) = w5 TS + f7 (x, 0) − w7 TS
(39)
f6 (x, 0) = w6 TS + f8 (x, 0) − w8 TS
(40)
Untuk kasus (b), kita harus menentuka syarat batas fungsi distribusi pada sisi selatan sedemikian sehingga perpindahan kalor pada arah sumbu-y sama dengan nol. Kita dapat mengaplikasikan persamaan (28) sehingga didapat hubungan, T (x, 0) = T (x, 1). Karena P P P T (x, y) = i fi (x, y) maka kita mendapat hubungan; i fi (x, 0) = i fi (x, 1). Masuk akal apabila kita mengambil kesimpulan bahwa, fi (x, 0) = fi (x, 1). Untuk fungsi distribusi yang tidak diketahui pada sisi-sisi lain dapat ditentukan dengan cara yang sama dengan kasus (a).
D.
Difusi 2D pada Channel
Pada kasus ini, kita meninjau suatu sistem yang terdiri dari bidang reservoir dan suatu bidang penampungan yang keduanya terhubung oleh suatu channel dengan lebar a seperti pada Gambar 8(a). Perlu diketahui meskipun T pada tulisan ini adalah temperatur, namun pada kenyataannya T bisa berupa fungsi skalar lainnya seperti molaritas gas, densitas elektron, dan zat-zat lain yang dapat mengalami difusi. Bagaimanapun, karena difusi kalor yang sedang kita tinjau, maka kita anggap saja bahwa T adalah temperatur dan bidang penyimpanan yang dimaksud adalah bidang penyimpanan kalor. Artinya, pada kasus ini, ktia mencoba memanaskan pelat di bagian kanan dengan cara menghubungkan pelat tersebut pada reservoir kalor melalui sebuah batang konduktor penghubung. Jika difusi gas yang sedang kita tinjau, maka channel dapat diartikan seba13
gai selang atau jika difusi elektron yang sedang kita tinjau, maka channel dapat diartikan sebagai kabel atau logam penghantar.
Gambar 8. (a)Kasus Difusi 2D dari reservoir menuju bidang penampungan melalui channel dengan lebar a dan (b) penyederhanannya
Karena reservoir berada dalam keadaan tunak (tidak terjadi perubahan apa-apa pada reservoir), maka sistem dapat disederhanakan menjadi seperti Gambar 8(b). Kita dapat memvariasikan nilai a dan melihat pengaruhnya pada distribusi temperatur di bidang penampungan. Syarat-syarat batas yang digunakan pada kasus ini dapat diturunkan dengan cara yang sama dengan pada kasus difusi 2D pada pelat persegi.
V.
HASIL DAN PEMBAHASAN A.
Difusi 1D, Syarat Batas Temperatur Konstan (Syarat batas Dirichlet)
Gambar 9 menunjukkan hasil simulasi difusi panas 1D dengan syarat batas temperatur konstan beserta perbandingan hasil secafa analitik dan secara numerik melalui metode finite difference (FDM). Secara analitik, distribusi temperatur untuk kasus difusi 1D, dengan syarat batas seperti pada kasus yang ditinjau di tulisan ini, pada waktu t yang jauh dari waktu steady state, memenuhi persamaan; T (x, t) = T (0, t) − erf
14
x √ 2 αt
! (41)
dimana; 2 erf(x) = √ π
Z
x
2
e−t dt
(42)
0
Gambar 9. Hasil simulasi LBM untuk difusi 1D dengan syarat batas temperatur konstan pada iterasi ke-100, 2000, dan 20000. Solusi analitik dan hasil simulasi dengan Finite Difference Method (FDM) disajikan pula sebagai pembanding.
Sementara saat waktu sudah mendekati waktu steady state, maka distribusi temperatur akan memenuhi persamaan linier, yang dapat diturunkan langsung dari persamaan difusi; T (x) = T (0) −
x 100
(43)
Dapat dilihat bahwa solusi yang didapatkan dari LBM dan FDM cocok dengan solusi analitik. 15
B.
Difusi 1D, Syarat Batas Adiabatik (Syarat batas Neumann)
Gambar 10 menunjukkan solusi hasil simulasi untuk kasus kedua. Terlihat bahwa FDM dan LBM menghasilkan solusi yang serupa.
Gambar 10. Hasil simulasi LBM dan FDM untuk difusi 1D dengan syarat batas adiabatik pada iterasi ke 100, 2000, dan 20000
C.
Difusi 2D Pada Plat Persegi
Gambar 11 menunjukkan distribusi temperatur di dalam plat setelah 10000 iterasi untuk kasus syarat batas temperatur konstan. Terlihat bahwa simulasi LBM dan FDM memberikan hasil yang identik. Perbandingan yang lebih jelas ditunjukkan pada Gambar 12 yang menyajikan profil temperatur sepanjang sumbu-x pada posisi y = Ly /2 serta sepanjang sumbu-y pada posisi x = Lx /2. 16
Gambar 11. Distribusi temperatur setelah 10000 iterasi, hasil simulasi (a) LBM dan (b) FDM untuk difusi 2D pada plat persegi dengan syarat batas temperatur konstan
Gambar 12. Profil temperatur sepanjang sumbu-x pada posisi y = Ly /2 dan sepanjang sumbu-y pada posisi x = Lx /2. LBM dan FDM memberikan hasil yang identik.
Untuk kasus dengan geometri yang sama namun menggunakan syarat batas Neumann (Adiabatik pada sisi selatan), distribusi temperatur dalam pelat ditunjukkan oleh Gambar 13. Dapat dilihat kembali bahwa metode LBM dan FDM memberikan hasil yang identik. 17
Gambar 13. Distribusi temperatur dalam pelat apabila bagian selatan pelat diberi isolator sehingga adiabatik. (a) LBM dan (b) FDM memberikan hasil yang identik.
D.
Difusi 2D pada Channel )
Pada kasus ini, disimulasikan difusi melalui 3 variasi lebar celah a yaitu a = 10, a = 20, dan a = 60. Hasil yang ditunjukkan oleh Gambar 14, Gambar 15, dan Gambar 16. Hasil menunjukkan bahwa semakin besar lebar channel maka proses difusi dari reservoir di sisi kiri menuju pelat di sisi kanan menjadi lebih cepat. Pada iterasi ke-20000, temperatur pelat di sisi kanan berada pada nilai sekitar 0.5 untuk a = 10, 0.65 untuk a = 20 dan 0.85 untuk a = 60.
VI.
KESIMPULAN
Pada tulisan ini telah ditunjukkan bahwa metode Lattice-Boltzmann dapat digunakan untuk menyelesaikan kasus difusi tanpa menyelesaikan secara langsung persamaan difusi melainkan menyelesaikan persamaan Transfer Boltzmann (Boltzmann Transport Equation) yang disesuaikan dengan kasus difusi. Penyelesaian dengan cara pandang yang berbeda ini terbukti sesuai dengan hasil analitik dan memberikan hasil yang identik dengan metode numerik lain yaitu metode beda hingga (Finite Difference Method ). Pada umumna, metode Lattice-Boltzmann digunakan dalam menyelesaikan kasus-kasus aliran fluida. Tulisan ini menunjukkan bahwa Lattice-Boltzmann juga dapat digunakan un18
tuk menyelesaikan persamaan difusi dan tidak menutup kemungkinan untuk aplikasi pada kasus perambatan gelombang (persamaan gelombang), baik gelombang mekanik, gelombang elektrobagnetik, atau penyelesaian persamaan gelombang pada fisika kuantum (Persamaan Schrodinger). Tantangannya adalah menentukan expresi dari fungsi distribusi equilibrium yang sesuai dengan sistem serta menentukan approksimasi yang sesuai untuk operator tumbukan dimana approksimasi BGKW mungkin tidak relevan, misalnya, untuk fenomena kuantum.
VII.
PERNYATAAN
Penulis mengucapkan terima kasih pada Prof. Umar Fauzi dari kelompok penelitian Fisika Batuan Departemen Fisika Institut Teknologi Bandung atas diskusi yang bermanfaat. Bagi pembaca, simulasi yang menunjukkan proses difusi termal pada kasus Difusi 2D pada channel dapat dilihat pada tautan berikut: https://youtu.be/vjY8eul1rPU. Kode pemrograman (MATLAB) dapat diunduh pada tautan berikut: https://drive.google. com/file/d/0BwQ0tdA_6p05ZjNxTW04eElSRVE/view?usp=sharing. Komentar, diskusi dan pertanyaan mengenai tulisan ini dapat dikirimkan melalui alamat e-mail ridho.akbar@s. itb.ac.id
[1] Mohamad, A. A.. 2011. Lattice Boltzmann Method, Fundamental and Engineering Application. London: Springer. [2] Wagner, Alexander. 2008. A Practical Introduction to Lattice-Boltzmann Method. Departemen Fisika, Universitas North Dakota, Amerika Serikat.
19
Gambar 14. Distribusi temperatur pada bidang untuk lebar channel a = 10. Simulasi untuk waktu t = 1, t = 5000, t = 10000, t = 20000
20
Gambar 15.
Distribusi temperatur pada bidang untuk lebar channel a = 20. Simulasi untuk
waktu t = 1, t = 5000, t = 10000, t = 20000
21
Gambar 16. Distribusi temperatur pada bidang untuk lebar channel a = 60. Simulasi untuk waktu t = 1, t = 5000, t = 10000, t = 20000
22
