시간간격을 조이는 음속은 그 24 m/s가 아니라, 충격파 뒤 물의 약 2530 m/s다.

5-Equation Model 직접 구현하기#

충격파 뒤 물의 압력은 $10^{9}$ , 앞의 물은 $10^{5}$ 다. 그 사이에 공기 버블이 있고, 체적분율 $\alpha_1$ 은 거의 0과 거의 1을 오간다. 1D 5-equation diffuse interface model은 이 점프를 식 다섯 개로 민다. 두 유체는 속도 하나, 압력 하나를 나눠 갖고, 체적분율만 따로 흐른다.

인공 확산 슬라이더를 0으로 내리면 체적분율 경계가 한 칸에서 끊긴다.

속도 u 1.00인공 확산 0.020초기 폭 0.20

5-방정식 모델의 핵심: 부피분율 α의 이류. upwind 만으로는 계면이 점차 번지므로 실제 코드는 THINC / interface sharpening을 추가한다.

슬라이더를 올리면 경계가 여러 칸으로 퍼진다. 그 칸의 $\alpha$ 가 혼합 압력의 분모에 그대로 들어간다.

지배 방정식 정리#

1D에서 보존 변수와 플럭스는 이렇다.

\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}

다섯 번째 줄의 플럭스는 0이다. $\alpha_1$ 은 비보존형이라 따로 민다.

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

한 속도와 한 압력으로 줄인 모델이 이 다섯 식이다. 각 상에 방정식을 되돌리는 쪽은 Baer–Nunziato 7-방정식이다.

상태방정식: Stiffened Gas EOS#

두 유체 모두 stiffened gas다.

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

혼합 압력은 이 규칙으로 구한다.

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}}

$p_k$ 는 $k$ 유체의 압력, $e_k$ 는 그 내부에너지, $\gamma_k$ 와 $p_{\infty,k}$ 는 EOS 상수다. 분모는 $\sum \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))

pressure는 위 혼합 규칙을 보존 변수로 다시 쓴 것이다. sound_speed의 $c_k^2=\gamma_k(p+p_{\infty,k})/\rho_k$ 가 순수 유체 음속이다.

HLLC Flux 계산#

보존 식 네 개의 수치 플럭스만 HLLC다. $\alpha$ 의 이류 속도 $S^\star$ 는 같이 돌려준다.

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)#

시간 전진은 전방 오일러다. $\Delta t$ 는 $\mathrm{CFL}\,\Delta x/\max(|u|+c)$ 이고, 기본 CFL은 0.5다.

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

$\alpha_1$ 갱신은 보존 플럭스가 아니다. 면 속도 uface의 부호로 풍상을 고른다.

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

위치 0.25보다 앞은 충격파를 지난 물이다. $\rho=1100$ , $u=50$ , $p=10^{9}$ . 0.4와 0.6 사이는 공기 버블, $\rho=1$ , $u=0$ , $p=10^{5}$ . 나머지는 주변 물, $\rho=1000$ , $u=0$ , $p=10^{5}$ 다.

# 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

compute_total_energy는 이 조각 안에 없다. 압력에서 에너지를 되돌리는 식은 EOS 쪽에 있다. $\rho e=(p+\gamma p_\infty)/(\gamma-1)$ 이다.

주의사항#

1. $\alpha$ 의 범위 제한#

$\alpha_1$ 이 정확히 0이나 1이면 보존량 $\alpha_1\rho_1$ 을 $\alpha_1$ 로 나누지 못한다. 코드는 $\alpha_1$ 을 $10^{-8}$ 과 $1-10^{-8}$ 사이에 자른다.

2. 혼합 음속: Wood의 공식#

혼합 칸의 음속은 두 유체보다 낮아질 수 있다. 직관과 반대지만, 식은 그렇게 말한다. CFL도 이 음속을 본다.

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

초기 조건의 순수 상태를 이 식에 넣으면 숫자가 나온다. 버블 공기는 $\rho=1$ , $p=10^{5}$ , $\gamma=1.4$ , $p_\infty=0$ 이다. $c^2=\gamma(p+p_\infty)/\rho$ 이면 $c=\sqrt{1.4\times 10^{5}}\approx 374\,\mathrm{m/s}$ 다. 주변 물은 $\rho=1000$ , $p=10^{5}$ , $\gamma=4.4$ , $p_\infty=6\times 10^{8}$ 이라 $c\approx 1625\,\mathrm{m/s}$ 다. 충격파 뒤 물은 $\rho=1100$ , $p=10^{9}$ 라 $c\approx 2530\,\mathrm{m/s}$ 다. $p+p_\infty$ 는 주변 물 대비 약 2.67배인데, 밀도는 $1100/1000=1.10$ 배뿐이다. $u=50$ 의 운동에너지는 $0.5\times 1100\times 50^{2}=1.375\times 10^{6}$ 이다. 같은 칸의 내부에너지 밀도 $\rho e\approx 1.07\times 10^{9}$ 의 약 0.13%다. 이 충격의 에너지는 속도보다 $p_\infty$ 쪽에 있다.

$\alpha_1=0.5$ 로 버블 공기와 주변 물을 한 압력 $10^{5}$ 에서 섞으면, Wood 식의 $1/(\rho c^2)$ 은 거의 공기 항이다. $\rho_1 c_1^2=1.4\times 10^{5}$ 이고 $\rho_2 c_2^2\approx 2.64\times 10^{9}$ 다. 혼합 밀도는 $0.5\times 1+0.5\times 1000=500.5$ 이고, $c_\mathrm{mix}\approx 24\,\mathrm{m/s}$ 다. 374도 아니고 1625도 아니다.

압력 분모도 같은 자리에서 민감하다. 순수 물의 $\sum \alpha/(\gamma-1)$ 은 $1/3.4\approx 0.294$ 다. 물에 $\alpha_1=0.01$ 만 새어도 분모는 $0.01/0.4+0.99/3.4\approx 0.316$ 으로 약 7.5% 커진다. 압력이 이 분모로 나뉘므로, 퍼진 계면의 $\alpha$ 오차는 압력 오차로 들어온다.

CFL의 $\Delta t=0.5\,\Delta x/\max(|u|+c)$ 이고, $\Delta x$ 는 linspace(0, 1, 1000)의 간격 $1/999$ 다. 제한은 혼합 칸이 아니다. 충격파 뒤는 $|u|+c\approx 50+2530=2580$ 이라 $\Delta t\approx 1.9\times 10^{-7}$ 이다. 음속 24 m/s인 칸만 보면 간격이 약 110배 크다. 혼합 음속이 낮아져도 시간간격은 그 칸이 조이지 않는다. 같은 $\Delta t$ 동안 계면의 음향이 거의 멈춰 보이는 쪽이 영향이다.

3. 고차 정확도로의 확장#

위 시간 적분은 1차다. 공간도 칸 값을 그대로 쓴다. 2차가 필요하면 MUSCL과 slope limiter, 그 위는 WENO다. 시간은 Runge–Kutta로 올린다. 차수를 올려도 $\alpha$ 의 비보존 갱신과 Wood 음속은 그대로 남는다.

결과 해석#

충격이 버블에 들어가면 뒤쪽으로 투과 충격파(transmitted shock)가 나가고, 왼쪽으로 반사 팽창파(reflected rarefaction)가 돌아온다. 버블은 급히 눌리며 안쪽 압력이 오른다. 제트(jet)는 2D나 3D에서 뒤쪽 벽에 생긴다. 수중 폭발과 충격파 쇄석술이 이 그림을 쓴다.

풍상 α만 밀면 접촉 압력이 먼저 흔들린다#

다섯 번째 플럭스는 0이고, $\alpha$ 는 면 속도의 부호로 풍상한다. 물에 $\alpha_1=0.01$ 만 새어도 압력 분모가 7.5% 움직이는 계산은, 계면 한 칸의 오차가 평형 압력을 건드린다는 뜻이다. 이 1차 갱신은 그 평형을 따로 강제하지 않는다. 비보존 항의 이산화가 압력을 그대로 두는지가 먼저 갈린다.

평형을 재는 문제는 아무 일도 없어야 하는 접촉이다. 거기서 압력이 어디로 새는지는 Abgrall 기준과 비보존 항이 묻는다.

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

5-Equation Model 직접 구현하기#

지배 방정식 정리#

상태방정식: Stiffened Gas EOS#

Python 구현#

초기 설정과 EOS#

HLLC Flux 계산#

시간 적분 (1st order)#

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

주의사항#

1. $\alpha$ 의 범위 제한#

2. 혼합 음속: Wood의 공식#

3. 고차 정확도로의 확장#

결과 해석#

풍상 α만 밀면 접촉 압력이 먼저 흔들린다#

관련

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

5-Equation Model 직접 구현하기#

지배 방정식 정리#

상태방정식: Stiffened Gas EOS#

Python 구현#

초기 설정과 EOS#

HLLC Flux 계산#

시간 적분 (1st order)#

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

주의사항#

1. α\alphaα의 범위 제한#

2. 혼합 음속: Wood의 공식#

3. 고차 정확도로의 확장#

결과 해석#

풍상 α만 밀면 접촉 압력이 먼저 흔들린다#

관련

1. $\alpha$ 의 범위 제한#