BAB III DESKRIPSI MODEL

BAB III DESKRIPSI MODEL

III.1 Konsep Model Dan Pendekatan Model yang dikembangkan merupakan model rambatan tsunami di perairan dangkal hingga ke darat. Berdasarkan definisi tsunami oleh Hamzah (2006), daerah ini dibatasi hingga kedalaman laut sekitar 50 meter. Parameter yang mempengaruhi rambatan gelombang adalah:

•

Tinggi, periode, panjang gelombang dan arah gelombang datang

•

Kontur dan bathimetri

•

Penutup lahan

•

Kedalaman/muka air normal

Pemodelan rambatan gelombang tsunami dilakukan dengan menggunakan dua persamaan pengatur, St.Venant dan Boussinesq bentuk standar. Batasan wet/dry hanya terdapat pada model St.Venant. Batasan wet/dry dimodelkan dengan memberikan nilai batas minimum kedalaman. Metode ini dipilih karena sederhana dan mudah diterapkan serta memberikan hasil yang baik pada pemodelan kasus dambreak oleh Tawatchai Tingsanchal (1999). Pada persamaan Boussinesq, penerapan batasan wet/dry ini tidak dapat diterapkan. Batasan wet/dry untuk persamaan Boussinesq dapat dilakukan dengan melakukan ekstrapolasi dari titiktitik basah. Akan tetapi, metode ini sulit diterapkan untuk kontur yang kompleks. Batasan wet/dry lain yang dapat digunakan untuk persamaan Boussinesq memerlukan modifikasi dari persamaan dan sulit dilakukan. Oleh karenanya, model run up hanya akan dimodelkan dengan persamaan St. Venant. Penyelesaian numerik untuk pemodelan rambatan tsunami seperti yang dapat dilihat pada flowchart berikut ini.

III-1

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

Gambar III-1. Flow Chart Penyelesaian Numerik

III-2

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

III.2 Input Data III.2.1 Topografi dan Bathimetri Kondisi topografi dan bathimetri dimodelkan dalam 1 layer (Z1). Untuk memberikan kedalaman air normal, maka diberikan layer kedua (Z2) yang memuat elevasi air normal dan topografi. Kedua layer tersebut dioverlay dengan ketentuan elevasi air kondisi awal adalah hasil dari layer Z1-Z2. Dengan demikian, maka pada nilai Z1-Z2 di darat akan bernilai 0, sedangkan di laut akan bernilai > 0.

III.2.2 Gelombang Datang Gelombang datang dimasukkan sebagai syarat batas pada pinggir domain laut. Input gelombang berupa perubahan elevasi muka air terhadap waktu. Diskritisasi waktu tidak terikat dengan interval waktu simulasi. Pada program, digunakan metode interpolasi linear sehingga jika diskritisasi gelombang datang memiliki interval waktu yang tidak sama dengan interval waktu simulasi, maka elevasi pada waktu yang diinginkan akan diinterpolasi secara otomatis.

III.2.3 Koefisien Manning Koefisien manning dimasukkan sebagai suatu nilai tetap bergantung pada kondisi penutup lahan yang dimodelkan. Untuk bathimetri laut, nilai koefisien manning ditetapkan 0.

III.3 Persamaan Pengatur III.3.1 St. Venant Pemodelan dilakukan dengan menggunakan persamaan pengatur st.venant dalam bentuk u,v,h. ∂h ∂ (uh ) ∂ (vh ) + + = 0 ...( III-1 ) ∂t ∂x ∂y

dimana :

III-3

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

u dan v = kecepatan arah x dan y h

= kedalaman air

Persamaan momentum 2D adalah : ∂ (h + z) ∂u ∂ ∂ + u u + v u + gh = − ghS fx ...( III-2 ) ∂t ∂x ∂y ∂x ∂ (h + z) ∂v ∂ ∂ + u v + v v + gh = − ghS fy ...( III-3 ) ∂t ∂x ∂y ∂y

dimana : Sfx, Sfy = kemiringan arah x dan y, dihitung dengan persamaan: Sfx

= n2 U x (U2+V2)1/2 / h4/3 ...( III-4 )

Sfx

= n2 V x (U2+V2)1/2 / h4/3 ...( III-5 )

Dengan n adalah koefisien manning.

III.3.2 Boussinesq Persamaan boussinesq yang akan digunakan adalah persamaan boussinesq oleh Nwogu (1993) bentuk standar.

⎧⎪⎛ zα 2 h 2 ⎞ ⎫⎪ h⎞ ⎛ − ⎟ h∇ ( ∇.u ) + ⎜ zα + ⎟ h∇ ( ∇.(hu ) ) ⎬ = 0 6 ⎠ 2⎠ ⎝ ⎪⎩⎝ 2 ⎪⎭

ηt + ∇. ⎡⎣( h + η ) u ⎤⎦ + ∇ ⎨⎜

...(

III-6 )

⎧z ⎫ ut + g ∇η + ( u.∇ ) u + zα ⎨ α ∇ ( ∇.ut ) + ∇ ( ∇.(hut ) ) ⎬ = 0 ⎩2 ⎭

...( III-7 )

Dimana: u = kecepatan (u,v) η = elevasi muka air h = kedalaman

∇ = (∂ / ∂x, ∂ / ∂y )

III-4

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

g = gravitasi Bentuk standar diperoleh dengan memasukkan nilai α = -1/3 (Zα/h = -1/2), yang memberikan nilai a1=0. a2=0, b1=1/6 dan b2 = -1/2

III.4 Skema Numerik III.4.1 Mc Cormack Persamaan pengatur St. Venant dikepingkan dengan metode Mc Cormack. Pada skema ini, nilai untuk time step berikutnya dihitung dengan dua langkah, predictor dan corrector. Pada langkah predictor, dilakukan pengepingan beda maju untuk menyelesaikan turunan dalam ruang dan waktu.

f i* − f in ∂f = ...( III-8 ) ∂t Δt f in+1 − f in ∂f = ...( III-9 ) ∂x Δx Pada langkah corrector, dilakukan pengepingan beda maju untuk menyelesaikan turunan dalam waktu, sedangkan untuk turunan dalam ruang dilakukan pengepingan beda mundur.

f i** − f i* ∂f = ...( III-10 ) ∂t Δt

f i* − f i*−1 ∂f = ...( III-11 ) ∂x Δx Nilai pada time step selanjutnya diperoleh sebagai berikut: f t+1 = (f t + f**)/2...( III-12 ) Penerapan pada persamaan pengatur menghasilkan Kontunuitas

III-5

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

PREDICTOR „ du = u(i+1,j) h(i+1,j)-u(i,j) h(i,j); „ dv = v(i,j+1) h(i,j+1)-v(i,j) h(i,j); „ hp(i,j) = h (i,j)-(dt/dx*du)-(dt/dy*dv);

CORRECTOR „ du = up(i,j) hp(i,j)-up(i-1,j) hp(i-1,j); „ dv = vp(i,j) hp(i,j)-vp(i,j-1) hp(i,j-1); „ hc(i,j) = h (i,j)-(dt/dx*du)-(dt/dy*dv);

FINAL „ h_final(i,j) = (h(i,j) + h_correct(i,j))/2

Momentum PREDICTOR „ dhfi= h(i+1,j)-h(i,j); „ dufi= u(i+1,j)-u(i,j); „ dufj= u(i,j+1)-u(i,j); „ dhfj= h(i,j+1)-h(i,j); „ dvfi= v(i+1,j)-v(i,j); „ dvfj= v(i,j+1)-v(i,j); „ up(i,j)=u(i,j)-(u(i,j)*dt/dx*dufi)-(v(i,j)*dt/dy*dufj)-(g*dt/dx*dhfi)-

g*dt*((sox(i,j))+sfx(i,j)); „ vp(i,j)=v(i,j)-(u(i,j)*dt/dx*dvfi)-(v(i,j)*dt/dy*dvfj)-(g*dt/dy*dhfj)-

g*dt*((soy (i,j))+sfy(i,j)); CORRECTOR „ dhfi= hp(i,j)-hp(i-1,j); „

dufi= up(i,j)-up(i-1,j); III-6

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

„

dufj= up(i,j)-up(i,j-1);

„

dhfj= hp(i,j)-hp(i,j-1);

„

dvfi= vp(i,j)-vp(i-1,j);

„

dvfj= vp(i,j)-vp(i,j-1);

„

uc(i,j)=up(i,j)-(up(i,j)*dt/dx*dufi)-(vp(i,j)*dt/dy*dufj)-(g*dt/dx*dhfi)g*dt*((sox (i,j))+sfx(i,j));

„

vc(i,j)=vp(i,j)-(up(i,j)*dt/dx*dvfi)-(vp(i,j)*dt/dy*dvfj)-(g*dt/dy*dhfj)g*dt*((soy (i,j))+sfy(i,j));

FINAL „ u_n(i,j) = (u(i,j) + u_c (i,j))/2 „ v_n(i,j) = (v(i,j) + v_c (i,j))/2

III.4.2 Adams Bashforth Persamaan pengatur Boussinesq dikepingkan dengan menggunakan skema Adam Bashford. Persamaan boussinesq memiliki orde yang lebih tinggi dibandingkan dengan persamaan St. Venant. Dalam langkah penyelesaian numerik, diperlukan adanya penyederhanaan. ηt = E (η , u , v) ...( III-13 ) U t = F (η , u , v) + [ F1 (v)]t ...( III-14 ) Vt = G (η , u , v) + [G1 (u )]t ...( III-15 )

Dimana: U (u ) = u + h[b1 + hu xx + b2 (hu ) xx ] ...( III-16 ) V (v) = v + h[b1 + hv yy + b2 ( hv) yy ] ...( III-17 )

Dengan mensubstitusikan U dan V maka diperoleh persamaan: E (η , u , v ) = −[( h + η )u ]x − [(h + η )v ] y ...( III-18 )

III-7

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

F (η , u , v) = − gη x − (uu x + vu y )

...( III-19 )

G (η , u , v) = − gη x − (vv y + uvx ) ...( III-20 )

F1 (v) = − h[(b1hvxy + b2 (hv) xy ] ...( III-21 ) G1 (u ) = − h[(b1hu xy + b2 (hu ) xy ] ...( III-22 )

Turunan pertama diselesaikan dengan menggunakan central difference scheme. f n − fi −n1, j ∂f ...( III-23 ) = i +1, j ∂x 2 Δx

Turunan pertama dalam arah y diperoleh dengan cara yang sama. Turunan kedua diselesaikan sebagai berikut:

f i +n1, j − 2 fi ,nj + f i −n1, j ∂f ...( III-24 ) = ∂x Δx 2

Turunan kedua dalam arah y diperoleh dengan cara yang sama. Turunan silang diselesaikan sebagai berikut:

fi +n1, j +1 + fi −n1, j −1 + fi −n1, j +1 + fi +n1, j −1 ∂f = ...( III-25 ) 4ΔxΔy ∂x

Pengepingan terhadap waktu dilakukan dengan metode Adams Bashforth ηin, +j 1 = ηin, j +

Δt [23Ein, j − 16in,−j1 + 5in,−j 2 ] ...( III-26 ) 12

„ hapiu(i,j)=(h(i,j)+pi(i,j))*u(i,j);

III-8

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

„ hapiv(i,j)=(h(i,j)+pi(i,j))*v(i,j); „ dqx=1/2/dx*(hapiv(i+1,j)-hapiv(i-1,j)); „ dqy=1/2/dy*(hapiv(i,j+1)-hapiv(i,j-1)); „ dqx=1/2/dx*(hapiv(i+1,j)-hapiv(i-1,j)); „ E(i,j)= -dqx-dqy; „ pip(i,j)=pi(i,j)+dt/12*(23*E(i,j)-16*EE(i,j)+5*EEE(i,j)); U in, +j 1 = U in, j + Vi ,nj+1 = Vi ,nj +

Δt [23Fi ,nj − 16 Fi ,nj−1 + 5Fi ,nj− 2 ] + 2 F1ni , j − 3F1in, −j 1 + 2 F1in, −j 2 ...( III-27 ) 12

Δt [23Gin, j − 16Gin, −j 1 + 5Gin, −j 2 ] + 2G1ni , j − 3G1ni ,−j 1 + 2G1ni ,−j 2 ...( III-28 ) 12

„ dux =1/2/dx*(u(i+1,j)-u(i-1,j)); „ duy =1/2/dy*(u(i,j-1)+u(i,j+1)); „ dvx =1/2/dx*(v(i+1,j)-v(i-1,j)); „ dvy =1/2/dy*(v(i,j+1)-v(i,j-1)); „ dpix =1/2/dx*(pi(i+1,j)-pi(i-1,j)); „ dpiy =1/2/dy*(pi(i,j+1)-pi(i,j-1)); „ duxy =1/4/dx/dy*(u(i+1,j+1)+u(i-1,j-1)+u(i-1,j+1)+u(i+1,j-1)); „ dvxy =1/4/dx/dy*(v(i+1,j+1)+v(i-1,j-1)+v(i-1,j+1)+v(i+1,j-1)); „ dhuxy=1/4/dx/dy*(qx(i+1,j+1)+qx(i-1,j-1)+qx(i-1,j+1)+qx(i+1,j-1)); „ dhvxy=1/4/dx/dy*(qy(i+1,j+1)+qy(i-1,j-1)+qy(i-1,j+1)+qy(i+1,j-1)); „ F(i,j) = -g*dpix-(u(i,j)*dux+v(i,j)*duy); „ F1(i,j) = -h(i,j)*(b1*h(i,j)+dvxy+b2*dhvxy); „ G(i,j) = -g*dpiy-(v(i,j)*dvy+u(i,j)*dvx); „ G1(i,j) = -h(i,j)*(b1*h(i,j)*duxy+b2*dhuxy);

III-9

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

Elevasi muka air dapat diselesaikan langsung. Demikian juga dengan variabel U dan V. Akan tetapi, untuk mendapatkan nilai u dan v, maka persamaan U dan V perlu diselesaikan secara implisit. Penyelesaian implisit dilakukan dengan melakukan diskritisasi dan menulis ulang persamaan: U (u ) = u + h[b1 + hu xx + b2 (hu ) xx ] ...( III-29 ) V (v) = v + h[b1 + hv yy + b2 ( hv) yy ] ...( III-30 )

Persamaan ditulis ulang ke dalam bentuk: Ai uin−+1,1 j + Bi uin, +j 1 + Ci uin++1,1 j + Di , j = 0

...( III-31 )

Dimana: „ B = -2*(hp(i,j)^2/dx^2*b1+b2/dx^2*hp(i,j)*hp(i,j))+1 „ A = hp(i,j)^2/dx^2*b1+b2/dx^2*hp(i,j)*hp(i-1,j) „ C = hp(i,j)^2/dx^2*b1+b2/dx^2*hp(i,j)*hp(i+1,j) „ D = U(i,j)

Sistem persamaan diatas pada tiap nilai j akan menghasilkan matriks tridiagonal berikut: ⎡B2 ⎢ ⎢ A3 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

C2 B3 A4

C3 B4 ....

C4 .... ....

.... .... AN − 2

.... B N −2 AN −1

⎤ ⎡ u 2 ⎤ ⎡ D2 ⎤ ⎥⎢ ⎥ ⎢ ⎥ ⎥ ⎢ u 3 ⎥ ⎢ D3 ⎥ ⎥ ⎢ u 4 ⎥ ⎢ D4 ⎥ ⎥⎢ ⎥ ⎢ ⎥ ⎥ ⎢ .... ⎥ = ⎢ .... ⎥ ⎥ ⎢ .... ⎥ ⎢ .... ⎥ ⎥⎢ ⎥ ⎢ ⎥ C N − 2 ⎥ ⎢u N − 2 ⎥ ⎢ D N − 2 ⎥ B N −1 ⎥⎦ ⎢⎣ u N −1 ⎥⎦ ⎢⎣ D N −1 ⎥⎦

Matrix diatas diselesaikan untuk memperoleh u pada tiap time step dengan menggunakan built in functioni yang ada di MATLAB (matrix LU, reff. MATLAB Help).

III-10

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

III.5 Initial dan Boundary Condition III.5.1 Initial Condition Sebagai kondisi awal, diasumsikan elevasi air laut adalah 0. hal ini berarti: h(x,y,0) = sealevel – bottom ...( III-32) u(x,y,0) = 0 ...( III-33) v(x,y,0) = 0 ...( III-34)

III.5.2 Syarat Batas Gelombang Gelombang yang datang diasumsikan sebagai gelombang tunggal dengan periode dan amplitudo yang ditentukan. Untuk kondisi batas tempat dimana gelombang datang, diasumsikan gelombang yang datang sebagai gelombang linear dengan tinggi dan kecepatan yang berubah terhadap waktu dengan persamaan: h=d+(A*-sin(2pi*(t/T))); ...( III-35) u= (g*h)0.5 ...( III-36) v= (g*h)0.5 ...( III-37) dimana: h

: ketinggian gelombang

d

: kedalaman air

g

: gravitasi

u

: vektor kecepatan arah x

v

: vektor kecepatan arah

A

: Amplitudo gelombang datang

T

: Periode gelombang datang

t

: waktu

Jika terdapat data hasil pengukuran, maka syarat batas gelombang dapat berupa fluktuasi elevasi muka air terhadap waktu/

III-11

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

Pada waktu (t) lebih besar dari periode gelombang (T), diasumsikan air kembali tenang pada posisi awal. Penerapan kondisi diatas dilakukan dengan memodelkan arah gelombang datang sebisa mungkin tegak lurus garis pantai dan berada pada satu sumbu saja sehingga kecepatan di sumbu yang lain adalah 0.

III.5.3 Syarat Batas Bebas Sedangkan untuk syarat batas di ujung-ujung lain menggunakan gabungan syarat batas dirichlet dan neumann. Syarat batas dirichlet (tetap) digunakan untuk kecepatan pada arah tegak lurus sumbu. Sedangkan syarat batas neumann digunakan untuk kecepatan searah sumbu dan kedalaman air.

Syarat batas untuk sumbu x: ∂h ∂ (uh ) + = 0 ...( III-38) ∂t ∂x ∂U ∂ ∂ (h + z ) + U U + gh = − ghS fx ...( III-39) ∂t ∂x ∂x

V (i,J) = 0 ...( III-40)

Syarat batas untuk sumbu y: ∂h ∂ (vh ) + = 0 ...( III-41) ∂t ∂y

U (I,j) = 0 ...( III-42) ∂V ∂ ∂ (h + z ) + V V + gh = − ghS fy ...( III-43) ∂t ∂y ∂y

III.5.4 Syarat Batas Dinding Syarat batas untuk memodelkan adanya dinding seperti pada kasus gelombang dalam flume ataupun kasus dimana gelombang membentur suatu dinding vertikal menggunakan syarat batas dinding. Kecepatan pada arah tegak lurus dinding III-12

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

diberi nilai 0 sedangkan ketinggian dan kecepatan pada arah paralel dengan dinding diberi syarat batas bebas.

III.6 Wet/Dry Condition Batasan wet dry condition dilakukan dengan memberikan suatu nilai minimum untuk kedalaman air. Jika nilai hasil perhitungan numerik h<=toler, maka node tersebut akan diberi nilai h,u, dan v = 0.

i+1,j i,j i+1,j h(i,j)<=toler

i+1,j

i,j

i+1,j

h(i,j)>toler

Gambar III-2. Wet/Dry Condition

III.7 Filter Untuk mengatasi ketidak stabilan akibat adanya Caustic Wave, digunakan filter numerik Hansen (1962). F(i,j)=C*F(i,j)+0.25*(1-C)* (F(i,j)+F(i+1,j)+F(i,j-1)+F(i,j+1)) ...( III-44) Dengan nilai C diambil sebesar 0.99. Filter diaplikasikan ke tiap parameter (h,u,v) disetiap titik basah untuk setiap time step.

III-13

Pemodelan Aliran Permukaan 2 D Pada Suatu Lahan Akibat Rambatan Tsunami

i,j+1 i-1,j

i,j

i+1,j

i,j-1 Gambar III-3. Filter Numerik

III-14