ISSN 0853-2982
Natakusumah, dkk.
Jurnal Teoretis dan Terapan Bidang Rekayasa Sipil
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode Volume Hingga-Cell Center Dantje K. Natakusumah Kelompok Keahlian Teknik Sumber Daya Air, Fakultas Teknik Sipil dan Lingkungan Institut Teknologi Bandung (ITB). Jl. Ganesa No. 10, Bandung, Indonesia. E-mail:
[email protected]
Muhammad Syahril Badri Kusuma Kelompok Keahlian Teknik Sumber Daya Air, Fakultas Teknik Sipil dan Lingkungan Institut Teknologi Bandung (ITB). Jl. Ganesa No. 10, Bandung, Indonesia. E-mail:
[email protected]
Mochamad Rizky Ramadhan Alumni Magister Teknik Sipil, Pengutamaan Teknik Sumberdaya Air Institut Teknologi Bandung (ITB). Jl. Ganesa No. 10, Bandung, Indonesia. E-mail:
[email protected] Abstrak Makalah ini menjelaskan model numerik untuk mensimulasikan fenomena banjir akibat keruntuhan bendungan. Fenomena ini perlu dikaji karena bencana keruntuhan bendungan dapat menimbulkan kerusakan dan menimbulkan korban jiwa. Model numerik yang dikembangkan adalah metode volume hingga explisit dengan skema sel terpusat yang bekerja pada suatu sistem grid yang tidak beraturan. Metode volume hingga yang digunakan ini, yang pertama kali dikembangkan untuk menyelesaikan persamaan-persamaan Euler Compresible yang selanjutnya dimodifikasi untuk menyelesaikan persamaan St. Vennant dua dimensi yang diturunkan dari persamaan hidrodinamika perairan dangkal. Pada model numerik ini, diskritisasi ruang diselesaikan menggunakan metode volume hingga dengan skema sel terpusat (cell-centered finite volume method) dan integrasi waktu diselesaikan menggunakan metode Runge-Kutta bertingkat. Untuk meredam osilasi numerik yang mungkin timbul akibat tidak adanya suku difusi digunakan disipasi numerik buatan. Sebagai kondisi batas dinding dipakai asumsi tidak ada aliran yang menembus dinding dan sebagai kondisi batas aliran digunakan metode karakteristik aliran. Studi banding hasil aplikasi model dilakukan untuk kasus-aliran akibat keruntuhan bendungan dan aliran permukaan. Perbandingan hasil dengan hasil eksperiment dan hasil model numerik lain ternyata memberikan kesesuaian hasil yang baik. Makalah ini juga membahas hasil aplikasi model untuk memprediksi rambatan aliran akibat kasus hipotetis bencana keruntuhan dari Bendungan Lawe-lawe. Kajian sepeti ini diperlukan untuk upaya mitigasi apabila terjadi keruntuhan bendungan tersebut. Kata-kata Kunci: Disipasi numerik buatan, Grid tak beraturan, Metode Runge-Kutta bertingkat, Metode volume hingga. Abstract This paper describes a numerical model to simulate the phenomenon of flooding due to dam collapse. This phenomenon needs to be studied because of catastrophic dam collapse could inflict damage and cause casualties. Numerical model developed is based on cell-centered finite volume method with explicit time step that works on an irregular grid system. The method used by this volume, which was first developed to solve the Euler equations Compressible which is further modified to solve the equation St. Vennant two-dimensional hydrodynamic equations derived from the shallow waters equations. In this numerical model, the discretization space is solved using the method of volume to cell centered scheme (cell-centered finite volume method) and the integration time is solved using the Runge-Kutta method stratified. To dampen numerical oscillations that may arise due to lack of diffusion terms was overcome using artificial numerical dissipation. For wall boundary conditions, no flow through the wall condition was implemented, and for flow boundary conditions, the method of flow characteristics was implemented. Comparative study of the results of the model application were done for case-flow due to the collapse of the dam and surface runoff. Comparison of the results with the results of experiments and other numerical model results shows a good results agreement. This paper also discusses the results of application of the model to predict the propagation of the flow due to the hypothetical case of a dam collapse disaster Lawe-lawe dam. A case study is required for mitigation in the event of collapse of Lawe-Lawe dam. Keywords: Dambreak, Finite volume, Flood propagation, Cell centered. Vol. 21 No. 1 April 2014
79
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
1. Pendahuluan Potensi Sumber daya air yang ada untuk pembangkit listrik, irigasi, air baku, pengendalian banjir, budidaya perikanan, wisata air dan konservasi. Usia rata-rata mayoritas bendungan yang ada pada saat ini telah mencapai lebih dari 40 tahun sehingga daerah tangkapan dan layananya telah mengalami perubahan tataguna lahan/lingkungan hidup yang mengakibatkan meningkatnya resiko bencana bila bendungan tersebut mengalami keruntuhan (dam break). Peningkatan resiko bencana akibat adanya peningkatan probabilitas terjadinya dambreak dan peningkatan kerawanan daerah genangan pada alur rambatan benjir pada saat terjadi dam break. Peningkatan probabilitas terjadinya dambreak terjadi selain karena adanya peningkatan beban direact runoff /banjir sebagai akibat perubahan tataguna lahan dan perubahan iklim dalam dua decade terakhir juga karena penurunan kapasitas bendungan akibat sedimentasi.
UH H 2 2 W = UH F = U H + 0.5 gH VH UVH VH 0 G= UVH S = gH S OX − S fx gH S OX − S fx V 2 H + 0.5 gH 2
( (
(2)
) )
dan
S fx =
UU C
2
dan S ox =
∂h ∂x
(3)
Persamaan (1) merupakan persamaan differensial orde ke-1 dan dapat disederhanakan menjadi: ∂W + ∇.H = S ∂t
(4)
dengan r r H = Fi + Gj
(5)
2.2 Metode numerik Peningkatan kerawanan daerah genangan terjadi karena adanya perubahan tataguna lahan daerah tersebut menjadi daerah pemukiman/lahan usaha yang lebih padat. (Kusuma et.al, 2008). Salah satu contoh dam break yang terjadi di Indonesia adalah Situ Gintung yang runtuh pada Maret 2009. Untuk mereduksi dampak akibat kegagalan bendungan ini maka dikembangkan suatu konsep metode numerik sebagai alat bantu mitigasi. Berbagai upaya pengembangan model matematik dua dimensi untuk rambatan banjir akibat dambreak aliran yang dilakukan berbasiskan metoda beda hingga memberikan hasil yang kurang memuaskan karena adanya permasalahan syarat batas geometri yang belum dapat diselesaikan secara akurat (Kusuma et.al, 2008). Dalam upaya menyelesaikan permasalahan bentuk syarat batas geometri yang tidak beraturan tersebutlah pengembangan model matematika berbasis volume hingga ini dilakukan.
Dalam metode volume hingga dengan skema cellcentered yang pertama dikembangkan oleh Jameson (1981), Persamaan (4) diintegral terhadap suatu bidang atau kontrol volume (Ω) dan menghasilkan:
∂ ∂t
∫∫W dΩ + ∫∫ ∇.H dΩ = ∫∫ S dΩ Ω
Ω
(6)
Ω
dengan menggunakan teorema divergensi Gauss, Persamaan (6) dikonversi menjadi integral permukaan sebagai berikut:
∂ ∂t
∫∫W dΩ +∫ H ⋅ n dΓ = ∫∫ S dΩ Ω
Γ
(7)
Ω
dimana n merupakan vektor normal dari permukaan Γ dan mempunyai nilai: ∂y r ∂x r (8) n= i − j ∂s ∂s
2. Uraian Penelitian
Sehingga Persamaan (6) dapat dituliskan kembali dalam bentuk:
2.1 Persamaan pengatur
∂ ∂t
Persamaan pengatur dalam penelitian ini adalah persamaan St. Vennant 2D yang didapat dari penurunan persamaan Hidrodinamika yang dirata-ratakan terhadap kedalaman atau dapat ditulis sebagai berikut (Ramadhan, 2013) : ∂W ∂F ∂G + + =S ∂t ∂x ∂y
dimana
80
Jurnal Teknik Sipil
(1)
∫∫ W dΩ + ∫ (F dy − G dx ) = ∫∫ S dΩ
Ωi , j
Γi , j
(9)
Ωi , j
2.3 Diskritisasi ruang model numerik Proses diskritisasi ruang FVM dimulai dengan membagi sebuah domain besar menjadi domain kecil dalam bentuk segitiga atau quadrilateral. Untuk ilustrasi, kita tinjau domain Ω yang dibagi menjadi volume sel Ω1, Ω2 dan Ω3 (Gambar 1). Dengan cara tersebut Persamaan (9) dapat ditulis kembali menjadi (Natakusumah, et.al., 2004):
Natakusumah, dkk.
∂ ∂t
∫∫W dΩ + ∫ Ω
H ⋅ n dΓ = ∫∫ S dΩ
ABCDE
(10)
Ω
Perhitungan komponen suku konvektif pada FVM sangat bergantung pada flux yang melewati volume sel dan dapat ditulis dengan:
C (Wk ) =
Np
∑ (Fi ∆yi − Gi ∆xi )
(17)
i =1
dimana Np adalah jumlah sisi yang membentuk volume sel k dan C(Wk) adalah operator konvektif yang merupakan pendekatan diskrit keseimbangan flux yang melewati semua sisi volume sel k. Dengan memperkenalkan flux kecepatan Qi untuk flux yang melewati sisi volume sel ke-i sebagai: Qi = ui ∆yi − vi ∆xi
Gambar 1. Domain Ω
Suku kedua Persamaan (10) merupakan suku konvektif yang didapat dengan menghitung flux yang melewati domain Ω dan ditulis sebagai: ∂ ∂t
∫∫ W dΩ + ∫ (F dy − G dx ) = ∫∫ S dΩ
Ωi , j
Γi , j
(11)
Ωi , j
dan DA
∫ (F dy − G dx ) ≈ ∑ (Fk ∆yk − Gk ∆xk )
(12)
maka operator konvektif untuk volume sel k dapat ditulis sebagai berikut:
Qi H i Np 1 C (Wk ) = ∑ QiU i H i + gH i2 ∆yi 2 i Q V H + 1 gH 2 ∆x i i i i i 2
Ak
Flux yang melewati domain Ω dapat dihitung dengan mudah dan hanya bergantung pada jumlah elemen disekitar domain Ω seperti ditunjukkan oleh Persamaan (12). Dengan demikian diskritisasi ruang FVM dapat dilakukan pada elemen dengan bentuk segitiga, segiempat, kurvalinier atau gabungan dari ketiganya. Variabel W dalam Persamaan (11) dianggap sebagai titik pada pusat volume sel Wi,j dan nilainya dianggap mewakili rata-rata nilai volume sel tersebut. Nilai dari variabel W dinyatakan sebagai berikut: Wi , j =
1 Ai , j
∫∫ W dΩ
(13)
Ωi, j
dimana Ai,j adalah luas volume sel Wi,j dan dihitung dengan rumus trapezoid rule:
Ai , j ≈
1 D (xk +1 + xk )( yk +1 + yk ) 2 k=A
∑
(14)
atau menggunakan rumus integral luas:
Ai , j ≈ −
1 D (x k ⋅ y k +1 − x k +1 ⋅ y k ) 2 k=A
∑
(15)
Besaran Ai,j dianggap tidak bergantung terhadap waktu (diasumsikan konstan). Dengan mensubtitusikan Persamaan (12) dan (13) ke dalam Persamaan (11) maka didapat diskritisasi ruang untuk FVM sebagai berikut: DA ∂ Ai, j Wi, j + (Fk ∆yk − Gk ∆xk ) = Ai, j ⋅ Si, j ∂t k=AB
∑
(16)
(19)
Sehingga Persamaan (16) dapat ditulis menjadi:
k = AB
Γi, j
(18)
∂ (Wk ) + C (Wk ) = Ak S k ∂
(20)
atau Qi Hi 0 Hk Np 1 2 ∂ Ak UkHk + QiUi Hi + gHi ∆yi = Ak gHk Soxk −Sfxk 2 ∂t i gHk Soyk −Sfyk 1 Vk Hk QiVi Hi + gHi2∆xi 2
∑
( (
) )
(21)
Persamaan (21) tidak memiliki suku viscous untuk meredam osilasi pada dirinya secara alami saat terjadi gelombang kejut. Untuk meredam osilasi yang timbul dan mendapatkan solusi yang tidak cacat dilakukan penambahan disipasi numerik buatan yang dapat ditulis sebagai berikut: Ak
∂ (W k ) + C (W k ) − D (W k ) = Ak S k ∂t
(22)
dimana D(Wk) adalah operator disipasi numerik buatan. Dalam penelitian ini digunakan operator disipasi numerik buatan yang dikembangkan oleh Jameson untuk persamaan Euler yang merupakan gabungan dari operator Laplacian: Np
∇ 2Wk = ∑ (Wi − Wk )
(23)
i =1
dan operator Biharmonik:
∇ 4Wk =
∑ (∇ 2Wi − ∇ 2Wk ) Np
(24)
i =1
Vol. 21 No. 1 April 2014
81
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
Operator disipasi numerik buatan D(Wk) didefinisikan sebagai berikut:
D(Wk ) = D 2 (Wk ) − D 4 (Wk )
∑
(26)
∑∈(2) ∆tikik (∇ 2Wi − Wk ) Np
A
i =1
dimana operator D2(Wk) dan D4(Wk) adalah operator Lapacian dan operator Biharmonik. Operator Biharmonik ditambahkan dalam domain aliran agar solusi perhitungan lebih halus. Pada saat terjadi gelombang kejut operator ini dimatikan dan digantikan oleh operator Laplacian yang dihidupkan untuk meredam osilasi. Pergantian (switching) kedua operator ini dilakukan dengan menggunakan sensor kejut berikut: Np
∈ik =∈(1)
∑ hi − hk i =1 Np
(27)
∑ hi + hk i =1
(
∈(2 ) = max 0,∈(2 ) − ∈ik
)
(28)
2.4 Diskritisasi waktu model numerik Proses diskritisasi waktu model numerik dilakukan dengan merubah Persamaan (22) menjadi persamaan diferensial biasa sebagai berikut: d (Wk ) = − R(Wk ) + S k k = 1,2,3,…,N dt
(29)
1 R(Wk ) = [C (Wk ) − D(Wk )] Ak
(30)
dimana N adalah jumlah volume sel. R(Wk) adalah total residu dari kesetimbangan flux konvektif yang melewati sisi volume sel yang melingkupi volume sel k. Perhitungan langkah waktu dilakukan dengan menggunakan metode eksplisit Runge-Kutta bertingkat dimana harga Wkn pada Persamaan (29) dihitung pada interval nDt dan (n+1)Dt untuk mendapatkan pendekatan harga Wkn+1.Dalam penelitian ini digunakan perhitungan metode Runge-Kutta dalam tiga tingkat sebagai berikut: W k(0 ) = W kn
[(
)
(
)]
[(
)
(
)]
[(
)
(
)]
W k(1 ) = W k(0 ) − α 1
∆tk C W k(0 ) − D W k(0 ) Ak
W k(2 ) = W k(0 ) − α 2
∆tk C W k(1 ) − D W k(0 ) Ak
W k(3 ) = W k(0 ) − α 3
∆tk C W k(2 ) − D W k(0 ) Ak
W kn + 1 = W k(3 )
82
Jurnal Teknik Sipil
α 2 = 0 .6
α 3 = 1.0
2.5 Kondisi awal dan kondisi batas model numerik
A D (Wk ) = ∈ik ik (Wi − Wk ) ∆t ik i =1 D 4 (Wk ) =
α1 = 0.6
(25)
Np
2
dengan koefisien α1, α2, dan α3 sebagai berikut:
(31)
Kondisi awal aliran model numerik adalah aliran diam dengan debit tetap disepanjang saluran sedangkan kondisi batas model numerik dibagi menjadi dua kondisi, yaitu kondisi batas dinding saluran (wall boundary conditions) dan kondisi batas aliran (flow boundary conditions). Syarat dari kondisi batas dinding saluran adalah tidak ada aliran yang menembus permukaan dinding atau flux kecepatan normal aliran yang melintasi permukaan dinding sama dengan nol yangdapat dinyatakan sebagai berikut Qw ⋅ n = 0
(32)
Untuk kondisi batas aliran baik pada aliran masuk (hulu) maupun aliran keluar (hilir) ditentukan dengan menggunakan metode karakteristik. 2.6 Perlakuan wet and dry Perlakuan wet and dry merupakan salah satu masalah tersulit dalam pemodelan aliran karena ada kriteria lain yang harus diperhatikan selain persamaan pengatur itu sendiri. Teknik yang digunakan dalam penelitian ini untuk memecahkan masalah tersebut adalah dengan menggunakan fungsi porositas (Casulli, 2008) berikut:
p ( x, y , z ) =
1 h ( x, y ) + z > 0 ( x, y ) ∈ Ω 0 otherwise
(33)
dimana evaluasi integral secara horizontal setiap volume sel pada z = ηin dinyatakan dengan: p (η i ) =
∫ p (x , y , x ) dΩ
(34)
Ω
Persamaan (34) menunjukkan saat p(η) = 0 maka volume sel dalam keadaan kering (dry) sebaliknya pada p(η) = 1 volume sel dalam keadaan basah (wet). Dengan demikian, total kedalaman H(x, y, z) pada setiap sel dinyatakan sebagai: z
H(x, y, z) = ∫ p(x, y, z) dz = max(Dmin, h(x, y) + z)
(35)
−∞
3. Hasil dan Analisis Algoritma numerik yang dibuat menggunakan metode volume hingga dengan skema sel terpusat ini akan diuji untuk kasus-kasus keruntuhan bendungan. Hal ini bertujuan untuk mengetahui apakah algoritma yang telah dibuat tersebut dapat digunakan untuke menyelesaikan permasalahan keruntuhan bendungan terutama pada kodisi dimana terjadi diskontinuitas.
Natakusumah, dkk.
3.1 Model saluran lurus dengan variasi kedalaman Kasus kesatu yang akan diuji adalah model saluran lurus yang memiliki variasi kedalaman dengan input pasang surut di salah satu ujungnya dan dinding di ujung lainya. Tujuan dari pengujian ini adalah untuk mengetahui kemampuan dalam mengatasi masalah kondisi wet and dry karena banyak skema numerik yang gagal menghadapi kondisi diskontinu akibat kedalaman aliran nol. Untuk uji kasus ini digunakan domain yang sama dengan uji kasus pertama tetapi saluran memiliki dasar yang miring. Dasar saluran memiliki kemiringan 0.01 dimana ujung sebelah kiri lebih rendah (ujung dengan masukan) dari ujung sebelah kanan (ujung dengan dinding). Kondisi awal untuk kasus ini adalah saluran memiliki kedalaman bervariasi sepanjang 4 km dan kedalaman air nol pada 1 km sisanya (Gambar 2). Amplitudo pasang surut yang digunakan sebesar 1.5 m dengan periode 12 jam dan simulasi dilakukan selama 48 jam dengan langkah waktu 1 detik. Kedalaman aliran hasil simulasi numerik ditunjukkan oleh Gambar 2. Hasil simulasi numerik di atas memperlihatkan Saat terjadi surut daerah kering pada saluran berada pada x = 2500 m sedangkan pada saat terjadi pasang sepanjang saluran terendam air. Hasil tersebut telah menunjukan bahwa model numerik sudah cukup baik dalam mengatasi kondisi wet and dry dengan memberikan nilai kedalaman minimum pada kondisi awal simulasi.
(a)
(b)
3.2 Model keruntuhan bendungan sebagian (partial dambreak) Kasus kedua yang akan diuji adalah runtuhnya sebagian dinding bendungan (partial dambreak) atau pintu air yang dibuka dengan cepat. Kondisi awal yang tidak kontinu memberikan kesulitan karena banyak skema mengalami kegagalan dalam kondisi ini. Penelitian sebelumnya dilakukan oleh, Loukili dan Soulaimani (2007) dan Tahershamsi dan Hessaroeyeh (2010) yang mengasumsikan sebuah kanal yang memiliki panjang 200 m dan lebar 200 m serta pemutusan bendungan yang tidak simetris dengan lebar 75 m dan tebal bendungan 10 m.
(c) Gambar 2. Kedalaman aliran hasil simulasi numerik sepanjang saluran pada (a) t = 3 jam, (b) t = 6 jam, (c) t = 9 jam
Dalam model ini digunakan kemiringan dasar saluran 0.01 dan koefisien Manning 0.025. Kondisi awal berupa kedalaman air yang berbeda antara hulu dan hilir bendungan sebesar 10 m di hulu dan 0.01 m di hilir. Domain komputasi yang digunakan dibagi menjadi 1656 titik simpul dan membentuk 3100 elemen segitiga dan simulasi perhitungan dilakukan selama 7.2 detik dengan selang waktu 0.01 detik. Gambar 3. Desain keruntuhan bendungan sebagian (Loukili dan Soulaimani, 2007). (Sumber: Tahershamsi dan Hessaroeyeh, 2010) Vol. 21 No. 1 April 2014
83
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
(b) Gambar 4. Domain dan mesh model keruntuhan bendungan sebagian
Gambar 5. Perbandingan elevasi air hasil simulasi model keruntuhan bendungan sebagian pada t = 7.2 detik (a) Keadaan dry bed dan (b) Keadaan wet bed. (Sumber: Loukili dan Soulaimani, 2007)
(a)
(a)
(b) Gambar 6. Perbandingan elevasi air dan kecepatan hasil simulasi model keruntuhan bendungan sebagian pada t = 7.2 detik (a) Keadaan dry bed dan (b) Keadaan wet bed
84
Jurnal Teknik Sipil
Natakusumah, dkk.
3.3 Model saluran lurus (triangular obstacle)
dengan
gangguan
Kasus ketiga yang akan diuji adalah model saluran yang memiliki gangguan berbentuk segitiga (Triangular Obstacle). Uji coba ini pertama kali dilakukan untuk melihat pengaruh aliran akibat kegagalan bendungan saat melewati sebuah gangguan. Penelitian saat itu pengukuran nilai tinggi muka air hanya dilakukan di beberapa titik yang diamati. Penelitian berikut (Soarez, 1999) melakukan dengan skala yang lebih kecil dibantu dengan kamera CCD berkecepatan tinggi untuk menangkap fenomena yang terjadi. Soarez menggunakan saluran sepanjang 5.6 m dan lebar 0.5 m dengan dinding yang terbuat dari kaca. Pada bagian ujung sebelah kiri terdapat reservoir sepanjang 2.39 m dengan tinggi air 0.111 m dalam keadaan diam sedangkan didepan reservoir terdapat sebuah gangguan dan saluran dalam keadaan kering. Pada bagian belakang gangguan terdapat sebuah kolam sedalam 0.025 m dan ujung sebelah kiri dibatasi dengan dinding. Kondisi awal dari penelitian Soarez ditunjukkan oleh Gambar 7.
Gambar 7. Set-up dan kondisi awal dari penelitian Soarez (dalam m). Sumber: Soarez and Zech 1999)
model ini dengan dasar saluran 0.33 m lebih tinggi dari dasar reservoir sehingga terdapat perbedaan vertikal pada pintu masuk saluran dan dipisahkan oleh penutup. Desain reservoir dan saluran pada penelitian CADAM ditunjukkan oleh Gambar 9. Kondisi awal kedalaman air pada reservoir sebesar 0.53 m dan saluran dianggap kering serta koefisien Manning yang digunakan adalah 0.0095. Domain simulasi dibuat sama dengan desain penelitian CADAM dan terdiri dari 1021 titik simpul yang membentuk 906 elemen segiempat (Gambar 10). Simulasi ini dilakukan selama 40 detik dan titik tinjau kedalaman aliran sama dengan eksperimen lab CADAM. Hasil simulasi numerik akan dibandingkan dengan hasil eksperimen lab dan hasilnya ditunjukkan oleh Gambar 11.
(a)
Pada simulasi model ini digunakan domain yang sama dengan dimensi saluran dari penelitian Soarez dan terbagi menjadi 1243 titik simpul dan membentuk 1120 elemen segiempat. Koefisien kekasaran dasar saluran sama dengan nol karena dasar dibuat dari bahan kaca. Kedalaman aliran hasil simulasi numerik ditunjukkan oleh Gambar 8. Dari perbandingan hasil simulasi numerik dan eksperimen laboratorium yang dilakukan Soarez dapat terlihat perubahan muka air saat melewati gangguan (obstacle) antara keduanya memiliki pola yang sama. Perbedaan pola perubahan muka air terjadi di bagian kolam dibelakang obstacle pada saat t = 3 detik dan t = 3.7 detik dimana simulasi numerik memberikan nilai yang lebih besar. Hal ini disebabkan kondisi batas dinding yang digunakan dalam simulasi numerik.
(b)
3.4 Model reservoir dengan saluran l-shape 90° Kasus keempat yang akan diuji adalah model reservoir dengan saluran berbentuk L (L-Shape) sudut 90o untuk melihat pengaruh bentuk saluran yang memiliki pembelokkan tajam terhadap pola aliran. Pada tahun 1997 tim CADAM melakukan eksperimen lab menggunakan (c) Vol. 21 No. 1 April 2014
85
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
(d) (a)
(e) Gambar 8. Kedalaman aliran hasil simulasi numerik sepanjang saluran pada (a) t = 1.8 detik, (b) t = 3 detik, (c) t = 3.7 detik, (d) t = 8.4 detik dan (e) t = 15.5 detik
(b)
(a)
(b) Gambar 9. Desain reservoir dan saluran pada penelitian CADAM (dalam m) (a) tampak atas dan (b) tampak samping
(c)
Gambar 10. Domain dan mesh model saluran dengan bentuk L (L-Shape) dengan sudut 90°
86
Jurnal Teknik Sipil
(d)
Natakusumah, dkk.
numerik dibandingkan dengan hasil eksperimen lab dan hasilnya ditunjukkan oleh Gambar 14.
(e) Gambar 11. Perbandingan kedalaman aliran pada titik (a) G1, (b) G2, (c) G3, (d) G4, (e) G5 dan (f) G6
Perbandingan antara hasil simulasi numerik dengan hasil eksperimen laboratorium pada titik G1 atau titik tinjau pada reservoir keduanya menunjukkan adanya kesamaan pola perubahan. Sedangkan pada titik G3 sampai G6, terdapat perbedaan nilai antara keduanya tetapi tetap memiliki pola yang sama. Perbedaan nilai terbesar terlihat pada titik G2 dan hal ini berkaitan dengan kestabilan perhitungan diawal simulasi numerik. Simulasi numerik ini juga dapat memperlihatkan dengan baik fenomena gelombang pantul akibat pengaruh pembelokkan saluran yang tajam. Gelombang pantul ini terjadi akibat air yang menabrak dinding saluran pada saat pembelokkan dan terlihat sebagai perubahan kenaikan elevasi pada titik tinjau G2, G3 dan G4. Pada t = 7 detik gelombang pantul melewati titik G4, kemudia pada t = 12 detik gelombang pantul sampai di titik G3 dan saat t = 14 detik sampai di titik G4. 3.5 Model reservoir dengan saluran L-Shape 45° Kasus kelima yang akan diuji adalah model reservoir dengan saluran berbentuk (L-Shape) sudut 45° untuk melihat efek damping pada sudut 45° saluran. Eksperimen laboratorium dilakukan pada tahun 1998 oleh UCL menggunakan sebuah reservoir berukuran panjang 2.44 m dan lebar 2.39 m yang dihubungkan dengan sebuah saluran yang memiliki belokan dengan sudut 45° dan dipisahkan oleh penutup. Desain reservoir dan saluran pada penelitian UCL ditunjukkan oleh Gambar 12. Pada eksperimen laboratorium ini tidak ada perbedaan tinggi dasar antara saluran dan reservoir. Kondisi awal kedalaman air resevoir sebesar 0.25 m dan saluran dianggap dalam keadaan kering serta koefisien Manning yang digunakan adalah 0.0095. Pada kasus ini domain simulasi sama dengan desain penelitian UCL dan terdiri dari 1099 titik simpul yang membentuk 971 elemen segiempat (Gambar 13). Simulasi dilakukan selama 40 detik dan titik tinjau kedalaman aliran sama dengan titik tinjau eksperimen lab UCL. Hasil simulasi
Perbandingan antara simulasi numerik dengan eksperimen laboratorium pada titik P1 atau titik tinjau pada reservoir keduanya sudah menunjukkan kesamaan pola perubahan. Sedangkan pada titik P2 dan P3, terdapat perbedaan nilai antara keduanya hal ini berkaitan dengan kestabilan perhitungan diawal simulasi numerik. Pada titik P4 antara hasil simulasi numerik dan eksperimen memiliki nilai dan pola yang sama hal ini menunjukkan pada titik simulasi memberikan perhitungan yang tepat. Sementara itu perbandingan hasil simulasi numerik dan eksperimen laboratorium pada titik P5 sampai P9 telah memperlihatkan adanya perbedaan nilai tetapi memiliki pola yang sama. Perbedaan ini terletak pada perbedaan fase pada perubahan kedalaman aliran. Hal ini berkaitan dengan waktu tiba aliran melewati titik-titik tersebut dimana simulasi numerik mempunyai fase yang lebih cepat. Hasil simulasi numerik ini juga dapat memperlihatkan dengan cukup baik fenomena efek damping dimana pada belokan saluran akan terjadi kenaikan kedalaman aliran dibagian luarnya. Hal ini ditunjukkan dengan adanya perbedaan tinggi antara titik P5, P6 dan P7 dimana titik P5 memiliki kedalaman aliran paling tinggi dan titik P7 memiliki kedalaman aliran yang paling rendah. 3.6 Model circular dambreak Kasus keenam yang akan diuji adalah model circular dambreak untuk melihat kemampuan solusi numeric dalam menggambarkan perambatan gelombang yang terjadi akibat perubahan muka air. Pada kasus ini simulasi numerik akan dibandingkan dengan hasil simulasi numerik menggunakan metode volume hingga dengan pendekatan Weighted Average Flux (WAF) yang dilakukan oleh Loukili dan Soulaimani (2007). Domain yang digunakan dalam melakukan simulasi ini adalah sebuah reservoir dengan ukuran panjang 40 m dan lebar 40 m seperti ditunjukkan oleh Gambar 15.
Gambar 12. Desain reservoir dan saluran pada penelitian UCL (dalam cm). Sumber: (Brufau and Garcia-Navarro, 2000)
Gambar 13. Domain dan mesh simulasi model saluran dengan bentuk L (L-Shape) dengan sudut 45° Vol. 21 No. 1 April 2014
87
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
(a)
(e)
(b)
(f)
(c)
(g)
(d)
88
Jurnal Teknik Sipil
(h)
Natakusumah, dkk.
(i) Gambar 14. Perbandingan kedalaman aliran pada titik (a) P1, (b) P2, (c) P3, (d) P4, (e) P5, (f) P6 (g) P7, (h) P8 dan (i) P9
Pada simulasi ini diasumsikan sebuah dam berbentuk lingkaran dengan jari-jari 2.5 m berada dibagian tengah domain. Kondisi awal muka air masing-masing di dalam dan luar dam adalah 2.5 m dan 0.5 m Perbandingan antara hasil simulasi numerik penelitian ini dengan simulasi penelitian yang dilakukan oleh Loukili dan Soulaimani (2007) dari model circular dambreak (Gambar 15). Berdasarlan grafik terlihat pada t = 0.4 detik gelombang kejut merambat kearah luar dam dan gelombang yang merambat kearah dalam serta mencapai bagian tengah dari dam. Pada t = 0.7 detik pemantulan gelombang rarefaction menyebabkan muka air dibagian tengah dam jatuh. Pada t = 1.4 detik pantulan gelombang rarefaction mengakibatkan muka air berada lebih rendah dari elevasi disekitar dam dan juga mulai membentuk gelombang kejut yang kedua. Setelah itu, pada t = 3.5 detik, kedua gelombang kejut masing masing merambat keluar dam dan pada t = 4.7 detik terjadi pantulan gelombang kejut dibagian tengah dan terus merambat keluar dari dam. Perbedaan hasil simulasi antara kedua model ini terletak pada kecepatan perubahan muka air dan osilasi gelombang kejut. Perbedaan perubahan air terlihat sejak t = 0.4 detik dimana tinggi muka air hasil simulasi model penelitian lebih rendah dari simulasi model WAF_Superbee. Walau terdapat perbedaan tetapi hasil simulasi numerik model dalam penelitian ini memiliki pola yang sama dengan hasil simulasi model WAF_Siperbee. Kesimpulan dari uji kasus circular dambreak ini adalah untuk menunjukkan model numerik pada penelitian ini cukup baik untuk mensimulasikan proses pembentukkan gelombang yang rumit.
Gambar 15. Domain dan mesh simulasi model circular dambreak
4. Studi Kasus Bendungan Lawe-lawe Salah satu tujuan dari penelitian ini adalah aplikasi model numerik pada keruntuhan bendungan dengan studi kasus untuk Bendungan Lawe-lawe. Simulasi numerik ini merupakan upaya mitigasi yang merupakan bagian dari perencanaan pembangunan Bendungan Lawe-lawe. Simulasi model ini dilakukan setelah upaya verifikasi model numerik terhadap hasil penelitian sebelumnya (refrensi) memberikan hasil yang cukup baik. Domain simulasi numerik digunakan ditunjukkan oleh Gambar 17. Bendungan ini direncakan akan memiliki tinggi 8 m sehingga domain yang digunakan pada simulasi ini adalah daerah rencana pembangunan bendung dengan kontur +13. Kondisi awal kedalaman air pada bendungan beragam di setiap titik di dalam domain yang kemudian dibuat agar rata pada keadaan 13 m sesuai batas domain. Sedangkan kondisi awal di depan bendung mempunyai kedalaman 0.01 m dan koefisien Manning yang digunakan sebesar 0.013. Hasil simulasi numerik memperlihatkan aliran air mencapai batas domain dalam waktu 300 detik atau 5 menit dan kedalaman aliran cenderung bertambah setelahnya. Di bagian dalam Bendungan dapat terlihat muka air perlahan mulai turun karena mengalir keluar dari Bendungan. Hasil simulasi masih kurang baik karena air cenderung menumpuk dan terus mengalir keluar dari domain simulasi.
Vol. 21 No. 1 April 2014
89
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
(a)
(d)
(b)
(e)
(c)
(f)
Gambar 16. Grafik perbandingan hasil simulasi model numerik WAF_Superbee dan model numerik penelitian pada (a) t = 0 detik, (b) t = 0.4 detik, (c) t = 0.7 detik, (d) t = 1.4 detik, (e) t = 3.5 detik dan (f) t = 4.7 detik
90
Jurnal Teknik Sipil
Natakusumah, dkk.
(a)
Gambar 17. Kontur ketinggian pada domain Bendungan Lawe-lawe
(b)
(c)
Gambar 18. Domain dan mesh simulasi pada Bendungan Lawe-lawe
(d) Vol. 21 No. 1 April 2014
91
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
(e)
(h)
(f)
(i)
(g)
(j)
Gambar 19. Kontur kedalaman aliran hasil simulasi numerik pada (a) t = 0 detik, (b) t = 60 detik,(c) t = 120 detik (d) t = 180 detik, (e) t = 240 detik,(f) t = 300 detik, (g) t = 360 detik, (h) t = 420 detik, (i) t = 480 detik, (j) t = 540 detik dan (k) t = 600 detik
92
Jurnal Teknik Sipil
Natakusumah, dkk.
5. Kesimpulan Berdasarkan hasi simulasi beberapa studi kasus yang telah dilakukan tersebut diatas dapat disimpulkan: 1. Metode numerik menggunakan metode elemen hingga dengan skema cell centered untuk diskritisasi ruang, metode Runge-Kutta tiga tingkat untuk diskritisasi waktu dan penambahan disipasi buatan memberikan hasil yang baik untuk penyelesain persamaan perairan dangkal. 2. Perhitungan secara eksplisit yang digunakan dalam penelitian ini sangat dipengaruhi oleh langkah waktu. Semakin besar langkah waktu yang digunakan semakin tidak stabil, sedangkan pada kondisi sebaliknya akan membutuhkan waktu simulasi yang lama. 3. Untuk penanganan kondisi wet and dry, tidak diperlukan aturan khusus untuk penentuan nilai minimum kondisi awal kedalaman air selama model numerik masih bisa melakukan perhitungan (tidak terjadi error). 4. Hasil simulasi menggunakan mesh berbentuk segiempat memberikan hasil yang lebih baik dibandingkan dengan mesh berbentuk segitiga dan ukuran mesh disesuaikan besar domain untuk mendapatkan hasil yang baik.
Loukili, Y., Soulaimani, A., 2007, Numerical Tracking of Shallow Water Waves by Unstructured Finite Volume WAF Approximation, International Journal for Computational Methods in Engineering Science and Mechanics, 8:1–14, 2007 Natakusumah D.K., Nuradil C, 2004, Simulasi Aliran di Perairan Dangkal dengan Menggunakan Metoda Volume Hingga pada Sistem Grid tak Beraturan, Jurnal Teknik Sipil, Volume 11 April 2004, No. 2. Ramadhan, M.R., 2013, Penerapan Metode Volume Hingga Dengan Skema Cell Centered untuk Mitigasi Bencana Banjir Akibat Keruntuhan Bendungan, Thesis Magister, Institut Teknologi Bandung, Indonesia Soarez, F.S., and Zech, Y., 1999, Effects of a Sharp Bend on Dam-break Flow, Proc., 28th IAHR Congress CD-ROM, Graz, Austria, August 1999. Tahershamsi, A., and Namin, M., 2010, Two Dimensional Modeling of Dam-break Flows, River Flow 2010 - Dittrich, Koll, Aberle & Geisenhainer (eds) - 2010 Bundesanstalt für Wasserbau ISBN 978-3-939230-00-7
6 Ucapan Terima Kasih Penulis menyampaikan ucapan terima kasih kepada Institut Teknologi Bandung sebagai sponsor penyedia dana penelitian melalui program ITB Innovation Research Program 2012.
Daftar Pustaka Brufau, P., Garcia, N.P., 2000., Two-Dimensional Dam Break Flow Simulation, International Journal for Numerical Methods in Fluids, Volume 33, Issue 1, pages 35–57, 15 May 2000. Casulli, V., 2008, A High Resolution Wetting and Drying Algorithm for Free Surfac Hydrodynamics, International Journal for Numerical Methods in Fluids. Jameson, A., Schmidt, W., and Turkel, E., 1981, Numerical Solution of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time Stepping Schemes, AIAA paper, 1981. Kusuma, M.S.B., Farid, M., and Bagus, M., 2008, Numerical Model Study of Two Dimension Flow Generated by a Dam Break, Proceeding of International Conferences on Earthquake Engineering and Disaster, ICEED. Vol. 21 No. 1 April 2014
93
Simulasi Numerik Perambatan Banjir Akibat Keruntuhan Bendungan dengan Metode...
94
Jurnal Teknik Sipil
