Abstrak

Paper ini memperkenalkan penggunaan Metoda Elemen Hingga Diskrit untuk menyelesaikan persamaan gelombang panjang. Metoda ini tidak memerlukan matrix dengan ukuran yang besar.

Pada metoda ini perhitungan variabel pada suatu elemen hanya memperhitungkan pengaruh variabel pada elemen yang berbatasan. Formulasi persamaan pengatur pada elemen hingga digunakan prosedur Galerkin,sedangkan penyelesaian differential waktu digunakan metoda prediktor-korektor dari Adams-Bashforth-Moulton.

Model dicoba untuk mensimulasikan propagasi gelombang panjang pada suatu kanal dimana solusi analitisnya sudah tersedia. Hasil model numeris terlihat sangat mendekati solusi model analitis.

Kata-kata Kunci : diskrit, elemen yang berbatasan.

1. Pendahuluan

Penyelesaian persamaan gelombang panjang atau persamaan sirkulasi arus dengan berbagai metoda numeris telah banyak dilakukan, mulai dengan metoda selisih hingga, metoda elemen hingga sampai dengan metoda boundary fitted. Berbagai metoda numeris tersebut menghasilkan persamaan simultan atau matrix dengan ukuran yang besar. Semakin komplex geometri daerah perhitungan, semakin kecil ukuran grid atau elemen, akan semakin besar ukuran matrixnya. Walaupun komputer jenis PC pada saat ini sudah mempunyai kapasitas yang sangat besar, tetapi masih sering dijumpai kesulitan karena besarnya ukuran matrix.

Pada metoda selisih hingga dengan prosedur explisit tidak terbentuk suatu persamaan simultan, tetapi metoda ini tidak banyak digunakan karena diperlukan pertambahan waktu (∆t) yang sangat kecil dalam menyelesaikan komponen diferensial terhadap waktu. Pada prosedur selisih hingga explisit ini, perhitungan pada suatu titik hanya melibatkan beberapa titik disekitarnya. Berdasarkan prosedur tersebut, maka pada metoda elemen hingga akan dapat juga dikerjakan dengan hanya memperhitungkan pengaruh elemen yang berbatasan.

Pada paper ini, metoda (DFEM) digunakan untuk menyelesaikan persamaan gelombang panjang yang juga merupakan persamaan sirkulasi arus. Penerapan model pada suatu kanal memberikan hasil yang mendekati hasil dari solusi analitis.

Catatan : Usulan makalah dikirimkan pada tanggal 9 September 2002 dan dinilai oleh peer reviewer pada tanggal 16 September 2002–7 Pebruari 2003. Revisi penulisan dilakukan antara tanggal 10 Pebruari 2003 hingga 27 Pebruari 2003.

¹⁾ Staf Pengajar Departemen Teknik Sipil, ITB.

2. Persamaan Pengatur

Persamaan pengatur gelombang panjang terdiri dari 3 buah persamaan yaitu sebagai berikut.

a. Persamaan kontinuitas

\[\frac{\partial \eta}{\partial t} + \frac{\partial (uH)}{\partial x} + \frac{\partial (uv)}{\partial y} = 0 \tag{1}\]

b. Persamaan momentum arah x

\[\frac{\partial u}{\partial x} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + g \frac{\partial \eta}{\partial x} = 0\] (2)

c. Persamaan momentum arah y

\[\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + g \frac{\partial \eta}{\partial x} = 0\] (3)

= fluktuasi muka air

= kecepatan arus pada arah x

= kecepatan arus pada arah y

\(= \eta + h\)

h = kedalaman perairan

g = kecepatan gravitasi

3. Formulasi Persamaan pada Elemen

Pada penelitian ini digunakan elemen segiempat kuadratis. Pemilihan elemen ini adalah agar diperoleh turunan pertama yang kontinu pada sisi ataupun titik penghubung antar elemen dan diharapkan dapat memperkecil jumlah elemen pada satu panjang gelombang, sedangkan pemilihan elemen segiempat adalah untuk mempermudah pengembangan metoda diskrit yang akan dijelaskan pada bagian berikutnya.

Gambar 1. Fluktuasi muka air dan kedalaman

Pada elemen segiempat kuadratis, terdapat sembilan titik simpul pada setiap elemen (gambar 2), dan variabel pada elemen dinyatakan sebagai berikut.

\[\left\{\mathbf{u}\right\} = \sum_{i=1}^{9} \mathbf{N}_{i} \mathbf{u}_{i} = \left\lfloor \mathbf{N} \right\rfloor \left\{\mathbf{u}\right\} \tag{3.a}\]

\[\{\mathbf{v}\} = \sum_{i=1}^{9} \mathbf{N}_{i} \mathbf{v}_{i} = \lfloor \mathbf{N} \rfloor \{\mathbf{v}\}\] (3.b)

\[\{\eta\} = \sum_{i=1}^{9} N_{i} \eta_{i} = \lfloor N \rfloor \{\eta\}\] (3.c)

\[\left\{h\right\} = \sum_{i=1}^{9} N_{i} h_{i} = \left\lfloor N \right\rfloor \left\{h\right\}\] (3.d)

dimana N adalah fungsi bentuk.

Gambar 2. Elemen segiempat kuadratis

Prosedur formulasi persamaan yang digunakan adalah prosedur Galerkin yaitu.

\[\int_{\Omega_{e}} \{N\} L(u, v, \eta) \partial \Omega_{e} = 0\] (4)

dimana L(u, v, η) adalah persamaan kontinuitas atau persamaan momentum arah x atau arah v.

a. Persamaan kontinuitas

\[\frac{\partial \eta}{\partial t} + \frac{\partial (uH)}{\partial x} + \frac{\partial (vH)}{\partial y} = 0\]

Prosedur Galerkin pada persamaan ini menghasilkan

\[\begin{split} &\int\limits_{\Omega e} \left\{ N \right\} \left( \frac{\partial \eta}{\partial t} + \frac{\partial \left( uH \right)}{\partial x} + \frac{\partial \left( vH \right)}{\partial y} \right) d\Omega e = 0 \\ &\int\limits_{\Omega e} \left\{ N \right\} \left( \frac{\partial \eta}{\partial t} d\Omega e = - \int\limits_{\Omega e} \left\{ N \right\} \left( \frac{\partial \left( uH \right)}{\partial x} + \frac{\partial \left( vH \right)}{\partial y} \right) d\Omega e \\ &\left[ M \right] \left\{ \frac{d\eta}{dt} \right\} = \left\{ E \right\} \end{split} \tag{5.a}\]

\[\{E\} = -\int_{\Omega} \{N\} \left( \frac{\partial (uH)}{\partial x} + \frac{\partial (vH)}{\partial y} \right) d\Omega e\] (5.b)

\[[M] = \int_{\Omega_c} \{ N \} \lfloor N \rfloor d\Omega e \qquad (5.c)\]

b. Persamaan momentum arah x

\[\frac{\partial u}{\partial x} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + g \frac{\partial \eta}{\partial x} = 0\] dengan mengerjakan prosedur Galerkin akan diperoleh

\[\left[M\right] \left\{ \frac{\partial \mathbf{u}}{\partial t} \right\} = \left\{ F \right\} \tag{6.a}\]

\[\{F\} = -\int_{\Omega} \{N\} \left( u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + g \frac{\partial \eta}{\partial x} \right) d\Omega e \qquad (6.b)\]

\[[M] = \int_{\Omega_e} \{N\} \lfloor N \rfloor d\Omega e\] (6.c)

c. Persamaan momentum arah y

Pengerjaan prosedur Galerkin pada persamaan momentum pada arah y, akan menghasilkan,

\[\left[\mathbf{M}\right] \left\{ \frac{\partial \mathbf{v}}{\partial \mathbf{t}} \right\} = \left\{ \mathbf{G} \right\} \tag{7.a}\]

\[\{G\} = -\int_{\Omega e} \{N\} \left( u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + g \frac{\partial h}{\partial y} \right) d\Omega e \qquad (7.b)\]

\[[M] = \int_{\Omega_e} \{N\} \lfloor N \rfloor d\Omega e\] (7.c)

4. Penyelesaian Differensial Waktu

Differensial waktu, diselesaikan berdasarkan metoda prediktor-korektor yang dikembangkan oleh Adams -Bashforth - Moulton, vaitu:

Prediktor:

\[\begin{split} \left[M\right] & \left\{\eta^{*}\right\} = \left[M\right] \left\{\eta\right\}^{t} + \frac{\Delta t}{12} \left(23\left\{E\right\}^{t} \\ & -16\left\{E\right\}^{t-\Delta t} + 5\left\{E\right\}^{t-2\Delta t}\right) \\ \left[M\right] & \left\{u^{*}\right\} = \left[M\right] \left\{u\right\}^{t} + \frac{\Delta t}{12} \left(23\left\{F\right\}^{t} \\ & -16\left\{F\right\}^{t-\Delta t} + 5\left\{F\right\}^{t-2\Delta t}\right) \\ \left[M\right] & \left\{v^{*}\right\} = \left[M\right] \left\{v\right\}^{t} + \frac{\Delta t}{12} \left(23\left\{G\right\}^{t} \\ & -16\left\{G\right\}^{t-\Delta t} + 5\left\{G\right\}^{t-2\Delta t}\right) \end{split}\]

\(\eta^*\), \(u^*\), \(v^*\) adalah harga pendekatan dari

\[\left\{\eta\right\}^{t+\Delta t}\] , \(\left\{u\right\}^{t+\Delta t}\) , \(\left\{v\right\}^{t+\Delta t}\)

Korektor:

\[\text{[rumus tidak dapat ditampilkan dengan baik — lihat PDF asli]}\]

\[[M] \{u\}^{t+\Delta t} = [M] \{u\}^{t} + \frac{\Delta t}{24} \left(9 \{F^{*}\}\right)^{t+\Delta t} + [M] \{F\}^{t} - 5\{F\}^{t-\Delta t} + [M] \{F\}^{t-\Delta t}\]

\[\begin{split} \big[ M \, \big] \big\{ v \big\}^{t + \Delta t} \, &= \big[ M \, \big] \big\{ v \big\}^t \, + \frac{\Delta t}{24} \bigg( 9 \Big\{ G^{\, \star} \Big\} \\ &+ 19 \, \big\{ F \big\}^t \, - 5 \big\{ G \big\}^{t - \Delta t} \, + \big\{ G \big\}^{t - 2 \Delta t} \, \bigg) \end{split}\]

Prosedur korektor ini diulang-ulang sampai diperoleh harga dengan perbedaan yang sangat kecil dengan hasil iterasi sebelumnya, tetapi pada umumnya cukup dengan 2-3 kali iterasi.

Pada prosedur tersebut terdapat perkalian [M ]{u }^t yang selain memakan waktu juga dapat menambah kesalahan pembulatan. Untuk menghindari perkalian tersebut dilakukan perubahan sebagai berikut.

Prediktor:

\[[M] \{d\eta\}^* = \frac{\Delta t}{12} \left( 23 \{E\}^t - 16 \{E\}^{t-\Delta t} + 5 \{E\}^{t-2\Delta t} \right)\]\[\{\eta\}^* = \{\eta\}^t + \{d\eta\}^*\]

Begitu juga dengan persamaan-persamaan momentum serta pada prosedur korektor dilakukan modifikasi yang sama.

5. Metoda Elemen Hingga Diskrit

Perhitungan dengan metoda elemen hingga biasanya dikerjakan sekaligus pada seluruh elemen pada seluruh daerah perhitungan, sehingga terbentuk matriks dengan ukuran yang sangat besar. Karena itu pada penelitian ini dikembangkan metoda elemen hingga diskrit. Pada metoda ini perhitungan dilakukan pada elemen demi elemen dengan memperhitungkan pengaruh elemen yang berbatasan baik berbatasan sisi, maupun hanya titik saja. Prosedur perhitungan ini sangat didukung oleh prosedur iterative dari Adams -Bashforth - Moulton, seperti yang telah dijelaskan pada bagian sebelumnya.

Gambar 3. Pembagian grid

Sebagai ilustrasi digunakan suatu daerah perhitungan dengan grid seperti pada gambar 3. Pada saat melakukan perhitungan pada elemen 1, hanya dilibatkan elemen 2. jadi perhitungan dilakukan dengan subdomain elemen 1 dan 2 (gambar 4.a.). Perhitungan juga akan menghasilkan harga variabel pada elemen ke 2, tetapi hasil pada elemen ini tidak digunakan. Selanjutnya perhitungan pada elemen ke 2 akan melibatkan elemen ke 1 dan ke 3, dengan subdomain seperti gambar 4.b. Perhitungan pada elemen ke 6, akan melibatkan elemen-elemen ke 5, 9, 10, 11, 7, 19, 18, dan 17, dengan hasil yang digunakan hanya pada elemen ke 6.

Gambar 4.a. Subdomain untuk perhitungan pada elemen ke 1

Gambar 4.b. Subdomain untuk perhitungan pada elemen ke 2

Gambar 4.c. Subdomain untuk perhitungan pada elemen ke 6

Dengan prosedur ini dapat dihindari pembentukan matrix global yang besar, sehingga perhitungan dapat dilakukan dengan jumlah elemen dan titik simpul dalam jumlah yang besar tanpa dibatasi oleh kapasitas komputer.

6. Syarat Batas

Penyelesaian persamaan gelombang panjang ini memerlukan beberapa syarat batas yaitu:

• Syarat batas gelombang datang, yaitu titik-titik dimana karakteristik gelombang diberikan (diketahui), baik fluktuasi muka air maupun kecepatan arusnya.
• Syarat batas permukaan benda padat, yaitu kecepatan pada arah normal permukaan benda padat adalah nol.

Gambar 5. Syarat batas permukaan benda padat

Gambar 6. Syarat batas pada permukaan benda padat pada titik diskontinu

Pada titik A yang diskontinu, besar sudut normal dihitung dengan prosedur sebagai berikut:

Pada titik yang sangat dekat dengan titik A dengan jarak mendekati nol, disebelah kiri dan kanannya

\[\begin{split} &u_{A}\cos\alpha_{1}+v_{A}\sin\alpha_{1}=0\\ &u_{A}\cos\alpha_{2}+v_{A}\sin\alpha_{2}=0\\ &u_{A}\left(\cos\alpha_{1}+\cos\alpha_{2}\right)+v_{A}\left(\sin\alpha_{1}+\sin\alpha_{2}\right)=0\\ &u_{A}=v_{A}\frac{\sin\alpha_{1}+\sin\alpha_{2}}{\cos\alpha_{1}+\cos\alpha_{2}} \end{split}\] pada titik A:

\[u = -v \tan \alpha_A\]

Maka

\[\tan \alpha_A = \frac{\sin \alpha_1 + \sin \alpha_2}{\cos \alpha_1 + \cos \alpha_2}\]

\[\alpha_{A} = a \tan \left( \frac{\sin \alpha_{1} + \sin \alpha_{2}}{\cos \alpha_{1} + \cos \alpha_{2}} \right)\]

Selanjutnya syarat batas permukaan padat pada titik A adalah:

\[u_A \cos \alpha_A + v_A \sin \alpha_A = 0\]

a. Pada kanal

Model dikerjakan pada suatu kanal dengan lebar 20 m, panjang 400 m dan kedalaman 1.0 m. Gelombang sinusoidal dengan periode 90 detik, amplitudo 0.05 m masuk kedalam kanal. Sebagai syarat batas, pada titik diambang kanal diberikan fluktuasi muka air

\[\eta = 0.05 \sin \frac{2\pi}{90} t\]

Ukuran element yang digunakan adalah 20 m, sehingga pada satu panjang gelombang terdapat 14 elemen, sedangkan pertambahan waktu digunakan sebesar 1.5 dt atau 60 kali perhitungan untuk satu perioda gelombang.

Solusi analitis permasalahan tersebut adalah:

\(\eta = 0.05 \sin(\sigma t - k x)\)

\[\sigma = \frac{2\pi}{T}\]

T = perioda gelombang

\[k = \frac{2\pi}{L} = \frac{2\pi}{CT}\] = bilangan gelombang

\(C = \sqrt{gh} = \text{kecepatan gelombang}\)

g = percepatan gravitasi

= kedalaman perairan

Gambar 7. Kanal untuk pengujian model (tidak berskala)

pada Kanal (tidak berskala)

100m

Gambar 9. Profil muka air pada detik ke 180

Pada gambar 9, terlihat bahwa solusi model numeris sangat dekat dengan solusi analitis. Terdapat sedikit perbedaan, hal ini dikarenakan solusi analitis tidak memperhitungkan kondisi ujung kanal yang tertutup.

Dari hasil pengujian pada kasus 1 ini, terlihat bahwa model memberikan hasil yang cukup baik. Ukuran elemen, dapat digunakan elemen yang cukup besar yaitu 1/14 panjang gelombang. Untuk gelombang pasang surut dengan panjang gelombang mencapai ratusan kilometer, maka dapat digunakan ukuran elemen mencapai 10-15 km, tergantung pada luas areal perhitungan. Interval waktu yang dapat digunakan juga cukup besar yaitu 1/60 perioda gelombang, untuk gelombang pasang surut dengan perioda 25 jam, maka dapat digunakan interval waktu sebesar 1500 detik.

b. Pada kolam

Kasus berikutnya yang ditinjau adalah kolam dengan bentuk dan ukuran seperti pada gambaran berikut. Kedalaman perairan adalah 1 m. gelombang yang masuk kekolam gelombang sinusoidal berperioda 180 dt, amplitudo 0.05 m.

Gambar 10. Bentuk dan ukuran kolam pada kasus 2 (tidak berskala)

Gambar 11. Grid-point untuk kasus 2 (tidak berskala)

Solusi analitis untuk kasus 2 ini, baik elevasi muka air maupun arusnya belum tersedia. Solusi analitis untuk kasus ini memang sulit untuk dibuat mengingat adanya perubahan lebar kanal yang tiba-tiba. Perubahan lebar tersebut akan menimbulkan juga peristiwa difraksi atau penyebaran energi gelombang pada arah tegak lurus arah gelombang.

Seperti terlihat pada pola arus hasil model, terlihat terjadinya penyebaran gelombang pada arah lebar kolam (Gambar 13) seperti yang diharapkan.

Pada perkembangan terakhir, analisis difraksi gelombang dikembangkan dengan menggunakan persamaan Boussinesq. Sedangkan persamaan gelombang panjang yang digunakan pada penelitian ini merupakan bagian dari persamaan Boussinesq, yaitu dengan membuang suku nonlinier dari persamaan Boussinesq akan diperoleh persamaan gelombang panjang.

Secara umum, model yang akan dikembangkan dapat menggambarkan fenomena difraksi gelombang dan juga dapat menggambarkan pola arus pada suatu basin perairan dengan baik.

Gambar 12. Pola arus pada detik ke 60, skala : 1 cm = 0.5 m/dt2

Gambar 13. Pola arus pada detik ke 120, skala : 1 cm = 0.5 m/dt

Kesimpulan

Dari hasil penelitian yang telah dilakukan tersebut, maka dapat disimpulkan bahwa :

• a. Metoda elemen hingga diskrit dapat digunakan untuk menyelesaikan persamaan gelombang panjang atau persamaan sirkulasi arus.
• b. Dengan metoda ini dapat digunakan ukuran elemen yang cukup besar yaitu 1/14 panjang gelombang.
• c. Metoda ini tidak memerlukan interval waktu yang sangat kecil, yaitu dapat digunakan interval waktu 1/60 perioda gelombang.
• d. Pada bagian depan gelombang yang pertama terlihat adanya nois. Pada penelitian ini belum didapatkan metoda untuk menghilangkan nois. Perlu dilakukan penelitian lebih lanjut untuk menghilangkan nois.
• e. Pada gambar 9, terlihat profil muka air yang bergelombang (tidak smooth), yaitu pada gelombang yang pertama. Hal ini dikarenakan pantulan gelombang nois. Bila gelombang pantul mencapai titik batas, akan menimbulkan ketidak stabilan pada model. Karena itu perlu digunakan suatu absorbing boundary dalam memodelkan perambatan gelombang didalam kanal tertutup atau kolam dimana pada penelitian ini belum digunakan.

Dafar Pustaka

• Becker, Eric B., 1981, "Finite Elements, An Introduction Volume I, Prentice-Hall.
• Chen, Y., Liu, P.L-F., 1994, "Modified Boussinesq equations and Associated Parabolic Models for Water Wave Propagation", Fluid Mechanic (1995), Vol. 288, PP. 351-381, Cambridge University Press.
• Dean, Robert G., 1984, "Water Wave Mechanics for Engineers and Scientist", Prentice-Hall.
• Herbich, John B., 2000, "Handbook of Coastal Engineering", McGraw-Hill.
• Pinder, George F., 1977., "Finite Element Simulation in Surface and Subsurface Hydrology",Academic Press.
• Stansa, Frank L., 1985., "Applied Finite Element Analysis for Engineers", CBS International Editions.

Zienkiewicz, O.C., 1983, "Finite Elements and Aproximation", John Wiley & Sons.