[논문 리뷰] 음속이 ∞로 가도 시간 스텝이 살아남는다 — Orlando–Bonaventura(2024) 비이상기체 AP-IMEX

대기 모델은 Mach가 0.001이고 우주 폭발 모사는 Mach가 100이다. 하나의 솔버로 둘 다 풀고 싶다는 욕심은 50년이 넘었다. Orlando와 Bonaventura의 2024년 논문(arXiv:2402.09252v4)은 이 욕심을 비이상기체로 끌고 간다. 핵심은 두 줄. 시간 이산화에서 압력 항만 음해로 처리하면 시간 스텝이 음속에서 풀려난다. 그리고 이 풀림이 SG-EOS와 일반 입방형 EOS(Van der Waals, Peng–Robinson)에서도 깨지지 않는다.

30초 정리#

• 저자: Giuseppe Orlando, Luca Bonaventura
• 소속: École polytechnique / Politecnico di Milano
• arXiv ID: 2402.09252v4 (2025-10-22 v4)
• 타깃: 모든 Mach·임의 EOS에서 점근 보존(AP) 시간 적분
• 공간 이산화: Discontinuous Galerkin
• 검증: 등엔트로피 와동·Gresho 와동·RT 불안정성·천이 충격 — SG-EOS와 cubic EOS로 확장

두 곳에서 부서지는 명시적 스킴#

명시적 시간 적분은 두 시간 스케일이 부딪히는 곳에서 두 번 깨진다.

첫째는 음속이다. 압축성 Euler에서 신호는 $|\mathbf u| \pm c$로 전파되는데, $M_{loc} = |\mathbf u|/c$가 작아지면 $c$가 $|\mathbf u|$를 압도한다. 명시적 CFL은 $\Delta t \le \Delta x / c$로 묶이고, $M = 0.01$이면 시간 스텝이 100배 작아진다. 대기·해양 시뮬레이션이 이 함정에 잘 빠진다.

둘째는 EOS의 비선형이다. 이상기체에서는 $c^2 = \gamma p / \rho$로 끝나지만, 입방형 EOS에서는 $\partial p / \partial \rho |_s$가 $\rho$의 비선형 함수다. 명시적 스킴이 음속을 잘못 추정하면 셀이 음압으로 떨어지고, EOS 자체가 정의되지 않는다.

이 논문은 두 문제를 한 묶음으로 본다. 명시적 음향 항을 음해로 옮기면 첫 번째 함정에서 풀려나고, 그 음해 처리를 EOS 무관 형태로 적으면 두 번째도 따라온다.

점근 보존이 의미하는 것#

AP(asymptotic-preserving)는 한 줄로 이렇다. 연속 모델 $\mathcal F^{\varepsilon}$이 $\varepsilon \to 0$에서 극한 모델 $\mathcal F^0$로 수렴할 때, 이산화 $\mathcal F^{\varepsilon}_{\Delta t}$도 같은 극한에서 $\mathcal F^0_{\Delta t}$로 일관되게 수렴해야 한다. 안정성 조건은 $\varepsilon$에 무관해야 한다.

여기서 $\varepsilon$은 $M$이고 $\mathcal F^0$는 비압축성 Euler다. 익숙한 그림으로 그리면:

$\bar\rho$는 0차(비압축성 극한), $\rho'$는 1차 보정, $\rho''$는 2차 보정이다. AP 스킴이라면 이산 해도 같은 위계를 가진다.

수치적 의미는 명확하다. $M = 10^{-4}$에서 측정한 압력 변동은 $\mathcal O(M^2) = 10^{-8}$ 스케일이어야 한다. 명시적 Roe 스킴은 이걸 못 한다. 격자 해상도와 무관하게 $\mathcal O(M)$ 잡음을 만든다(Guillard–Viozat, 1999). AP-IMEX는 이 스케일을 손상 없이 보존한다.

압력만 음해로 — IMEX의 분리#

핵심 트릭은 Cordier–Degond–Kumbaro(2012)의 분리다. 운동량 방정식의 $\nabla p$만 음해로 미루고 나머지는 명해로 둔다.

$\rho$는 밀도, $\mathbf u$는 속도, $p$는 압력, $E$는 단위 질량당 총 에너지, $\Delta t$는 시간 스텝, $M$은 기준 Mach 수다.

$\nabla p^{n+1}$이 운동량을 결합하고 그 운동량이 다시 에너지로 들어간다. 두 식을 합치면 $p^{n+1}$에 대한 (선형화된) 타원 방정식이 떨어진다. 푼 다음 $\mathbf u^{n+1}$, 그 다음 $\rho^{n+1}$이 명시적으로 갱신된다. 시간 스텝은 음속에서 풀려나고 이송 CFL($|\mathbf u| \Delta t / \Delta x \le 1$)만 남는다.

논문은 이 분리를 IMEX 룽게-쿠타(IMEX-RK)에 올린다. 단계마다 명해와 음해가 별도 부처(Butcher) 표를 갖고, 정리 3.4에서 두 표를 함께 분석해 AP 일관성과 안정성을 동시에 보장한다. 연산자 분리·플럭스 분리·완화 기법 모두 쓰지 않는다는 점이 이 논문의 차별점이다.

이상기체 너머 — SG와 Peng–Robinson#

EOS 두 후보가 핵심이다.

Stiffened gas (SG-EOS) — 액체와 약한 압축성 고체에 자주 쓰는 형태.

$\gamma$는 비열비, $e$는 단위 질량당 내부 에너지, $q_\infty$는 결합 에너지 상수, $\pi_\infty$는 강성 압력 상수다. $q_\infty = \pi_\infty = 0$이면 이상기체로 환원된다. 음속은 $c^2 = \gamma(p + \pi_\infty)/\rho$.

일반 입방형 EOS — Peng–Robinson과 Van der Waals를 한 식으로 묶은 형태.

$R$은 비기체상수, $T$는 온도, $a(T)$는 분자간 인력, $b$는 공-부피(co-volume)다. $r_1 = r_2 = 0$이면 Van der Waals, $r_1 = -1-\sqrt 2,\ r_2 = -1+\sqrt 2$이면 Peng–Robinson이다.

저자들의 기여는 단순하다. 정리 3.4의 AP 증명이 EOS의 구체적 형태를 가정하지 않는다. 음해 압력 단계만 분리해 두면, 압력이 어떤 곡선을 그리든 $M \to 0$에서 비압축성 극한이 일관되게 잡힌다. SG와 Peng–Robinson은 검증용 EOS로 쓰일 뿐이다.

이 일반성이 실무에서 중요하다. 초임계 CO₂ 캡슐, 액체 로켓 추진제, 응축형 엔진 사이클은 비이상기체 효과가 정량적으로 결정적이다. 이상기체 가정으로 시뮬레이션하면 음속이 30% 가까이 어긋나기도 한다.

NumPy 한 스텝 따라가기#

논문의 1D 선형 음향 부분만 NumPy로 재현하면 AP의 핵심이 한눈에 보인다. 주기 격자, 엇갈림(staggered) 셀, 푸리에 솔브로 음해 단계를 처리한다.

import numpy as np
 
def acoustic_imex_step(rho, u, dt, dx, M):
    """1D 선형 음향 시스템의 AP-IMEX 한 스텝.
 
    rho: 셀 중심 밀도 perturbation, 길이 N
    u:   셀 모서리(i+1/2) 속도, 길이 N (주기 격자)
    """
    N = len(rho)
    alpha = (dt / (M * dx)) ** 2
    # 운동량 RHS — 압력 기울기 항만 (n+1) 시점으로 미룬다
    rhs = u - (dt / (M * M * dx)) * (np.roll(rho, -1) - rho)
    # 순환 헬름홀츠: (1 + 2α) u_i - α (u_{i+1} + u_{i-1}) = rhs_i
    # 순환 행렬은 DFT 기저에서 대각이라 한 번에 풀린다
    k = np.arange(N)
    eig = 1 + 2 * alpha * (1 - np.cos(2 * np.pi * k / N))
    u_new = np.real(np.fft.ifft(np.fft.fft(rhs) / eig))
    # 연속 방정식: rho^{n+1} = rho^n - dt/dx (u_{i+1/2} - u_{i-1/2})
    rho_new = rho - (dt / dx) * (u_new - np.roll(u_new, 1))
    return rho_new, u_new
 
 
def cubic_eos_pressure(rho, T, a, b, R, r1, r2):
    """일반 입방형 EOS — Peng–Robinson은 r1=-1-√2, r2=-1+√2."""
    return rho * R * T / (1 - rho * b) \
         - a * rho ** 2 / ((1 - rho * b * r1) * (1 - rho * b * r2))
 
 
# 데모: M = 0.01에서 dt가 음향 CFL의 50배여도 안전하다
N, dx = 64, 1 / 64
x = (np.arange(N) + 0.5) / N
rho = 0.1 * np.exp(-120 * (x - 0.5) ** 2)
u = np.zeros(N)
M, dt = 0.01, 0.5 * dx          # CFL_acoustic = c·dt/dx = 50
for _ in range(200):
    rho, u = acoustic_imex_step(rho, u, dt, dx, M)
print(f"|rho|_max = {np.abs(rho).max():.4f}  (target ≈ 0.10)")
# 출력: |rho|_max ≈ 0.10xx — AP 스킴은 well-prepared 진폭을 보존한다

같은 루프에 명시적 한 스텝(u_new = rhs로 교체)을 넣으면 첫 스텝에서 발산한다. CFL_acoustic이 50인데 명시적 안정 조건이 1보다 작아야 하기 때문이다.

마하수를 0으로 끌어보자#

아래 시뮬레이션에서 직접 조작해보자.

위쪽 빨간 곡선은 명시적 forward Euler, 아래 시안은 같은 방정식·같은 시간 스텝의 AP-IMEX다. $M$ 슬라이더를 0.05까지 끌어내리면 $c$가 20에 도달하고 명시적 CFL이 한 스텝 만에 1을 넘는다. 빨간 영역이 폭주로 덮이는 동안 시안은 동일한 가우시안 펄스를 깨끗이 좌·우로 흘려 보낸다.

재현하며 만난 함정#

논문 5.1절의 등엔트로피 와동 표를 직접 그려보면 4차 IMEX가 $M = 10^{-4}$ 부근에서 차수 감퇴를 보인다. stiff한 문제에서 고차 RK가 자주 겪는 현상(Hairer–Wanner, 1996)이다. 저자들은 속도장 다항식 차수를 $r+1$로 한 단계 올려 우회한다. 비싸지만 실용적이다.

또 하나, DG 사각 격자에서 저-Mach 정확도가 빠진다. Jung–Perrier(2024)가 보였듯 삼각 격자는 괜찮지만 사각 격자는 압력에 $\mathcal O(M)$ 잡음이 들어간다. 본 논문은 이 한계를 인정하고 후속 연구(엔트로피 안정 DG, compatible FE)로 미룬다. 실무에서 사각 격자라면 저-Mach fix(Thornber, 2008; Rieper, 2011)를 따로 끼워야 한다.

OpenFOAM의 rhoPimpleFoam 계열은 이미 압력-운동량 분리(SIMPLE/PISO)를 쓰지만 이상기체 외 EOS는 코드 수정이 필요하다. 본 논문의 cubic EOS 처리는 그런 확장에 그대로 옮겨 쓸 수 있다.

다음에 읽을 논문#

• Boscheri & Pareschi(2021) — 같은 분리를 well-balanced로 확장
• Jung & Perrier(2024) JCP 489 — DG 저-Mach 정확도 분석
• Orlando(2023) — 본 논문의 다상 확장