[논문 리뷰] 음속이 ∞로 가도 시간 스텝이 살아남는다 — 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)에서도 그 풀림은 유지된다.

$M$ 슬라이더를 0.05까지 내리면 음속 $c$가 20이 되고, 같은 시간 스텝에서 명시적 CFL만 한 스텝에 1을 넘긴다.

위 빨간 곡선은 명시적 forward Euler, 아래 시안은 같은 방정식과 같은 $\Delta t$의 AP-IMEX다. 빨간 영역이 폭주로 덮일 때 시안은 가우시안 펄스를 좌우로 흘린다.

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 룽게-쿠타(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로 돌리면, 음해 한 스텝이 무슨 일을 하는지 숫자로 보인다. 주기 격자, 엇갈림(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 진폭을 보존한다

주석의 CFL_acoustic = 50은 데모 숫자에서 바로 나온다. $M = 0.01$이면 무차원 음속은 $1/M = 100$이다. $dt = 0.5\,dx$이므로 $c\,dt/dx = 100 \times 0.5 = 50$이다. 양해 쪽은 이 값이 1 미만이어야 한다. u_new = rhs로 바꾸면 첫 스텝에서 터진다. 이유가 이 50이다.

$\alpha$ 식도 코드에 있다. $\alpha = (dt/(M\,dx))^2$에 $dt/dx = 0.5$, $M = 0.01$을 넣으면 $\alpha = 50^2 = 2500$이다. 고유값은 $1 + 2\alpha(1 - \cos(2\pi k/N))$이다. $k = 0$은 $\cos = 1$이라 고유값이 1이다. 평균 밀도는 건드리지 않는다. $N = 64$에서 가장 짧은 모드는 $k = 32$다. $\cos(\pi) = -1$이므로 고유값은 $1 + 4\alpha = 10001$이다. DFT 한 번에 그 모드를 10001로 나눈다. 초기 펄스는 진폭 0.1이고, 폭은 $\exp(-120(x-0.5)^2)$다. 200 스텝을 돌려도 print가 가리키는 목표는 $|\rho|_{\max} \approx 0.10$이다. 음향 CFL 50이 well-prepared 진폭을 깎지 못한 셈이다.

$M = 10^{-4}$로 가면 압력 변동은 $\mathcal{O}(M^2) = 10^{-8}$이어야 한다. 같은 자리에서 명시적 Roe는 $\mathcal{O}(M) = 10^{-4}$를 남긴다. 두 숫자 사이는 네 자리다.

같은 루프에 명시적 한 스텝(u_new = rhs로 교체)을 넣으면 첫 스텝에서 발산한다. CFL_acoustic이 50인데 명시적 안정 조건이 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 처리는 그런 확장에 그대로 옮겨 쓸 수 있다.

사각 격자를 남길지#

정리 3.4는 EOS 곡선을 가정하지 않는다. 사각 DG는 그 정리 밖에 있다. Jung–Perrier(2024)의 $\mathcal{O}(M)$ 압력 잡음은 시간 적분을 AP로 바꿔도 남는다. 격자를 삼각으로 못 바꾸면 Thornber(2008)나 Rieper(2011)를 플럭스에 얹는다. 엔트로피 안정 DG는 논문이 뒤로 민 선택이다.

음속 30%는 다른 자리다. 입방형 EOS를 이상기체로 바꾸면 $c$가 그 정도로 어긋날 수 있고, 셀이 음압이면 EOS 평가가 멈춘다. 시간 스텝을 음속에서 빼는 다른 절단은 Peluchon의 음향-수송 IMEX에 있다. 그 절단도 사각 격자의 잡음은 못 지운다. 데모의 200 스텝과 $|\rho|_{\max} \approx 0.10$은 주기 1D 음향에서만 읽은 숫자다.