[논문 리뷰] Poisson 없이 두 상이 산다 — 일반 압력 방정식과 MSTACS의 완전 명시적 결합

Poisson을 hyperbolic GPE로 대체해 두 상 유동을 완전 명시적으로 푼다

LBM, ACM, EDAC, GPE — 압력 Poisson을 안 풀고 비압축 유동을 모사하는 네 갈래의 길이 있다. Bodhanwalla·Raghunathan·Sudhakar(2025)는 그중 마지막을 골라 두 상까지 끌어왔다. 게다가 MSTACS라는 sharper algebraic VOF를 operator-split과 묶어 선형 솔버 호출 0회의 완전 명시적 솔버를 만들었다. 1000:1 밀도비에서도 도는데, 비밀은 압력 방정식을 hyperbolic으로 두고 음속을 인공으로 박아 넣은 것이다.

논문 정보#

저자: H. Bodhanwalla, D. Raghunathan, Y. Sudhakar
소속: IIT Goa, School of Mechanical Sciences
저널: Computers & Fluids (2025)
제목: A general pressure equation based method for incompressible two-phase flows
키워드: general pressure equation, VOF, MSTACS, operator-split, weakly compressible

네 가지 "Poisson 버리기"의 비교 매트릭스#

방법	압력 식의 성격	두 상 적용	메모리	병렬 효율
LBM	kinetic 분포함수 합	대 밀도비에서 불안정	큼 (Q19/Q27)	우수
ACM	인공 hyperbolic ( $\partial_t p + \beta \nabla \cdot \mathbf{u} = 0$ )	dual-time 필요	보통	보통
EDAC	음향파 감쇠 parabolic	계면에서 switch parameter 필요	보통	우수
GPE	열역학 유도 hyperbolic + diffusion	본 논문이 최초 직접 적용	보통	우수

다른 셋이 비슷한 출발점에서 갈라졌다면, GPE는 **상태방정식 가정 $\gamma = \mathrm{Pr}$ **에서 유도된다는 점이 다르다. 인공 변수가 아니라 thermodynamic limit이 만들어 준 식이다.

GPE — 압력에 음속을 빌려주는 식#

운동량은 평범한 Navier–Stokes 그대로. 변한 것은 압력 식이다.

\rho \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}\right) = -\nabla p + \nabla \cdot \left[\mu(\nabla\mathbf{u} + \nabla\mathbf{u}^\top)\right] + \mathbf{F}

\frac{\partial p}{\partial t} + \rho c_s^2 (\nabla \cdot \mathbf{u}) = \frac{1}{\rho}\nabla \cdot (\mu \nabla p)

$c_s$ 는 인공 음속, $\rho$ 는 혼합 밀도( $\rho = C\rho_1 + (1-C)\rho_2$ , $C$ 는 부피분율), $\mu$ 는 혼합 점성. 우변의 압력 diffusion은 GPE 유도 과정에서 자연스럽게 떨어진 항이지 EDAC처럼 사람이 추가한 항이 아니다.

비차원 형태로 쓰면 압축성 효과는 인공 Mach 수 $\mathrm{Ma} = U/c_s$ 로 정량화된다.

\frac{\partial p}{\partial t} + \frac{\rho^*}{\mathrm{Ma}^2}(\nabla \cdot \mathbf{u}) = \frac{1}{\rho^*\mathrm{Re}}\nabla \cdot (\mu^* \nabla p)

$\mathrm{Ma}$ 가 작을수록 압력이 빨리 전파돼 비압축 극한에 가까워지지만, 안정성을 위해 시간 스텝은 $\Delta t \le \Delta x / c_s$ 로 줄어든다. 음속을 빌려 온 대가다.

논문 표 5에 따르면 $\Delta t_{\mathrm{GPE}} \approx \mathrm{Ma}\,\Delta t_{\mathrm{INS}}$ . 직렬로 돌리면 INS Poisson 솔버보다 느리지만, Poisson 호출이 없어 HPC 확장성에서 역전된다는 게 저자들의 주장이다.

MSTACS — donor–acceptor에 cos⁴θ를 끼우다#

VOF의 부피분율 전달은 $\partial_t C + \nabla \cdot (\mathbf{u} C) = C(\nabla \cdot \mathbf{u})$ . 압축성 보정 항 $C(\nabla \cdot \mathbf{u})$ 가 우변에 붙은 것에 주목하자. 약압축 framework가 만든 GPE 호환 형태다.

핵심은 면 값 $C_{\text{face}}$ 를 계산하는 방법. donor–acceptor 형식으로

C_{\text{face}} = (1 - \beta) C_D + \beta C_A

여기서 $\beta$ 는 정규화 변수 $\widetilde{C}_D = (C_D - C_U)/(C_A - C_U)$ 로 결정된다. MSTACS는 정규화된 면 값을 두 스킴의 혼합으로 둔다.

\widetilde{C}_{\text{face}} = \gamma \, \widetilde{C}_{\text{face}}^{\mathrm{CDS}} + (1 - \gamma) \, \widetilde{C}_{\text{face}}^{\mathrm{HR}}

CDS (compressive): 계면을 1셀 안으로 압축하지만 평면이 아니면 주름을 만든다.
HR (high resolution): 평탄해서 안전하지만 계면을 부풀린다.
블렌드 $\gamma = \cos^4\theta$ : $\theta$ 는 계면 법선 $\hat{\mathbf{n}}$ 과 donor–acceptor 방향 $\hat{\mathbf{d}}$ 의 각도. 계면이 면에 수직이면 $\theta=0$ , $\gamma=1$ → CDS. 계면이 면과 평행하면 $\theta=\pi/2$ , $\gamma=0$ → HR.

자연 발생한 게이트다. 평탄한 면에선 sharp CDS를 쓰고, 비스듬한 면에선 안전한 HR로 부드럽게 넘어간다.

연산자 분리 + SSP-RK3#

논문의 흐름도를 정리하면:

$C^{(n+1)}$ = OS-MSTACS( $C^{(n)}, \mathbf{u}^{(n)}$ ). x→y→z 순차 sweep, 매 step 마다 sweep 순서 교체.
$\rho^{(n+1)}, \mu^{(n+1)}$ = 혼합 규칙으로부터 업데이트.
$\kappa^{(n+1)}$ = height function으로 곡률 계산.
SSP-RK3 3단계로 $\mathbf{u}^{(n+1)}, p^{(n+1)}$ 계산. 각 단계에서 운동량 → GPE 순서.

SSP-RK3는

\Psi^{(1)} = \Psi^{(n)} + \Delta t L(\Psi^{(n)})

\Psi^{(2)} = \tfrac{3}{4}\Psi^{(n)} + \tfrac{1}{4}\Psi^{(1)} + \tfrac{1}{4}\Delta t L(\Psi^{(1)})

\Psi^{(n+1)} = \tfrac{1}{3}\Psi^{(n)} + \tfrac{2}{3}\Psi^{(2)} + \tfrac{2}{3}\Delta t L(\Psi^{(2)})

전체 알고리즘에서 선형 솔버 호출은 없다. 그래서 GPU에 그대로 옮길 수 있고, MPI 통신은 stencil halo만 보내면 된다.

직접 구현 — 1D 음향 펄스#

GPE의 작동 원리를 가장 작은 코드로 확인해보자. 점성을 끄면 1D 선형 음향계가 된다.

import numpy as np
 
def gpe_step_1d(p, u, cs, nu, dx, dt):
    """1D linear acoustics + GPE pressure-diffusion, periodic BC.
    p: 압력 fluctuation, u: 속도, cs: 인공 음속, nu: GPE 확산 계수
    """
    # 주기 경계 패딩
    pL = np.roll(p, 1); pR = np.roll(p, -1)
    uL = np.roll(u, 1); uR = np.roll(u, -1)
    # Rusanov 플럭스 (선형 시스템, eigenvalue ±cs)
    fp = 0.5*cs*cs*(uL + u) - 0.5*cs*(p - pL)   # face i-1/2
    fu = 0.5*(pL + p) - 0.5*cs*(u - uL)
    fp_r = 0.5*cs*cs*(u + uR) - 0.5*cs*(pR - p) # face i+1/2
    fu_r = 0.5*(p + pR) - 0.5*cs*(uR - u)
    # GPE 확산 항 (중심차분)
    diff = nu * (pL - 2*p + pR) / (dx*dx)
    p_new = p - dt/dx*(fp_r - fp) + dt*diff
    u_new = u - dt/dx*(fu_r - fu)
    return p_new, u_new
 
def run_demo(cs=5.0, nu=0.0, T=0.2, N=160):
    dx = 1.0/N
    dt = 0.4*dx/cs                # 음향 CFL
    x = (np.arange(N) + 0.5)*dx
    p = 0.4*np.exp(-((x - 0.5)/0.05)**2)
    u = np.zeros_like(x)
    steps = int(T/dt)
    for _ in range(steps):
        p, u = gpe_step_1d(p, u, cs, nu, dx, dt)
    return x, p, u
 
x, p, u = run_demo()                       # cs=5, no diffusion
x2, p2, u2 = run_demo(cs=15.0, nu=0.0)     # cs=15 → 음향파 3배 속도, dt 1/3로 감소

$c_s$ 를 5에서 15로 키우면 음향파는 3배 빨리 퍼지고, 안정성을 위한 $\Delta t$ 는 1/3로 줄어든다. 같은 물리 시간을 시뮬레이션하려면 3배 더 많은 step. GPE의 본질적 거래다.

인터랙티브 — $c_s$ 가 시간 스텝을 어떻게 잡아먹는가#

아래 시뮬레이션에서 직접 조작해보자. 슬라이더로 인공 음속과 압력 diffusion을 바꾸면 가운데 가우시안 펄스가 좌우 두 음향파로 갈라지며 전파한다.

artificial sound speed c_s: 5.0

pressure diffusion ν: 0.010

Acoustic Δt = 0.4·Δx/c_s ≈ 0.00050 | sim t = 0.000

관찰 포인트:

$c_s$ 를 5 → 20으로 올리면 음향파가 4배 빨리 도메인을 휩쓴다. 허용 $\Delta t$ 는 4배 줄어든다.
$\nu$ 를 0.05까지 올리면 음향파의 등마루가 빠르게 깎인다. EDAC의 인공 점성 항과 정확히 같은 역할.
$\nu = 0$ 이면 펄스가 주기 경계를 돌아 영구히 살아남는다. 비압축 한계에 가까울수록 좋아 보이지만, 거친 흐름에서 음파가 안 사라지면 spurious oscillation의 씨앗이 된다.

인터랙티브 — MSTACS의 블렌드를 풀어보기#

CDS, HR, MSTACS 세 면 보간 스킴이 같은 top-hat을 어떻게 다르게 옮기는지 직접 비교해보자.

Courant Cou_adv: 0.30

interface angle θ: 30° → γ = cos⁴θ = 0.563

steps: 120

Top-hat advected by U=0.5 with periodic BC. CDS = compressive (sharp but can wrinkle), HR = high-resolution (smoother), MSTACS = γ·CDS + (1−γ)·HR.

관찰 포인트:

CDS only: 모서리를 칼처럼 유지하지만 Courant ≤ 1/3 밖에선 살짝 튀어오른다.
HR only: 안정하지만 step이 점점 부드러워진다. 점성 같은 거짓 확산.
MSTACS with $\theta = 0$ °: $\gamma = 1$ → CDS와 동일.
MSTACS with $\theta = 90$ °: $\gamma = 0$ → HR과 동일.
실제 2D/3D에선 매 면마다 $\theta$ 가 다르게 잡혀, 평면 부분만 sharp하고 경사 부분은 안전한 HR로 풀린다.

검증 결과 — 2D RT, 댐 붕괴, 기포 상승, 3D RT#

케이스	$\rho_1/\rho_2$	비교 대상	정량 결과
2D Rayleigh–Taylor	3:1	He & Doolen (1999)	spike·bubble 위치 1% 이내
댐 붕괴 (broken dam)	1000:1	Martin & Moyce (1952) 실험	전면 진입 길이 5% 이내
단일 기포 상승	1000:1, $Eo=10$	Hysing 벤치마크 (2009)	terminal velocity 3% 이내
3D Rayleigh–Taylor	3:1	Tryggvason (1988)	spike 모양 정성적 일치

OS-MSTACS 자체의 부피분율 보존 오차는 $E_G = 6.57 \times 10^{-8}$ , 직전 OS-CICSAM의 $E_G = 8.18 \times 10^{-4}$ 보다 4 자릿수 개선됐다. 부피분율 redistribution이 mass conservation을 자릿수 정확도로 보장한 결과다.

비판 포인트 — 논문이 숨긴 비용#

GPE 그 자체가 가져온 부담이 분명히 있다.

시간 스텝의 음속 제약: $\Delta t \le \Delta x / c_s$ . $c_s$ 를 너무 작게 잡으면 인공 압축성 효과가 비물리적으로 새고, 너무 크게 잡으면 시간이 1/Ma 배로 늘어난다. 본문은 $\mathrm{Couacs} \le 1.25$ 까지 안전하다고만 한다.
GPE 확산항의 물리적 정당화: 우변의 $\frac{1}{\rho}\nabla \cdot (\mu \nabla p)$ 는 단상에서는 thermodynamic 유도가 깔끔하지만, 두 상에서 $\rho$ 가 계면을 가로질러 1000배 점프할 때 이 확산 계수가 무엇을 의미하는지 본문은 모호하다.
표면장력의 well-balanced 미보장: CSF + height function 곡률은 잘 알려진 표준 조합이지만 spurious current가 cs와 함께 어떻게 스케일하는지 정량적 분석이 없다. 정적 droplet 테스트만 정성 확인.
시리얼 비용 vs HPC 약속의 갭: $\Delta t_{\mathrm{GPE}} \approx \mathrm{Ma}\,\Delta t_{\mathrm{INS}}$ 이므로 1코어에서는 늘 손해. HPC 비교는 future work로 미뤘다.

OpenFOAM의 interFoam은 PIMPLE + MULES 조합으로 검증·튜닝이 충분히 쌓였다. 본 논문의 가장 큰 매력은 거기서 못 가는 GPU + iteration-free 영역인데, GPU 구현은 공개되지 않았다.

재현 가능성 점수#

알고리즘 의사코드 (Section 3.3): ★★★★★ — 8단계로 빈틈없이.
MSTACS 식 (13)–(15): ★★★★☆ — Cou_adv 분기에 typo 의심 (식 (13)이 두 번 등장). 결국 식 자체는 Anghan et al. 2018 원본을 봐야.
코드 공개: ★★☆☆☆ — 없음. 1D는 200줄, 3D OS-MSTACS는 1500줄 정도 추정.
벤치마크 데이터: ★★★★☆ — 4개 케이스 모두 정량 비교, 격자 수렴 표 포함.

일요일 오후, OpenFOAM 없이 가벼운 두 상 솔버를 처음부터 짜고 싶다면 이 논문이 의외로 친절한 출발점이다. Poisson 솔버 코드를 안 쓰는 대신 음속과 부피분율 redistribution 알고리즘과 친해질 것.

다음에 읽을 논문#

Toutant (2017) — GPE의 원조 유도, 단상.
Saincher & Sriram (2022) — operator-split 알고리즘 원조.
Kajzer & Pozorski (2021) — EDAC 기반 두 상 변종, 비교 대상.
Yang & Aoki (2021) — projection-based pressure evolution for VOF.