5-Equation Model 직접 구현하기

이론만으로는 감이 오지 않는 법입니다. 이 글에서는 1D 5-equation diffuse interface model을 Python으로 직접 구현하고, 충격파가 기체 버블을 관통하는 문제를 시뮬레이션합니다.

지배 방정식 정리

1D 5-equation model의 보존 변수와 플럭스:

\mathbf{U} = \begin{bmatrix} \alpha_1 \rho_1 \\ \alpha_2 \rho_2 \\ \rho u \\ E \\ \alpha_1 \end{bmatrix}, \quad \mathbf{F} = \begin{bmatrix} \alpha_1 \rho_1 u \\ \alpha_2 \rho_2 u \\ \rho u^2 + p \\ (E+p)u \\ 0 \end{bmatrix}

5번째 방정식 ( $\alpha_1$ )은 비보존형이므로 별도 처리합니다:

\frac{\partial \alpha_1}{\partial t} + u \frac{\partial \alpha_1}{\partial x} = 0

상태방정식: Stiffened Gas EOS

두 유체 모두 stiffened gas EOS를 사용합니다:

p_k = (\gamma_k - 1)\rho_k e_k - \gamma_k p_{\infty,k}

혼합 규칙(mixture rule)으로 전체 압력을 구합니다:

p = \frac{\rho e - \sum_k \alpha_k \frac{\gamma_k p_{\infty,k}}{\gamma_k - 1}}{\sum_k \frac{\alpha_k}{\gamma_k - 1}}

Python 구현

초기 설정과 EOS

import numpy as np
import matplotlib.pyplot as plt
 
# Domain
N = 1000
x = np.linspace(0, 1, N)
dx = x[1] - x[0]
 
# EOS parameters
# Fluid 1: Air (ideal gas)
gamma1, pinf1 = 1.4, 0.0
# Fluid 2: Water (stiffened gas)
gamma2, pinf2 = 4.4, 6e8
 
def pressure(alpha1, arho1, arho2, E, u):
    """Compute mixture pressure from conserved variables."""
    alpha2 = 1.0 - alpha1
    rho = arho1 + arho2
    ke = 0.5 * rho * u**2
    rhoe = E - ke
 
    num = rhoe - alpha1*gamma1*pinf1/(gamma1-1) - alpha2*gamma2*pinf2/(gamma2-1)
    den = alpha1/(gamma1-1) + alpha2/(gamma2-1)
    return num / den
 
def sound_speed(alpha1, arho1, arho2, p):
    """Compute mixture sound speed (Wood's formula)."""
    alpha2 = 1.0 - alpha1
    rho = arho1 + arho2
    rho1 = arho1 / np.maximum(alpha1, 1e-12)
    rho2 = arho2 / np.maximum(alpha2, 1e-12)
 
    c1sq = gamma1 * (p + pinf1) / rho1
    c2sq = gamma2 * (p + pinf2) / rho2
 
    # Wood's formula
    inv_rhoc2 = alpha1/(rho1*c1sq) + alpha2/(rho2*c2sq)
    return np.sqrt(1.0 / (rho * inv_rhoc2))

HLLC Flux 계산

def hllc_flux(UL, UR, alpha1L, alpha1R):
    """HLLC numerical flux for the 4 conservation equations."""
    arho1L, arho2L, momL, EL = UL
    arho1R, arho2R, momR, ER = UR
 
    rhoL = arho1L + arho2L
    rhoR = arho1R + arho2R
    uL = momL / rhoL
    uR = momR / rhoR
    pL = pressure(alpha1L, arho1L, arho2L, EL, uL)
    pR = pressure(alpha1R, arho1R, arho2R, ER, uR)
    cL = sound_speed(alpha1L, arho1L, arho2L, pL)
    cR = sound_speed(alpha1R, arho1R, arho2R, pR)
 
    # Wave speed estimates
    SL = min(uL - cL, uR - cR)
    SR = max(uL + cL, uR + cR)
    Sstar = (pR - pL + rhoL*uL*(SL-uL) - rhoR*uR*(SR-uR)) \
          / (rhoL*(SL-uL) - rhoR*(SR-uR))
 
    FL = np.array([arho1L*uL, arho2L*uL, momL*uL + pL, (EL+pL)*uL])
    FR = np.array([arho1R*uR, arho2R*uR, momR*uR + pR, (ER+pR)*uR])
 
    if SL >= 0:
        return FL, Sstar
    elif Sstar >= 0:
        coeff = rhoL * (SL - uL) / (SL - Sstar)
        UstarL = coeff * np.array([
            arho1L/rhoL,
            arho2L/rhoL,
            Sstar,
            EL/rhoL + (Sstar - uL)*(Sstar + pL/(rhoL*(SL-uL)))
        ])
        return FL + SL*(UstarL - UL), Sstar
    elif SR > 0:
        coeff = rhoR * (SR - uR) / (SR - Sstar)
        UstarR = coeff * np.array([
            arho1R/rhoR,
            arho2R/rhoR,
            Sstar,
            ER/rhoR + (Sstar - uR)*(Sstar + pR/(rhoR*(SR-uR)))
        ])
        return FR + SR*(UstarR - UR), Sstar
    else:
        return FR, Sstar

시간 적분 (1st order)

def solve(U, alpha1, dx, t_end, cfl=0.5):
    """Forward Euler time integration."""
    t = 0.0
    while t < t_end:
        rho = U[0] + U[1]
        u = U[2] / rho
        p = pressure(alpha1, U[0], U[1], U[3], u)
        c = sound_speed(alpha1, U[0], U[1], p)
 
        dt = cfl * dx / np.max(np.abs(u) + c)
        if t + dt > t_end:
            dt = t_end - t
 
        # Compute fluxes at cell interfaces
        flux = np.zeros((4, N+1))
        uface = np.zeros(N+1)
 
        for i in range(N+1):
            iL = max(i-1, 0)
            iR = min(i, N-1)
            UL_local = U[:4, iL]
            UR_local = U[:4, iR]
            flux[:, i], uface[i] = hllc_flux(
                UL_local, UR_local, alpha1[iL], alpha1[iR]
            )
 
        # Update conserved variables
        for k in range(4):
            U[k, 1:-1] -= dt/dx * (flux[k, 2:-1] - flux[k, 1:-2])
 
        # Update alpha1 (non-conservative, upwind)
        for i in range(1, N-1):
            if uface[i] >= 0:
                alpha1[i] -= dt/dx * uface[i] * (alpha1[i] - alpha1[i-1])
            else:
                alpha1[i] -= dt/dx * uface[i] * (alpha1[i+1] - alpha1[i])
 
        alpha1 = np.clip(alpha1, 1e-8, 1-1e-8)
        t += dt
 
    return U, alpha1

초기 조건: 충격파-기체 버블

# Initial condition: shock hitting an air bubble in water
U = np.zeros((4, N))
alpha1 = np.zeros(N)
 
for i in range(N):
    if x[i] < 0.25:
        # Post-shock water (high pressure)
        rho2 = 1100.0; u_val = 50.0; p_val = 1e9
        alpha1[i] = 1e-8
        arho1 = alpha1[i] * 1.0  # negligible air
        arho2 = (1 - alpha1[i]) * rho2
    elif 0.4 < x[i] < 0.6:
        # Air bubble
        rho1 = 1.0; u_val = 0.0; p_val = 1e5
        alpha1[i] = 1 - 1e-8
        arho1 = alpha1[i] * rho1
        arho2 = (1 - alpha1[i]) * 1000.0
    else:
        # Ambient water
        rho2 = 1000.0; u_val = 0.0; p_val = 1e5
        alpha1[i] = 1e-8
        arho1 = alpha1[i] * 1.0
        arho2 = (1 - alpha1[i]) * rho2
 
    rho = arho1 + arho2
    E = compute_total_energy(alpha1[i], arho1, arho2, p_val, u_val)
 
    U[0, i] = arho1
    U[1, i] = arho2
    U[2, i] = rho * u_val
    U[3, i] = E

주의사항

1. $\alpha$ 의 범위 제한

$\alpha_1$ 이 정확히 0이나 1이 되면 EOS 계산에서 division by zero가 발생합니다. 항상 $\epsilon < \alpha_1 < 1 - \epsilon$ ( $\epsilon \sim 10^{-8}$ )으로 clipping합니다.

2. 혼합 음속: Wood의 공식

혼합 영역에서의 음속은 직관과 다르게 매우 낮아질 수 있습니다. 이는 물리적으로 올바른 현상이며, CFL 조건에 영향을 줍니다.

\frac{1}{\rho c^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}

$\alpha_1 = 0.5$ 근처에서 $c_{mix}$ 는 두 유체의 음속보다 훨씬 작아집니다.

3. 고차 정확도로의 확장

위 코드는 1차 정확도입니다. 실제 연구에서는:

MUSCL 재구성 + slope limiter로 2차
WENO 재구성으로 5차 이상
Runge-Kutta 시간 적분으로 시간 정확도 향상

결과 해석

시뮬레이션을 실행하면 충격파가 기체 버블을 압축하면서:

**투과 충격파(transmitted shock)**가 버블 뒤쪽으로 진행
**반사 팽창파(reflected rarefaction)**가 좌측으로 이동
버블이 급격히 압축되며 내부 압력 상승
후방 벽에서 제트(jet) 형성 (2D/3D에서)

이 현상은 수중 폭발, 충격파 쇄석술 등에서 핵심적인 물리입니다.

다음 단계

다음 글에서는 수치해석의 기본 개념들을 게임처럼 배워보는 인터랙티브 콘텐츠를 준비했습니다. 수치 안정성, CFL 조건 등을 직접 체험해볼 수 있습니다.

5-Equation Model 구현 가이드: Python으로 만드는 다상유동 솔버