Bab V Prosedur Numerik
Pada bab ini, metode numerik digunakan untuk menghitung medan kecepatan, yakni dengan menghitung batas dan domain integral. Tensor tegangan tak Newton melalui persamaan Maxwell Linear dan perubahan (Evolution) bentuk batas dari setiap permukaan S (l) , l = 1, ..., K melalui persamaan kinematik. V.1
Integral Batas
Pada bab ini, kita akan membangun proses numerik dari persamaan integral batas melalui proses diskritisasi batas dan diskritisasi domain integral. Diskritisasi batas domain S yang mulus (smooth boundary) atau permukaan (interface) fluida menjadi N b segmen Si , i = 1, ..., N b sedemikian sehingga membentuk himpunan titik diskritisasi. Penghubung antar dua titik diskritisasi dinamakan elemen batas (Boundary Elements) dilambangkan dengan Sei , i = 1, ..., N b, dan titik-titik diskritisasi dinamakan titik-titik ekstrim (Extreme Points) atau node pada elemen batas. Pada daerah domain didiskritisasi menjadi N r sel integral. Hasil jumlahan kecepatan ui pada masing-masing elemen batas menghampiri nilai kecepatan pada batas S dan hampiran batas b e Se = ∪N I=1 Si . Selisih antara nilai kecepatan pada batas S dengan nilai hampiran kecepatan Sei dinamakan dengan error diskritisasi. Elemen batas dipilih sedemikian sehingga menghasilkan error dikritisasi yang minimum.
Melalui persamaan integral batas pada bab V (l) (l)
cik ui (x) −
R ∂Ω(l)
(l)
(l)
Kijk (r)ui (y)nj (y)dΩ =
R
(l) J (r)∂j ∂t τ ij (y)dS ΩR(l) ik (l) + η1 ∂S (l) Jik (r)πij (y)nj (l) (y)dS 1 η
(5.1)
dengan
x ∈ Ω; δij , 1/2δij , x ∈ S dan mulus di x; cik (x) = 0, x ∈ ∂Ω dan xdi luar domain
34 dari kondisi dinamik (5.23), diperoleh persamaan integral batas R R (l) (l) (l) (l) (l) cik ui (x) − ∂Ω(l) Kijk (r)ui (y)nj (y)dΩ = η1 Ω(l) Jik (r)∂j ∂t τ ij (y)dS ³ ´ R − ση ∂S (l) Jik (r) R11 + R12 ni (l) (y)dS (5.2) (l) dan πij = −P δij + λ(∂i uj + ∂j ui ). Tensor tegangan pada fluida Non-Newton τijd ditentukan dari persamaan Maxwell linear ³ ´ ¢ (l) ¡ (l) (l) ∂ (5.3) τ ij = η ∂i uj + ∂j ui 1 + De ∂t Pada interface awal S 0 diberi gangguan ¡ ¢ f (x1 , x2 , t) = ad 1 + ²e−ikx1
(5.4)
dengan ² sebagai parameter kecil dan tensor tegangan Non-Newton awal, kita asumsikan sebagai distribusi tegangan isotropik. τ ij (0) = Qδij
(5.5)
Dari persamaan (4.34) dan persamaan (5.1), diperoleh bentuk benang yang baru pada waktu t(l) . Nilai kecepatan u(l) dari persamaan (5.1) digunakan pada persamaan (4.34) untuk mendapatkan tensor tegangan tak Newton baru pada waktu t(l) . Iterasi ini berulang hingga mencapai energi minimum, yakni benang sudah terdeformasi menjadi satu tetesan (droplet). Algoritma 5.1 Proses iterasi ini dituliskan dalam algoritma berikut: Step 1 Menentukan bentuk permukaan benang awal dengan menyelesaikan persamaan (5.1) dengan memanfaatkan posisi awal x dan tensor tegangan awal (5.4) Step 2 Menghitung tensor tegangan Non-Newton τijd pada tl dengan menyelesaikan persamaan (5.2) untuk interface S l Step 3 Menghitung kecepatan pada ti dengan menyelesaikan persamaan (5.1) Step 4 Mengulang step 1 sampai step 3 Kecepatan untuk fluida Newton dapat ditentukan melalui kekontinuan kecepatan, yakni ud = uc , dengan algoritma sebagai berikut: Algoritma 5.2 Proses iterasi ini dituliskan dalam algoritma berikut: Step 1 Mengambil bentuk permukaan benang awal dengan menyelesaikan persamaan (5.1) dengan memanfaatkan kondisi awal (5.3) dan tensor tegangan awal (5.4) τijd pada ti dengan menyelesaikan Step 2 Menghitung tensor ´ ³ tegangan Newton (c)
(l)
(l)
persamaan τij = η ∂i uj + ∂j ui
untuk interface S (l) (ti )
35 Step 3 Menghitung kecepatan pada ti dengan menyelesaikan persamaan (5.1) Step 4 Lakukan ulang step 1 sampai step 3. Algoritma ini digambarkan sebagai berikut:
Input : Q, eps, etac, etad,sigma, E, a, L, x_{0}, y_{0}, tau11,tau12,tau21, tau22
Kecepatan pada batas {u_bts}
Kecepatan pada domain {udom}
Tensor tegangan baru melalui persamaan Maxwell Linear
Penentuan posisi baru melalui kondisi kinematik
Keluaran : u_bts, udom, tau11,tau12,tau21,tau22,x_baru
akhir
V.2
Integral Domain
Proses deformasi benang yang bergerak pada sistem koordinat Ox1 x2 ditransformasi ke sistem koordinat polar. Daerah domain didiskritisasi menjadi N r internal sel. Masing-masing internal sel berbentuk segitiga. Hal-hal ini dapat terlihat pada gambar berikut:
Tensor tegangan tak Newton pada tiap-tiap titik integrasi di daerah domain dihitung dengan menggunakan metode Gauss Legendre 7 titik. Tensor tegangan
36 untuk waktu t berikutnya diformulasikan µ (l) (l) uj (n+1)−uj (n) (l+1) ηht + τij = De (l) (l) xi (n+1)−xi (n) ¡ ¢ (l) ht + 1 − De τij
(l)
(l)
(l)
(l)
¶
ui (n+1)−ui (n)
(5.6)
xj (n+1)−xj (n)
dan (l+1) 1 τ xj (i+1)−xj (i) ij
³
=
uj (n+1)−uj (n) (ηht )2 De2 (xj (i+1)−xj (i)) xi (n+1)−xi (n) ht (1− De ) (l) + xj (i+1)−xj (i) τij
+
ui (n+1)−ui (n) xj (n+1)−xj (n)
´ (5.7)
dan disubstitusi ke integral domain Z (τij ∂j Jik ) dΩ Ω
Perhitungan numerik untuk persamaan integral batas (5.5) terbagi menjadi dua bagian: a. Elemen batas dan internal cell tidak mengandung titik asal x (Elemen atau cell regular) Pada kasus ini, jarak antara titik asal x dengan titik-titik hasil diskritisasi y lebih besar dari nol kx − yk > 0 sedemikian sehingga singularitas kernel berada di luar domain integral. Simulasi numerik untuk kasus ini menggunakan Quadratur Gauss. Pada elemen batas menggunakan Quadratur Gauss Legendre 12 titik. Pada inner domain menggunakan modifikasi Quadratur Gauss Legendre. b. Elemen batas dan internal cell mengandung titik asal x Elemen atau cell singular Pada kasus ini, jarak antara titik asal x dengan titik-titik hasil diskritisasi y sama dengan nol kx − yk = 0 sedemikian sehingga elemen batas 1 mengandung singularitas kernel Jik (r). Kernel ∂j Jik (r) = krk juga pada singular pada integral domain. Ambil sembarang internal sel pada domain, dengan titik masing-masing a = (u1 , v1 ), b = (u2 , v2 ), c = (u3 , v3 ). Titik-titik ini ditransformasi ke internal sel baru x1 = (0, 0) = (y1 , y2 ), x2 = (1, 0) = (y1 , y2 ), x3 = (0, 1) = (y1 , y2 ) sedemikian sehingga diperoleh y1 = (u3 − u1 )y1 + (u2 − u1 )y2 + u1 y2 = (v3 − v1 )y1 + (v2 − v1 )y2 + v1 dengan Jacobian à Jac =
∂x∗ ∂y1 ∂y ∗ ∂y1
∂x∗ ∂y2 ∂y ∗ ∂y2
!
(5.8) (5.9)
37 Selanjutnya, pada masing-masing titik-titik segitiga terdapat tensor tegangan A, B, C, yang masing-masing dinyatakan µ A=
¶
A11 A12 A21 A22
µ ;B =
B11 B12 B21 B22
¶
µ ;C =
C11 C12 C21 C22
¶ ;
dan mengalami proses transformasi sedenikian sehingga y1 A + y2 B + (1 − y1 − y2 )C
(5.10)
yang mana x1 , x2 , x3 menyatakan titik-titik segitiga (vertices). Dengan demikian integral domain dinyatakan dengan Z
Z
Z
1
1−u
τij ∂j Jik dΩ = Ω
τij ∂j Jik dudv u=0
(5.11)
v=1
dengan 1 Jik = 4π
Z
1
Z
µ
1−u
τij ∂j 0
v=1
ri rk −δik ln |r| + 2 |r|
¶ dudv
Misal f ((y1 , y2 ) ; (α, β)) = ln
p
(y1 − α)2 + (y2 − β)2
dengan ∂f y1 − α = p ∂y1 (y1 − α)2 + (y2 − β)2 ∂f y2 − β = p ∂y2 (y1 − α)2 + (y2 − β)2 dan (y1 − α)2 (y1 − α)2 + (y2 − β)2 (y1 − α)(y2 − β) g12 ((y1 , y2 ) ; (α, β)) = (y1 − α)2 + (y2 − β)2 (y1 − α)(y2 − β) g21 ((y1 , y2 ) ; (α, β)) = (y1 − α)2 + (y2 − β)2 (y2 − β)2 g22 ((y1 , y2 ) ; (α, β)) = (y1 − α)2 + (y2 − β)2 g11 ((y1 , y2 ) ; (α, β)) =
38 dan turunan dari fungsi g, yakni: ∂g11 ∂y1 ∂g11 ∂y2 ∂g12 ∂y1 ∂g12 ∂y2 ∂g21 ∂y1 ∂g21 ∂y2 ∂g22 ∂y1 ∂g22 ∂y2
= = = = = = = =
2(y1 − α) − (y1 − α)2 + (y2 − β)2 2(y1 − α)2 (y2 − β) ((y1 − α)2 + (y2 − β)2 )2 (y2 − β) − (y1 − α)2 + (y2 − β)2 (y1 − α) − (y1 − α)2 + (y2 − β)2 (y2 − β) − (y1 − α)2 + (y2 − β)2 (y1 − α) − (y1 − α)2 + (y2 − β)2 2(y1 − α)2 (y2 − β) ((y1 − α)2 + (y2 − β)2 )2 2(y2 − β) − (y1 − α)2 + (y2 − β)2
2(y1 − α)3 ((y1 − α)2 + (y2 − β)2 )2
2(y1 − α)2 (y2 − β) ((y1 − α)2 + (y2 − β)2 )2 2(y2 − β)2 (y1 − α) ((y1 − α)2 + (y2 − β)2 )2 2(y1 − α)2 (y2 − β) ((y1 − α)2 + (y2 − β)2 )2 2(y2 − β)2 (y1 − α) ((y1 − α)2 + (y2 − β)2 )2
2(y2 − β)3 ((y1 − α)2 + (y2 − β)2 )2
Jika (u1 , v1 ) = (α, β), maka (y1 − α)2 = ((u3 − u1 ) y1 + (u2 − u1 ) y2 )2 (y2 − β)2 = ((v3 − v1 ) y1 + (v2 − v1 ) y2 )2 Misal : u3 − u1 = a; u2 − u1 = b; v3 − v1 = c; v2 − v1 = d; Dengan demikian penentuan integral domain Z 1Z 1 Z 1 Z 1−u τ ij ∂j Jik dudv = (y1 A + y2 B + (1 − y1 − y2 )C) 0 1 0 1 Ã ! ay1 + by2 2(ay1 + by2 ) p + (ay1 + by2 )2 + (cy1 + dy2 )2 (ay1 + by2 )2 + (cy1 + dy2 )2 µ ¶ 2(ay1 + by2 )3 − (y1 A + y2 B + (1 − y1 − y2 )C) dy1 dy2 ((ay1 + by2 )2 + (cy1 + dy2 )2 )2 V.3
Integral Waktu
Pada bagian ini, kita akan menentukan tensor tegangan tak Newton dan bentuk permukaan domain (interface) S (l) untuk waktu t berikutnya. Bentuk permukaan domain S (l) baru ditentukan melalui perhitungan posisi titik-titik diskritisasi batas domain S (l) , yakni perhitungan integral waktu pada kondisi kinematik (2.24) dengan menggunakan metode skema Euler Forward. ¡ nM n ¢ nM (5.12) · n nj , n = 1, 2, ..., N b xnM j (ti+1 ) = xj (ti ) + ∆t u
39 V.4
Hasil Numerik
Telah diperoleh persamaan Stokes nonhomogen (l)
(l)
η∂jj ui − ∂i P (l) = De∂j ∂t τij
dalam Ω
dengan kondisi awal x2 (0) x2 (ad ) x1 (0) xL
= = = =
0; ¡ ¢ ad 1 + ²e−ikx1 0; L;
dan kondisi batas [|τ ij tj |] = 0 [|τ ij nj |] = −σ dxi = ui dt
µ
kondisi dinamik pada batasS ¶ 1 1 + ni kondisi dinamik pada S R1 R2 kondisi kinematik padaS
Berdasarkan persamaan integral batas pada bab 4, diperoleh Z Z 1 cik ui (x) − Kijk (r)ui (y)nj (y)dS − Jik (r)κni (y)dS λ S S Z 1 NN τij (r)∂j Jik (y)dΩ = λ Ω Melalui hasil numerik dari persamaan integral batas, proses deformasi fluida tak Newton menjadi droplet dinyatakan sebagai perubahan bentuk permukaan interface (perubahan x2 ) pada setiap iterasi waktu t Berdasarkan gambar di atas, pembentukan tetesan (droplet) terjadi pada iterasi terakhir. Tetesan (droplet) terlihat pada saat bentuk interface mancapai nilai minimum. Namun, proses numerik pada subbab (V.2) dan subbab (V.3) memiliki kekurangan, yakni adanya kernel ∂j Jik (r) pada integral domain yang mengandung singularitas r . Oleh karen itu, kita tidak memperoleh hasil yang maksimal. Perubahan bentuk yang telah diperoleh yang digambarkan sebagai berikut Pada gambar di atas nmpak jelas adanya singularitas pada titik x1 (L) = 1.
40
Profil interface fluida, dt = 0.0005 0.06
0.055
0.05
0.045
0.04
0.035
0.03
0.025
0
0.2
0.4
0.6
0.8
1
Profil interface fluida, De = 0, dt = 10−6 0.012
0.0115
0.011
0.0105
0.01
0.0095
0.009
0
0.05
0.1
0.15
0.2
0 1 2 3 4 5 6 7 8 9
41
T a b e l d a ta : P o s is i y b a r u e ta d = 0 .0 1 E = 4 Q = 0 .1 - 9 7 .8 5 1 - 5 8 .4 2 5 - 1 2 .4 7 8 - 8 9 .0 7 7 - 5 4 .0 3 8 - 1 2 .2 6 2 - 8 0 .8 3 1 - 4 9 .9 1 5 - 1 2 .0 6 3 - 7 4 .8 4 8 - 4 6 .9 2 4 - 1 1 .9 2 4 - 7 2 .6 5 0 - 4 5 .8 2 5 - 1 1 .8 8 4 - 7 3 .0 9 5 - 4 6 .0 4 7 - 1 1 .9 1 4 - 7 6 .4 0 4 - 4 7 .7 0 2 - 1 2 .0 1 8 - 8 3 .6 6 1 - 5 1 .3 3 1 - 1 2 .2 2 4 - 9 4 .1 2 4 - 5 6 .5 6 2 - 1 2 .5 1 2 -1 1 0 - 6 4 .4 0 9 - 1 2 .9 3 7 - 4 3 .9 0 5 - 3 1 .4 5 2 - 1 1 .1 5 4 7 0 .8 7 5 3 5 .4 3 7 0 .1 9 1 3 1 1 8 5 .6 5 1 0 .2 3 3 1 .4 9 8 1 3 5 .4 1 7 7 .2 0 5 1 3 .5 1 4 1 3 1 .2 7 7 5 .1 3 4 1 3 .3 1 6 1 4 8 .1 2 8 3 .5 5 8 1 3 .7 4 5 1 7 6 .6 9 9 7 .8 4 7 1 4 .5 0 9 2 1 2 .1 5 1 1 .5 5 7 1 5 .4 6 9 2 5 3 .2 6 1 3 .6 1 3 1 6 .5 7 8 2 9 3 .4 9 1 5 .6 2 4 1 .7 6 6 3 2 5 .4 2 1 7 .2 2 1 1 8 .5 0 8 3 3 8 .5 3 1 7 .8 7 6 1 8 .8 2 4 6 6 .6 3 4 4 2 .8 1 7 1 1 .5 5 3 - 1 5 .6 0 8 - 0 .7 8 0 4 2 0
Iterasi ke 1, De = 4 , t = 0.02 3
2
1
0
−1
−2
−3 0.5
0.6
0.7
0.8
0.9
1
