두 상에 각자의 방정식을 준다 — Baer–Nunziato 7-방정식 모델

비평형 7-방정식 모델의 구조, 비보존 항, 완화로 가는 모델 계층

1986년, Mark Baer와 Nunziato는 고체 추진제 속에서 불꽃이 어떻게 폭굉으로 돌변하는지를 풀려 했다. 화약은 미세한 알갱이의 더미였다. 그 틈으로 뜨거운 가스가 스며들며 알갱이를 태운다. 같은 점에 고체와 가스가 서로 다른 압력, 서로 다른 속도로 공존했다. 단일 압력·단일 속도를 가정하는 기존 혼합 모델로는 이 순간을 담을 수 없었다. 그래서 두 사람은 두 상에 각자의 보존 법칙을 통째로 주는 모델을 세웠다. 이 글은 그 7-방정식 모델의 구조를 따라가며, 왜 방정식이 일곱 개나 필요한지, 비보존 항이 왜 골칫거리인지, 그리고 완화(relaxation)가 어떻게 더 단순한 모델들을 낳는지 살펴본다. 마지막엔 압력 완화 연산자를 직접 적분해 두 상이 평형으로 가는 과정을 재현한다.

화약 알갱이에서 시작된 일곱 방정식#

확산 계면(diffuse interface, 두 상의 경계를 유한 두께로 번지게 보는 관점) 방법에는 갈래가 있다. 가장 단순한 4-방정식은 두 상이 압력·속도·온도까지 모두 같다고 본다. Baer–Nunziato(BN)는 정반대 끝에 선다. 아무 평형도 가정하지 않는다. 각 상은 자기만의 밀도, 속도, 압력을 가진다. 거기에 한 상이 차지하는 부피 비율, 즉 부피분율(volume fraction) $\alpha$ 의 수송식이 한 줄 더 붙는다.

왜 이 일반성이 중요할까. 비평형을 허용하면 모델이 **무조건 쌍곡형(hyperbolic)**이 된다. 단일 압력을 강제한 6-방정식 모델은 특정 조건에서 쌍곡성이 깨져 음의 음속 제곱 같은 비물리적 해를 낳는다. BN은 그 함정을 피한다. 대신 대가를 치른다. 방정식이 일곱 개로 늘고, 표준 보존형으로 쓸 수 없는 항이 끼어든다.

일곱 개의 방정식을 펼치면#

두 상을 $k = 1, 2$ 로 쓰자. 먼저 부피분율 수송식이다.

\frac{\partial \alpha_1}{\partial t} + u_I \frac{\partial \alpha_1}{\partial x} = \mu\,(p_1 - p_2)

좌변은 부피분율이 계면 속도 $u_I$ 로 실려 가는 것을, 우변은 압력 완화를 뜻한다. $\mu$ 는 완화율(클수록 빨리 평형). 나머지 여섯은 각 상의 질량·운동량·에너지 보존이다.

\frac{\partial (\alpha_k \rho_k)}{\partial t} + \frac{\partial (\alpha_k \rho_k u_k)}{\partial x} = 0

\frac{\partial (\alpha_k \rho_k u_k)}{\partial t} + \frac{\partial (\alpha_k \rho_k u_k^2 + \alpha_k p_k)}{\partial x} = p_I \frac{\partial \alpha_k}{\partial x} + \lambda\,(u_{k^*} - u_k)

\frac{\partial (\alpha_k E_k)}{\partial t} + \frac{\partial \big(\alpha_k (E_k + p_k)\,u_k\big)}{\partial x} = p_I u_I \frac{\partial \alpha_k}{\partial x} + \lambda u_I (u_{k^*} - u_k)

여기서 $\rho_k$ 는 상밀도, $u_k$ 는 상속도, $p_k$ 는 상압력, $E_k = \rho_k e_k + \tfrac12 \rho_k u_k^2$ 는 상의 총에너지 밀도, $k^*$ 는 상대 상이다. $\lambda$ 는 속도 완화율이다. 각 상은 자기 상태방정식으로 닫힌다. 추진제든 캐비테이션이든, 흔히 강성 기체(stiffened gas)를 쓴다.

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

$\gamma_k$ 는 비열비, $p_{\infty,k}$ 는 분자 간 반발을 나타내는 압력 상수다. 물처럼 거의 비압축인 액체는 $p_\infty$ 가 수 GPa로 크다.

비보존 항과 계면 변수 — 진짜 어려운 부분#

위 운동량·에너지 식의 우변에 $\partial \alpha_k / \partial x$ 가 있다. 이 항은 발산형으로 묶이지 않는다. 보존형 $\partial_t U + \partial_x F(U) = 0$ 으로 쓸 수 없다는 뜻이다. 표준 Godunov류 유한체적법은 이런 비보존(non-conservative) 항 앞에서 흔들린다. 불연속을 가로지르며 $\alpha$ 가 점프할 때 어떤 값으로 곱해야 하는지가 애매하기 때문이다. 잘못 이산화하면 균일한 압력·속도 장에서도 계면을 가로질러 가짜 진동이 튄다.

이 함정을 거르는 시금석이 Abgrall의 판정 기준이다. 압력과 속도가 처음에 균일하면, 부피분율이 어떻게 분포하든 그 균일함은 유지되어야 한다. 두 상이 같은 압력·속도로 흐르면 계면은 그저 수송될 뿐, 새 파동을 만들면 안 된다. 비보존 항의 이산화는 이 조건을 정확히 만족하도록 설계된다.

계면 변수 $u_I$ , $p_I$ 도 모델링의 산물이다. 가장 흔한 Saurel–Abgrall 선택은 질량가중 속도와 부피가중 압력이다.

u_I = \frac{\alpha_1 \rho_1 u_1 + \alpha_2 \rho_2 u_2}{\alpha_1 \rho_1 + \alpha_2 \rho_2}, \qquad p_I = \alpha_1 p_1 + \alpha_2 p_2

이 선택은 모델을 대칭으로 만들어, 모든 셀에서 같은 방정식·같은 수치 기법을 쓸 수 있게 한다.

완화: 평형으로 가는 길#

우변의 두 소스, $\mu(p_1-p_2)$ 와 $\lambda(u_{k^*}-u_k)$ 가 완화 연산자다. $\mu$ , $\lambda$ 를 크게 키우면 두 상은 거의 즉시 같은 압력·속도로 끌려간다. 이것이 BN의 묘미다. 완화율을 무한대로 보내는 극한이 곧 더 단순한 평형 모델이 된다.

아래 시뮬레이션에서 직접 조작해보자. 두 상이 서로 다른 압력·속도로 출발하고, 완화율 $\mu$ · $\lambda$ 가 둘을 평형으로 끌어당긴다.

압력 완화율 μ: 18속도 완화율 λ: 14초기 압력차 Δp: 1.6초기 속도차 Δu: 0.8

종료 시 압력차 0.000 bar · 속도차 0.000 m/s

μ·λ를 0으로 내리면 두 상은 영영 평형에 도달하지 못한다 — 이것이 7-방정식의 “비평형” 상태다.

완화율을 0으로 내리면 두 곡선은 영영 만나지 못한다. 이것이 7-방정식이 허용하는 비평형 상태다. $\mu$ 를 키울수록 압력 곡선이 더 가파르게 한 점으로 모인다. 압력과 속도는 서로 독립적으로 완화된다는 점도 눈여겨보자.

모델 계층: 7 → 6 → 5 → 4#

완화를 순간(instantaneous)으로 가정하면, 즉 평형을 강제하면 방정식 하나가 대수적 제약으로 바뀌며 사라진다. 속도를 평형으로 묶으면 6-방정식, 거기에 압력까지 묶으면 Kapila의 5-방정식, 온도까지 묶으면 4-방정식이다. BN은 이 계층의 꼭대기다.

아래에서 각 모델을 눌러, 무엇을 평형으로 가정하고 무엇을 잃는지 비교해보자.

→

Baer–Nunziato (full)

없음 — 압력·속도·온도 모두 독립

미지수: α, (ρ, u, p)₁, (ρ, u, p)₂

장점
가장 일반적. 무조건 쌍곡형(hyperbolic). 상마다 독립 EOS.

대가
비보존 항 + 두 속도장. 계면 변수 PI·uI 모델링이 까다롭다.

평형을 하나씩 더 부과할수록 방정식은 줄지만 — 일반성도 함께 줄어든다.

5-방정식(Kapila)이 계면 포착에 가장 널리 쓰이지만, 혼합 음속이 부피분율에 대해 비단조(Wood 곡선)라 부피분율 소스가 뻣뻣하다. BN은 그 뻣뻣함을 완화 소스로 떼어내 다룬다는 점에서 수치적으로 더 유연하다.

Python — 압력 완화 연산자를 적분한다#

완화의 핵심만 떼어내 보자. 두 상의 부분질량과 내부에너지는 보존하고, 부피분율 $\alpha$ 만 $d\alpha/dt = \mu(p_1 - p_2)$ 로 움직인다. 강성 기체 압력이 $\alpha$ 에 따라 어떻게 변하며 평형에 도달하는지 적분한다.

import numpy as np
 
# 강성 기체(stiffened gas) 두 상의 '압력 완화' 연산자.
# 각 상의 부분질량 m_k 와 부분 내부에너지 E_k 는 보존,
# 부피분율 alpha 만 완화 ODE  da/dt = mu (p1 - p2) 로 변한다.
 
def phase_pressures(alpha, mass, energy, gamma, p_inf):
    frac = np.array([alpha, 1.0 - alpha])
    rho  = mass / frac          # 상밀도 = 부분질량 / 부피분율
    e    = energy / mass        # 비내부에너지
    return (gamma - 1.0) * rho * e - gamma * p_inf
 
def relax_pressure(alpha0, mass, energy, gamma, p_inf,
                   mu=4.0e-8, dt=1.0e-3, nmax=200000, tol=1.0):
    alpha, hist = alpha0, []
    for n in range(nmax):
        p = phase_pressures(alpha, mass, energy, gamma, p_inf)
        hist.append((n * dt, alpha, p[0], p[1]))
        gap = p[0] - p[1]
        if abs(gap) < tol:                       # 압력차가 사라지면 종료
            break
        alpha += dt * mu * gap                    # da/dt = mu (p1 - p2)
        alpha = min(max(alpha, 1e-4), 1.0 - 1e-4)
    return alpha, np.array(hist)
 
# 물 같은 상(1, 강성 기체) + 공기 같은 상(2)
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])      # Pa
mass   = np.array([480.0, 0.55])     # 부분질량  (kg/m^3)
energy = np.array([3.93e8, 1.10e6])  # 부분 내부에너지 (J/m^3)
 
alpha_eq, hist = relax_pressure(0.50, mass, energy, gamma, p_inf)
p0  = phase_pressures(0.50, mass, energy, gamma, p_inf)
peq = phase_pressures(alpha_eq, mass, energy, gamma, p_inf)
print(f"초기: alpha=0.500  p1={p0[0]/1e5:8.3f} bar  p2={p0[1]/1e5:7.3f} bar")
print(f"평형: alpha={alpha_eq:.3f}  p1={peq[0]/1e5:7.3f} bar  p2={peq[1]/1e5:7.3f} bar")
print(f"반복 {len(hist)}회,  최종 압력차 {abs(peq[0]-peq[1]):.3e} Pa")

실행 결과는 이렇다.

초기: alpha=0.500  p1= 324.000 bar  p2=  8.800 bar
평형: alpha=0.506  p1=  8.906 bar  p2=  8.906 bar
반복 75회,  최종 압력차 9.012e-01 Pa

눈여겨볼 대목은 부피분율이 0.500에서 0.506으로 0.6%만 움직였는데 액체 압력이 324 bar에서 8.9 bar로 무너졌다는 점이다. 강성 기체로 모델한 물은 거의 비압축이라, 미세한 부피 변화가 거대한 압력 변화를 부른다. 완화 평형이 가스 압력 근처에 자리 잡는 이유도 같다. 뻣뻣한 액체는 거의 팽창하지 않아도 압력을 맞출 수 있기 때문이다.

기억할 점#

7-방정식 = 비평형. 두 상에 각자의 압력·속도·에너지를 준 대가로 무조건 쌍곡형을 얻는다.
비보존 항이 핵심 난점. $p_I \partial_x \alpha$ 는 보존형으로 안 묶인다. Abgrall 판정 기준(균일 압력·속도 보존)이 이산화의 시금석이다.
완화가 계층을 만든다. $\mu, \lambda \to \infty$ 극한이 6 → 5 → 4-방정식 모델이다. BN은 그 꼭대기에서 가장 일반적인 그림을 준다.