1 Pendahuluan

Pada perencanaan jaringan tata air, masalah ekspansi (pelebaran) saluran seringkali tidak dapat dihindari sebagai akibat adanya kebutuhan penyesuaian dimensi saluran terhadap perubahan debit, persilangan bangunan dengan saluran, kondisi topografi, dll.

Penelitian mengenai struktur aliran pada saluran ekspansi telah banyak dilakukan, namun untuk aliran dengan fluida cair, kebanyakan dari penelitian tersebut bersifat eksperimental dan fokus pada masalah gerusan dan aliran di luar daerah sirkulasi, sementara itu untuk struktur aliran pada zona resirkulasi banyak dilakukan terutama untuk aliran dengan fluida udara. Uraian pada paragraf di bawah ini membahas rangkuman dari beberapa hasil penelitian mengenai struktur aliran turbulen akibat penurunan dasar saluran yang sebagian besar berlaku untuk kasus fluida udara.

Pada batas awal sebuah saluran ekspansi, massa fluida akan lepas landas dari dinding membentuk free mixing layer di bagian hilir dan akan menyentuh dinding kembali pada titik yang disebut reattachment point setelah menempuh jarak tertentu yang disebut reattachment length [1]. Pada saluran yang mengalami ekspansi sebesar H, posisi reattachment point berfluktuasi pada daerah sepanjang 2 H dan panjang reattachment length bervariasi antara 4-10 H, bergantung pada turbulensi aliran pada bagian hulu ekspansi [2]. Sepanjang jarak reattachment length tersebut aliran mengalami gradient tekanan positif dan tidak stabil karena membentuk zona resirkulasi pada daerah sekitar dinding. Pada zona ini, arah aliran berlawanan dengan arah aliran utama. Kondisi ini mengakibatkan aliran tidak stabil dan bersifat turbulen [3-5]. Pada zona tersebutlah struktur coherent turbulen akan terbentuk secara periodik melalui apa yang disebut fenomena bursting [2, 3, 5]. Pada zona resirkulasi kecepatan fluktuasi menjadi besar dan dapat melebihi kecepatan rata-rata aliran. Pada titik reattachment akan terjadi pembelahan vorteks yang berlawanan arah, satu vorteks menuju zona resirkulasi dan satu vorteks mengalir menuju hilir [2, 6]. Pada zona resirkulasi, vorteks mempunyai kemampuan untuk menstimulir pembentukan vorteks yang berlawanan arah sebagai konsekuensi dari kekekalan massa [2]. Pada prinsipnya, proses yang sama dengan karakteristik yang berbeda akan dialami oleh sebuah massa fluida cair bila dialirkan melalui saluran tersebut di atas. Perbedaan karakteristik ini terutama diakibatkan adanya pengaruh gravitasi pada fluida cair. Perbedaan yang dapat dengan mudah di identifikasi antara lain adalah bilangan Reynold karakteristik turbulen, panjang

reattachment lenght, bursting karakteristik dan ketebalan lapisan mixing layer [7-10].

Metoda penyelesaian masalah terbulen tersebut diatas telah banyak dikembangkan, antara lain adalah metoda Zero equation model, One equation model, Reynolds stress model dan Two equation \(\kappa\)-\(\varepsilon\) model. Namun demikian, karena tingkat akurasi yang diberikan cukup baik, metoda two equation \(\kappa\)-\(\varepsilon\) model lebih banyak dikembangkan, seperti yang akan disajikan/dibahas pada makalah ini.

2 Persamaan Pengatur

Persamaan pengatur aliran ini diturunkan berdasarkan persamaan Navier-Stokes untuk aliran turbulen 2 dimensi tak mampu mampat dalam bentuk depthaveraged velocity. Beberapa anggapan yang dipakai dalam melakukan penurunan persamaan pengatur tersebut adalah sbb:

• 1. Fluida tak mampu mampat (incompressible)
• 2. Kecepatan aliran yang ditinjau adalah kecepatan rata-rata
• 3. Aliran adalah dua dimensi (arah –x dan arah-y)
• 4. Kemiringan dasar saluran relatif kecil (sin~tangen~kemiringan saluran)
• 5. Distribusi tekanan fluida bersifat hidrostatis (viscous stress diabaikan)
• 6. Pengaruh gaya putaran bumi (efek coriolis) diabaikan
• 7. Konstituent tercampur merata (well mixed).

2.1 Persamaan Kontinuitas & Momentum

Dengan menerapkan prinsip dekomposisi Reynold terhadap besaran aliran turbulen, persamaan rata-rata kontinyuitas dan momentum untuk aliran tiga dimensi dapat dituliskan secara berturutan dalam bentuk sbb.:

\[\frac{\partial \overline{\mathbf{u}_{i}}}{\partial \mathbf{x}_{i}} = 0 \tag{1}\]

\[\frac{\partial \overline{u_i}}{\partial t} + \overline{u_j} \frac{\partial \overline{u_i}}{\partial x_j} = g_i - \frac{1}{\rho} \frac{\partial \overline{p}}{\partial x_i} + \frac{1}{\rho} \frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right) - \rho \overline{u_i' u_j'} \right]\](2)

Integrasi persamaan (1) dan (2) terhadap kedalaman dengan menggunakan metode "Leibnitz Rule" akan memberikan persamaan rata-rata kontinyuitas dan momentum dalam bentuk Depth-Averaged Velocity untuk aliran 2 dimensi tak mampu mampat sbb.:

Persamaan Kontinuitas

\[\left[\frac{\partial H}{\partial t}\right] + \left[\frac{\partial \left(UH\right)}{\partial x}\right] + \left[\frac{\partial \left(VH\right)}{\partial y}\right] = 0 \tag{3}\]

Persamaan Momentum Arah x

\[\frac{\partial \left(hU\right)}{\partial t} + \frac{\partial}{\partial x} \left(hU^2 + \frac{gh^2}{2}\right) + \frac{\partial \left(hUV\right)}{\partial y} = ghS_{0x} - \frac{1}{\rho}\tau_{bx} + \frac{1}{\rho}\frac{\partial \left(h\tau_{xx}\right)}{\partial x} + \frac{1}{\rho}\frac{\partial \left(h\tau_{xy}\right)}{\partial y} + \frac{1}{\rho}\tau_{Sx}\] \[\tag{4}\]

Arah y

\[\frac{\partial \left(hV\right)}{\partial t} + \frac{\partial}{\partial y} \left(hV^2 + \frac{gh^2}{2}\right) + \frac{\partial \left(hUV\right)}{\partial x} = ghS_{0y} - \frac{1}{\rho}\tau_{by} + \frac{1}{\rho}\frac{\partial \left(h\tau_{xy}\right)}{\partial y} + \frac{1}{\rho}\frac{\partial \left(h\tau_{yy}\right)}{\partial y} + \frac{1}{\rho}\tau_{sy}\] (5)

dimana

\[\begin{split} U &= \frac{1}{h + \eta} \int_{-h}^{\eta} u dz & V &= \frac{1}{h + \eta} \int_{-h}^{\eta} v dz \\ \tau_{xx} &= 2\mu \frac{\partial \overline{u}}{\partial x} - \rho \overline{u} \overline{u} \\ \tau_{yy} &= 2\mu \frac{\partial \overline{v}}{\partial y} - \rho \overline{v} \overline{v} \\ \tau_{xy} &= \mu \left( \frac{\partial \overline{u}}{\partial y} + \frac{\partial \overline{v}}{\partial x} \right) - \rho \overline{u} \overline{v} \end{split}\]

\(\tau_{bx}\) dan \(\tau_{sx}\) masing-masing merupakan tegangan dasar saluran dan tegangan permukaan air, sedangkan \(S_{0x}\) dan \(S_{0y}\) masing-masing merupakan kemiringan enerji pada arah x dan y.

2.2 Persamaan Model Turbulen κ- ε

Persamaan (2) sangat kompleks dan untuk penyelesaiannya membutuhkan Closure Problem dari suku "Reynold Stress" yang dalam hal ini dapat didekati dengan menggunakan pendekatan Boussinesq Eddy Viscosity sbb.:

\[\overline{-u_{i}'u_{j}'} = \hat{v}_{t} \left[ \frac{\partial (U_{i})}{\partial x_{i}} + \frac{\partial (U_{j})}{\partial x_{i}} \right] - \frac{2}{3} \hat{k}h \delta_{ij}\] (6)

dimana

\(\delta_{ii}\) = delta Kronecker

\[\delta_{ij} \begin{cases} i = j, \delta = 1 \\ i \neq j, \delta = 0 \end{cases}\] (7)

Dari analisis dimensional, didapatkan bahwa besaran eddy viscosity \((v_t)\) sebanding dengan karakteristik skala kecepatan (v) dan skala panjang \((\ell)\), yaitu \(v_t \approx v\ell\).

Model turbulen κ-ε dikembangkan untuk menyelesaikan persamaan Reynold Stress dengan menggunakan 2 persamaan tambahan yaitu Turbulence Kinetic-Energy Equation dan Turbulence Energy Dissipation Rate Equation. Persamaan Enerji Kinetik Turbulen tersebut dapat dituliskan sbb.:

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

Pada persamaan (8) tersebut, suku (i), (ii), (iii), (iv) dan (v) adalah suku-suku yang menunjukkan laju perubahan energi turbulen, convective diffusion energi turbulen, produksi energi turbulen, kerja turbulen stresses dan turbulent viscous dissipation. Bila v diasumsikan konstan maka persamaan (8) menjadi:

\[\frac{\partial \mathbf{k}}{\partial t} + \overline{\mathbf{u}}_{\mathbf{j}} \frac{\partial \mathbf{k}}{\partial x_{\mathbf{j}}} = -\frac{\partial}{\partial x_{\mathbf{i}}} \left[ \overline{\mathbf{u}_{\mathbf{i}}' \left( \frac{1}{2} \mathbf{u}_{\mathbf{j}}' \mathbf{u}_{\mathbf{j}}' + \frac{\mathbf{p}'}{\rho} \right)} \right] - \overline{\mathbf{u}_{\mathbf{i}}' \mathbf{u}_{\mathbf{j}}'} \frac{\partial \overline{\mathbf{u}}_{\mathbf{j}}}{\partial x_{\mathbf{i}}} + v \frac{\partial^{2} \mathbf{k}}{\partial x_{\mathbf{j}} \partial x_{\mathbf{j}}} - v \frac{\overline{\partial \mathbf{u}_{\mathbf{i}}'} \partial \mathbf{u}_{\mathbf{i}}'}{\partial x_{\mathbf{j}}} \frac{\partial \mathbf{u}_{\mathbf{i}}'}{\partial x_{\mathbf{j}}}\](9)

Launder & Spalding [11] menuliskan persamaan (9) diatas dalam bentuk:

\[\frac{\partial \mathbf{k}}{\partial \mathbf{t}} + \overline{\mathbf{u}}_{j} \frac{\partial \mathbf{k}}{\partial \mathbf{x}_{j}} = -\frac{\partial}{\partial \mathbf{x}_{j}} \left( \frac{\mathbf{v}_{t}}{\mathbf{\sigma}_{k}} \frac{\partial \mathbf{k}}{\partial \mathbf{x}_{j}} \right) + \mathbf{v}_{t} \left( \frac{\partial \overline{\mathbf{u}}_{i}}{\partial \mathbf{x}_{j}} + \frac{\partial \overline{\mathbf{u}}_{j}}{\partial \mathbf{x}_{i}} \right) \frac{\partial \overline{\mathbf{u}}_{i}}{\partial \mathbf{x}_{j}} - C_{D} \frac{\mathbf{k}^{3/2}}{\ell}\](10)

Sementara itu, persamaan Turbulence Energy Dissipation Rate (ε) dapat dituliskan sbb.:

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

dimana suku (i), (ii)+(iii), (iv), (v), (vi)+(vii)+(viii) adalah suku-suku yang menunjukkan laju perubahan \(\varepsilon\), transport energi kinetis oleh interaksi dengan gerakan rata-rata (generation of \(\varepsilon\) by mean flow), transfer energi kinetis oleh efek kecepatan (generation of \(\varepsilon\) by self stretching of vortex tube), viscous destruction dan difusi. Persamaan (11) disederhanakan oleh Launder and Spalding [11] menjadi sebagai berikut:

\[\frac{\partial \varepsilon}{\partial t} + \overline{u}_{j} \frac{\partial \varepsilon}{\partial x_{j}} = -\frac{\partial}{\partial x_{j}} \left( \frac{v_{t}}{\sigma_{\varepsilon}} \frac{\partial \varepsilon}{\partial x_{j}} \right) + C_{1\varepsilon} \frac{\varepsilon}{k} P - C_{2\varepsilon} \frac{\varepsilon^{2}}{k}\] (12)

dimana.

\[P = v_t \left( \frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i} \right) \frac{\partial \overline{u}_i}{\partial x_j}\] (13)

Persamaan tersebut diatas hanya dapat digunakan untuk menyelesaikan closure-problem dari persamaan reynolds stress dalam bentuk depth averaged setelah dimodifikasi dalam bentuk depth integrated. Dalam hal ini, Chapman and Kuo [12] telah memodifikasi model Rastogi and Rodi [13] dalam bentuk persamaan yang lebih konsisten dengan persamaan kontinuitas depth averaged dan persamaan momentum depth averaged sehingga diperoleh hubungan depth-integrated Reynolds stresses dengan depth integrated strain rates sebagai berikut:

\[-\frac{1}{h} \int_{z_{b}}^{h+z_{b}} \overline{u'_{i}u'_{j}} dz = \hat{v} \left[ \frac{\partial (v U_{i})}{\partial x_{j}} + \frac{\partial (v U_{j})}{\partial x_{i}} \right] - \frac{2}{3} \hat{k}h \delta_{ij}\] (14)

dimana \(\delta_{ij}\) adalah delta kronecker dan \(\hat{V_t}\) adalah Viskositas Turbulen untuk depth averaged yang dapat dituliskan sebagai berikut:

\[\hat{\mathcal{V}}_{t} = C_{\mu} \frac{\hat{k}^{2}}{\hat{\varepsilon}} \tag{15}\]

Pemodelan pada makalah ini menggunakan model Turbulen κ-ε yang diperoleh dari pendekatan Chapman and Kuo [12] tersebut di atas sehingga diperoleh persamaan energi kinetis turbulen dan laju disipasi energi turbulen yang berlaku untuk persamaan depth averaged sbb.:

• Persamaan K

\[\frac{\partial(h\hat{k})}{\partial t} + \frac{\partial(h \cup \hat{k})}{\partial x} + \frac{\partial(h \cup \hat{k})}{\partial y} = \frac{\partial}{\partial x} \left[ \frac{\hat{v}_t}{\sigma_k} \frac{\partial(h\hat{k})}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{\hat{v}_t}{\sigma_k} \frac{\partial(h\hat{k})}{\partial y} \right] + p_h + p_k - \hat{\epsilon}h (16)\]

Persamaan ε

\[\frac{\partial(h\hat{\varepsilon})}{\partial t} + \frac{\partial(h\,U\,\hat{\varepsilon})}{\partial x} + \frac{\partial(h\,V\,\hat{\varepsilon})}{\partial y} = \frac{\partial}{\partial y} \left[ \frac{\hat{v}_{t}}{\sigma_{\varepsilon}} \frac{\partial(h\hat{\varepsilon})}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{\hat{v}_{t}}{\sigma_{\varepsilon}} \frac{\partial(h\hat{\varepsilon})}{\partial y} \right] + \frac{\hat{\varepsilon}}{\hat{k}} (C_{1}p_{h} - C_{2}\hat{\varepsilon}h) + p_{\varepsilon} \tag{17}\]

dimana.

\[P_{h} = \frac{\hat{v}_{t}}{h} \left\{ 2 \left[ \frac{\partial (h U)}{\partial x} \right]^{2} + 2 \left[ \frac{\partial (h V)}{\partial y} \right]^{2} + \left[ \frac{\partial (h U)}{\partial y} + \frac{\partial (h V)}{\partial x} \right]^{2} \right\}\](18)

\[P_{k} = \frac{g}{C^{2}} q^{3}, P_{\varepsilon} = \frac{C_{2} C_{\mu}^{1/2} g^{5/4} q^{4}}{h D^{1/2} C^{5/2}}, q = \sqrt{U^{2} + V^{2}}\] (19)

dimana besaran koefisien-koefisien pada persamaan tersebut adalah : \(C_{\mu}\)= 0.09, \(C_1\)= 1.44, \(C_2\) = 1.92, \(\sigma_k\) = 1.0, \(\sigma_\epsilon\) = 1.3 dan D = 0.075

Secara lengkap persamaan hidrodinamik dan \(\kappa\)-\(\epsilon\) diatas dapat dituliskan kembali dalam bentuk sebagai berikut:

\[\frac{\partial}{\partial t} \begin{bmatrix} H \\ HU \\ HV \\ H\hat{k} \\ H\hat{\epsilon} \end{bmatrix} + \frac{\partial}{\partial x} \begin{bmatrix} UH \\ UVH \\ UVH \\ UH \hat{k} \\ H\hat{\epsilon} \end{bmatrix} + \frac{\partial}{\partial y} \begin{bmatrix} VH \\ UVH \\ HV \\ HV \hat{\epsilon} \end{bmatrix} = \frac{\partial}{\partial x} \begin{bmatrix} 0 \\ 2v_t \frac{\partial HU}{\partial x} - \frac{2}{3}H\hat{k} \\ \hat{v}_t \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial y} \right) \\ \hat{v}_t \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial y} \right) \\ \frac{\hat{v}_t}{\sigma_k} \frac{\partial (H\hat{k})}{\partial y} \\ \frac{\hat{v}_t}{\sigma_k} \frac{\partial (H\hat{k})}{\partial y} \end{bmatrix} + \frac{\partial}{\partial y} \begin{bmatrix} \frac{1}{y} \left[ \frac{\partial (hU)}{\partial y} \right]^2 + 2 \left[ \frac{\partial (hV)}{\partial y} \right]^2 + 2 \left[ \frac{\partial (hV)}{\partial y} \right]^2 \right] \\ \frac{1}{y} \left[ \frac{gHS_{ox}}{\partial x} - \frac{gV\sqrt{U^2 + V^2}}{2V^2} + \frac{p_aC^*W_xW}{\rho} \right] \\ \frac{\hat{v}_t}{\hat{v}_t} \frac{\partial (H\hat{k})}{\partial y} \\ \frac{\hat{v}_t}{\sigma_k} \frac{\partial (H\hat{k})}{\partial y} \end{bmatrix} + \frac{\partial}{\partial y} \begin{bmatrix} \frac{1}{y} \left[ \frac{\partial (hU)}{\partial y} + \frac{\partial (hV)}{\partial y} \right]^2 + 2 \left[ \frac{\partial (hV)}{\partial y} \right]^2 + 2 \left[ \frac{\partial (hV)}{\partial y} \right]^2 \\ \frac{1}{k} \left[ \frac{g}{v} \left[ \frac{1}{v} \left[ \frac{1}{v} \left[ \frac{\partial (hU)}{\partial y} + \frac{\partial (hV)}{\partial y} \right]^2 \right] - C_2\hat{\epsilon}h \right] \\ \frac{1}{k} \left[ \frac{C_2C_{\mu}^{1/2}g^{5/4} \left( \sqrt{U^2 + V^2} \right)^4}{hD^{1/2}C^{5/2}} \right] \end{bmatrix}\] \[(20)\]

2.3 Persamaan Transport Konstituent Kualitas Air

Persamaan transport konveksi-difusi, yang penurunannya didasarkan pada penggunaan hukum kekekalan massa di ruang tilik, dapat dituliskan sbb.:

\[\frac{\partial}{\partial t} (\Phi) + \frac{\partial}{\partial x} (U\Phi) + \frac{\partial}{\partial x} (V\Phi) = \frac{\partial}{\partial x} \left( D_x \frac{\partial \Phi}{\partial x} \right) + \frac{\partial}{\partial x} \left( D_y \frac{\partial \Phi}{\partial y} \right)\](21)

2.4 Syarat Batas

Pada kondisi batas perairan diterapkan syarat batas Inward Difference, sedangkan pada dinding dianggap kecepatan aliran sama dengan nol.

3 Penyelesaian Numerik

Metoda penyelsaian numerik dari persamaan yang digunakan adalah sbb.:

• 1. Persamaan hidrodinamik diselesaikan dengan Skema Mac Cormack-Splitting.
• 2. Persamaan κ-ε diselesaikan dengan Skema QUICKEST pada suku konveksi, Skema Central Difference pada suku difusi, dan Skema Euler pada suku reaksi.
• 3. Persamaan kualitas air/material transport (salinitas) dengan skema QUICKEST.

3.1 Penyelesaian Numerik Hidrodinamika dengan Teknik Mac Cormack-Splitting

Secara lengkap skema splitting untuk persamaan Hidrodinamik adalah:

\[F^{n+2} = \left\lceil \left( L_x L_y L_{xx} L_{yy} L_s \right) \bullet \left( L_s L_{yy} L_{xx} L_y L_x \right) \right\rceil F^n\]

Suku L_x, L_y, L_xx, L_yy dan L_s diselesaikan dengan metode Mac Cormack sbb.:

1. Penyelesaian numerik persamaan diferensial orde satu \(L_x\)

Predictor

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n} - \frac{\Delta t}{\Delta x} \left\{ \begin{bmatrix} UH \\ U^{2}H \\ UVH \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} UH \\ U^{2}H \\ UVH \end{bmatrix}_{i-1,j}^{n} \right\} - gH_{i,j}^{n} \frac{\Delta t}{\Delta x} \left\{ \begin{bmatrix} 0 \\ H+z \\ 0 \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} 0 \\ H+z \\ 0 \end{bmatrix}_{i-1,j}^{n} \right\} + \frac{\Delta t}{\Delta x} \left\{ \begin{bmatrix} 0 \\ -\frac{2}{3}H\hat{k} \\ 0 \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} 0 \\ -\frac{2}{3}H\hat{k} \\ 0 \end{bmatrix}_{i-1,j}^{n} \right\} \tag{22}\]

Corrector

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} - \frac{\Delta t}{\Delta x} \left\{ \begin{bmatrix} UH \\ U^{2}H \\ UVH \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} UH \\ U^{2}H \\ UVH \end{bmatrix}_{i,j}^{*} \right\} - gH_{i,j}^{n} \frac{\Delta t}{\Delta x} \left\{ \begin{bmatrix} 0 \\ H+z \\ 0 \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} 0 \\ H+z \\ 0 \end{bmatrix}_{i,j}^{*} \right\} + \frac{\Delta t}{\Delta x} \left\{ \begin{bmatrix} 0 \\ -\frac{2}{3}H\hat{k} \\ 0 \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} 0 \\ -\frac{2}{3}H\hat{k} \\ 0 \end{bmatrix}_{i,j}^{*} \right\} \tag{23}\]

<u>n+1</u>

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n+1} = \frac{1}{2} \left\{ \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} + \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} \right\}\](24)

2. Penyelesaian numerik persamaan diferensial orde satu L_v

Predictor

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n} - \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} UH \\ UVH \\ V^{2}H \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} UH \\ UVH \\ V^{2}H \end{bmatrix}_{i-1,j}^{n} \right\} - gH_{i,j}^{n} \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} 0 \\ 0 \\ z+H \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} 0 \\ 0 \\ z+H \end{bmatrix}_{i-1,j}^{n} \right\} + \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} 0 \\ 0 \\ 0 \\ -\frac{2}{3}H\hat{k} \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} 0 \\ 0 \\ -\frac{2}{3}H\hat{k} \end{bmatrix}_{i-1,j}^{n} \right\} (25)\]

Corrector

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} - \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} UH \\ UVH \\ V^{2}H \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} UH \\ UVH \\ V^{2}H \end{bmatrix}_{i,j}^{*} \right\} - gH_{i,j}^{n} \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} 0 \\ 0 \\ z+H \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} 0 \\ 0 \\ z+H \end{bmatrix}_{i,j}^{*} \right\} + \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} 0 \\ 0 \\ -\frac{2}{3}H\hat{k} \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} 0 \\ 0 \\ -\frac{2}{3}H\hat{k} \end{bmatrix}_{i,j}^{*} \right\} \tag{26}\]

n+1

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n+1} = \frac{1}{2} \left\{ \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} + \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} \right\}\] (27)

3. Penyelesaian numerik persamaan diferensial orde dua \(L_{xx}\)

Predictor

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n} + \frac{\Delta t}{\Delta x} \begin{bmatrix} 0 \\ 2\nu_{t} \frac{\partial HU}{\partial x} \\ \hat{\nu}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} 0 \\ 2\nu_{t} \frac{\partial HU}{\partial x} \\ \hat{\nu}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i-1,j}^{n}\](28)

Corrector

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} + \frac{\Delta t}{\Delta y} \begin{bmatrix} 0 \\ 2\nu_{t} \frac{\partial HU}{\partial x} \\ \hat{\nu}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} 0 \\ 2\nu_{t} \frac{\partial HU}{\partial x} \\ \hat{\nu}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i,j}^{*}\] \[(29)\]

n+1

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n+1} = \frac{1}{2} \left\{ \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} + \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} \right\}\] (30)

4. Penyelesaian numerik persamaan diferensial orde satu L_vv

Predictor

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n} + \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} 0 \\ \hat{v}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} 0 \\ \hat{v}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i,j}^{n} - \begin{bmatrix} 0 \\ \hat{v}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i-1,j}^{n} \right\}\] \[(31)\]

<u>Corrector</u>

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} + \frac{\Delta t}{\Delta y} \left\{ \begin{bmatrix} 0 \\ \hat{v}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i+1,j}^{*} - \begin{bmatrix} 0 \\ \hat{v}_{t} \left( \frac{\partial HU}{\partial y} + \frac{\partial HV}{\partial x} \right) \end{bmatrix}_{i,j}^{*} \right\}\](32)

n+1

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n+1} = \frac{1}{2} \left\{ \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^* + \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n+1} \right\}\] (33)

5. Penyelesaian numerik persamaan reaksi L_s

Predictor

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n} + \Delta t \begin{bmatrix} 0 \\ -g \frac{U\sqrt{U^{2} + V^{2}}}{C^{2}H} + + \frac{\rho_{a}C * W_{x}W}{\rho} \\ -g \frac{V\sqrt{U^{2} + V^{2}}}{C^{2}H} + \frac{\rho_{a}C * W_{y}W}{\rho} \end{bmatrix}_{i,j}^{n}\](34)

Corrector

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} = \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{*} + \Delta t \begin{bmatrix} 0 \\ -g \frac{U\sqrt{U^{2} + V^{2}}}{C^{2}H} + \frac{\rho_{a}C^{*}W_{x}W}{\rho} \\ -g \frac{V\sqrt{U^{2} + V^{2}}}{C^{2}H} + \frac{\rho_{a}C^{*}W_{y}W}{\rho} \end{bmatrix}_{i,j}^{*}\] (35)

n+1

\[\begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{n+1} = \frac{1}{2} \left\{ \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^* + \begin{bmatrix} H \\ UH \\ VH \end{bmatrix}_{i,j}^{**} \right\}\] (36)

3.2 Penyelesaian Numerik Suku Turbulen κ-ε dengan Teknik Splitting

1. Penyelesaian numerik persamaan konveksi \(\kappa\)-\(\epsilon\) \(L_{xke}\)

Konveksi \(\kappa\)-\(\epsilon\) dalam arah X diselesaikan dengan menggunakan skema QUICKEST:

\[\left[h\hat{k}\right]_{i,j}^{n+1} = \left[h\hat{k}\right]_{i,j}^{n} - \frac{\Delta t}{\Delta x} \left(F\left(h\hat{k}\right)_{i+\frac{1}{2},j} - F\left(h\hat{k}\right)_{i-\frac{1}{2},j}\right)\]

dimana:

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

2. <u>Penyelesaian numerik persamaan konveksi κ-ε L_yke</u>

Konveksi \(\kappa\)-\(\epsilon\) dalam arah Y diselesaikan dengan menggunakan skema QUICKEST:

\[\left[h\hat{k}\right]_{i,j}^{n+1} = \left[h\hat{k}\right]_{i,j}^{n} - \frac{\Delta t}{\Delta v} \left(F\left(h\hat{k}\right)_{i,j+\frac{1}{2}} - F\left(h\hat{k}\right)_{i,j-\frac{1}{2}}\right)\] (3.1)

dimana:

\[\begin{split} F\left(h\hat{k}\right)_{i,j+\frac{1}{2}} &= V_{i,j+\frac{1}{2}} * \Phi_{i,j+\frac{1}{2}} \\ V_{i,j+\frac{1}{2}} &= \frac{V_{i,j} + V_{i,j+1}}{2} \\ \left(h\hat{k}\right)_{i,j+\frac{1}{2}} &= \left[ \varphi_{lin} - Cr^R Grad_{i,j}^R \frac{\Delta y}{2} + \left\{ \frac{\alpha_R}{2} - \frac{1}{6} \left( 1 - \left(Cr^R \right)^2 \right) \right\} Curv_{i,j}^R \Delta y^2 \right] \\ \left(h\hat{k}\right)_{lin} &= \frac{\left(h\hat{k}\right)_{i,j} + \left(h\hat{k}\right)_{i,j+1}}{2} \\ Cr^R &= \frac{V_{i,j+\frac{1}{2}} \Delta t}{\Delta y} \\ Grad_{i,j}^R &= \frac{\varphi_{i,j+1} - \varphi_{i,j}}{\Delta y} \\ \alpha^R &= \frac{D_{i,j+\frac{1}{2}} \Delta t}{\left(\Delta y\right)^2} = \frac{D_{i,j} + D_{i,j+1}}{2} * \frac{\Delta t}{\left(\Delta y\right)^2} \\ Curv_{i,j}^R &= \frac{\left(h\hat{k}\right)_{i,j+1}^n - 2\left(h\hat{k}\right)_{i,j}^n + \left(h\hat{k}\right)_{i,j-1}^n}{\Delta y^2} \text{ jika } V_{i,j+\frac{1}{2}} > 0 \\ &= \frac{\left(h\hat{k}\right)_{i,j+2}^n - 2\left(h\hat{k}\right)_{i,j+1}^n + \left(h\hat{k}\right)_{i,j}^n}{\Delta y^2} \text{ jika } V_{i,j+\frac{1}{2}} > 0 \end{split}\]

\[\begin{split} F\left(h\hat{k}\right)_{i,j-\frac{1}{2}} &= V_{i,j-\frac{1}{2}} *\Phi_{i-\frac{1}{2},j} \\ V_{i,j-\frac{1}{2}} &= \frac{V_{i,j} + V_{i,j-1}}{2} \\ \left(h\hat{k}\right)_{i,j-\frac{1}{2}} &= \left[\left(h\hat{k}\right)_{lin} - Cr^L Grad_{i,j}^L \frac{\Delta y}{2} + \left\{\frac{\alpha_L}{2} - \frac{1}{6}\left(1 - \left(Cr^L\right)^2\right)\right\} Curv_{i,j}^L \Delta y^2\right] \\ \left(h\hat{k}\right)_{lin} &= \frac{\left(h\hat{k}\right)_{i,j} + \left(h\hat{k}\right)_{i,j-1}}{2} \\ Cr^L &= \frac{V_{i,j-\frac{1}{2}} \Delta t}{\Delta y} \\ Grad_{i,j}^L &= \frac{\varphi_{i,j} - \varphi_{i,j-1}}{\Delta y} \\ \alpha^L &= \frac{D_{i,j-\frac{1}{2}} \Delta t}{\left(\Delta x\right)y^2} = \frac{D_{i,j} + D_{i,j-1}}{2} *\frac{\Delta t}{\left(\Delta y\right)^2} \\ Curv_{i,j}^L &= \frac{\left(h\hat{k}\right)_{i,j}^n - 2\left(h\hat{k}\right)_{i,j-1}^n + \left(h\hat{k}\right)_{i,j-2}^n}{\Delta y^2} \; jika \; V_{i,j-\frac{1}{2}} > 0 \\ &= \frac{\left(h\hat{k}\right)_{i,j+1}^n - 2\left(h\hat{k}\right)_{i,j}^n + \left(h\hat{k}\right)_{i,j-1}^n}{\Delta y^2} \; jika \; V_{i,j-\frac{1}{2}} < 0 \end{split}\]

3. Penyelesaian numerik persamaan difusi \(\kappa\)-\(\epsilon\) \(L_{xxke}\)

Difusi κ-ε dalam arah X diselesaikan dengan menggunakan Central Scheme:

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

dimana:

\[\begin{bmatrix} \frac{\partial (H\hat{k})}{\partial x} \\ \frac{\partial (H\hat{\epsilon})}{\partial x} \end{bmatrix}_{i,j}^{n} = \begin{bmatrix} \frac{(H\hat{k})_{i+1,j} - (H\hat{k})_{i-1,j}}{\Delta x} \\ \frac{(H\hat{\epsilon})_{i+1,j} - (H\hat{\epsilon})_{i-1,j}}{\Delta x} \end{bmatrix}^{n}\]

4. Penyelesaian numerik persamaan difusi κ-ε L_yyke

Difusi κ-ε dalam arah Y diselesaikan dengan menggunakan Central Scheme:

\[\begin{bmatrix} H\hat{k} \\ H\hat{\epsilon} \end{bmatrix} = \begin{bmatrix} \frac{\hat{v}_t}{\sigma_k} \\ \frac{\hat{v}_t}{\sigma_{\epsilon}} \end{bmatrix} \frac{\partial}{\partial y} \begin{bmatrix} \frac{\partial(H\hat{k})}{\partial y} \\ \frac{\partial(H\hat{\epsilon})}{\partial y} \end{bmatrix} + \begin{bmatrix} \frac{\partial(H\hat{k})}{\partial y} \\ \frac{\partial(H\hat{\epsilon})}{\partial y} \end{bmatrix} \frac{\partial}{\partial y} \begin{bmatrix} \frac{\hat{v}_t}{\sigma_k} \\ \frac{\hat{v}_t}{\sigma_{\epsilon}} \end{bmatrix}\] (38)

\[\begin{bmatrix} H\hat{k} \\ H\hat{\epsilon} \end{bmatrix}_{i,j}^{n+1} = \begin{bmatrix} H\hat{k} \\ H\hat{\epsilon} \end{bmatrix}_{i,j}^{n} + \begin{bmatrix} \frac{\hat{v}_{t}}{\sigma_{k}} \\ \frac{\hat{v}_{t}}{\sigma_{\epsilon}} \end{bmatrix}_{i,j}^{n} \underbrace{\frac{\Delta t}{\Delta y}} \left\{ \begin{bmatrix} \frac{\partial (H\hat{k})}{\partial y} \\ \frac{\partial (H\hat{\epsilon})}{\partial y} \end{bmatrix}_{i,j+1}^{n} - \begin{bmatrix} \frac{\partial (H\hat{k})}{\partial y} \\ \frac{\partial (H\hat{\epsilon})}{\partial y} \end{bmatrix}_{i,j-1}^{n} \right\} \\ + \begin{bmatrix} \frac{\partial (H\hat{k})}{\partial y} \\ \frac{\partial (H\hat{\epsilon})}{\partial x} \end{bmatrix}_{i,j}^{n} \underbrace{\frac{\Delta t}{\Delta y}} \left\{ \begin{bmatrix} \frac{\hat{v}_{t}}{\sigma_{k}} \\ \frac{\hat{v}_{t}}{\sigma_{\epsilon}} \end{bmatrix}_{i,j+1}^{n} - \begin{bmatrix} \frac{\hat{v}_{t}}{\sigma_{k}} \end{bmatrix}_{i,j-1}^{n} \right\} \tag{39}\]

dimana.

\[\begin{bmatrix} \frac{\partial (H\hat{k})}{\partial y} \\ \frac{\partial (H\hat{\epsilon})}{\partial y} \end{bmatrix}_{i,i}^{n} = \begin{bmatrix} \frac{(H\hat{k})_{i,j+1} - (H\hat{k})_{i,j-1}}{\Delta y} \\ \frac{(H\hat{\epsilon})_{i,j+1} - (H\hat{\epsilon})_{i,j-1}}{\Delta y} \end{bmatrix}^{n}\]

5. Penyelesaian numerik persamaan reaksi κ-ε L_ske

Reaksi κ-ε diselesaikan dengan menggunakan Euler Scheme:

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

3.3 Penyelesaian Persamaan Transport

Bentuk skema numerik QUICKEST untuk persamaan transport di atas adalah:

\[\Phi_{i,j}^{n+1} = \Phi_{i,j}^{n} - \frac{\Delta t}{\Delta x} \left( F(\Phi)_{i+\frac{1}{2},j} - F(\Phi)_{i-\frac{1}{2},j} \right) - \frac{\Delta t}{\Delta y} \left( F(\Phi)_{i,j+\frac{1}{2}} - F(\Phi)_{i,j-\frac{1}{2}} \right) \\ + \frac{\Delta t}{\Delta x} \left( D_{x} F(\Phi)_{i+\frac{1}{2},j} - D_{x} F(\Phi)_{i-\frac{1}{2},j} \right) + \frac{\Delta t}{\Delta y} \left( D_{y} F(\Phi)_{i,j+\frac{1}{2}} - D_{y} F(\Phi)_{i,j-\frac{1}{2}} \right) \tag{41}\]

Fluks konvektif akibat U

\[\begin{split} F\left(\Phi\right)_{i+\frac{1}{2},j} &= flukscx_{i,j}^{R} \\ &= U_{i+\frac{1}{2},j} * \Phi_{i+\frac{1}{2},j} \\ U_{i+\frac{1}{2},j} &= \frac{U_{i,j} + U_{i+1,j}}{2} \\ \Phi_{i+\frac{1}{2},j} &= \left[ \varphi_{lin} - Cr^{R} \, Grad_{i,j}^{R} \, \frac{\Delta x}{2} + \left\{ \frac{\alpha_{R}}{2} - \frac{1}{6} \left( 1 - \left( Cr^{R} \, \right)^{2} \, \right) \right\} Curv_{i,j}^{R} \Delta x^{2} \, \right] \\ \varphi_{lin} &= \frac{\varphi_{i,j} + \varphi_{i+1,j}}{2} \\ Cr^{R} &= \frac{U_{i+\frac{1}{2},j}}{\Delta x} \\ Grad_{i,j}^{R} &= \frac{\varphi_{i+1,j} - \varphi_{i,j}}{\Delta x} \\ \alpha^{R} &= \frac{D_{i+\frac{1}{2},j}}{(\Delta x)^{2}} = \frac{D_{i,j} + D_{i+1,j}}{2} * \frac{\Delta t}{(\Delta x)^{2}} \\ Curv_{i,j}^{R} &= \frac{\varphi_{i+1,j}^{n} - 2\varphi_{i,j}^{n} + \varphi_{i-1,j}^{n}}{\Delta x^{2}} \quad jika \, U_{1+\frac{1}{2},j} > 0 \\ &= \frac{\varphi_{i+2,j}^{n} - 2\varphi_{i+1,j}^{n} + \varphi_{i,j}^{n}}{\Delta x^{2}} \quad jika \, U_{1+\frac{1}{2},j} < 0 \end{split}\]

\[\begin{split} F\left(\Phi\right)_{i-\frac{1}{2},j} = & = \text{flukscx}_{i,j}^L \\ &= U_{i-\frac{1}{2},j} * \Phi_{i-\frac{1}{2},j} \\ U_{i-\frac{1}{2},j} = & = \frac{U_{i,j} + U_{i-l,j}}{2} \\ \Phi_{i-\frac{1}{2},j} = & = \left[ \phi_{lin} - Cr^L Grad_{i,j}^L \frac{\Delta x}{2} + \left\{ \frac{\alpha_L}{2} - \frac{1}{6} \left( 1 - \left( Cr^L \right)^2 \right) \right\} Curv_{i,j}^L \Delta x^2 \right] \\ \phi_{lin} = & = \frac{\phi_{i,j} + \phi_{i-l,j}}{2} \\ Cr^L = & = \frac{U_{i-\frac{1}{2},j}}{\Delta x} \\ Grad_{i,j}^L = & = \frac{\Phi_{i,j} - \Phi_{i-l,j}}{\Delta x} \\ \alpha^L = & = \frac{D_{i-\frac{1}{2},j}}{(\Delta x)^2} = \frac{D_{i,j} + D_{i-l,j}}{2} * \frac{\Delta t}{(\Delta x)^2} \\ Curv_{i,j}^L = & = \frac{\phi_{i,j}^n - 2\phi_{i,j-1}^n + \phi_{i-2,j}^n}{\Delta x^2} \quad \text{jika } U_{l+\frac{1}{2},j} > 0 \\ = & = \frac{\phi_{i+l,j}^n - 2\phi_{i,j}^n + \phi_{i-l,j}^n}{\Delta x^2} \quad \text{jika } U_{l+\frac{1}{2},j} < 0 \end{split}\]

3.3.1 Fluks Difusif Akibat U

\[\begin{split} &D_{x}F\left(\Phi\right)_{i+\frac{1}{2},j}=&fluksdx_{i,j}^{R}\\ &D\hat{F}\left(\Phi\right)_{i+\frac{1}{2},j}=\Delta x.\alpha^{R}\cdot\frac{\partial\Phi^{*}}{\partial x}\bigg|_{i+\frac{1}{2},j}\\ &\left.\frac{\partial\Phi^{*}}{\partial x}\bigg|_{i+\frac{1}{2},j}=&Grad_{i,j}^{R}-Cr_{i,j}^{R}Curv_{i,j}^{R}\frac{\Delta x}{2}\\ &\alpha^{R}=\frac{D_{i+\frac{1}{2},j}^{\Delta t}}{\left(\Delta x\right)^{2}}=\frac{D_{i,j}+D_{i+1,j}}{2}*\frac{\Delta t}{\left(\Delta x\right)^{2}}\\ &Grad_{i,j}^{R}=\frac{\Phi_{i+1,j}-\Phi_{i,j}}{\Delta x}\\ &Cr_{i,j}^{R}=\frac{U_{i+\frac{1}{2},j}^{\Delta t}}{\Delta x}\\ &Curv_{i,j}^{R}=\frac{\Phi_{i+1,j}^{n}-2\Phi_{i,j}^{n}+\Phi_{i-1,j}^{n}}{\Delta x^{2}}\quad jika\;U_{i+\frac{1}{2},j}>0\\ &=\frac{\Phi_{i+2,j}^{n}-2\Phi_{i+1,j}^{n}+\Phi_{i,j}^{n}}{\Delta x^{2}}\quad jika\;U_{i+\frac{1}{2},j}<0 \end{split}\]

\[\begin{split} D_x F\left(\Phi\right)_{i-1/2,j} &= fluks dx_{i,j}^L \\ D\hat{F}\left(\Phi\right)_{i-1/2,j} &= \Delta x.\alpha^R \cdot \frac{\partial \Phi^*}{\partial x} \bigg|_{i-1/2,j} \\ &\frac{\partial \Phi^*}{\partial x} \bigg|_{i-1/2,j} &= Grad_{i,j}^L - Cr_{i,j}^L Curv_{i,j}^L \frac{\Delta x}{2} \\ &\alpha^L = \frac{D_{i-\frac{1}{2},j}^{-1} \Delta t}{\left(\Delta x\right)^2} = \frac{D_{i,j} + D_{i-1,j}}{2} * \frac{\Delta t}{\left(\Delta x\right)^2} \\ Grad_{i,j}^L &= \frac{\Phi_{i,j}^{-1} - \Phi_{i-1,j}}{\Delta x} \\ Cr_{i,j}^L &= \frac{U_{i-\frac{1}{2},j}^{-1} \Delta t}{\Delta x} \\ Curv_{i,j}^L &= \frac{\varphi_{i,j}^n - 2\varphi_{i-1,j}^n + \varphi_{i-2,j}^n}{\Delta x^2} \quad jika \ U_{1+\frac{1}{2},j} > 0 \\ &= \frac{\varphi_{i+1,j}^n - 2\varphi_{i,j}^n + \varphi_{i-1,j}^n}{\Delta x^2} \quad jika \ U_{1+\frac{1}{2},j} < 0 \end{split}\]

3.3.2 Fluks Konvektif Akibat V

\[\begin{split} F\left(\Phi\right)_{i,j+\frac{1}{2}} &= flukscy_{i,j}^{R} \\ &= V_{i,j+\frac{1}{2}} * \Phi_{i,j+\frac{1}{2}} \\ V_{i,j+\frac{1}{2}} &= \frac{V_{i,j} + V_{i,j+1}}{2} \\ \Phi_{i,j+\frac{1}{2}} &= \left[ \varphi_{lin} - Cr^{R} Grad_{i,j}^{R} \frac{\Delta y}{2} + \left\{ \frac{\alpha_{R}}{2} - \frac{1}{6} \left( 1 - \left( Cr^{R} \right)^{2} \right) \right\} Curv_{i,j}^{R} \Delta y^{2} \right] \\ \varphi_{lin} &= \frac{\varphi_{i,j} + \varphi_{i,j+1}}{2} \\ Cr^{R} &= \frac{V_{i,j+\frac{1}{2}} \Delta t}{\Delta y} \\ Grad_{i,j}^{R} &= \frac{\varphi_{i,j+1} - \varphi_{i,j}}{\Delta y} \\ \alpha^{R} &= \frac{D_{i,j+\frac{1}{2}} \Delta t}{\left( \Delta y \right)^{2}} = \frac{D_{i,j} + D_{i,j+1}}{2} * \frac{\Delta t}{\left( \Delta y \right)^{2}} \\ Curv_{i,j}^{R} &= \frac{\varphi_{i,j+1}^{n} - 2\varphi_{i,j}^{n} + \varphi_{i,j-1}^{n}}{\Delta y^{2}} \quad jika \ V_{i,j+\frac{1}{2}} > 0 \\ &= \frac{\varphi_{i,j+2}^{n} - 2\varphi_{i,j+1}^{n} + \varphi_{i,j}^{n}}{\Delta y^{2}} \quad jika \ V_{i,j+\frac{1}{2}} < 0 \end{split}\]

\[\begin{split} F\left(\Phi\right)_{i,j-\frac{1}{2}} &= flukscy_{i,j}^L \\ &= V_{i,j-\frac{1}{2}} * \Phi_{i-\frac{1}{2},j} \\ V_{i,j-\frac{1}{2}} &= \frac{V_{i,j} + V_{i,j-l}}{2} \\ \Phi_{i,j-\frac{1}{2}} &= \left[ \varphi_{lin} - Cr^L Grad_{i,j}^L \frac{\Delta y}{2} + \left\{ \frac{\alpha_L}{2} - \frac{1}{6} \left( 1 - \left( Cr^L \right)^2 \right) \right\} Curv_{i,j}^L \Delta y^2 \right] \\ \varphi_{lin} &= \frac{\varphi_{i,j} + \varphi_{i,j-l}}{2} \\ Cr^L &= \frac{V_{i,j-\frac{1}{2}} \Delta t}{\Delta y} \\ Grad_{i,j}^L &= \frac{\varphi_{i,j} - \varphi_{i,j-l}}{\Delta y} \\ \alpha^L &= \frac{D_{i,j-\frac{1}{2}} \Delta t}{(\Delta x) y^2} = \frac{D_{i,j} + D_{i,j-l}}{2} * \frac{\Delta t}{(\Delta y)^2} \\ Curv_{i,j}^L &= \frac{\varphi_{i,j}^n - 2\varphi_{i,j-1}^n + \varphi_{i,j-2}^n}{\Delta y^2} \quad jika \ V_{i,j-\frac{1}{2}} > 0 \\ &= \frac{\varphi_{i,j+l}^n - 2\varphi_{i,j}^n + \varphi_{i,j-l}^n}{\Delta y^2} \quad jika \ V_{i,j-\frac{1}{2}} < 0 \end{split}\]

3.3.3 Fluks Difusif Akibat V

\[\begin{split} D_{y}F\left(\Phi\right)_{i,j+\frac{1}{2}} &= fluksdy_{i,j}^{R} \\ D\hat{F}\left(\Phi\right)_{i,j+\frac{1}{2}} &= \Delta y.\alpha^{R} \cdot \frac{\partial \Phi}{\partial y} \bigg|_{i,j+\frac{1}{2}} \\ &\frac{\partial \Phi}{\partial y} \bigg|_{i,j+\frac{1}{2}} = Grad_{i,j}^{R} - Cr_{i,j}^{R}Curv_{i,j}^{R} \frac{\Delta y}{2} \\ &\alpha^{R} = \frac{D_{i,j+\frac{1}{2}}\Delta t}{\left(\Delta y\right)^{2}} = \frac{D_{i,j} + D_{i,j+1}}{2} * \frac{\Delta t}{\left(\Delta y\right)^{2}} \\ Grad_{i,j}^{R} &= \frac{\Phi_{i,j+1} - \Phi_{i,j}}{\Delta y} \\ &Cr_{i,j}^{R} &= \frac{U_{i,j+\frac{1}{2}}\Delta t}{\Delta y} \\ Curv_{i,j}^{R} &= \frac{\Phi_{i,j+1}^{n} - 2\Phi_{i,j}^{n} + \Phi_{i,j-1}^{n}}{\Delta y^{2}} \quad jika\ V_{i,j+\frac{1}{2}} > 0 \\ &= \frac{\Phi_{i,j+2}^{n} - 2\Phi_{i,j+1}^{n} + \Phi_{i,j}^{n}}{\Delta y^{2}} \quad jika\ V_{i,j+\frac{1}{2}} < 0 \end{split}\]

\[\begin{split} D_{y}F\left(\Phi\right)_{i,j-\frac{1}{2}} &= fluksdy_{i,j}^{L} \\ D\hat{F}\left(\Phi\right)_{i,j-\frac{1}{2}} &= \Delta y.\alpha^{R} \cdot \frac{\partial \Phi *}{\partial y} \bigg|_{i,j-\frac{1}{2}} \\ &\frac{\partial \Phi *}{\partial x} \bigg|_{i,j-\frac{1}{2}} &= Grad_{i,j}^{L} - Cr_{i,j}^{L} Curv_{i,j}^{L} \cdot \frac{\Delta y}{2} \\ &\alpha^{L} = \frac{D_{i,j-\frac{1}{2}}\Delta t}{\left(\Delta y\right)^{2}} = \frac{D_{i,j} + D_{i,j-1}}{2} * \frac{\Delta t}{\left(\Delta y\right)^{2}} \\ Grad_{i,j}^{L} &= \frac{\Phi_{i,j} - \Phi_{i,j-1}}{\Delta x} \\ &Cr_{i,j}^{L} &= \frac{U_{i,j-\frac{1}{2}}\Delta t}{\Delta y} \\ Curv_{i,j}^{L} &= \frac{\varphi_{i,j}^{n} - 2\varphi_{i,j-1}^{n} + \varphi_{i,j-2}^{n}}{\Delta y^{2}} \quad jika \ V_{i,j-\frac{1}{2}} > 0 \\ &= \frac{\varphi_{i,j+1}^{n} - 2\varphi_{i,j}^{n} + \varphi_{i,j-1}^{n}}{\Delta y^{2}} \quad jika \ V_{i,j-\frac{1}{2}} < 0 \end{split}\]

4 Hasil Pemodelan

Kalibrasi dan verifikasi dari model ini telah dilakukan dengan membandingkan hasil pemodelan dengan solusi analitis untuk kasus saluran lurus dan hasil uji lapangan untuk kasus aliran terbuka pada sebuah kolam. Pada proses tersebut juga dilakukan studi banding dengan hasil pemodelan berdasarkan model non turbulen [9]. Pada paper hanya akan disajikan komparasi hasil simulasi model turbulen κ-ε yang dikembangkan dengan hasil model non turbulen untuk kasus saluran lurus dan saluran ekspansi satu sisi (lihat Tabel 1). Kasus saluran lurus dilakukan untuk melihat efektifitas suku difusi dan adveksi dalam model. Kasus ekspansi satu sisi dilakukan untuk melihat peninkatan akurasi dari model dalam melingkup fenomena pembangkitan vorteks pada zona sirkulasi. Sketsa masingmasing kasus dapat dilihat pada uraian dalam subbab dibawah ini.

Tabel 1 Data-data masukan test.

4.1.1 Kasus Saluran Lurus

Gambar 1 Model Saluran Lurus.

Gambar 2 Vektor Kecepatan.

Pada kasus saluran lurus ini terlihat pergerakan salinitas sesuai dengan arah kecepatan, dengan kecepatan yang sama. Tampak bahwa pergerakan salinitas dipengaruhi oleh konveksi dan difusi, yang ditandai dengan pergerakan puncak dan penyebaran (distribusi) disertai penurunan puncak, dimana pergerakan puncak memiliki kecepatan yang sama dengan kecepatan aliran.

Gambar 3 Grafik Pergerakan Salinitas.

4.1.2 Kasus Ekspansi Satu Sisi

Gambar 4 Model saluran ekspansi tiba-tiba.

Gambar 5 Vektor kecepatan pada jam ke-7 untuk kasus non turbulen (a) dan turbulen (b).

Gambar 6 Kontur Kecepatan pada jam ke 7 untuk kasus non turbulen (a) dan turbulen (b).

Gambar 7 Kontur Salinitas pada jam ke -7 untuk kasus non turbulen (a) dan kasus turbulen (b).

Pada saluran ekspansi terlihat turbulensi terbentuk akibat ekspansi tiba-tiba pada dinding batas saluran. Pada model turbulen pola turbulensi menguat, terdapat peningkatan kecepatan, serta pembentukan pusaran-pusaran baru yang ukurannya lebih kecil terlihat lebih jelas dan stabil dibandingkan dengan hasil model non turbulen. Produksi turbulensi yang sangat tinggi (fenomena bursting) terjadi secara repetitif dan berfluktuasi pada zona resirkulasi dengan frekuensi tertentu. Fenomena bursting ini ditandai dengan terbentuknya vorteks dan terjadinya proses resirkulasi. Pusaran-pusaran yang besar terbentuk akibat kondisi batas (dalam hal ini dinding), sedangkan pusaran-pusaran yang kecil terbentuk oleh viscous force [14].

Pada aliran turbulen pemisahan struktur pola arus/kecepatan di tengah dan bagian sayap saluran terlihat sangat signifikan, begitu pula pusaran arus menjadi bentuk pola arus yang dominan di bagian sayap saluran. Pada aliran non turbulen pola pusaran arus di bagian sayap terlihat lebih lemah. Pola aliran sangat berpengaruh pada kecepatan penyebaran salinitas. Adanya turbulensi menyebabkan salinitas lebih cepat terdistribusi ke seluruh bagian saluran, jika dibandingkan dengan penyebaran salinitas pada aliran non turbulen. Hal ini disebabkan sifat difusif aliran turbulen yang mempercepat proses pencampuran (mixing) berbagai kuantitas, termasuk massa (konsentrasi).

5 Kesimpulan

Pada kasus ekspansi tiba-tiba, model turbulen Kappa-Epsilon mampu memberikan pola arus yang lebih baik dibandingkan model non turbulen. Hal ini terlihat dengan munculnya fenomena bursting dan resirkulasi yang tidak dapat dihasilkan oleh model non turbulen. Selain itu, model turbulen κ-ε telah dapat memberikan gambaran pola aliran yang lebih stabil pada aliran turbulen.

Selain itu model turbulen κ-ε juga memberikan pola penyebaran salinitas (konvektif dan difusif) yang lebih baik dan akurat dibandingkan model non turbulen.

Acknowledgement

Penelitian ini dilakukan dengan biaya DIKTI melalui program hibah Bersaing XII dan Hibah Pasca Sarjana P2M DIKTI 2007.