Pemodelan Matematis & Algoritma Komputasi untuk Project Teknik
Kumpulan Pemodelan Matematis dan Algoritma Komputasi Project Komputasi Teknik
Dokumen ini merangkum pemodelan matematis dan algoritma komputasi yang dapat digunakan untuk Project Komputasi Teknik, terutama yang mencakup:
- CFD 2D airfoil NACA 4112 dengan validasi terhadap NACA 4412,
- Metode Finite Volume (FVM) untuk konduksi panas 1D,
- Implementasi numerik seperti SIMPLE/SIMPLEC/PISO, TDMA, GCI, dan validasi,
- Simulasi 1D di OpenModelica sebagai contoh implementasi komputasi.
1. Peta Model dan Algoritma
| Komponen | Model Matematis | Algoritma/Numerik | Output |
|---|---|---|---|
| Aliran fluida umum | Persamaan konservasi massa, momentum, energi | CFD, FVM | Medan aliran, tekanan, kecepatan |
| Aliran turbulen | RANS + model turbulensi SST (k-\omega) | Solver pressure-based, SIMPLEC | (C_L), (C_D), distribusi (C_p) |
| Airfoil NACA 4112 | Geometri NACA 4-digit + RANS steady-incompressible | Mesh C, validasi, GCI | Karakteristik aerodinamika |
| Konduksi panas 1D | Persamaan difusi 1D | FVM, Central Differencing, TDMA | Distribusi temperatur |
| Simulasi 1D | Persamaan diskrit keseimbangan energi | OpenModelica, for-loop | Temperatur tiap sel |
| Validasi | Error relatif, GCI | Mesh convergence | Tingkat akurasi model |
2. Pemodelan Matematis Umum CFD
2.1 Persamaan Konservasi
Persamaan kontinuitas:
[
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0
]
Persamaan momentum NavierโStokes:
[
\frac{\partial (\rho \mathbf{u})}{\partial t}
- \nabla \cdot (\rho \mathbf{u}\mathbf{u})
= -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g}
]
Persamaan energi:
[
\frac{\partial (\rho E)}{\partial t}
- \nabla \cdot (\rho \mathbf{u} E)
= \nabla \cdot (k \nabla T) + \Phi + S_E
]
Untuk aliran steady-incompressible yang digunakan pada simulasi airfoil:
[
\nabla \cdot \mathbf{u} = 0
]
[
\nabla \cdot (\mathbf{u}\mathbf{u})
= -\frac{1}{\rho}\nabla p + \nabla \cdot \left(\nu_{\text{eff}} \nabla \mathbf{u}\right)
]
dengan (\nu_{\text{eff}} = \nu + \nu_t).
2.2 Persamaan Transport Umum
Persamaan transport skalar (\phi):
[
\underbrace{\frac{\partial (\rho \phi)}{\partial t}}{\text{transien}}
+
\underbrace{\nabla \cdot (\rho \mathbf{u}\phi)}{\text{konvektif}}
\underbrace{\nabla \cdot (\Gamma \nabla \phi)}{\text{difusif}}
+
\underbrace{S}{\text{source}}
]
Empat terminologi utama:
- Transien: perubahan terhadap waktu.
- Konvektif: perpindahan akibat aliran fluida.
- Difusif: perpindahan akibat gradien.
- Source: sumber/sink eksternal.
Untuk kasus steady, transien nol. Untuk benda padat tanpa aliran, konvektif nol. Untuk konduksi 1D tanpa sumber, source nol.
2.3 Model Turbulensi SST (k-\omega)
Model SST (k-\omega) dipilih karena mampu menangani adverse pressure gradient dan separasi lapisan batas. Bentuk umum:
[
\frac{\partial (\rho k)}{\partial t}
- \nabla \cdot (\rho \mathbf{u} k) \nabla \cdot (\Gamma_k \nabla k)
- P_k – \beta^* \rho k \omega
]
[
\frac{\partial (\rho \omega)}{\partial t}
- \nabla \cdot (\rho \mathbf{u} \omega) \nabla \cdot (\Gamma_\omega \nabla \omega)
- \frac{\alpha P_k}{\nu_t}
- \beta \rho \omega^2
- S_\omega
]
Konstanta closure mengikuti rekomendasi Menter. Model ini menggabungkan keunggulan (k-\omega) di dekat dinding dan (k-\varepsilon) di daerah free-shear.
3. Model 1D Konduksi Panas dan Diskretisasi FVM
3.1 Persamaan Dasar
Untuk konduksi panas 1D steady tanpa sumber:
[
\frac{d}{dx}\left(\Gamma \frac{d\phi}{dx}\right) + S = 0
]
Dengan (\phi = T) dan (\Gamma = k):
[
\frac{d}{dx}\left(k \frac{dT}{dx}\right) + S = 0
]
Jika (k) konstan dan (S=0):
[
\frac{d^2T}{dx^2} = 0
]
3.2 Diskretisasi Finite Volume
Integrasi persamaan pada volume kontrol:
[
\int_{CV} \frac{d}{dx}\left(\Gamma \frac{d\phi}{dx}\right)dV
+
\int_{CV} S,dV = 0
]
Bentuk diskrit umum:
[
a_P \phi_P = a_W \phi_W + a_E \phi_E + b
]
Untuk difusi murni 1D dengan grid seragam:
[
a_E = \frac{\Gamma_e A}{\delta x_{PE}}, \quad
a_W = \frac{\Gamma_w A}{\delta x_{WP}}
]
[
a_P = a_E + a_W – S_P \Delta V
]
[
b = S_u \Delta V
]
Untuk material homogen dan grid seragam:
[
a_E = a_W = \frac{kA}{\Delta x}
]
[
a_P = a_E + a_W
]
Source term dilinearisasi:
[
S = S_u + S_P \phi_P
]
dengan (S_P < 0) untuk stabilitas numerik.
3.3 Penanganan Kondisi Batas Temperatur Tetap
Jika sel batas berbatasan langsung dengan dinding bertemperatur (T_B), konduktansi dari node ke batas adalah:
[
G_B = \frac{kA}{\Delta x/2} = \frac{2kA}{\Delta x}
]
Kontribusi kondisi batas masuk ke suku sumber:
[
b \leftarrow b + G_B T_B
]
dan menambah koefisien diagonal:
[
a_P \leftarrow a_P + G_B
]
Untuk batas konveksi:
[
q” = h(T_B – T_P)
]
Maka:
[
S_u = hA T_B, \quad S_P = -hA
]
3.4 Algoritma FVM 1D
Pseudocode:
1. Definisikan domain, properti material, kondisi batas, jumlah sel N.
2. Hitung lebar sel ฮx = L/N.
3. Inisialisasi koefisien aW, aE, aP, b untuk setiap sel.
4. Untuk setiap sel i:
a. Jika sel batas kiri/kanan:
terapkan kontribusi kondisi batas ke aP dan b.
b. Jika sel internal:
aW[i] = kA/ฮx
aE[i] = kA/ฮx
aP[i] = aW[i] + aE[i] - Sp[i]
b[i] = Su[i]
5. Susun sistem persamaan linear:
aP[i] T[i] = aW[i] T[i-1] + aE[i] T[i+1] + b[i]
6. Selesaikan dengan TDMA atau metode matriks.
7. Jika nonlinear/iteratif, ulangi sampai konvergen.
8. Post-processing: plot distribusi temperatur.3.5 Implementasi OpenModelica
Model latihan2:
model latihan2
parameter Real K=1000;
parameter Real HL=10;
parameter Real TP=10;
parameter Real TB=100;
parameter Real L=0.5;
parameter Real A=0.05;
parameter Integer N=5;
parameter Real deltaX=L/N;
Real T[N];
Real Tw[N], Te[N];
Real Aw[N], Ae[N], Ap[N], Sp[N], Su[N];
equation
// Sel batas kiri i=1
Aw[1] = HL*A;
Ae[1] = K*A/deltaX;
Sp[1] = -HL*A;
Su[1] = HL*A*TP;
Ap[1] = Aw[1] + Ae[1] - Sp[1];
Tw[1] = TP;
Te[1] = T[2];
// Sel internal
for i in 2:N-1 loop
Aw[i] = K*A/deltaX;
Ae[i] = K*A/deltaX;
Sp[i] = 0;
Su[i] = 0;
Ap[i] = Aw[i] + Ae[i] - Sp[i];
Tw[i] = T[i-1];
Te[i] = T[i+1];
end for;
// Sel batas kanan i=N
Aw[N] = K*A/deltaX;
Ae[N] = HL*A;
Sp[N] = -HL*A;
Su[N] = HL*A*TB;
Ap[N] = Aw[N] + Ae[N] - Sp[N];
Tw[N] = T[N-1];
Te[N] = TB;
// Persamaan keseimbangan energi
for i in 1:N loop
Ap[i]*T[i] = Aw[i]*Tw[i] + Ae[i]*Te[i] + Su[i];
end for;
end latihan2;4. Pemodelan Matematis Airfoil NACA 4112
4.1 Geometri NACA 4-Digit
NACA 4112:
- (m = 0.04) โ camber maksimum 4% chord,
- (p = 0.10) โ posisi camber maksimum 10% chord,
- (t = 0.12) โ ketebalan maksimum 12% chord.
Distribusi camber:
[
y_c(x) =
\begin{cases}
\dfrac{m}{p^2}\left(2px – x^2\right), & 0 \le x \le p \[8pt]
\dfrac{m}{(1-p)^2}\left((1-2p) + 2px – x^2\right), & p < x \le 1
\end{cases}
]
Distribusi ketebalan:
[
y_t(x) = 5t\left(
0.2969\sqrt{x}
-0.1260x
-0.3516x^2
+0.2843x^3
-0.1015x^4
\right)
]
Kemiringan camber:
[
\theta = \tan^{-1}\left(\frac{dy_c}{dx}\right)
]
Koordinat permukaan atas:
[
x_u = x – y_t \sin\theta, \quad
y_u = y_c + y_t \cos\theta
]
Koordinat permukaan bawah:
[
x_l = x + y_t \sin\theta, \quad
y_l = y_c – y_t \cos\theta
]
4.2 Model Aliran
Simulasi menggunakan:
- RANS steady-incompressible,
- model turbulensi SST (k-\omega),
- aliran 2D,
- fluida Newtonian,
- dinding no-slip, smooth, adiabatic,
- tanpa perpindahan panas,
- turbulen penuh tanpa model transisi.
Bilangan Reynolds:
[
Re = \frac{\rho U_\infty c}{\mu}
]
Kecepatan inlet:
[
U_\infty = \frac{Re \mu}{\rho c}
]
Untuk sudut serang (\alpha):
[
u_{\text{inlet}} = U_\infty \cos\alpha, \quad
v_{\text{inlet}} = U_\infty \sin\alpha
]
4.3 Domain dan Mesh
- Chord (c = 1) m.
- Topologi C-mesh.
- Inlet 10c, outlet 20c, batas atas/bawah 10c.
- Mesh quadrilateral terstruktur.
- Inflation layer 20โ25, growth ratio 1.10โ1.15.
- Tiga level mesh: coarse, medium, fine.
- Uji mesh convergence dengan GCI.
- Target (y^+) sesuai model dinding/turbulensi.
4.4 Kondisi Batas
| Batas | Jenis | Nilai |
|---|---|---|
| Inlet | Velocity inlet | (U_\infty), (\alpha), turbulensi 1%, viscosity ratio 10 |
| Outlet | Pressure outlet | 0 Pa gauge |
| Airfoil | Wall | No-slip, smooth, adiabatic |
| Atas/bawah | Farfield | Symmetry atau velocity inlet |
4.5 Metode Numerik
- Solver pressure-based.
- Pressure-velocity coupling: SIMPLEC.
- Skema momentum dan turbulensi: second-order upwind.
- Gradien: least squares cell-based.
- Residual kontinuitas, momentum, (k), (\omega) di bawah ambang konvergensi.
- Monitor (C_L), (C_D), (C_M) stabil: perubahan <0,1% dalam 200 iterasi, minimal 1000 iterasi.
4.6 Hasil Validasi dan Prediksi
Validasi terhadap NACA 4412:
| Parameter | Eksperimen | CFD SST (k-\omega) | Error relatif |
|---|---|---|---|
| (C_L) | 0,84 | 0,86 | 2,4% |
| (C_D) | 0,0105 | 0,0118 | 12,4% |
| (C_p) permukaan atas | โ | โ | 2,6% |
Mesh convergence:
| Mesh | (C_L) | (C_D) | GCI |
|---|---|---|---|
| Coarse | 0,845 | 0,0112 | โ |
| Medium | 0,857 | 0,0117 | 3,5% |
| Fine | 0,861 | 0,0119 | 1,7% |
Prediksi NACA 4112:
[
C_L \approx 0,94, \quad
C_D \approx 0,0126, \quad
\frac{C_L}{C_D} \approx 74,6
]
Karakteristik:
- NACA 4112 menghasilkan (C_L) lebih tinggi daripada NACA 4412.
- Suction peak lebih tajam dan lebih maju karena camber maksimum pada 10% chord.
- Pressure recovery lebih cepat.
- Risiko adverse pressure gradient dan separasi lebih tinggi pada sudut serang besar.
- Aliran masih attached pada kasus yang disimulasikan.
5. Algoritma Komputasi CFD 2D
5.1 Alur Pre-processing
1. Bangkitkan geometri airfoil NACA 4112.
2. Buat domain komputasi C-mesh.
3. Lakukan mesh generation dengan refinement near-wall.
4. Tetapkan kondisi batas.
5. Hitung Uโ dari Re.
6. Definisikan model turbulensi SST k-ฯ.
7. Atur skema numerik dan kriteria konvergensi.5.2 Alur Solver
1. Inisialisasi medan u, v, p, k, ฯ.
2. Hitung flux dan koefisien diskrit untuk persamaan momentum.
3. Selesaikan persamaan momentum.
4. Selesaikan persamaan tekanan/koreksi tekanan.
5. Koreksi kecepatan dan tekanan.
6. Selesaikan persamaan transport k dan ฯ.
7. Perbarui properti turbulensi.
8. Periksa residual dan monitor C_L, C_D, C_M.
9. Jika belum konvergen, ulangi langkah 2โ8.
10. Post-processing: kontur C_p, streamline, vektor kecepatan.5.3 Algoritma SIMPLE/SIMPLEC/PISO
- SIMPLE: iteratif, memprediksi tekanan (p^*), menyelesaikan momentum, menghitung koreksi tekanan (p’), lalu mengoreksi (p), (u), (v).
- SIMPLEC: versi SIMPLE dengan koreksi tekanan-kecepatan yang lebih konsisten; digunakan pada proyek NACA 4112.
- PISO: cocok untuk aliran unsteady atau kopling tekanan-kecepatan yang lebih kompleks dengan splitting operator.
5.4 TDMA
Sistem persamaan tri-diagonal:
[
a_i \phi_i = b_i \phi_{i+1} + c_i \phi_{i-1} + d_i
]
Transformasi:
[
\phi_i = P_i \phi_{i+1} + Q_i
]
dengan:
[
P_i = \frac{b_i}{a_i – c_i P_{i-1}}
]
[
Q_i = \frac{d_i + c_i Q_{i-1}}{a_i – c_i P_{i-1}}
]
Lalu back-substitution:
[
\phi_i = P_i \phi_{i+1} + Q_i
]
TDMA efisien untuk sistem 1D dan dapat diterapkan secara line-by-line pada 2D/3D.
5.5 Validasi dan GCI
Error relatif:
[
\varepsilon = \frac{|X_{\text{CFD}} – X_{\text{exp}}|}{|X_{\text{exp}}|} \times 100%
]
GCI fine mesh:
[
p = \frac{\ln\left|\frac{f_3 – f_2}{f_2 – f_1}\right|}{\ln r}
]
[
GCI_{\text{fine}} =
F_s \frac{\left|\frac{f_2 – f_1}{f_1}\right|}{r^p – 1}
\times 100%
]
dengan (F_s = 1{,}25) dan (r) adalah rasio refinement. Target GCI < 2%. Proyek NACA 4112 memperoleh GCI fine sekitar 1,7%, sehingga mesh fine dianggap memadai.
8. Takeaway Utama
Pemodelan matematis project ini bertumpu pada:
- persamaan konservasi massa, momentum, dan energi,
- persamaan transport umum,
- RANS dengan model turbulensi SST (k-\omega),
- geometri NACA 4-digit,
- diskretisasi FVM.
Algoritma komputasinya mencakup:
- pembentukan mesh dan kondisi batas,
- diskretisasi volume hingga,
- penyelesaian sistem linear dengan TDMA,
- kopling tekanan-kecepatan SIMPLE/SIMPLEC/PISO,
- uji konvergensi mesh dengan GCI,
- validasi terhadap data eksperimen,
- implementasi 1D di OpenModelica.
Dengan demikian, project ini menunjukkan jembatan antara teori analitik, eksperimen, dan komputasi numerik dalam menyelesaikan masalah fluida dan perpindahan panas.
