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.
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 vang 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.
2. Uraian Penelitian
2.1 Persamaan pengatur
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):
\[\frac{\partial W}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = S \tag{1}\]
dimana
\[W = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix} F = \begin{bmatrix} UH \\ U^2H + 0.5gH^2 \\ UVH \end{bmatrix}\] \[G = \begin{bmatrix} VH \\ UVH \\ V^2H + 0.5gH^2 \end{bmatrix} S = \begin{bmatrix} 0 \\ gH(S_{OX} - S_{fx}) \\ gH(S_{OX} - S_{fx}) \end{bmatrix}\] (2)
\[S_{fx} = \frac{|U|U}{C^2} \quad \text{dan} \quad S_{ox} = \frac{\partial h}{\partial x}\] (3)
Persamaan (1) merupakan persamaan differensial orde ke-1 dan dapat disederhanakan menjadi:
\[\frac{\partial W}{\partial t} + \nabla . H = S \tag{4}\]
\[H = F\vec{i} + G\vec{j} \tag{5}\]
2.2 Metode numerik
Dalam metode volume hingga dengan skema cellcentered yang pertama dikembangkan oleh Jameson (1981), Persamaan (4) diintegral terhadap suatu bidang atau kontrol volume (\(\Omega\)) dan menghasilkan:
\[\frac{\partial}{\partial t} \iint_{\Omega} W \, d\Omega + \iint_{\Omega} \nabla . H \, d\Omega = \iint_{\Omega} S \, d\Omega \tag{6}\] dengan menggunakan teorema divergensi Gauss, Persamaan (6) dikonversi menjadi integral permukaan sebagai berikut:
\[\frac{\partial}{\partial t} \iint_{\Omega} W \, d\Omega + \oint_{\Gamma} H \cdot n \, d\Gamma = \iint_{\Omega} S \, d\Omega \tag{7}\]
dimana n merupakan vektor normal dari permukaan Γ dan mempunyai nilai:
\[n = \frac{\partial y}{\partial s}\vec{i} - \frac{\partial x}{\partial s}\vec{j} \tag{8}\]
Sehingga Persamaan (6) dapat dituliskan kembali dalam bentuk:
\[\frac{\partial}{\partial t} \iint_{\Omega_{i,j}} W \, d\Omega + \oint_{\Gamma_{i,j}} (F \, dy - G \, dx) = \iint_{\Omega_{i,j}} S \, d\Omega \tag{9}\]
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 \(\Omega\) yang dibagi menjadi volume sel \(\Omega\)1, \(\Omega\)2 dan \(\Omega\)3 (Gambar 1). Dengan cara tersebut Persamaan (9) dapat ditulis kembali menjadi (Natakusumah, et.al., 2004):
\[\frac{\partial}{\partial t} \iint_{\Omega} W \, d\Omega + \oint_{ABCDE} H \cdot n \, d\Gamma = \iint_{\Omega} S \, d\Omega \tag{10}\]
Gambar 1. Domain Ω
Suku kedua Persamaan (10) merupakan suku konvektif yang didapat dengan menghitung flux yang melewati domain \(\Omega\) dan ditulis sebagai:
\[\frac{\partial}{\partial t} \iint_{\Omega_{l,i}} W \, d\Omega + \oint_{\Gamma_{l,i}} (F \, dy - G \, dx) = \iint_{\Omega_{l,i}} S \, d\Omega \tag{11}\] dan
\[\oint_{\Gamma_{i,i}} (F \, dy - G \, dx) \approx \sum_{k=AB}^{DA} (F_k \Delta y_k - G_k \Delta x_k) \tag{12}\]
Flux yang melewati domain \(\Omega\) dapat dihitung dengan mudah dan hanya bergantung pada jumlah elemen disekitar domain \(\Omega\) 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 \(W_{i,j}\) dan nilainya dianggap mewakili rata-rata nilai volume sel tersebut. Nilai dari variabel W dinyatakan sebagai berikut:
\[W_{i,j} = \frac{1}{A_{i,j}} \iint_{\Omega_{i,j}} W \, d\Omega \tag{13}\]
dimana \(A_{i,j}\) adalah luas volume sel \(W_{i,j}\) dan dihitung dengan rumus trapezoid rule:
\[A_{i,j} \approx \frac{1}{2} \sum_{k=A}^{D} (x_{k+1} + x_k) (y_{k+1} + y_k)\] (14)
atau menggunakan rumus integral luas:
\[A_{i,j} \approx -\frac{1}{2} \sum_{k=A}^{D} (x_k \cdot y_{k+1} - x_{k+1} \cdot y_k)\] (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:
\[A_{i,j} \frac{\partial}{\partial t} W_{i,j} + \sum_{k=AB}^{DA} (F_k \Delta y_k - G_k \Delta x_k) = A_{i,j} \cdot S_{i,j}\] (16)
Perhitungan komponen suku konvektif pada FVM sangat bergantung pada flux yang melewati volume sel dan dapat ditulis dengan:
\[C(W_k) = \sum_{i=1}^{Np} (F_i \Delta y_i - G_i \Delta x_i)\] (17)
dimana \(N_p\) adalah jumlah sisi yang membentuk volume sel k dan \(C(W_k)\) adalah operator konvektif yang merupakan pendekatan diskrit keseimbangan flux yang melewati semua sisi volume sel k. Dengan memperkenalkan flux kecepatan \(Q_i\) untuk flux yang melewati sisi volume sel ke-i sebagai:
\[Q_i = u_i \, \Delta y_i - v_i \, \Delta x_i \tag{18}\]
maka operator konvektif untuk volume sel k dapat ditulis sebagai berikut:
\[C(W_{k}) = \sum_{i}^{Np} \begin{bmatrix} Q_{i}H_{i} \\ Q_{i}U_{i}H_{i} + \frac{1}{2}gH_{i}^{2}\Delta y_{i} \\ Q_{i}V_{i}H_{i} + \frac{1}{2}gH_{i}^{2}\Delta x_{i} \end{bmatrix}\](19)
Sehingga Persamaan (16) dapat ditulis menjadi:
\[A_k \frac{\partial}{\partial} (W_k) + C(W_k) = A_k S_k \tag{20}\] atau
\[A_{k} \frac{\partial}{\partial t} \begin{bmatrix} H_{k} \\ U_{k}H_{k} \\ V_{k}H_{k} \end{bmatrix} + \sum_{i}^{N_{p}} \begin{bmatrix} QH_{i} \\ QU_{i}H_{i} + \frac{1}{2}gH_{i}^{2}\Delta y_{i} \\ QV_{i}H_{i} + \frac{1}{2}gH_{i}^{2}\Delta x_{i} \end{bmatrix} = A_{k} \begin{bmatrix} 0 \\ gH_{k}(S_{o_{x}} - S_{f_{x_{k}}}) \\ gH_{k}(S_{o_{x}} - S_{f_{y_{k}}}) \end{bmatrix}\] \[(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:
\[A_k \frac{\partial}{\partial t} (W_k) + C(W_k) - D(W_k) = A_k S_k\] (22)
dimana \(D(W_k)\) 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:
\[\nabla^2 W_k = \sum_{i=1}^{Np} (W_i - W_k)\] (23)
dan operator Biharmonik:
\[\nabla^4 W_k = \sum_{i=1}^{Np} \left( \nabla^2 W_i - \nabla^2 W_k \right) \tag{24}\]
Operator disipasi numerik buatan \(D(W_k)\) didefinisikan sebagai berikut:
\[D(W_k) = D^2(W_k) - D^4(W_k)\] (25)
\[D^{2}(W_{k}) = \sum_{i=1}^{Np} \in_{ik} \frac{A_{ik}}{\Delta t_{ik}} (W_{i} - W_{k})\] \[D^{4}(W_{k}) = \sum_{i=1}^{Np} \in^{(2)} \frac{A_{ik}}{\Delta t_{ik}} (\nabla^{2}W_{i} - W_{k})\] (26)
dimana operator \(D^2(W_k)\) dan \(D^4(W_k)\) 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
\[\epsilon_{ik} = \epsilon^{(1)} \frac{\sum_{i=1}^{Np} |h_i - h_k|}{\sum_{i=1}^{Np} |h_i + h_k|} \tag{27}\]
\[\in^{(2)} = \max\left(0, \in^{(2)} - \in_{ik}\right) \tag{28}\]
2.4 Diskritisasi waktu model numerik
Proses diskritisasi waktu model numerik dilakukan dengan merubah Persamaan (22) menjadi persamaan diferensial biasa sebagai berikut:
\[\frac{d}{dt}(W_k) = -R(W_k) + S_k \quad k = 1, 2, 3, ..., N\] (29)
\[R(W_k) = \frac{1}{A_k} [C(W_k) - D(W_k)]\] (30)
dimana N adalah jumlah volume sel. \(R(W_k)\) adalah total residu dari kesetimbangan flux konvektif vang melewati sisi volume sel yang melingkupi volume sel k.
Perhitungan langkah waktu dilakukan menggunakan metode eksplisit Runge-Kutta bertingkat dimana harga \(W_k^n\) pada Persamaan (29) dihitung pada interval nDt dan (n+1)Dt untuk mendapatkan pendekatan harga \(W_k^{n+1}\). Dalam penelitian ini digunakan perhitungan metode Runge-Kutta dalam tiga tingkat sebagai berikut:
\[W_{k}^{(0)} = W_{k}^{n}\] \[W_{k}^{(1)} = W_{k}^{(0)} - \alpha_{1} \frac{\Delta t_{k}}{A_{k}} \left[ C\left(W_{k}^{(0)}\right) - D\left(W_{k}^{(0)}\right) \right]\] \[W_{k}^{(2)} = W_{k}^{(0)} - \alpha_{2} \frac{\Delta t_{k}}{A_{k}} \left[ C\left(W_{k}^{(1)}\right) - D\left(W_{k}^{(0)}\right) \right]\] \[W_{k}^{(3)} = W_{k}^{(0)} - \alpha_{3} \frac{\Delta t_{k}}{A_{k}} \left[ C\left(W_{k}^{(2)}\right) - D\left(W_{k}^{(0)}\right) \right]\] \[W_{k}^{n+1} = W_{k}^{(3)}\] \[(31)\] dengan koefisien \(\alpha_1\), \(\alpha_2\), dan \(\alpha_3\) sebagai berikut:
\[\alpha_1 = 0.6\] \(\alpha_2 = 0.6\) \(\alpha_3 = 1.0\)
2.5 Kondisi awal dan kondisi batas model numerik
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
\[Q_w \cdot n = 0 \tag{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) = \begin{cases} 1 & h(x, y) + z > 0 \\ 0 & otherwise \end{cases} (x, y) \in \Omega\] (33)
dimana evaluasi integral secara horizontal setiap volume sel pada \(z = \eta\)in dinyatakan dengan:
\[p(\eta_i) = \int_{\Omega} p(x, y, x) d\Omega\] (34)
Persamaan (34) menunjukkan saat \(p(\eta) = 0\) maka volume sel dalam keadaan kering (dry) sebaliknya pada \(p(\eta) = 1\) volume sel dalam keadaan basah (wet). Dengan demikian, total kedalaman H(x, y, z) pada setiap sel dinyatakan sebagai:
\[H(x, y, z) = \int_{-\infty}^{z} p(x, y, z) dz = \max(D_{\min}, 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 vang telah dibuat tersebut dapat digunakan untuke menyelesaikan permasalahan keruntuhan bendungan terutama pada kodisi dimana terjadi diskontinuitas.
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.
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.
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 2. Kedalaman aliran hasil simulasi numerik sepanjang saluran pada (a) t = 3 jam, (b) t = 6 jam, (c) t = 9 jam
Gambar 3. Desain keruntuhan bendungan sebagian (Loukili dan Soulaimani, 2007). (Sumber: Tahershamsi dan Hessaroeyeh, 2010)
Gambar 4. Domain dan mesh model keruntuhan bendungan sebagian
(a)
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)
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
3.3 Model saluran lurus dengan gangguan (triangular obstacle)
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)
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.
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
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 CA-DAM. Hasil simulasi numerik akan dibandingkan dengan hasil eksperimen lab dan hasilnya ditunjukkan oleh Gambar 11.
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
Gambar 9. Desain reservoir dan saluran pada penelitian CADAM (dalam m) (a) tampak atas dan (b) tampak samping
Gambar 10. Domain dan mesh model saluran dengan bentuk L (L-Shape) dengan sudut 90°
(d)
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
numerik dibandingkan dengan hasil eksperimen lab dan hasilnya ditunjukkan oleh Gambar 14.
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°
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.
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
Gambar 17. Kontur ketinggian pada domain Bendungan Lawe-lawe
Gambar 18. Domain dan mesh simulasi pada Bendungan Lawe-lawe
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
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.
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.