11
PEMBAHASAN Pada karya ilmiah ini, persamaan Boltzmann yang akan dicari solusinya adalah persamaan Boltzmann spasial homogen, yaitu persamaan Boltzmann dengan ∂f ∂x bernilai nol, dituliskan: ∂f ∂t =
∫ ∫ g (cos θ )[ f (v ) f (w ) − f (v ) f (w )]dwde. '
2
S R
'
(24)
3
Ruas kiri persamaan Boltzmann mengandung bentuk diferensial, sehingga solusi masalah tersebut bergantung pada nilai awal yang dipilih. Selanjutnya akan dicari solusi eksak dan solusi numerik dari persamaan Boltzmann dengan menggunakan nilai awal Bobylev. Solusi numerik akan dicari dengan simulasi menggunakan Metode DSMC satu dimensi.
Solusi Eksak Persamaan Boltzman dengan Nilai Awal Bobylev
Misalkan dipilih nilai awal Bobylev berikut: 2 2 f (v ,0) = f 0 (v ) = ( A0 + B0 v ) exp⎛⎜ − C 0 v ⎞⎟ , ⎠ ⎝
(25)
dengan A0 , B0 ≥ 0 dan C 0 > 0 . Selanjutnya, akan dicari nilai dan hubungan antara koefisien A0 , B0 dan C 0 dengan memanfaatkan sifat-sifat sebagai berikut: 1.
Karena f 0 (v ) merupakan fungsi kepekatan peluang, maka: ∫ f 0 (v )dv = 1 .
2.
3 1 v f 0 (v )dv , Dari Persamaan (20), T = ∫ 2 d 2
3.
β 0 = 2TC 0 − 1 .
2
Maka, dari sifat 1,2 dan 3 (lihat Lampiran 4) diperoleh: 3
⎡ (1 + β 0 ) 2 ⎤ ⎡ (1 + β 0 ) ⎤ 2 ⎧ ⎡ (1 + β 0 ) 2 3 ⎤ ⎫ f (v ,0) = F (v , β 0 ) = ⎢ v v ⎥ 1 exp β + − ⎨ ⎬ 0 ⎢− 2T ⎥ ⎢ 2T 2 ⎥⎦ ⎭ ⎣ ⎦ ⎣ 2πT ⎦ ⎩ ⎣ (26)
12
dengan
β 0 = 2TC 0 − 1 ,
parameter
0 ≤ β0 ≤
2 3
menunjukkan
simpangan
kesetimbangan dari sebaran awal. Selanjutnya, akan dicari solusi spasial homogen Persamaan Boltzmann, yaitu fungsi kepekatan peluang partikel pada saat t, yaitu f (v, t ) = F (v , β (t )) , dengan β (t ) merupakan simpangan dari sebaran awal pada saat t dan β (0) = β 0 . Pertama, akan dihitung nilai integral dari fungsi kerapatan partikel pada saat t = 0 terhadap v. Dengan melakukan substitusi koordinat bola (lihat Lampiran 4), diperoleh:
∫
π
f (0, v ) dv =
3
2 5
2C 0 2
R3
(2 A0 C0 + 3B0 ) .
(27)
Dengan demikian, fungsi kerapatan partikel pada waktu t didefinisikan sebagai:
∫
d=
R
f (t , v )dv =
π
3
2
5
(2 AC + 3B ) .
(28)
2C 2
3
Akan dicari nilai A, B, C yang memenuhi persamaan di atas. Nilai tersebut dapat diperoleh dari dua fungsi lainnya, yaitu dari aliran impuls dan temperatur. Dari Persamaan (18), aliran impuls didefinisikan sebagai:
M=
∫ vv R
T
f (t , v )dv ,
3
sehingga M (aliran impuls) dapat dinyatakan dalam bentuk: M =
π
3
4
2
(2 AC + 5B ) I 7
3,
(29)
C 2
Bukti persamaan ini dapat dilihat pada Lampiran 6. Di sisi lain, temperatur (T) didefinisikan sebagai (lihat lampiran 7): 3
1 π 2 (2 AC + 5B ) T= trM = , 7 3d 4d C 2 sedemikian sehingga:
(30)
13
3
dC 2 ⎛ 5 ⎞ A= − 3TC ⎟ , ⎜ 3 ⎝2 ⎠ 2
(31)
π
dan 5
B=
dC 2
π
3
(2TC − 1) .
(32)
2
Dengan demikian, fungsi kerapatan partikel pada saat t, f (t , v ) , dapat dinyatakan sebagai: f (t , v ) =
5 7 d ⎛ 5 3 ⎞ ⎜⎜ 2T v 2 C 2 − ⎛⎜ 3T + v 2 ⎞⎟C 2 + C 2 ⎟⎟ exp⎛⎜ − C v 2 ⎞⎟ , 3 ⎝ ⎠ ⎝ ⎠ 2 ⎠ π 2⎝
(33)
sedemikian sehingga, fungsi sebaran partikel pada saat t sangat bergantung pada kecepatan partikel tersebut. Ruas kiri persamaan Boltzmann (24) merupakan turunan fungsi sebaran partikel terhadap waktu t, dengan f (t , v ) pada persamaan (33) diperoleh: 1 ⎡ 4 2 15 ⎤ 2 C 2 (1 − 2CT )⎢C 2 v − 5C v + ⎥ exp⎛⎜ − C v ⎞⎟ ∂C ∂t . 3 ⎝ ⎠ 4⎦ ⎣ 2
d
∂f ∂t =
π
(34)
Selanjutnya, akan dicari nilai ruas kanan dari persamaan Boltzmann secara bertahap sebagai berikut: I[ f , f ] =
∫ ∫ g (cos θ )[ f (v ) f (w ) − f (v ) f (w )]dwde. '
'
S 2 R3
2 ⎛ 2 2 2⎞ 2 2 f (v ' ) f ( w ' ) − f (v ) f ( w ) = B 2 ⎜⎜ v ' w ' − v w ⎟⎟ exp⎛⎜ − C ⎛⎜ v + w ⎞⎟ ⎞⎟ ⎠⎠ ⎝ ⎝ ⎝ ⎠
= B 2 ⎡ U ,V ⎢⎣
2
2 2 2 − v U T eeT U ⎤ exp⎛⎜ − C ⎛⎜ v + w ⎞⎟ ⎞⎟. ⎥⎦ ⎝ ⎠⎠ ⎝
(35)
Dengan memanfaatkan substitusi dan penyederhanaaan pada persamaan (8), (9) dan (11), maka: T ∫ g (cosθ )ee de = µI + S
dengan:
2
UU T v
2
(η − 3µ ),
(36)
14
π
π
0
0
v+w U= , µ = π ∫ g (cos θ )sin 3 θdθ ,η = 2π ∫ g (cos θ )sin θdθ . 2 Ruas kanan persamaan Boltzmann dapat dituliskan kembali dalam bentuk (lihat Lampiran 11): 3
B 2π 2 µ ⎡ 2 4 2 15 ⎤ 2 C v − 5C v + ⎥ exp⎛⎜ − C v ⎞⎟. I( f , f ) = ⎢ 7 ⎠ ⎝ 4⎦ ⎣ 2C 2
(37)
Dengan mempergunakan bentuk baru ruas kanan (37) dan kiri persamaan Boltzmann (34), diperoleh: 2 1 d ⎡ 4 2 15 ⎤ −C v ∂C ∂t C 2 (1 − 2CT )⎢C 2 v − 5C v + ⎥ e 3 4⎦ ⎣ 2 π 3 B 2π 2 µ = 7 C 2
⎡ 2 4 2 15 ⎤ 2⎞ ⎛ ⎢⎣C v − 5C v + 4 ⎥⎦ exp⎜⎝ − C v ⎟⎠,
∂C ∂t = −
dµ C (2CT − 1) . 2
(38)
Substitusi β = 2TC − 1 dan α =
dµ , menghasilkan: 2
∂β = −α (β + 1)β . ∂t
(39)
Persamaan (39) merupakan persamaan diferensial biasa, sehingga solusinya adalah:
β (t ) =
Ae −αt (1 − Ae −αt ) , A konstanta.
(40)
Dengan melakukan substitusi t = 0, pada Persamaan (40), diperoleh:
A=
β0 β (0) = , (1 + β (0)) (1 + β 0 )
(41)
sehingga,
β (t ) =
β 0 e −αt
1 + β 0 (1 − e −αt )
.
(42)
Jadi, solusi eksak persamaan Boltzmann homogen spasial mengambil bentuk yang sama dengan fungsi awal f (v , β 0 ) dengan menggantikan β 0 dengan β (t ) , yaitu:
15
f (v , t ) = F (v , β (t )) 3
⎛ (1 + β (t )) 2 ⎞ (43) ⎡ (1 + β (t )) 2 3 ⎤ ⎫ ⎡ (1 + β (t )) ⎤ 2 ⎧ ( ) 1 exp β + − = d⎢ t v v ⎟, ⎜− ⎨ ⎬ ⎢⎣ 2T 2 ⎥⎦ ⎭ 2T ⎠ ⎝ ⎣ 2πT ⎥⎦ ⎩ dengan β (t ) dari Persamaan (42) serta d > 0, T > 0, 0 < β t <
2 . 3
Besaran Makroskopik Gas
Jika fungsi sebaran partikel pada saat t diketahui, maka besaran makroskopik dari masalah nilai awal Bobylev dapat dihitung, sebagai berikut: 1. Fungsi Densitas d (t , x ) =
∫ f (t , x, v )dv
R3
⎡ (1+ β (t )) v 2T
3
− ⎡ (1 + β (t )) 2 3 ⎤ ⎫ ⎢⎣ ⎡ (1 + β (t )) ⎤ 2 ⎧ ( ) β 1 v t e + − = ∫ d⎢ ⎨ ⎬ ⎢ 2T 2πT ⎥⎦ ⎩ 2 ⎥⎦ ⎭ ⎣ 3 ⎣
2⎤
⎥ ⎦ dv
R
= 1.
2. Kecepatan Rata-rata: V (t , x ) =
1 d
1 = d
∫ vf (v, t )dv
R3
3
⎛ (1 + β (t )) 2 ⎞ ⎡ (1 + β (t )) ⎤ 2 ⎧ ⎡ (1 + β (t )) 2 3 ⎤ ⎫ ∫ vd ⎢⎣ 2πT ⎥⎦ ⎨⎩1 + β (t )⎢⎣ 2T v − 2 ⎥⎦ ⎬⎭ exp⎜⎝ − 2T v ⎟⎠dv 3
R
= 0. 3. Aliran Impuls M (t , x ) =
T ∫ vv f (t , x, v )dv
R3
=
3 T ⎡ (1 + β (t )) ⎤ 2 ⎧
∫ vv d ⎢⎣ 3
R
= dTI 3 .
4. Aliran Energi
2πT
⎥⎦
⎡ (1+ β (t )) v 2T
− ⎡ (1 + β (t )) 2 3 ⎤ ⎫ ⎢⎣ v − ⎥ ⎬e ⎨1 + β (t )⎢ 2 ⎦⎭ ⎣ 2T ⎩
2⎤
⎥ ⎦ dv
16
r (t , x ) =
1 2
1 = 2
2
∫v v R
f (t , x, v )dv
3
∫ vv
R3
3 2 ⎡ (1 + β (t )) ⎤ 2 ⎧
d⎢ ⎣
⎡ (1 + β (t )) 2 3 ⎤ ⎫ ⎛ (1 + β (t )) 2 ⎞ v − ⎥ ⎬ exp⎜ − v ⎟dv ⎨1 + β (t )⎢ 2T 2 ⎦⎭ ⎣ 2T ⎝ ⎠ ⎩
⎥⎦
2πT
= 0. 5. Energi E (t , x ) =
1 2
1 = 2 =
∫v R
∫
R3
2
f (t , x , v )dv
3
3 2 ⎡ (1 + β (t )) ⎤ 2 ⎧
v d⎢ ⎣
2πT
⎥ ⎦
⎛ (1 + β (t )) 2 ⎞ ⎡ (1 + β (t )) 2 3 ⎤ ⎫ v − ⎥ ⎬ exp⎜ − v ⎟ dv ⎨1 + β (t )⎢ 2 ⎦⎭ 2T ⎝ ⎠ ⎣ 2T ⎩
3 dT . 2
Simulasi dengan menggunakan metode DSMC satu dimensi
Metode DSMC merupakan salah satu metode yang dapat digunakan untuk mensimulasikan mekanisme tumbukan secara langsung. Inti dari metode ini adalah membuat representasi sederhana mengenai sebaran awal partikel, gerak, tumbukan dan pemberian indeks terhadap setiap partikel. Seperti program simulasi yang lainnya, program ini juga mengalami beberapa penyederhanaan, antara lain pada jumlah partikel yang dijadikan subjek pengamatan dan pada dimensi posisi yang digunakan. Posisi partikel diperhatikan hanya berdasarkan sumbu x saja. Misalkan ruang yang dipergunakan sebagai sistem ada pada
0 ≤ x ≤ 1 , 0 ≤ y ≤ 1 dan 0 ≤ z ≤ 1 . Untuk memudahkan proses inisialisasi posisi awal partikel, sumbu x dibagi menjadi beberapa sel dan setiap selnya dibagi lagi menjadi beberapa sub sel (Bird 1994). Proses simulasi dilakukan dengan menggunakan MATLAB 7.0. Program utama diberi nama NSBIC.m, yaitu program untuk menguji prosedur tumbukan pada gas seragam sederhana. Program NSBIC.m dibuat dengan algoritma sebagai berikut: 1. Menentukan nilai awal variabel dan aliran partikel pada t = 0. 2. Menentukan nilai awal variabel sampel.
17
3. Mendeskripsikan pergerakan sejumlah partikel pada selang waktu tertentu. 4. Menentukan urutan partikel dalam sel dan sub sel. 5. Menghitung banyaknya tumbukan selama selang waktu tertentu. 6. Menentukan contoh aliran partikel. 7. Menampilkan hasil. Sistem dibagi menjadi beberapa sel dan subsel. Pada simulasi kali ini, sistem akan dibagi menjadi 50 sel dan 400 sub sel dengan jumlah partikel maksimal 1000 partikel. Secara keseluruhan, program cukup besar sehingga perlu dipecah menjadi beberapa subroutine. Pada program NSBIC.m, terdapat 8 buah subroutine yaitu sebagai berikut: 1. DATAOS.m Subroutine ini berisi data awal yang berkaitan dengan sifat-sifat fisis partikel gas seperti kerapatan, suhu, banyak partikel sebenarnya yang disimulasikan oleh partikel simulasi, interval waktu (time step), jumlah subsel pada masing-masing sel, massa serta diameter partikel, tetapan kekentalan, serta scattering parameter. 2. INITOS4.m Subroutine ini berisi nilai variabel awal dan sebaran partikel pada saat
t = 0 . Nilai awal yang didefinisikan antara lain adalah konstanta Boltzmann, k = 1.3806e − 23 , collision cross section, informasi geometri setiap sel dan sub sel (termasuk nomor sel dan sub sel), serta kecepatan awal masingmasing partikel. 3. SAMPIOS.m SAMPIOS.m merupakan sub routine yang berisi inisialisasi seluruh variabel sampling, antara lain: banyaknya tumbukan pada t = 0, jumlah sampel, banyaknya partikel yang berpindah posisi, serta banyaknya partikel yang terseleksi untuk bertumbukan dan terpisah lagi. 4. MOVEOS.m MOVEOS.m merepresentasikan gerak perpindahan partikel dari satu posisi ke posisi lainnya selama selang waktu tertentu berdasarkan posisi di sumbu x, melakukan pendataan terhadap sel sebelum dan setelah tumbukan, dengan asumsi tumbukan yang terjadi antar partikel dan dengan dinding
18
pembatas merupakan tumbukan lenting sempurna, sedemikian sehingga gerakan pantulnya mengikuti sifat pantulan cermin. 5. INDEXS.m Subroutine INDEX.m mengatur penomoran partikel berdasarkan susunan sel dan sub selnya. 6. COLLS3.m COLLS3.m merupakan subroutine yang mensimulasikan tumbukan antara dua partikel, yaitu mengatur partikel partikel yang akan bertumbukan serta menghitung kecepatan relatif partikel, sudut elevasi, azimuth, sudut defleksi, serta kecepatan partikel setelah tumbukan. 7. SAMPLEOS.m SAMPLEOS.m melakukan sample terhadap partikel dalam aliran. 8. OUTOS.m OUTOS.m bertugas menampilkan hasil output pada setiap langkah waktu tertentu secara terus menerus. Hasil dari simulasi selama selang waktu t tertentu menghasilkan pola sebaran sebagai berikut:
Gambar 1 Kurva sebaran kecepatan partikel hasil solusi eksak dan solusi numerik dengan nilai awal Bobylev.
19
Kurva merah menunjukkan sebaran kecepatan yang diperoleh melalui solusi eksak. Kurva biru menunjukkan solusi hasil simulasi. Karena fungsi sebaran partikel terhadap kecepatan dan waktu dari hasil simulasi telah diketahui, maka dihitung nilai beberapa besaran makroskopik, antara lain: 1. Fungsi Kepekatan Peluang Dengan menggunakan algoritma: Density = sum(fPV)/1000, Diperoleh Density = 1, yang sesuai dengan hasil yang diperoleh secara eksak. Keterangan: fPV = frekuensi speed partikel. 1000 = jumlah partikel. 2. Kecepatan Rata-rata Dari hasil simulasi diperoleh: Vrata = [0.2407, 0.5707, 0.0044] yang hasilnya mendekati nol, sesuai dengan hasil yang diperoleh secara eksak. Dibawah ini adalah kurva yang menunjukkan sebaran partikel terhadap masing-masing komponen kecepatan partikel.
Gambar 2 Kurva sebaran komponen kecepatan x hasil simulasi.
20
Gambar 3 Kurva sebaran komponen kecepatan y hasil simulasi.
Gambar 4 Kurva sebaran komponen kecepatan z hasil simulasi.
