[논문 리뷰] 등속으로 흐르는 물-공기 계면에서 압력이 8억 Pa 튀었다 — 압력기반 BN 모델의 비교란 조건

균일한 흐름에 계면 하나, 압력은 8억 Pa가 튀었다#

속도 1 m/s, 압력 1 bar로 균일한 1D 관이 있다. 그 안에 물 덩어리 하나가 공기 사이에 떠서 흘러간다. 물리적으로는 아무 일도 없어야 한다. 그런데 질량 플럭스만 2차로 올리고 체적분율은 1차로 남겼더니, 120스텝 뒤 계면 셀의 압력이 $8\times10^{8}$ Pa 어긋났다. 이 글은 Re와 Abgrall의 2022년 논문(IJNMF, A pressure-based method for weakly compressible two-phase flows under a Baer–Nunziato type model)을 Q&A 다섯 문항으로 읽는다. 읽고 나면 다상 코드의 비보존항을 어떻게 차분해야 이 숫자가 0으로 돌아오는지 알게 된다.

Q1. 이 논문은 무엇을 풀려고 했나?#

대상은 약압축성(weakly compressible, 마하수가 작아 밀도 변화가 작은) 2상 유동이다. 탱크 안 슬로싱이나 벽 근처의 갇힌 기포가 예다. 액체는 거의 비압축이지만 기포가 섞인 혼합물은 강하게 압축된다. 그래서 두 상의 압축성을 모두 남겨야 한다.

모델은 Saurel–Abgrall의 Baer–Nunziato형(BN, 각 상이 자기 압력·속도·온도를 갖는 7방정식) 모델이다. 논문은 압력·속도 이완(relaxation, 두 상의 값을 평형으로 끌어당기는 원천항) 계수를 무한대가 아닌 유한값으로 둔다. 상태방정식(EOS)은 stiffened gas와 Peng–Robinson 등 일반형을 허용한다.

밀도기반 솔버는 저마하에서 예조건(preconditioning)이 필요하다. 그런데 2상 혼합물은 음속이 하나로 정해지지 않는다. 그래서 저자들은 압력기반 정식화를 택했다. Bijl–Wesseling 스케일링으로 압력을 두 변수로 나누고, 음향 항을 암시적으로 푼다. 그 결과 음향 CFL이 10을 넘어도 안정했다. 엇갈림 격자(staggered grid)로 체커보드 진동도 막았다. Toro–Vázquez 분할이 압력 부분계로 강성을 옮기는 방식과 같은 문제의식이다.

Q2. 비교란 조건이란 무엇인가?#

비교란 조건(non-disturbance condition)은 Abgrall이 1996년에 정리한 요구다. 압력과 속도가 균일한 상태는 계면이 지나가도 균일하게 남아야 한다. 연속 방정식에서는 자명하다. 이산 방정식에서는 자명하지 않다.

BN 모델의 상 $\sigma$ 운동량 식에는 비보존항 $P_I\,\partial\alpha_\sigma/\partial x$가 있다.

$\alpha_\sigma$는 체적분율, $P_I$는 계면 압력이다. $P$가 균일하면 좌변 압력항 $P\,\partial\alpha/\partial x$와 우변이 정확히 상쇄된다. 이산화에서도 상쇄되려면 두 항이 같은 차분이어야 한다. 논문의 식 (49)가 그 답이다.

$k$는 셀 $i$와 $i+1$ 사이의 엇갈림 셀이다. 압력 구배는 $-[(\alpha P)_{i+1}-(\alpha P)_i]$로 차분한다. 같은 두 셀의 같은 차이이므로 균일 $P$에서 합이 0이다.

아래 시뮬레이션에서 $H_P$의 차분을 직접 바꿔보자.

averaged(셀 중심 중앙차분 기울기를 면으로 평균)를 고르면 계면 양옆 면에 빨간 힘이 생기고 물의 속도가 1에서 벗어난다. 계면 폭을 1칸에서 6칸으로 넓히면 그 힘이 줄지만 사라지지 않는다. 식 (49)는 폭과 압력에 상관없이 초록 막대, 즉 0을 유지한다.

Q3. 체적분율 식은 왜 질량 플럭스를 따라야 하나?#

같은 논리가 체적분율 식 $\partial_t\alpha_\sigma + u_I\,\partial_x\alpha_\sigma = \dots$에도 적용된다. 상 밀도는 $\rho_\sigma = (\alpha\rho)_\sigma/\alpha_\sigma$로 복원된다. 분자와 분모가 서로 다른 확산을 받으면 균일했던 $\rho_\sigma$가 계면 셀에서 바뀐다.

음속 $c_\sigma$ 근처로 선형화한 EOS를 쓰면 혼합 압력 오차가 한 줄로 나온다.

$\rho_{0,\sigma}$는 기준 밀도다. 괄호 안은 "질량 식이 옮긴 양"과 "체적분율 식이 옮긴 양"의 차이다. 물은 $c^2 = 2.25\times10^6$ m²/s²이므로, 1 kg/m³의 불일치도 2.25 MPa가 된다.

그래서 논문은 식 (45)에서 $u_I\,\partial_x\alpha$를 질량 식의 Rusanov 플럭스(국소 최대 파속으로 확산을 주는 1차 플럭스)와 같은 꼴로 쓴다.

압력기반이라 파속 $S$에 음속이 들어가지 않는다. 질량 식의 Rusanov 플럭스도 $S=\max|u|$를 쓴다. 두 식이 같은 $S$를 공유하니 균일 $\rho_\sigma$가 보존된다. 이는 $\hat F$가 진짜 면 플럭스가 아니라 "질량 플럭스를 흉내 낸 연산자"라는 뜻이다.

Q4. 코드로 확인하면 어떤가?#

물(1000 kg/m³, 1500 m/s)과 공기(1.2 kg/m³, 340 m/s)의 계면을 $u=1$로 0.3만큼 옮긴다. 체적분율과 두 상의 질량을 네 가지 조합으로 푼다. 이어서 균일 $P$에서 운동량 면 잔차 $R_k$를 두 가지 $H_P$로 잰다.

import numpy as np
 
# 균일한 u, p 속을 흐르는 1D 물-공기 계면 (균일하게 남아야 한다)
N, L, u, CFL = 200, 1.0, 1.0, 0.5
dx = L / N
dt = CFL * dx / u
x = (np.arange(N) + 0.5) * dx
rho0 = np.array([1000.0, 1.2])     # 물, 공기 기준 밀도 [kg/m^3]
c0 = np.array([1500.0, 340.0])     # 상별 음속 [m/s]
p0, eps = 1.0e5, 1e-6
 
def minmod_slope(q):
    a = q - np.roll(q, 1)
    b = np.roll(q, -1) - q
    return np.where(a * b > 0, np.sign(a) * np.minimum(abs(a), abs(b)), 0.0)
 
def face_flux(q, order):
    # S=|u|인 Rusanov 플럭스 (압력기반: 음속 없음), u>0
    qL = q + 0.5 * minmod_slope(q) if order == 2 else q
    return u * qL
 
def central_flux(q):
    return 0.5 * u * (q + np.roll(q, -1))
 
def mixture_dp(a1, m1, m2):
    a = np.array([a1, 1 - a1])
    rho = np.array([m1, m2]) / a
    p = p0 + (c0**2)[:, None] * (rho - rho0[:, None])
    return np.max(abs((a * p).sum(0) - p0))
 
def march(alpha_rule, mass_order, t_end=0.3):
    a1 = np.where(abs(x - 0.3) < 0.15, 1 - eps, eps)
    m1, m2 = rho0[0] * a1, rho0[1] * (1 - a1)
    n = int(round(t_end / dt))
    for _ in range(n):
        Fa = central_flux(a1) if alpha_rule == "central" else face_flux(a1, alpha_rule)
        F1, F2 = face_flux(m1, mass_order), face_flux(m2, mass_order)
        a1 = np.clip(a1 - dt / dx * (Fa - np.roll(Fa, 1)), eps, 1 - eps)
        m1 = m1 - dt / dx * (F1 - np.roll(F1, 1))
        m2 = m2 - dt / dx * (F2 - np.roll(F2, 1))
    return n, mixture_dp(a1, m1, m2)
 
cases = [("A  alpha 1st / mass 1st", 1, 1),
         ("B  alpha 1st / mass 2nd", 1, 2),
         ("C  alpha 2nd / mass 2nd", 2, 2),
         ("D  alpha central / mass 1st", "central", 1)]
print("uniform u=1, p=1e5 Pa, N=200, CFL=0.5")
for name, ar, mo in cases:
    n, dp = march(ar, mo)
    print(f"{name:30s} steps={n:4d}  max|p_mix - p0| = {dp:10.3e} Pa")
 
def hp_residual(a, P, rule):
    # 엇갈림 면 k(셀 i와 i+1 사이)의 운동량 원천
    grad = -(np.roll(a, -1) * P - a * P)          # -[(aP)_{i+1} - (aP)_i]
    if rule == "jump":
        HP = P * (np.roll(a, -1) - a)             # 논문 식 (49)
    else:
        g = 0.5 * (np.roll(a, -1) - np.roll(a, 1))   # 셀 중심 중앙차분 기울기
        HP = P * 0.5 * (g + np.roll(g, -1))          # 면으로 평균
    return grad + HP
 
print()
print("momentum source on faces, uniform P=1e5 Pa, water-rich faces (alpha>0.5)")
for w in (1, 3, 6):
    aw = 0.5 * (1 + np.tanh((0.15 - abs(x - 0.3)) / (w * dx)))
    aw = np.clip(aw, eps, 1 - eps)
    af = 0.5 * (aw + np.roll(aw, -1))
    for rule in ("jump", "averaged"):
        R = hp_residual(aw, p0, rule)
        du = dt * R / (dx * af * rho0[0])
        wet = af > 0.5
        print(f"width={w}  H_P={rule:8s}  max|R| = {np.max(abs(R)):9.3e} Pa"
              f"  du_water/step = {np.max(abs(du[wet])):9.3e} m/s")

uniform u=1, p=1e5 Pa, N=200, CFL=0.5
A  alpha 1st / mass 1st        steps= 120  max|p_mix - p0| =  5.116e-07 Pa
B  alpha 1st / mass 2nd        steps= 120  max|p_mix - p0| =  8.024e+08 Pa
C  alpha 2nd / mass 2nd        steps= 120  max|p_mix - p0| =  1.769e-05 Pa
D  alpha central / mass 1st    steps= 120  max|p_mix - p0| =  1.665e+09 Pa
 
momentum source on faces, uniform P=1e5 Pa, water-rich faces (alpha>0.5)
width=1  H_P=jump      max|R| = 1.319e-11 Pa  du_water/step = 6.777e-15 m/s
width=1  H_P=averaged  max|R| = 1.203e+04 Pa  du_water/step = 1.874e+00 m/s
width=3  H_P=jump      max|R| = 8.185e-12 Pa  du_water/step = 4.670e-15 m/s
width=3  H_P=averaged  max|R| = 8.326e+02 Pa  du_water/step = 4.142e-01 m/s
width=6  H_P=jump      max|R| = 1.283e-11 Pa  du_water/step = 6.439e-15 m/s
width=6  H_P=averaged  max|R| = 1.126e+02 Pa  du_water/step = 8.670e-02 m/s

위 코드를 브라우저에서 직접 돌려보자.

기본값에서 120스텝이 끝나면 네 숫자가 출력 블록과 같아진다. c_water를 1500에서 300으로 내리면 B와 D의 오차가 약 25배 줄어든다. 오차가 $c^2$에 비례한다는 Q3의 식 그대로다. A와 C는 어떤 값에서도 반올림 수준에 머문다.

A와 C는 각각 $5\times10^{-7}$, $2\times10^{-5}$ Pa다. 차수가 아니라 두 식의 차분이 같은가가 결과를 갈랐다. B는 질량만 MUSCL로 올렸다. 계면 셀에서 $(\alpha\rho)_1$은 덜 번지고 $\alpha_1$은 더 번진다. 공기 쪽 셀에 물 질량이 체적보다 적게 남아 $\rho_1$이 기준값에서 크게 어긋난다. D의 중앙차분은 $\alpha$에 진동을 만들고, 그 진동이 그대로 압력이 된다.

운동량 쪽에서는 식 (49)가 $10^{-11}$ Pa로 기계 정밀도다. 평균 기울기는 1칸 계면에서 면 잔차 $1.2\times10^4$ Pa를 낸다. 물 속도가 한 스텝에 1.87 m/s 바뀐다. 흐름 속도 자체가 1 m/s이니 첫 스텝에 흐름이 뒤집힌다.

Q5. 그럼 2차 정확도는 포기해야 하나?#

아니다. C가 보여주듯 두 식에 같은 재구성을 쓰면 된다. minmod는 1차 동차 함수다. 그래서 $(\alpha\rho)_\sigma = \rho_0\alpha_\sigma$이면 기울기도 $\rho_0$배가 되어 비율이 유지된다. 비선형 리미터를 상별 질량에 따로 걸면 이 성질이 깨질 수 있다. 그때는 $\alpha$와 $\rho_\sigma$ 같은 원시변수를 재구성하고 곱으로 질량 플럭스를 만든다.

주의할 점이 두 가지 있다. 첫째, 이 실험은 $u$가 상수라 $u_I\,\partial_x\alpha$와 $\partial_x(\alpha u)$가 같다. 실제 코드에서는 $u_I$가 셀마다 다르다. 그래서 논문은 $(u_I)_i$를 면이 아닌 셀 값으로 쓴다. 균일 $\alpha$에서 $H_u=0$이 되게 하려는 선택이다. 둘째, 논문은 이 이산화가 "강건하지만 스킴 의존적"이라고 밝힌다. 질량 플럭스를 바꾸면 $H_u$와 $H_P$를 다시 유도해야 한다. 압축형 리미터처럼 체적분율 쪽 재구성만 바꾸는 시도는 정확히 이 지점에서 위험하다.

논문의 다른 결과도 짧게 정리한다. 물-알루미늄, 물-공기 충격관과 CO₂ 증기-액체 혼합물 문제에서 기준해와 몇 % 안에서 맞았다. 충격 부근의 차이는 비보존 모델에서 예상되는 수준이다. 이완 계수를 유한하게 두면 두 상의 압력·속도가 평형으로 가는 과정이 해에 그대로 보인다.

다음 리미터를 바꾸기 전에 돌려볼 등속 계면 테스트#

균일한 $u$, $p$ 속에 날카로운 계면 하나를 놓고 100스텝 진행한다. 혼합 압력의 최대 편차를 출력한다. 기계 정밀도가 아니면 $\alpha$ 식과 질량 식의 차분이 어긋난 것이다. 운동량은 첫 스텝의 면 잔차 $R_k$만 보면 된다. 플럭스나 리미터를 바꿀 때마다 이 테스트를 다시 돌린다. 몇 줄짜리 테스트지만 8억 Pa짜리 버그를 첫날 잡는다.