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

1986년, Baer와 Nunziato는 고체 추진제 속 불꽃이 어떻게 폭굉으로 넘어가는지를 풀려 했다. 화약은 미세한 알갱이 더미였다. 틈으로 뜨거운 가스가 스며 알갱이를 태운다. 같은 점에 고체와 가스가 서로 다른 압력, 서로 다른 속도로 있었다. 압력 하나, 속도 하나를 가정하는 혼합 모델로는 그 순간이 안 들어갔다. 그래서 두 상에 보존 법칙을 따로 줬다.

완화율 μ와 λ를 0에 두면 두 압력 곡선이 만나지 않고, 값을 올리면 한 점으로 끌려온다.

압력 곡선과 속도 곡선이 따로 붙는지 본다. 한쪽 슬라이더만 올려도 다른 곡선은 제자리인지 확인하면 된다.

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

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

이 일반성이 쌍곡성을 지킨다. 비평형을 허용하면 모델은 무조건 쌍곡형이다. 압력을 하나로 묶은 6-방정식은 특정 조건에서 쌍곡성이 깨지고, 음속 제곱이 음수가 되는 해도 낸다. BN은 그 함정을 피한다. 대가가 있다. 식이 일곱이 되고, 표준 보존형으로 못 쓰는 항이 끼어든다.

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

두 상을 $k=1,2$로 쓴다. 먼저 부피분율이다.

좌변은 부피분율이 계면 속도 $u_I$로 실려 가는 것이다. 우변은 압력 완화다. $\mu$가 클수록 평형이 빠르다. 나머지 여섯은 각 상의 질량, 운동량, 에너지다.

$\rho_k$는 상 밀도, $u_k$는 상 속도, $p_k$는 상 압력이다. $E_k=\rho_k e_k+\tfrac12\rho_k u_k^2$는 상의 총에너지 밀도다. $k^*$는 상대 상이다. $\lambda$는 속도 완화율이다. 각 상은 자기 상태방정식으로 닫힌다. 흔히 강성 기체다.

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

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

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

시금석은 Abgrall 기준이다. 압력과 속도가 처음에 균일하면, 부피분율이 어떻게 깔려 있든 그 균일함은 유지되어야 한다. 같은 압력, 같은 속도로 흐르면 계면은 수송될 뿐, 새 파동을 만들면 안 된다. 균일한데 압력이 21% 튄 기록은 비보존 항을 잘못 곱하면 그 기준이 바로 깨진다는 예다.

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

이 선택은 모델을 대칭으로 만든다. 모든 셀에서 같은 식, 같은 수치 기법을 쓸 수 있다.

완화: 평형으로 가는 길#

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

앞 화면에서 완화율을 0으로 내리면 두 곡선은 만나지 않는다. 7-방정식이 허용하는 비평형이다. $\mu$를 키울수록 압력 곡선이 가파르게 한 점으로 모인다. 압력과 속도는 따로 완화된다.

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

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

5-방정식이 계면 포착에 가장 많이 쓰인다. 혼합 음속은 부피분율에 대해 비단조인 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

모델을 누르면 무엇을 평형으로 묶었는지, 음속 곡선이 어디서 비단조로 꺾이는지가 바뀐다. 5식으로 내리는 순간 부피분율 소스가 뻣뻣해지는 자리가 그 화면이다.

출력의 324 bar는 두 큰 항의 차다. $\alpha=0.5$에서 액체 밀도는 480을 0.5로 나눈 960이다. 비내부에너지는 $3.93\times10^{8}/480=8.1875\times10^{5}$다. $(\gamma-1)\rho e$는 $2.6724\times10^{9}$ Pa다. $\gamma p_\infty$는 $4.4\times6\times10^{8}=2.64\times10^{9}$ Pa다. 빼면 $3.24\times10^{7}$ Pa, 324 bar다. 각 항의 1.2%만 남는다.

기울기도 그 식에서 나온다. $\alpha=0.5$에서 액체 압력의 $\alpha$ 미분은 약 $-5.34\times10^{9}$ Pa다. 부피분율 0.001은 액체 압력 약 53 bar다. 인쇄된 이동은 0.500에서 0.506, 차이 0.006이다. 기체 압력은 $4.4\times10^{5}/(1-\alpha)$라 0.5에서 8.8 bar, 0.506 근처에서 8.9 bar다. 평형은 기체가 서 있는 자리로 액체가 내려가는 쪽이다.

첫 스텝의 크기는 코드에 적힌 계수다. $\mu=4\times10^{-8}$, 시간 간격은 $10^{-3}$, 초기 차이는 324 bar에서 8.8 bar를 뺀 $3.152\times10^{7}$ Pa다. 부피분율 증분은 그 셋을 곱한 약 $1.26\times10^{-3}$이다. 인쇄된 총 이동 0.006의 약 5분의 1이다. 75스텝은 그다음부터 차이가 줄어 잘게 잘린 결과다. 허용오차 1 Pa는 마지막에나 의미가 있다.

부피분율은 0.500에서 0.506으로 0.6%만 움직였는데, 액체 압력은 324 bar에서 8.9 bar로 무너졌다. 강성 기체로 둔 물은 거의 안 눌려서, 미세한 부피 변화가 거대한 압력 변화를 부른다. 완화 평형이 가스 압력 근처에 앉는 이유도 같다. 뻣뻣한 액체는 거의 안 팽창해도 압력을 맞출 수 있다.

첫 스텝이 부피분율 0.001을 먹는 완화#

첫 스텝이 부피분율을 이미 $1.26\times10^{-3}$ 움직인다. 액체 압력은 그 이동만으로 수십 bar가 빠진다. 비보존 항의 $p_I$는 이 붕괴를 계면 기울기에 곱한다. 완화와 플럭스를 같은 단계에 두면, 가짜 파의 크기가 이 첫 스텝에서 정해진다.

이류와 압력을 플럭스에서 완전히 갈라 두면 그 곱이 어디로 가는가. Baer–Nunziato 플럭스를 이류와 압력으로 가른 TV 분리는 그 분리가 비보존 항까지 덮는지 아직 열어 둔다. 324 bar가 한 스텝에 무너질 때, 갈라 둔 압력 플럭스는 계면을 조용히 두는가.