Metode Elemen Hingga untuk Penyelesaian Persamaan Aliran Turbulen k-ε PDF

Preview

Summary

Jurnal Natur Indonesia 6(1): 61-66 (2003) ISSN 1410-9379 Metode Elemen Hingga Aliran Turbulen k-ε 61 Metode Elemen Hingga untuk Penyelesaian Persamaan Aliran Turbulen k-εε Sigit Sutikno1, Adam Pamudji Rahardjo2 1 Laboratorium Hidraulika & Pengembangan Sumber Daya Air, Jurusan Teknik Sipil, FT, U...

Description

Accelerat ing t he world's research.

Metode Elemen Hingga untuk Penyelesaian Persamaan Aliran Turbulen k-ε Sigit Sutikno Jurnal Natur Indonesia

Cite this paper

Downloaded from Academia.edu 

Get the citation in MLA, APA, or Chicago styles

Related papers

Download a PDF Pack of t he best relat ed papers 

laporan mekt an Syarifudin Bahri

PENGGUNAAN MET ODE PN DAN SN UNT UK PERPINDAHAN KALOR RADIASI PADA REAKT OR NUKLIR T E… Alexander Agung 01-Lampiran-I Silabus-SAP KL S1 28Feb2013 Gest y Set iawat y

Jurnal Natur Indonesia 6(1): 61-66 (2003) ISSN 1410-9379

Metode Elemen Hingga Aliran Turbulen k-ε

61

Metode Elemen Hingga untuk Penyelesaian Persamaan Aliran Turbulen k-εε Sigit Sutikno1, Adam Pamudji Rahardjo2 1

Laboratorium Hidraulika & Pengembangan Sumber Daya Air, Jurusan Teknik Sipil, FT, Universitas Riau, Pekanbaru 28293 2 Laboratorium Komputasi, Jurusan Teknik Sipil, FT, Universitas Gadjah Mada, Yogyakarta, Diterima 26-02-2003

Disetujui 04-09-2003

ABSTRACT The phenomenon of turbulence flow such as erosion, water pollution, diffusion, and transport of material is affected by the rate of turbulence. In fact, to measure turbulence phenomenon accurately is difficult and expensive. A threedimensional turbulence flow has been developed in this research. Turbulence is modelled by the standard k-εε model with wall function. The numerical discretization has been done by finite element method (FEM) with PetrovGalerkin formulation. The numerical model could be used to simulate turbulence flow in domain of arbitrary shape. The verification process result shows that the developed numerical model has a good agreement with experiment data. The numerical model also able to use to simulat flow in domain of arbitrary shape. Keywords: Finite Element Method, Petrov-Galerkin, turbulence

PENDAHULUAN Penyelesaian numeris untuk pembuatan model dalam bidang rekayasa hidraulika ada beberapa macam yaitu metode karakteristik, metode beda hingga dan metode elemen hingga. metode elemen hingga (finite element method) dalam bidang hidraulika relatif masih baru. Dibanding metode-metode yang lain, metode elemen hingga lebih mudah dikembangkan menjadi alat yang akrab dengan pengguna dengan bekal pengetahuan yang relatif sedikit (user friendly). Selain itu dengan metode elemen hingga hitungan numerik dapat lebih mudah dilakukan untuk geometri yang rumit. Karakteristik aliran turbulen sangat penting dalam bidang rekayasa hidraulika. Fenomena-fenomena aliran seperti erosi, pencemaran air, dan angkutan suatu kuantitas sangat dipengaruhi oleh tingkat turbulensi aliran. Pemakaian metode elemen hingga untuk menyelesaikan persamaan aliran turbulen k-ε masih sangat jarang dilakukan. Beberapa peneliti yang telah melakukannya adalah Parparis (2000, komunikasi pribadi), Lew et al, (2000) dan Codina & Soto (1999). Mengingat masih sangat terbatasnya peneliti-peneliti yang menggunakan metode elemen hingga, penelitian ini diharapkan mampu mengeksplorasi lebih lanjut fenomena-fenomena yang terjadi jika persamaanpersamaan aliran turbulen diselesaikan dengan metode elemen hingga. Metode elemen hingga membagi daerah yang ditinjau dalam pias-pias kecil yang disebut elemen. Persamaan-persamaan yang merupakan proses fisik

diberlakukan pada elemen-elemen tersebut sehingga diperoleh rumusan dalam bentuk hubungan nilai-nilai yang dicari diantara elemen-elemen. Pemilihan bentuk elemen sesuai dengan filosofi metode elemen hingga yaitu menyederhanakan bentuk rumit daerah yang sedang ditinjau sehingga permasalahan dapat dipecahkan. Proses-proses yang terlibat di dalam metode elemen hingga adalah interpolasi, integrasi dan fungsi pembobotan. Metode elemen hingga menggunakan metode interpolasi untuk merumuskan persamaan diskret dari suatu proses fisik. Perumusan persamaan diskret ini melibatkan proses integrasi dari seluruh nilai-nilai dari domain yang ditinjau. Proses integrasi ini dibawa ke perumusan yang mengandung nilai-nilai di titik-titik tertentu, yaitu titik-titik nodal (nodes). Interpolasi dilakukan untuk mendapatkan nilai fungsi pendekat hˆ di suatu tempat dalam koordinat Cartesian dari nilai-nilai h di titik-titik sudut elemen yang bersangkutan, yaitu hi , dan dapat dituliskan sebagai berikut hˆ =

n

∑hN i

(1) i

i= 1

dengan hˆ adalah nilai fungsi yang ditaksir, hi nilai fungsi di titik nodal i dalam elemen yang ditinjau, Ni fungsi bentuk (shape function), dan n adalah jumlah titik nodal dalam satu elemen. Sedangkan interpolasi untuk turunan dari fungsi hˆ untuk arah x, y, dan z adalah  ∂ hˆ    =  ∂x 

n

∑ hˆ i =1

i

∂ N i  ∂ hˆ  ;  = ∂x  ∂y 

n

∑ hˆ i= 1

i

∂ N i  ∂ hˆ  ;  = ∂y  ∂z 

n

∑ hˆ i= 1

i

∂Ni ∂z

(2a)

62

Jurnal Natur Indonesia 6(1): 61-66 (2003)

 ∂ 2 hˆ   2 =  ∂x 

n

∑ hˆ i =1

i

∂ 2N i  ∂ 2hˆ  ; = ∂x 2  ∂y 2 

n

∑ hˆ i =1

i

∂ 2N i  ∂ 2hˆ  ; = ∂y 2  ∂z 2 

Sutikno & Pamudji.

n

∑ hˆ i =1

i

∂ 2N i ∂z 2

(2b)

Pada kasus dimana suatu fungsi berubah terhadap waktu, maka interpolasi untuk fungsi tersebut adalah n  ∂hˆ   ∂h   Ni   = ∑  ∂ t   i =1  ∂ t 

(3)

Fungsi bentuk (shape function). Fungsi bentuk (shape function) atau fungsi dasar (basic function) suatu titik nodal dalam interpolasi pada metode elemen hingga mempunyai sifat khusus, yaitu mempunyai nilai satu pada titik nodal tersebut dan mempunyai nilai nol pada titik nodal yang lain dalam elemen yang sama. Untuk menyederhanakan bentuk persamaan fungsi dasar maka dipakai sistem koordinat lokal, dimana tiap elemen untuk tiap-tiap arah masingmasing mempunyai nilai posisi antara –1 dan 1. Bentuk persamaan fungsi bentuk dalam sumbu koordinat lokal memberikan keuntungan dalam proses integrasi secara numerik karena hitungannya jauh lebih sederhana. Metode Sisa Berbobot. Proses penaksiran atau pendekatan suatu nilai fungsi dengan menggunakan teknik interpolasi seperti diuraikan di atas memberikan hasil penaksiran yang berbeda dengan penyelesaian eksak. Beda tersebut disebut juga sebagai kesalahan (error) atau sisa (residu) R. Fungsi kesalahan atau sisa R dinyatakan dalam bentuk R(x,y,z) = h(x,y,z) - hˆ (x,y,z) (4) dengan h(x, y, z ) adalah fungsi eksak dan hˆ x, y, z adalah fungsi pendekat. Bermula dari ide meminimumkan kesalahan tersebut secara keseluruhan dalam daerah yang dihitung, Metode Sisa Berbobot (weighted residual method) membentuk suatu formulasi dengan membuat integrasi perkalian antara fungsi kesalahan dan suatu fungsi pembobot pada seluruh domain hitungan sehingga sama dengan nol W(x,y,z)R(x,y,z)dxdydz = 0 (5)

∫

Ω

dengan W(x,y,z) adalah fungsi pembobot. Berdasarkan pemilihan fungsi pembobot yang dipakai dalam metode Elemen Hingga, dikenal dua metode yaitu metode Bubnov-Galerkin atau metode Galerkin standar dan metode Petrov-Galerkin. Pada metode Bubnov-Galerkin digunakan fungsi pembobot W yang sama dengan fungsi dasar N (basis function) yang digunakan dalam proses interpolasi (W = N). Pada penelitian ini dipakai fungsi pembobot berdasarkan formulasi Petrov-Galerkin yang diberikan dalam persamaan berikut ini

α h U m ∂N k ~ W k = Nk + α W k = Nk + 2 U ∂x m

(6)

dengan indeks m = 1,2,3 untuk arah x,y,z; dan indeks k sesuai dengan urutan nomor titik nodal dalam elemen dan α adalah koefisien upwinding. Transformasi Koordinat. Proses integrasi pada Metode Sisa Berbobot jika dilakukan dalam koordinat global akan sangat rumit bila dibandingkan integrasi dalam koordinat lokal. Untuk itu maka diperlukan transformasi fungsi diskret dari koordinat global x, y, z ke dalam koordinat lokal ξ, η, dan ζ yang masing-masing berkorespondensi satu-satu dengan x, y, dan z seperti ditunjukkan pada Gambar 1.

+1

+1 0

transformasi z

0 -1

y -1 -1

x

+1

0

(a)

(b)

Gambar 1. Transformasi koordinat global (a) ke koordinat lokal (b).

Jika proses transformasi ditulis secara keseluruhan dengan menggunakan koordinat lokal yang telah dinormalisasi, maka dapat dituliskan sebagai berikut

∫ G dV ≈ v

+1+1+1

∫ ∫ ∫ Gˆ (ξ, η, ζ )d ξ d η d ζ

(7)

−1 −1 −1

dengan fungsi G tergantung dari fungsi bentuk N atau turunannya dalam koordinat global dan G adalah fungsi hasil transformasi dari G pada koordinat lokal dikalikan dengan determinan matriks Jacobian J. Integrasi Numeris. Integrasi secara numeris dalam penelitian ini menggunakan metode GaussLegendre quadrature, yaitu metode integrasi numeris yang memanfaatkan titik-titik Gauss (Gauss points) yang masing-masing telah mempunyai nilai posisi dalam koordinat lokal dan faktor bobot (weighting factor) tertentu. Apabila suatu fungsi yang didekati pada koordinat lokal telah diketahui maka proses integrasi dengan metode Gauss-Legendre adalah sebagai berikut +1 +1 +1

∫ ∫ ∫ g( ξ, η, ζ ) d ξ d η d ζ ≈ ∑ ∑ ∑ w w w g(ξ , η , ζ ) NGP 1 NGP 2 NGP 3

i

−1 −1 −1

i =1

j =1

j

k

i

j

k

(8)

k =1

Nilai NGP1, NGP3, dan NGP3 pada Persamaan (8) masing-masing adalah jumlah titik Gauss pada arah ξ,

Metode Elemen Hingga Aliran Turbulen k-ε

η, dan ζ. Pada penelitian ini ditetapkan jumlah titik Gauss yang sama untuk tiap arah pada koordinat lokal, yaitu 3 titik, sehingga dalam satu elemen terdapat 3x3x3 atau 27 titik Gauss. Penetapan jumlah titik Gauss tersebut berdasarkan pada biaya komputasi yang relatif rendah dan tingkat ketelitian yang cukup tinggi (Sadtopo 2001). Pada jumlah titik Gauss yang lebih dari 3 titik untuk masing-masing arahnya, perbedaan akurasinya dengan jumlah titik Gauss 3 titik relatif sangat kecil, sehingga jumlah 3 titik Gauss untuk masing-masing arah adalah kondisi yang paling optimal. Persamaan Model Aliran Turbulen. Persamaan model aliran turbulen terdiri atas persamaan-persamaan kontinuitas, momentum dan persamaan k-ε. Persamaan-persamaan tersebut berturut-turut adalah (Lew et al, 2000), ∇ .u = 0

(10)

∂k + u.∇ k − ∇.( D k ∇ k ) + ε = Fk ∂t

(11a)

ε2 ∂ε + u .∇ ε − ∇ .( D k ∇ ε ) + c 2 = Fε ∂t k

(11b)

Dk = Fk =

υt +ν σk νt ∇u + ∇uT 2

νt = cµ

k2 ε

Pemakaian Model Fungsi Dinding (Wall Function Model) dapat mengatasi permasalahan tersebut. Aplikasi model ini selain grid hitungan di dekat dinding solid tidak harus rapat, juga akurasi hasil masih cukup baik. Pada penelitian ini kondisi batas pada dinding solid dipakai model Fungsi Dinding (Wall Function Model) jenis dua lapis ( a Two-Layer Model). Jika suatu titik berada pada daerah Fully Turbulent maka pada daerah tersebut berlaku hubungan berikut, U p+ =

Dε = 2

υt +ν σε

dan

Fε =

(11c)

(

2

Kp =

u∗ Cµ

2

(11d)

dengan Du = ν +νt, adalah koefisien difusi, u komponen kecepatan, p tekanan statis sesaat, ν viskositas molekuler, νt viskositas turbulen, F gaya medan (body forces), k energi kinetik turbulen, dan ε adalah kehilangan energi turbulen (dissipation of turbulence energy). Nilai-nilai konstanta persamaan tersebut yang direkomendasikan oleh Launder & Spalding (Rodi 1993) adalah c1ε1.44, C2ε 1.93, σk 1.0, σε 1.3 dan cµ 0,09. Kondisi batas pada dinding solid. Struktur turbulen pada daerah dekat dinding solid merupakan fenomena yang sangat kompleks, karena pada daerah tersebut variasi nilai properti turbulen relatif besar. Salah satu cara untuk melakukan penanganan pada daerah tersebut adalah dengan memakai model k-ε untuk angka Reynolds rendah (Low-Reynolds-Number k-ε) (Peng et al, 1997). Pemakaian Low-Reynolds-Number k-ε, walau memberikan akurasi yang baik namun model ini terbatas penggunaannya pada angka Reynolds yang rendah dan grid hitungan pada daerah dekat dinding harus cukup rapat, sehingga model ini pada kasus tertentu tidak ekonomis lagi.

εp =

u∗ κ y

3

(13)

p

Energi kinetik turbulen dan kehilangan energi kinetik turbulen pada titik P dihitung dengan menggunakan persamaan-persamaan sebagai berikut, Kp =

(11e)

(12)

Selanjutnya jika suatu titik berada pada daerah laminer sublayer maka pada daerah tersebut berlaku hubungan sebagai berikut, (14) U p+ = y p+

2

c 1k ∇u + ∇uT 2

)

1 ln y p+ E k

dengan, k=0.42, U p + =U p /u * , y p + =y p .u */½, dan u * kecepatan gesek. Energi kinetik turbulen dan kehilangan energi kinetik turbulen dihitung dengan menggunakan persamaan-persamaan sebagai berikut,

(9)

∂u 1 + u.∇ u = − ∇ p + ∇.(D u ∇ u ) + F ∂t ρ

63

u* Cµ

 yp   yt

  

3

εp =

2u ∗ C 0µ . 25 y 1

(15)

Kondisi batas pada muka air. Kondisi batas untuk permukaan air yang berhubungan langsung dengan udara adalah,

[

Cf kf Cµ εf = κYf

]

3 2

(16)

dengan εf dan kf (=u*3/Cm0.5), masing-masing adalah energi disipasi dan energi kinetik turbulen, Yf jarak terdekat terhadap muka air, dan C f adalah 0.164 (Kironoto 1992).

METODE Formulasi Model Elemen Hingga. Agar persamaan k-ε menjadi linier, maka diperkenalkan suatu koefisien, yaitu koefisien reaksi (reaction coefficients) (Lew et al, 2000) g³k=e/k sehingga Persamaan (11a) dan (11b) menjadi persamaan yang identik dan dapat ditulis dalam bentuk umun berikut, ∂Φ + u.∇ Φ − ∇ (D Φ ∇ Φ ) + γ Φ Φ = S Φ ∂t

(17)

Selanjutnya bentuk persamaan diskret dari Persamaan (17) adalah sebagai berikut, R=

ˆ ∂Φ ˆ − ∇ (D Φ ∇ Φ ˆ ) + γΦΦ ˆ = SΦ + u.∇ Φ ∂t

(18)

Jurnal Natur Indonesia 6(1): 61-66 (2003)

64

Sutikno & Pamudji.

Bentuk persamaan tersebut bisa dijabarkan lagi dalam bentuk persamaan interpolasi. Interpolasi tersebut adalah usaha untuk mendapatkan Φ ˆ disuatu tempat dalam koordinat (x,y,z) dari nilai-nilai Φ ˆ di titiktitik sudut elemen yang bersangkutan, sehingga persamaannya menjadi ˆ  ∂Φ ∂ 2Nl ˆ ∂N ˆ ˆ lNt − SΦ  + uilNt k Φ Φl + γ Φ Φ R = Nl  k − ΓΦ ∂xi2 ∂xi  ∂t l

(19)

Bila metode sisa berbobot diaplikasikan, dengan integrasi perkalian antara fungsi kesalahan dan suatu fungsi pembobot adalah sama dengan nol, maka diperoleh persamaan berikut ini,   ∂Φ ˆ  ∂N ˆ ∂ Nk ˆ ˆ lNl − SΦ dΩ = 0  + uilNl k Φ Φl + γ Φ Φ ∫ Wj Nl  k − ΓΦ ∂x i ∂x i2   ∂t l  2

(20)

dengan [L] dan [N] adalah suatu matriks bentukan baru. Dalam persamaan matriks tersebut, hanya vektor {Φˆ n+1 } yang tidak diketahui, sehingga persamaan tersebut bisa diselesaikan dengan penyelesaian persamaan matrik biasa. Penyusunan Algoritma dan program komputer. Penyusunan algoritma program komputer dilakukan dalam bahasa Fortran dan mengikuti bagan alur seperti ditunjukkan pada gambar berikut ini. Mulai Baca data-data input Tetapkan kondisi awal, Hitung fungsi bentuk Ni dan turunannya

Ω

Pemakaian fungsi pembobot pada masingmasing suku tidak sama. Fungsi pembobot berdasarkan metode Petrov-Galerkin hanya diberlakukan pada suku-suku konveksi dan difusi saja, sedangkan pada suku-suku yang lain dipakai fungsi pembobot berdasarkan metode standar Galerkin. Selanjutnya setelah dilakukan integrasi pada persamaan tersebut, maka diperoleh formulasi berikut ini.  ˆ

− Φˆ n   + [K ] Φˆ + γ Φ [M ] Φˆ − {S Φ } ∆t 

n +1

[M ] Φ 

{}

{}

(21)

Matriks [M] disebut sebagai matriks massa (mass matrix). Nilai matriks [M] ini tidak berubah dalam tiap langkah waktu, sehingga matriks [M] dihitung hanya pada langkah waktu pertama saja untuk kemudian digunakan pada langkah-langkah waktu selanjutnya. Sedangkan matriks [K] adalah matriks kekakuan (stiffness matrix) untuk suku konveksi dan suku difusi. Matriks ini mengandung interpolasi dari fungsi kecepatan u sehingga komponennya tergantung dari nilai u. Berhubung nilai kecepatan selalu berubah terhadap waktu, maka matriks kekakuan ini dihitung pada setiap langkah waktu pula. Persamaan (21) bisa diekspresikan dalam bentuk deferensi hingga (finite difference) sebagai berikut (Desai 1979), − Φˆ n  n +1 + (1 − θ ) Φˆ n  + [K ] θ Φˆ ∆t   γ Φ [M ] θ Φˆ n + 1 + ( 1 − θ ) Φˆ n = {S Φ }  ˆ

[M ] Φ

n+1

{{

{{

}

{ }}

}

{ }}+

(22)

dengan θ adalah suatu besaran skalar. Jika θ = 1, maka dipakai skema implisit, jika θ = 0.5 dipakai skema CrankNicholson dan jika θ = 0 skemanya adalah eksplisit. Dengan memisahkan variabel berdasarkan langkah waktu, Persamaan (22) dapat ditulis lagi dalam bentuk persamaan baru sebagai berikut ini. [L] Φˆ n+1 + [N] Φˆ n = {SΦ }∆t (23)

{ }

{ }

Integrasi numeris dan merakit mass matrix global M

Integrasi numeris dan merakit stiffness matrix global K. Selesaikan Persamaan NAVIER-STOKES Pembatasan Viskositas Turbulen Hitung Energi Kinetik turbulen

Pembatasan Koefisien Hitung Energi Dissipasi Hitung Viskositas Turbulen

Tidak

t=Nt ?

Ya Selesai

Gambar 2. Bagan alur penyelesaian persamaan aliran turbulen.

HASIL DAN PEMBAHASAN Verifikasi Model Numerik. Untuk menguji model numerik, dilakukan perbandingan antara hasil hitungan model numerik dengan hasil eksperimental dan dengan model-model numerik yang lain. Hasil model numerik lain yang dipakai untuk perbandingan adalah hasil dari model elemen hingga persamaan Navier-Stokes (Sadtopo 2001) dimana pada penyelesaian persamaan ini parameter viskositas turbulen belum dimodelkan. Data pembanding keluaran numerik yang lain adalah data dari keluaran numerik yang dilakukan oleh Ghia et al, (Sadtopo 2001). Sedangkan data pengukuran yang dipakai adalah data pengukuran yang dilakukan oleh Koseff & Street (Sadtopo 2001) yaitu berupa data profil kecepatan pada aliran dalam kotak untuk angka Reynolds Re = 3200. Model numerik yang telah dibuat dicoba diaplikasikan pada problem lid-driven cavity flow. Pada kasus ini problem diselesaikan pada kondisi aliran

Metode Elemen Hingga Aliran Turbulen k-ε

permanen, yaitu dengan cara melakukan pengulangan hitungan numerik sampai pada langkah waktu yang banyak. Aliran pada kotak ini mempunyai kondisi batas kecepatan konstan pada sisi atas geometrinya, sedangkan sisi-sisi lainnya berupa dinding. Pada dinding diterapkan no-slip boundary conditions, dengan vektor kecepatan u, v, dan w (arah x, y, dan z) masingmasing adalah nol. Domain pada kasus ini yaitu Ω = {(x,y) | 0 ≤ x ≤ 1, ≤ z ≤ 1}, dengan jumlah elemen yang digunakan yaitu 12 x 12 x 1. Pada z = 1 vektor kecepatan u = 1, sedangkan v = w = 0. Kondisi batas tekanan p = 0 diterapkan pada titik tengah dari sisi bawah geometri. Sedangkan untuk kuantitas turbulennya, pada dinding-dinding solid diberlakukan kondisi batas untuk solid wall boundary condition yaitu dengan mengaplikasikan wall function. Pada sisi atas geometri, yaitu pada z = 1 diberlakukan kondisi batas tipe moving wall boundary, dimana prosedur hitungannya sama dengan hitungan pada solid wall boundary condition, bedanya terletak pada besarn...