[논문 리뷰] 두 압력으로 풀고 한 압력으로 합치다 — Saurel(2009) 압축성 다상유동 이완법

물·공기 한 칸에서 Wood 음속은 약 24 m/s까지 꺼지고, 임피던스 평균은 공기 압력 쪽에 붙는다.

미국 아이다호 국립연구소의 핵연료 안전 시뮬레이션이 멈췄다. 액체 나트륨 사이에 갇힌 가스 거품 한 칸의 음속이 갑자기 25 m/s까지 곤두박질친 탓이다. 평형 압력을 가정한 5-방정식 다상유동 모델의 한계였다. Saurel과 동료는 2009년에 한 발 물러선다. 두 압력을 따로 풀고, 스텝 끝에서만 붙인다.

$\rho_k$ 와 $c_k$ 슬라이더를 움직이면 Wood 곡선과 frozen 곡선이 벌어지는 폭이 바뀐다.

ρ₁ = 1000 kg/m³c₁ = 1500 m/sρ₂ = 10.0 kg/m³c₂ = 340 m/scursor α₁ = 0.500log scale on y-axis

Wood (equilibrium) sound speed dips far below either pure-fluid value near intermediate α; frozen speed stays inside the [c₂, c₁] band.

물 1500 m/s, 공기 340 m/s 근처에서 Wood 쪽만 바닥으로 꺼진다. frozen 쪽은 두 순수 음속 사이에 머문다.

한 줄 요지#

평형 5-방정식이 음속에서 무너진다. 6-방정식으로 한 단계 물러나 풀고, 강체 이완으로 다시 합쳐라. 인터페이스 압력은 임피던스 평균 $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1+Z_2)$ 한 줄이 모든 것을 받쳐준다.

Wood 음속이 모래주머니다#

Kapila(2001)의 5-방정식 모델은 압력 평형을 가정한다. 두 상이 한 셀 안에서 같은 압력을 갖는다는 가정이다. 깔끔해 보이지만 음속이 이상해진다. 혼합 음속 $c_{eq}$ 는 Wood(1930) 공식

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

을 따르는데, 여기서 $\alpha_k$ 는 상 $k$ 의 부피분율, $\rho_k$ 는 밀도, $c_k$ 는 순수 음속이다. 문제는 이 곡선이 단조하지 않다는 점이다. 물(1500 m/s)과 공기(340 m/s)가 섞인 셀의 평형 음속은 두 끝값보다 한참 낮은 약 24 m/s까지 내려간다.

수치 확산 영역에서는 두 끝점 사이의 가상 혼합이 잠깐 만들어지고, 음파는 그 안에서 갑자기 느려진다. 음파의 도달 시간이 어긋나 충격파 추적이 망가진다. 게다가 $\alpha$ 가 0이나 1에 닿을 때 양의 부피분율을 보장하기도 까다롭다.

코드에 적힌 칸을 Wood 식에 넣으면 본문의 24 m/s가 다시 나온다. $\alpha_1 = 0.7$ , $\rho_1 = 1000$ , $c_1 = 1500$ 이면 $\rho_1 c_1^2 = 2.25\times 10^9$ 이다. 공기 쪽 $\alpha_2 = 0.3$ , $\rho_2 = 1$ , $c_2 = 340$ 이라 $\rho_2 c_2^2 = 1.156\times 10^5$ 다. 오른쪽 항은 공기가 $2.6\times 10^{-6}$ 근처, 물은 $3\times 10^{-10}$ 이라 거의 안 보인다. 혼합 밀도는 $0.7\times 1000 + 0.3 = 700.3$ 이다. $1/(\rho c_{eq}^2)$ 이 공기 항에 먹히면 $c_{eq}$ 는 약 24 m/s다. 같은 칸의 frozen 식 $c_f^2 = Y_1 c_1^2 + Y_2 c_2^2$ 는 반대다. 질량의 거의 전부가 물이라 $c_f$ 는 1500 m/s 옆에 붙는다. 한 칸에서 두 음속이 예순 배쯤 벌어진다. 이 꺼짐은 물과 공기를 반씩 섞은 음속이 더 파고든다.

우회로 — 압력을 일부러 다르게 둔다#

Saurel은 한 단계 뒤로 물러난다. 압력 평형 가정을 깨고 두 상이 서로 다른 압력 $p_1$ , $p_2$ 를 가질 수 있게 한다. 7-방정식 Baer–Nunziato 모델에서 속도 이완 극한만 먼저 취한 결과가 6-방정식 모델이다.

\partial_t \alpha_1 + u\,\partial_x \alpha_1 = \mu(p_1 - p_2)

\partial_t (\alpha_k \rho_k) + \partial_x(\alpha_k \rho_k u) = 0

\partial_t (\rho u) + \partial_x(\rho u^2 + \alpha_1 p_1 + \alpha_2 p_2) = 0

\partial_t (\alpha_k \rho_k e_k) + \partial_x(\alpha_k \rho_k e_k u) + \alpha_k p_k\,\partial_x u = \pm p_I \mu (p_1 - p_2)

여기서 $\mu$ 는 압력 이완율(stiffness), $p_I$ 는 인터페이스 압력, $e_k$ 는 상별 내부에너지다. 이 모델의 음속은

c_f^2 = Y_1 c_1^2 + Y_2 c_2^2

로 단조하다. $Y_k = \alpha_k \rho_k / \rho$ 는 질량분율. Wood가 만든 모래주머니가 사라진다.

인터페이스 압력 $\hat{p}_I$ — 임피던스 가중평균#

내부에너지 식의 우변에 들어가는 $p_I$ 는 7-방정식 모델의 비대칭 항을 속도 평형 극한에서 정리한 결과로,

p_I = \frac{Z_2 p_1 + Z_1 p_2}{Z_1 + Z_2}, \quad Z_k = \rho_k c_k

라는 깔끔한 형태로 떨어진다. $Z_k$ 는 음향 임피던스. 더 단단한 상이 더 큰 가중치를 갖는다. 직관적으로 보면 두 상이 한 인터페이스를 통해 일을 주고받을 때, 단단한 쪽이 더 적게 움직이고 부드러운 쪽이 더 많이 움직이는 비율이다. 이 한 줄이 에너지 보존과 엔트로피 부등식을 동시에 보장한다.

이완 단계 — 단일 변수 방정식 한 줄#

수치적으로 $\mu \to \infty$ 극한을 직접 풀지는 않는다. 대신 시간 스텝 안에서 hyperbolic 부분을 먼저 푼 뒤, 이완 항만 남은 ODE를 적분한다. 이완 단계에서는 $u$ , $\rho E$ , $\alpha_k \rho_k$ 가 보존되고 두 상의 비체적 $v_k = 1/\rho_k$ 만 변한다. Stiffened gas EOS

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

를 쓰면 적분식이 닫힌 형태로 떨어진다.

v_k(p) = v_k^0\,\frac{p^0 + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}{p + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}

포화 조건 $\sum_k (\alpha\rho)_k v_k(p) = 1$ 만 풀면 평형 압력 $p$ 가 한 번에 나온다. 이분법 한 번이면 끝. 이 게 Saurel 알고리즘의 핵심이다.

재초기화 — 보존 압력으로 다시 정렬#

이완으로 얻은 $p$ , $\alpha_k$ 가 각 상의 EOS에는 일치하지만, 혼합 보존 에너지와는 어긋날 수 있다. 충격파 영역에서 비보존 내부에너지 식의 누적 오차 때문이다. 그래서 Saurel은 한 번 더 손질한다. 보존된 혼합 총에너지로부터 혼합 EOS

p(\rho, e, \alpha_1, \alpha_2) = \frac{\rho e - \big(\frac{\alpha_1 \gamma_1 p_{\infty,1}}{\gamma_1 - 1} + \frac{\alpha_2 \gamma_2 p_{\infty,2}}{\gamma_2 - 1}\big)}{\frac{\alpha_1}{\gamma_1 - 1} + \frac{\alpha_2}{\gamma_2 - 1}}

로 압력을 다시 계산하고, 그 압력으로 각 상의 $e_k$ 를 재초기화한다. 비보존 모델로 시작했지만 충격파 위치에서는 진짜 보존 모델처럼 행동한다. 비싼 추가 비용은 없다 — 셀당 한 번의 EOS 평가뿐.

코드: 한 셀에서 두 압력이 만나는 순간#

이완 단계만 NumPy 30줄로 따라가 보자. Stiffened gas + 이분법으로 단일 방정식을 푸는 코어 루프다.

import numpy as np
 
def stiffened_gas_volume(p, p0, v0, gamma, p_inf, p_hat_I):
    num = p0 + gamma * p_inf + (gamma - 1.0) * p_hat_I
    den = p   + gamma * p_inf + (gamma - 1.0) * p_hat_I
    return v0 * num / den
 
def saturation_residual(p, alpha_rho, p0, v0, gamma, p_inf, p_hat_I):
    total = 0.0
    for k in range(2):
        v_k = stiffened_gas_volume(p, p0[k], v0[k], gamma[k], p_inf[k], p_hat_I)
        total += alpha_rho[k] * v_k
    return total - 1.0
 
def stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z, tol=1e-9):
    p_hat_I = (Z[1] * p0[0] + Z[0] * p0[1]) / (Z[0] + Z[1])
    lo, hi = 1.0, 1.0e10
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        r = saturation_residual(mid, alpha_rho, p0, v0, gamma, p_inf, p_hat_I)
        if r > 0:
            lo = mid       # 너무 낮은 압력 → 부피 합 > 1
        else:
            hi = mid
        if hi - lo < tol * mid:
            break
    return 0.5 * (lo + hi), p_hat_I
 
# 물(α=0.7) + 공기(α=0.3), 두 압력이 다름
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])
rho0   = np.array([1000.0, 1.0])
alpha0 = np.array([0.7, 0.3])
p0     = np.array([5.0e5, 1.0e5])
Z      = np.array([1500.0 * 1000.0, 340.0 * 1.0])
 
alpha_rho = alpha0 * rho0
v0        = 1.0 / rho0
 
p_eq, p_hat_I = stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z)
print(f"임피던스 평균 ^pI = {p_hat_I:.3e} Pa")
print(f"이완 후 평형압  p  = {p_eq:.3e} Pa")

물의 임피던스( $Z_1 = 1.5\times 10^6$ )가 공기( $Z_2 = 340$ )보다 약 4400배 크다. 식은 $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1 + Z_2)$ 다. 큰 가중치는 작은 쪽 압력에 붙는다. $p_1 = 5\times 10^5$ , $p_2 = 1\times 10^5$ 를 넣으면 분자의 $Z_1 p_2$ 가 $1.5\times 10^{11}$ 이고 $Z_2 p_1$ 은 $1.7\times 10^8$ 에 그친다. 그래서 $\hat{p}_I$ 는 물 압력 옆이 아니라 공기 압력 $1\times 10^5$ 바로 옆, 약 $1.001\times 10^5$ Pa에 선다. 단단한 상이 계면을 독점하는 그림이 아니다. 이완 뒤의 공통 압력은 그 다음, 포화 조건의 이분법이 따로 찾는다.

관찰: 두 압력이 만나는 시간#

코드 한 셀로는 동적 감각이 안 잡힌다. 1D 격자에서 단순화한 ODE 모델을 풀어보자. 셀마다 $p_1, p_2$ 가 임피던스 평균 $\hat{p}_I$ 를 향해 지수적으로 수렴한다 — $\mu \cdot \Delta t$ 가 클수록 빨리.

Z₁ = 1.50 MPa·s/mZ₂ = 3.40 kPa·s/mμ·Δt = 0.050p₁(left) = 5.00 barp₂(left) = 1.00 barp₁(right) = 1.00 barp₂(right) = 5.00 bar

Stiff relaxation drives both p₁ and p₂ toward the impedance-weighted average p̂ᴵ = (Z₂p₁ + Z₁p₂)/(Z₁+Z₂). Increase μ·Δt for faster pull-in.

$Z_1$ 을 키우면 $\hat{p}_I$ 가 더 빠르게 phase 1 압력 쪽으로 끌려간다. 임피던스 비율이 그대로 가중치다. $\mu \cdot \Delta t$ 를 1에 가깝게 두면 한 스텝 안에서 사실상 평형으로 떨어진다 — 이게 stiff 이완의 의미다. Saurel 논문은 이 절차를 매 시간 스텝 끝에 무조건 한 번씩 적용해 5-방정식 모델의 강체 극한을 회복한다.

한계 — 7-방정식의 그림자#

6-방정식은 만능이 아니다. 속도 평형은 여전히 가정이다. 다공성 매체나 분말 폭약처럼 두 상의 속도가 크게 다른 경우는 7-방정식 Baer–Nunziato로 돌아가야 한다. 충격파 강도가 매우 클 때 비보존 내부에너지 식의 누적 오차도 무시할 수 없어, Petitpas의 인공 열교환 보정을 추가한다(논문 5장). 그리고 압력 이완은 $\mu \to \infty$ 를 가정하지만 실제 인터페이스 응답 시간은 유한하다 — 표면장력이나 질량전달이 들어오면 별도 모델이 필요하다.

그럼에도 6-방정식 + 이완은 압축성 다상 시뮬레이션에서 가장 널리 쓰이는 토대다. AP-IMEX 시간 적분이나 MUSCL-THINC-BVD 재구성이 모두 이 6-방정식 위에서 돈다.

5-방정식을 그 칸에 둘 것인가#

Wood 음속이 24 m/s까지 꺼진 칸을 5-방정식에 맡기면 음파 도달이 밀리고 충격 위치가 흔들린다. 그 칸은 6-방정식으로 물러나 두 압력을 따로 두고, 스텝 끝에서만 임피던스 평균으로 붙인다.

속도까지 어긋나면 이 타협은 끝이다. 다공 매질이나 분말 폭약은 7-방정식으로 돌아간다. 식을 더 줄여 계면을 닫을지는 4-방정식 평형이 따로 묻는다.