Pemodelan Aliran Fluida 2-D Pada Kasus Aliran Permukaan Menggunakan Metode Beda Hingga. 2-D Fluid Surface Flow Modeling using Finite-Difference Method

Jurnal Matematika & Sains, Desember 2013, Vol. 18 Nomor 3

Pemodelan Aliran Fluida 2-D Pada Kasus Aliran Permukaan Menggunakan Metode Beda Hingga Cahyo Aji Hapsoro dan Wahyu Srigutomo Kelompok Keilmuan Fisika Bumi dan Sistem Kompleks Fakultas Matematika dan Ilmu Pengetahuan Alam Institut Teknologi Bandung, Bandung e-mail: [email protected] Diterima 8 Juni 2013, disetujui untuk dipublikasikan 18 Agustus 2013 Abstrak Persamaan dasar aliran fluida yang disebut persamaan Navier-Stokes merupakan persamaan diferensial parsial non linier yang kompleks. Untuk menyelesaikan dan memodelkan aliran fluida perlu dilakukan pendekatan numerik, salah satunya dengan metode beda hingga. Penyelesaian persamaan Navier-Stokes dilakukan dengan meninjau beberapa asumsi penyederhanaan yaitu: fluida bersifat tak termampatkan, parameter aliran bergantung pada arah spasial x dan y, serta semua variabel dianggap sebagai fungsi periodik. Pemodelan numerik dilakukan untuk menghitung ketidakstabilan Kelvin-Helmholtz lapisan campuran, evolusi struktur vorteks dan dipol vorteks. Ketidakstabilan Kelvin-Helmholtz divariasikan dengan nilai panjang gelombang gangguan sebesar  = 0.5Lx dan  = 0.25Lx. Bilangan Reynolds (Re) divariasikan dengan nilai 1000, 3000, dan 5000. Hasil pemodelan menunjukkan bahwa untuk ketiga Re tersebut, aliran fluida bersifat laminar, kritis, dan turbulen. Hal ini terindikasi oleh arah medan vortisitas serta distribusi massa jenis fluida. Semakin besar kecepatan aliran maka sifat aliran akan menjadi semakin acak atau turbulen. Kata kunci: Aliran fluida, Persamaan Navier-Stokes, Bilangan Reynolds, Vortisitas, Ketidakstabilan KelvinHelmholtz.

2-D Fluid Surface Flow Modeling using Finite-Difference Method Abstract Navier-Stokes equation is a complex non-linear partial differential equation describing a fluid flow. A numerical method called finite difference method is frequently used. In solving and modeling the fluid flow. Several simplication assumptions are incorporated in solving numerically Navier-Stokes equation which are: the fluid is incompressible, fluid flow parameters depend on its positions x and y, and all variables are considered as periodic functions. In this paper numerical calculation has been carried out to model the Kelvin-Helmholtz instability of mixed layer, evolution of vortex structure and vortex dipole. The calculation is done by varying perturbation wavelength  = 0.5Lx and  = 0.25Lx. Reynolds number (Re) is varied at 1000, 3000, and 5000. The results show that for the three values of Re, the properties of the flows are laminar, critical, and turbulent, respectively as indicated by the vorticity direction and distribution of fluid density. The larger value of fluid velocity, the more random and turbulent the fluid is. Keywords: Fluid flow, Navier-Stokes equation, Reynolds number, Vorticity, Kelvin-Helmholtz instability. bahwa perubahan momentum partikel-partikel fluida bergantung pada gaya gesekan (viskositas) yang bekerja pada fluida. Dengan kata lain persamaan Navier-Stokes menjelaskan kesetimbangan gaya-gaya yang bekerja pada fluida. Persamaan Navier-Stokes merupakan persamaan diferensial parsial non-linier yang kompleks. Keadaan ini menyebabkan penyelesaian persamaan tersebut tidak dapat dilakukan secara analitik. Oleh karena itu, perlu dilakukan pendekatan numerik untuk menyelesaikannya (Batty, 2010). Chen dan Zhang (2009) menjelaskan beberapa pendekatan numerik yang dilakukan selama ini antara lain menggunakan metode beda hingga, elemen hingga, volume hingga, Lattice Gas Automata

1. Pendahuluan Aliran fluida merupakan fenomena alam yang sering diamati di dunia ini. Aliran angin di atmosfer, aliran air di sungai, aliran lahar pada letusan gunung berapi, pergerakan awan dan asap yang membentuk suatu pola adalah beberapa contoh aliran fluida yang dapat teramati secara makroskopis. Sedangkan aliran darah dalam pembuluh merupakan aliran fluida yang teramati dalam skala mikroskopis. Dalam skala makroskopis, dinamika fluida dapat dijelaskan dengan persamaan Navier-Stokes (Guo dkk., 2006). Persamaan Navier-Stokes merupakan serangkaian persamaan terkopel yang menjelaskan pergerakan suatu fluida baik cairan maupun gas. Persamaan-persamaan ini menyatakan 81

82

Jurnal Matematika & Sains, Desember 2013, Vol. 18 Nomor 3

(LGA), Latice Bolzzmann Method (LBM), Cellular Automata (CA), dan sebagainya. Chen dan Zhang memodelkan secara numerik reservoar anisotropi, dinding horisontal, dan aliran multifase dengan memperhitungkan gaya gravitasi dan faktor lapisan. Dalam penelitiannya, Chen dan Zhang (2009) menggunakan beberapa metode, yaitu beda hingga, volume hingga, dan beda hingga campuran. Hasil penelitian tersebut dapat memodelkan beberapa parameter dalam sistem reservoar sumur. Dalam artikel ini, pendekatan numerik yang dilakukan menggunakan metode beda hingga (finite difference). Metode beda hingga adalah metode numerik yang menggunakan pendekatan solusi persamaan diferensial untuk menentukan fungsi diskrit yang memenuhi keterkaitan antara turunan fungsi dalam ruang dan/atau waktu dengan syarat batas sepanjang tepi domain (LeVeque, 2007). Keunggulan metode beda hingga dalam penyelesaian numerik persamaan Navier-Stokes adalah pada fitur formulasi tekanan dalam aliran (Rusli dkk., 2011). Dalam meninjau aliran fluida, sering digunakan konsep diskrit meskipun aliran itu sendiri dianggap sebagai sesuatu yang kontinu (Kambe, 2006). 2. Landasan Teori Prinsip mendasar dalam analisis aliran fluida termaktub dalam persamaan kontinuitas. Persamaan kontinuitas dibentuk dari prinsip kekekalan massa fluida. Jumlah aliran massa fluida atau massa per satuan waktu yang keluar dari suatu elemen volume pada aliran fluida sama dengan pengurangan massa pada volume tersebut. Tinjau elemen segiempat fluida, dengan sisi masing-masing pada sumbu x dan y adalah dx dan dy, serta ketebalan b, seperti ditunjukkan pada Gambar 1. Ketebalan b diukur secara tegak lurus terhadap bidang kertas. Kecepatan dalam arah x dan y secara berurutan adalah u dan v. Untuk arah x, massa fluida yang tersimpan dalam elemen fluida per satuan waktu dapat diperoleh dengan mengurangi laju aliran massa yang masuk dengan laju aliran massa yang keluar dari elemen fluida, seperti yang diungkapkan dalam persamaan:





  u 



x

  u  x

 dx  bdy 

(1)

bdxdy

Dengan cara yang sama, massa fluida yang tersimpan dalam elemen fluida per satuan waktu dalam arah y adalah:

   v  y

bdxdy

sebesar    bdxdy  persatuan waktu berdasarkan t

elemen fluida yang tersimpan ini. Berikutnya diperoleh persamaan sebagai berikut:



  u  x

  v

bdxdy 

y

bdxdy 

   bdxdy  t

atau

    u     v    0 t x y

y  

dy

    y

ρu

(2)

dy

u 

  u  x

dx

ρv  dx

O

2.1 Persamaan kontinuitas

 ubdy    u 

Massa elemen fluida (ρbdxdy) dapat meningkat

x

Gambar 1. Elemen segiempat fluida dengan sisi dx, dy, dan tebal b. Persamaan (1) disebut Persamaan Kontinuitas. Persamaan ini diaplikasikan untuk aliran fluida termampatkan tidak tunak. Untuk kasus aliran tunak, suku pertama dari persamaan tersebut harus sama dengan nol   t  0  . Untuk fluida tak termampatkan nilai ρ konstan, sehingga diperoleh persamaan sebagai berikut:

u v  0 x y

(3)

Persamaan (3) di atas dapat diaplikasikan untuk kedua fluida, tunak dan tidak tunak. Untuk kasus sebuah aliran simetrik aksial, maka Persamaan (3) menggunakan koordinat silinder, sehingga menjadi: u 1   rv   0 x r r

(4)

Persamaan kontinuitas merupakan persamaan yang tidak bergantung apakah fluida tersebut viskos atau tidak. Persamaan (4) dapat diaplikasikan untuk sebuah fluida ideal. 2.2 Persamaan Navier-Stokes Persamaan Navier-Stokes merupakan persamaan fluida yang didasarkan pada kekekalan massa (kontinuitas), kekekalan momentum, dan kekekalan energi. Persamaan Navier-Stokes

Hapsoro dan Srigutomo, Pemodelan Aliran Fluida 2-D Pada Kasus Aliran Permukaan Menggunakan .................. 83 diungkapkan dalam suku inersia, body force, tekanan, dan viskos seperti pada persamaan:

 u u u  u v   X   y  t x    2u  2u  p   2  2   x x y  



 v v v  Y u v    t x y  suku body  force  

(5)



p y 

suku tekanan

u

u

suku viskos

x  v

u

y

 , u

v

x  v

v

y



disebut

dengan percepatan konvektif. Pada persamaan Navier-Stokes di atas, dengan mengasumsikan suku body force sama dengan nol, dan dilakukan eliminasi suku tekanan dengan diferensiasi parsial, serta mensubstitusi persamaan v u  vortisitas   , maka diperoleh persamaan x y sebagai berikut:          u v     2  2 x y  y  t  x



2

(7)

Yilmaz (1982) menyatakan bahwa bilangan Reynolds tersebut akan mempengaruhi penurunan tekanan, sehingga menyebabkan aliran menjadi turbulen. Apabila Re<2000 maka aliran akan bersifat laminar dan akan menjadi turbulen jika Re>4000. Sedangkan apabila 2000
 2 y 2v   2  2   x y    

dengan ρ adalah kerapatan fluida, μ adalah viskositas fluida, X dan Y masing-masing adalah gaya yang bekerja pada arah sumbu x dan y, p adalah tekanan dalam fluida, serta u dan v masing-masing adalah kecepatan gerak fluida pada arah x dan y. Pada suku inersia, laju perubahan kecepatan terhadap posisi dan persamaan

 vd  d   v

3. Metode Pemodelan

suku inersia



Re 

2

  

(6)

dengan ζ merupakan vortisitas fluida. Untuk aliran ideal,  = 0, persamaan tersebut menyatakan bahwa vortisitas tidak berubah dalam proses aliran ideal. Pernyataan ini disebut dengan Teori Vorteks Helmholtz. 2.3 Bilangan Reynolds Transisi aliran fluida merupakan fenomena kompleks yang hanya dapat dikarakterisasi dengan nilai rata-rata beberapa parameter yang disebut dengan Bilangan Reynolds (Ding dan Kawahara, 1998). Dalam mekanika fluida, bilangan Reynolds merupakan perbandingan antara gaya inersia (ρvd) terhadap gaya viskositas () yang mengkuantifikasi hubungan gaya tersebut dengan suatu kondisi aliran tertentu. Bilangan ini mengidentifikasi jenis aliran yang berbeda, yaitu laminar dan turbulen. Nama dari Reynolds diambil dari Osborne Reynolds (18421912) yang mengusulkannya pada tahun 1883. Bilangan Reynolds dinyatakan dengan (Nakayama dan Boucher, 1999):

Penyelesaian persamaan Navier-Stokes dilakukan dengan meninjau beberapa asumsi di antaranya adalah fluida bersifat tak termampatkan (ρ=ρ0), parameter aliran bergantung pada variabel x dan y (2-D), serta semua variabel dianggap sebagai fungsi periodik (Danaila dkk., 2006). 3.2 Fluida tak termampatkan Matyka (2003), Borggaard dkk. (2006), serta Colella dan Puckett (1994) menyatakan bahwa untuk fluida tak termampatkan, persamaan Navier-Stokes tereduksi menjadi persamaan dengan orde rendah, dan dapat diubah ke dalam bentuk tak berdimensi:    u 1 2 (8)    u   u     u t Re Persamaan (7) dapat dijabarkan menjadi u u 2 uv p 1   2u  2u          t x y x Re  x 2 y 2    v uv v 2 p 1   2 v  2 v          t x y y Re  x 2 y 2  

(9)

3.3 Domain Komputasi, Staggered Grids, dan Syarat batas Penyelesaian persamaan Navier-Stokes secara numerik dapat disederhanakan dengan meninjau domain segiempat, yang mempunyai sisi Lx  Ly, dan dikenai syarat batas periodik di semua posisi. Titiktitik pada setiap solusi akan dikomputasikan dalam domain grid segiempat. Pertama, didefinisikan grid primer dengan mengambil titik-titik nx sepanjang x dan ny sepanjang y, seperti pada Persamaan (10) dan (11) berikut ini:

xc  i    i  1  x,  x  yc  j    j  1  y,  y 

Lx , i  1, , nx (10) nx  1 Ly ny 1

, j  1, , n y

(11)

Sistematika pendefinisian sub-domain komputasi yang digunakan dalam skema pemodelan ditampilkan pada Gambar 2.

84

Jurnal Matematika & Sains, Desember 2013, Vol. 18 Nomor 3 3.6 Kondisi awal

Ly

3.6.1 Ketidakstabilan Kelvin-Helmholtz

Y yc (j+1)

u(i,j)

ym (j)

p(i,j)

yc (j)

v(i,j) Lx X

Xc (i+1)

xc (i)

0

xm (i)

0

Gambar 2. Domain komputasi, staggered grids, dan kondisi syarat batas (digambar berdasarkan Danaila dkk., 2006). 3.4 Diskritisasi beda hingga Diskritisasi suku viskos pada persamaan Navier-Stokes dilakukan dengan metode beda hingga terpusat (Skogqvist, 1999):

H un  i , j   

u 2 uv  i, j    i, j  x y

(12)

H vn  i, j   

uv v 2  i, j    i, j  x y

(13)

Persamaan (12) menunjukkan diskritisasi suku viskositas untuk kecepatan fluida pada arah sumbu x sedangkan Persamaan (13) untuk arah sumbu y. Nilai pertambahan waktu dibatasi oleh syarat Courant-Friedrichs-Lewy (CFL) yang diungkapkan dengan: dt 

cfl  u v max    x y

(14)

  

Nilai pertambahan waktu di mana cfl < 1 disebut konstanta pengontrol pertambahan waktu. 3.5 Visualisasi aliran Visualisasi aliran dilakukan dengan menghitung medan vortisitas ω yang dikomputasikan di titik ((xc)(i), yc(j)) sebagai berikut:

  i, j  



v  i , j   v im  i  , j

x



u  i , j   u i , jm  j 





(15)

y Isokontur vortisitas yang dihasilkan dari pemodelan merepresentasikan karakteristik struktur vorteks dalam aliran.

Pada eksperimen pemodelan ini, ditentukan domain komputasi berbentuk segiempat dengan panjang sisi untuk arah x adalah Lx dan untuk arah y adalah Ly. Pada saat kondisi awal, diberikan dua buah lapisan fluida dengan massa jenis yang berbeda. Fluida dengan massa jenis yang lebih besar dilapiskan di atas fluida dengan massa jenis yang lebih kecil. Selanjutnya, pada lapisan tersebut diberikan perturbasi atau gangguan awal untuk mengetahui respon keadaan dari lapisan tersebut. Menurut Strubelj dan Tiselj (2005), ketidakstabilan Kelvin-Helmholtz merupakan salah satu ketidakstabilan sistem lapisan fluida yang mempengaruhi antarmuka lapisan. Chen dan Forbes (2011) mengungkapkan bahwa ketika perturbasi diberikan pada dua lapisan fluida terkopel yang mengalir secara paralel satu sama lain, maka akan terjadi pola ketidakstabilan pada antarmuka kedua lapisan. Hal ini terjadi karena adanya pergeseran lapisan cairan dengan perbedaan massa jenis (Danaila dkk., 2006). Simulasi pertama dilakukan dengan memberikan nilai panjang gelombang gangguan (λ) sebesar 0.5 Lx. Variasi iterasi waktu dan besarnya nilai bilangan Reynolds juga dilakukan untuk mengetahui pengaruhnya dalam perilaku medan vortisitas dan distribusi massa jenis pada aliran fluida. Variasi bilangan Reynolds dilakukan dengan nilai Re = 1000 (aliran laminar), Re = 3000 (aliran kritis), dan Re = 5000 (aliran turbulen). Simulasi kedua dilakukan dengan variasi panjang gelombang gangguan (λ) sebesar 0.25 Lx. Kemudian untuk variasi iterasi waktu dan besarnya bilangan Reynlods juga dilakukan seperti pada simulasi yang pertama. Bentuk antarmuka di antara dua cairan pada saat dikenai gangguan ditunjukkan pada Gambar 3. Pada Gambar 3.a, besarnya vortisitas dalam aliran ditunjukkan oleh kumpulan warna merah muda pada bagian atas dan biru muda pada bagian bawah. Dalam aliran tersebut, kedua bagian memiliki nilai vortisitas yang sama besar, namun arah rotasi berlawanan. Besarnya vortisitas dinyatakan dalam satuan radian per sekon. Pada Gambar 3.b menunjukkan proses terjadinya distribusi massa jenis setelah lapisan tersebut diberikan gangguan. Simulasi numerik untuk aliran tersebut dimulai dari kondisi awal yang menggunakan kecepatan awal yang sesuai dengan lapisan geser membentuk kontur jet sebagai berikut:



v  x, y   0, u  x, y   u1  y  1  u2  x 



(16)

u1 adalah kecepatan rata-rata: u1  y  

U0 2

   1  tanh  1 P 1  Ly / 2  y j  2  Rj    

    (17)    

Hapsoro dan Srigutomo, Pemodelan Aliran Fluida 2-D Pada Kasus Aliran Permukaan Menggunakan .................. 85 u 2 adalah gangguan yang menimbulkan ketidakstabilan Kelvin-Helmholtz:  x  u2  x   Ax sin  2   x  

massa jenis pada aliran fluida. Medan kecepatan dari vorteks dipol ditunjukkan pada Gambar 4.

(18)

(a)

(a)

(b) Gambar 3. Evolusi ketidakstabilan KelvinHelmholtz. Gangguan lapisan geser membentuk kontur jet (a), gulungan vortisitas Kelvin (cat eyes) (b). Kontur warna kuning dan biru muda pada Gambar 3.a menunjukkan nilai kecepatan rotasi dari vortex yang besarnya sama namun arahnya berlawanan. Kontur warna dari biru tua hingga merah tua pada Gambar 3.b menunjukkan kontras nilai densitas fluida yang semakin tinggi. Warna kuning dan biru muda pada gambar tersebut terbentuk pada batas dua fluida. Hal ini sebagai akibat dari terjadinya pencampuran dua fluida dengan densitas yang berbeda. 3.6.2 Evolusi dipol vorteks Pada eksperimen ini, akan dimodelkan pasangan vortisitas dengan arah rotasi yang saling berlawanan, yang menimbulkan dipol vorteks. Selain itu, dipol simetris untuk beberapa vorteks yang mempunyai besar vortisitas sama dalam satu aliran fluida juga akan dimodelkan. Pemodelan tersebut dilakukan dengan variasi iterasi waktu dan nilai bilangan Reynolds, yaitu dengan nilai Re = 1000 (aliran laminar), Re = 3000 (aliran kritis), dan Re = 5000 (aliran turbulen) untuk mengetahui respon medan vortisitas dan distribusi

(b) Gambar 4. Medan kecepatan. Vorteks tunggal (a), dipol vorteks, (b). Gambar 4.a menunjukkan keadaan struktur vorteks tunggal pada saat kondisi awal atau iterasi sama dengan nol. Sedangkan Gambar 4.b menunjukkan dua vorteks yang sama besar dan bergerak dengan arah yang berlawanan sehingga membentuk pasangan vorteks. Penyelesaian dipol vorteks secara numerik dilakukan dengan memasukkan medan kecepatan pada vorteks. Keadaan masing-masing vorteks ditentukan oleh keadaan pusat rotasinya, yang secara fisis didefinisikan oleh parameter-parameter yang mempengaruhinya seperti posisi (xv, yv), ukuran lv, dan intensitas 0. Parameter-parameter tersebut membangun suatu fungsi aliran yang secara analitik dinyatakan dengan 

  x, y    0 exp    

 x  xv 2   y  yv 2  lv2

 

(19)

Fungsi aliran pada persamaan di atas digunakan untuk menurunkan komponen kecepatan berikut:

86

Jurnal Matematika & Sains, Desember 2013, Vol. 18 Nomor 3

 y  yv   x , y    2   y lv2    x  xv   x , y   2 v    x lv2 

u

(20)

4. Hasil dan Pembahasan 4.1 Ketidakstabilan Kelvin-Helmholtz Dari pemodelan yang telah dilakukan, bentuk vorteks Kelvin cat eyes terjadi secara bertahap seiring dengan waktu iterasi dalam dua daerah geser aliran jet. Distribusi spasial dari jet tersebut dibentuk oleh panjang gelombang gangguan awal (λ). Pada lapisan dua fluida dengan massa jenis berbeda sesuai dengan yang telah dimodelkan, apabila diberikan gangguan dengan panjang gelombang tertentu, maka cairan dengan massa jenis lebih besar akan mulai mengalir dalam cairan yang mempunyai massa jenis lebih kecil dan akan mendorong cairan tersebut menuju tepi aliran. Pola seperti ini disebut Kelvin cat eyes. Pemodelan vortisitas yang telah dilakukan tersebut menunjukkan perilaku medan vortisitas. Secara fisis, respon yang teramati pada lapisan fluida tersebut adalah terjadinya proses pencampuran dua fluida dengan massa jenis berbeda. Medan vortisitas yang direpresentasikan oleh garis-garis panah quiver menunjukkan bahwa besarnya medan vortisitas proporsional dengan peregangan elemen-elemen garis vorteks. Simulasi secara periodik menjelaskan bahwa evolusi struktur vorteks dalam aliran jet atau lapisan geser sesuai dengan hasil eksperimen. Selain itu, pada eksperimen ini juga diamati bagaimana pengaruh vortisitas terhadap distribusi massa jenis selama proses pencampuran dua fluida. Hasil pemodelan menunjukkan bahwa distribusi massa jenis terjadi pada daerah pertemuan antara fluida dengan massa jenis yang berbeda. Pada bagian tengah fluida yang memiliki massa jenis lebih besar, distribusi massa jenis hampir tidak terjadi. Hal ini terjadi karena fluida yang memiliki viskositas lebih besar cenderung akan mempertahankan tegangan gesernya (shear stress) dari adanya perturbasi. Perturbasi atau gangguan pada lapisan fluida itu sendiri dipengaruhi gaya tekan, gaya tegang permukaan, dan gravitasi. Besarnya perturbasi sebanding dengan besarnya gaya tekan, dan berbanding terbalik adanya gaya tegang permukaan dan gaya gravitasi. Apabila gaya tekan tersebut lebih besar daripada total gaya tegang permukaan dan gaya gravitasi, maka akan menimbulkan kecepatan lapisan geser yang menyebabkan terjadinya ketidakstabilan Kelvin-Helmholtz (Kobayashi dkk., 2008). Pada variasi nilai bilangan Reynolds, terlihat bahwa nilai Re = 1000 seperti pada Gambar 5.a dan 5.b merepresentasikan sifat laminar pada aliran fluida. Hal ini ditunjukkan dengan garis-garis quiver yang cenderung berarah sejajar mengikuti arah aliran fluida. Untuk pemodelan dengan nilai Re = 3000 seperti pada Gambar 5.c dan 5.d diperoleh bahwa

garis-garis quiver mulai memperlihatkan perubahan sifat dari laminar menjadi acak. Peristiwa perubahan sifat partikel-partikel pada fluida seperti ini disebut dengan keadaan kritis. Sedangkan untuk nilai Re = 5000 seperti pada Gambar 5.e dan 5.f, diperoleh bahwa fluida bersifat turbulen. Hal ini juga ditunjukkan oleh perilaku garis-garis quiver yang semakin acak ketika diberikan bilangan Reynolds yang besar. Dengan demikian, besarnya bilangan Reynolds mempengaruhi perilaku partikel-partikel fluida yang akhirnya mempengaruhi sifat aliran fluida tersebut. Pemodelan kedua yang direpresentasikan oleh Gambar 6 menunjukkan fenomena fisis yang sama dengan Gambar 5. Untuk nilai Re = 1000 seperti pada Gambar 6.a dan 6.b, masing-masing menunjukkan perilaku medan vortisitas dan distribusi massa jenis fluida. Medan vortisitas yang dihasilkan cenderung berarah sejajar dan megikuti arah aliran. Pada Gambar 6.c dan 6.d yang diberikan nilai Re = 3000 mulai terlihat adanya sifat keacakan pada aliran yang ditunjukkan dengan arah garis-garis medan vortisitas. Sedangkan untuk nilai Re = 5000 yang dimodelkan oleh Gambar 6.e dan 6.f terlihat sifat aliran yang acak, dengan arah garis-garis medan vortisitas yang terlihat tidak beraturan pada tepi aliran. Untuk aliran dengan panjang gelombang perturbasi λ = 0.5 Lx secara umum terlihat perbedaan perilaku medan vortisitas yang signifikan untuk masing-masing nilai Re. Namun, untuk aliran dengan panjang gelombang perturbasi λ = 0.25 Lx, perbedaan perilaku medan vortisitas tidak terlalu signifikan untuk masing-masing nilai Re. Hal ini sesuai dengan eksperimen bahwa panjang gelombang perturbasi atau rasio λ/Lx yang semakin kecil akan menghasilkan jumlah vorteks yang semakin banyak dan letaknya saling berdekatan. Pada saat iterasi mencapat nilai tertentu, vorteks-vorteks tersebut dapat saling berhimpit dan menimbulkan perturbasi antar vorteks. Apabila perturbasi terus terjadi akan menyebabkan kecepatan angular menjadi semakin kecil dan menuju nol sehingga akan bergabung dengan arus aliran yang didominasi oleh kecepatan translasi. 4.2 Evolusi dipol vorteks tunggal Dinamika vorteks (vortisitas) berkaitan dengan perilaku aliran yang berotasi dalam fluida pada saat efek viskos (gesekan internal) relatif lemah, dan dapat diabaikan (Moffat, 2010). Vortisitas tersebut ditimbulkan oleh ketidakstabilan KelvinHelmholtz yang berotasi dengan arah yang sama. Pada evolusi dipol vorteks, dipol merambat secara harmonik sepanjang sumbu horizontal menuju syarat batas pada sisi kanan. Dipol vorteks ini merupakan struktur stabil yang merambat sepanjang sumbu simetrisnya dengan kecepatan translasi yang dibangkitkan oleh mekanisme induksi diri.

Hapsoro dan Srigutomo, Pemodelan Aliran Fluida 2-D Pada Kasus Aliran Permukaan Menggunakan .................. 87

(a)

(c)

(b)

(d)

(e) (f) Gambar 5. Evolusi ketidakstabilan Kelvin-Helmholtz untuk λ/Lx = 0.5, vortisitas (a); dan distribusi massa jenis untuk Re = 1000, (b); vortisitas (c) dan distribusi massa jenis untuk Re = 3000, (d), vortisitas (e) dan distribusi massa jenis untuk Re = 5000 (f).

88

Jurnal Matematika & Sains, Desember 2013, Vol. 18 Nomor 3

(a)

(b)

(c)

(d)

(e)

(f)

Gambar 6. Evolusi ketidakstabilan Kelvin-Helmholtz untuk λ/Lx = 0.25 : vortisitas (a); dan distribusi massa jenis untuk Re = 1000 (b); vortisitas (c) dan distribusi massa jenis untuk Re = 3000 (d); vortisitas (e); dan distribusi massa jenis untuk Re = 5000 (f).

Hapsoro dan Srigutomo, Pemodelan Aliran Fluida 2-D Pada Kasus Aliran Permukaan Menggunakan .................. 89

(a)

(b)

(c)

(d)

(e)

(f)

Gambar 7. Evolusi dipol vorteks tunggal, vortisitas (a) dan distribusi massa jenis untuk Re = 1000 (b), vortisitas (c) dan distribusi massa jenis untuk Re = 3000 (d), vortisitas (e) dan distribusi massa jenis untuk Re = 5000 (f).

90

Jurnal Matematika & Sains, Desember 2013, Vol. 18 Nomor 3

(a)

(b)

(c)

(d)

(e)

(f)

Gambar 8. Evolusi multi dipol vorteks identik, vortisitas (a) dan distribusi massa jenis untuk Re = 1000 (b), vortisitas (c) dan distribusi massa jenis untuk Re = 3000 (d), vortisitas (e) dan distribusi massa jenis untuk Re = 5000 (f).

Hapsoro dan Srigutomo, Pemodelan Aliran Fluida 2-D Pada Kasus Aliran Permukaan Menggunakan .................. 91

Kecepatan tersebut diinduksi oleh dipol yang memicu gerak skalar pasif, yang pada awalnya dalam keadaan diam dengan bentuk seperti jamur. Hasil pemodelan evolusi dipol vorteks tunggal ditampilkan pada Gambar 7 yang menunjukkan evolusi struktur vorteks dan pola distribusi massa jenis. Seperti pada pemodelan yang telah dilakukan untuk ketidakstabilan Kelvin-Helmholtz dengan nilai Re = 1000, medan vortisitas pada Gambar 7.a yang direpresentasikan oleh garis-garis panah quiver menunjukkan bahwa besarnya medan vortisitas proporsional dengan peregangan elemen-elemen garis vorteks. Pada Gambar 7.b menunjukkan bahwa distribusi massa jenis terjadi pada daerah pertemuan antara fluida dengan massa jenis yang berbeda. Pada bagian tengah fluida yang memiliki massa jenis lebih besar, distribusi massa jenis hampir tidak terjadi. Hal ini terjadi karena fluida dengan viskositas lebih besar cenderung akan mempertahankan tegangan geser pada fluida tersebut dari perturbasi. Kedua gambar tersebut masing-masing memperlihatkan keteraturan pola rotasi medan vortisitas dan distribusi massa jenis untuk lapisan campuran yang menunjukkan sifat aliran laminar. Pada Gambar 7.c dan 7.d dengan nilai Re = 3000 diperoleh bahwa perilaku medan vortisitas mulai bersifat acak, yang direpresentasikan oleh arah garisgaris quiver. Sedangkan untuk nilai Re = 5000 yang ditampilkan pada Gambar 7.e dan 7.d terlihat perilaku medan vortisitas yang acak dan aliran bersifat turbulen. Peristiwa ini terjadi karena bilangan Reynolds akan mempengaruhi penurunan tekanan sehingga menyebabkan aliran menjadi turbulen, sesuai dengan hasil ekperimen yang dilakukan oleh Yilmaz (1982). 4.3 Evolusi multi dipol vorteks identik Evolusi multi dipol vorteks merupakan dinamika yang tersusun lebih dari dua vorteks. Evolusi ini menunjukkan interaksi antara empat dipol pada posisi kiri dan kanan lapisan dengan massa jenis besar dan memiliki intensitas yang sama namun berlawanan arah. Empat pasang vorteks tersebut merupakan perturbasi yang mempengaruhi distribusi massa jenis lapisan campuran. Simulasi ditampilkan untuk nilai integrasi waktu yang lebih besar untuk melihat interaksi vorteks dengan lapisan campuran. Hal ini timbul sebagai akibat dari periodesitas, dipol meninggalkan domain dan masuk kembali melewati syarat batas yang berlawanan. Evolusi seperti ini disebut dengan evolusi kontinu. Pada Gambar 8.a dan 8.b dengan nilai Re = 1000 memperlihatkan keteraturan pola vorteks yang menunjukkan perilaku medan vortisitas. Peristiwa ini merepresentasikan aliran yang bersifat laminar. Untuk nilai Re = 3000, sifat medan vortisitas yang mulai terlihat acak. Seperti yang ditampilkan pada Gambar 8.c dan 8.d, yang menunjukkan bahwa aliran bersifat kritis atau mengalami transisi dari keadaan

laminar menjadi turbulen. Sedangkan untuk nilai Re = 5000 terlihat jelas perilaku medan vortisitas yang menunjukkan arah medan vortisitas yang acak dalam dan merepresentasikan sifat turbulen pada aliran tersebut. 5. Kesimpulan Ketidakstabilan Kelvin-Helmholtz terjadi apabila gaya tekan pada fluida lebih besar daripada total gaya tegang permukaan dan gaya gravitasi. Perilaku medan vortisitas dari hasil pemodelan tersebut menunjukkan bahwa besarnya medan vortisitas proporsional dengan peregangan elemenelemen garis vorteks. Distribusi massa jenis terjadi pada daerah pertemuan antara fluida dengan massa jenis yang berbeda. Pada bagian tengah fluida yang memiliki massa jenis lebih besar, distribusi massa jenis hampir tidak terjadi. Hal ini terjadi karena fluida dengan viskositas lebih besar cenderung akan mempertahankan tegangan geser pada fluida tersebut dari perturbasi. Akibat adanya tegangan geser pada aliran, vortisitas akan semakin mengecil seiring dengan waktu fluida mengalir. Variasi bilangan Reynolds juga menunjukkan adanya sifat keacakan atau turbulen pada fluida untuk bilangan Reynolds yang semakin besar. Panjang gelombang perturbasi atau rasio λ/Lx yang semakin kecil akan menghasilkan jumlah vorteks yang semakin banyak dan letaknya saling berdekatan. Pada saat iterasi mencapat nilai tertentu, vorteks-vorteks tersebut dapat saling berhimpit dan menimbulkan perturbasi antar vorteks. Apabila perturbasi terus terjadi akan menyebabkan kecepatan angular menjadi semakin kecil dan menuju nol sehingga akan bergabung dengan arus aliran yang didominasi oleh kecepatan translasi. Daftar Pustaka Batty, C., 2010. Simulating Viscous Incompressible Fluids with Embedded Boundary Finite Difference Methods. Ph.D. Thesis, University of British Columbia. Borggaard, J., A. Hay, and D. Pelletier, 2006, Interval-Based Reduced Order Models for Unsteady Fluid Flow, Int J. Num Annal Modell., 4, 353-367. Chen, M. J. and L. K. Forbes, 2011, Accurate Methods for Computing Inviscid and Viscous Kelvin-Helmholtz Instability, J Comput Phys., 230, 1499-1515. Chen, Z. and Y. Zhang, 2009, Well Flow Models for Various Numerical Methods. Int J. Num Annal Modell. 6:3, 375-388. Colella, P. and E. G. Puckett, 1994, Modern Numerical Methods of Fluid Flow. University of California. Danaila, I., J. Pascal, S. M. Kaber, and M. Postel, 2006, An Introduction to Scientific

92

Jurnal Matematika & Sains, Desember 2013, Vol. 18 Nomor 3

Computing, Twelve Computational Projects Solved with MATLAB. Springer, Paris. Ding, Y. and M. Kawahara, 1998, Linear Stability of Incompressible Fluid Flow in Cavity Using Finite Element Method, Int J Numer Meth Fl, 27, 139-157. Guo, Z., T. S. Zhao, Y. Shi, 2006, Generalized Hydrodynamics Model for Fluid Flows, From Nanoscale to Macroscale, J Phys Fluid, 18, 067107. Kambe, T., 2006. Elementary Fluid Mechanics. World Scientific Publishing, Singapore. Kobayashi, Y., M. Kato, K.T.A. Nakamura, T.K.M. Nakamura, and M. Fujimoto, 2008, The Structure of Kelvin-Helmholtz Vortices with Super-Sonic Flow, Advances in Space Research, 41, 1325-1330. LeVeque, R. J., 2007, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, SIAM, Philadelphia. Matyka, M., 2003, Solution to two-dimensional Incompressible Navier-Stokes Equations with SIMPLE, SIMPLER and VorticityStream Function Approaches, University of Wroclaw in Poland.

Moffat, H. K., 2010, A Brief Introduction to Vortex Dynamics and Turbulence, Lecture Notes Series, Institute for Mathematical Science, National University of Singapore, 21, 1-27. Nakayama, Y., and R. F. Boucher, 1999, Introduction to Fluid Mechanics. Yokendo, Tokyo. Rusli, N., E. H. Kasiman, A. K. B. Hong, A. Y. M. Yassin, and N. Amin, 2011, Numerical Computation of a Two-Dimensional NavierStokes Equation Using an Improved Finite Difference Method, Journal of Mathematics UTM, 27, 1-9. Skogqvist, P., 1999, An Adaptive Finite Difference Method for the Simulation of Combustible Flows, Thesis, Department of Numerical Analysis and Computing Science, Stockholm University. Strubelj, L. and I. Tiselj, 2005, CFD simulation of Kelvin-Helmholtz instability, International Conference Nuclear Energy for New Europe 2005, September 5-8, 2005, Bled, Slovenia, Proceedings, 002.1- 002.10. Yilmaz, T., 1982, Numerical Solution of NavierStokes Equation for Laminar Fluid Flow in Rows of Plates in Staggered Arrengement, Int. J Heat Fluid Fl, 3:4, 201-206.