[論文レビュー] 音速が発散しても時間ステップは生き残る — Orlando

圧力だけを陰的に置くと、M=0.01、dt=0.5 dx の音響 CFL は 50 のままで、200 ステップ後の密度振幅の目標は約 0.10 です。

大気モデルはマッハ 0.001 です。宇宙の爆発模擬はマッハ 100 です。一つのソルバで両方を解きたいという願いは半世紀続いています。Orlando と Bonaventura の 2024 年論文(arXiv:2402.09252v4)は、その願いを非理想気体まで押し出します。時間離散で圧力項だけを陰的に置けば、時間ステップは音速から外れます。SG-EOS でも一般三次状態方程式(Van der Waals, Peng–Robinson)でも、その外れ方は残ります。

$M$ スライダーを 0.05 まで下げると音速 $c$ は 20 になり、同じ時間ステップでは陽的 CFL だけが 1 ステップで 1 を超えます。

Mach number M = 0.20 (c = 5.00)CFL_ac = 2.00 (explicit unstable)

drag M low → c grows → explicit dt-cap shrinks; IMEX keeps marching

上の赤い曲線は陽的 forward Euler、下のシアンは同じ方程式・同じ $\Delta t$ の AP-IMEX です。赤い領域が発散で覆われるとき、シアンはガウシアンパルスを左右へ流します。

30 秒サマリー#

著者: Giuseppe Orlando, Luca Bonaventura
所属: École polytechnique / Politecnico di Milano
arXiv: 2402.09252v4 (2025-10-22 v4)
目的: 任意のマッハ数・任意の EOS で漸近保存(AP)時間積分を実現する
空間離散化: Discontinuous Galerkin
検証: 等エントロピー渦・Gresho 渦・RT 不安定性・遷音速衝撃 — SG-EOS と三次 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 分の 1 です。大気・海洋の計算がよく落ちる場所です。

次は 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 です。馴染みの形に書くと次です。

\rho(\mathbf x, t) = \bar\rho(\mathbf x, t) + M\rho'(\mathbf x, t) + M^2 \rho''(\mathbf x, t) + \mathcal O(M^3)

$\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^{n+1} = \rho^n - \Delta t\, \nabla\!\cdot\!(\rho \mathbf u)^{n}

(\rho\mathbf u)^{n+1} = (\rho\mathbf u)^n - \Delta t\, \nabla\!\cdot\!(\rho\mathbf u\otimes\mathbf u)^n - \frac{\Delta t}{M^2}\nabla p^{n+1}

(\rho E)^{n+1} = (\rho E)^n - \Delta t\, \nabla\!\cdot\!\big[(\rho E + p)\mathbf u\big]^{n+1}

$\rho$ は密度、 $\mathbf u$ は速度、 $p$ は圧力、 $E$ は単位質量当たり総エネルギー、 $\Delta t$ は時間ステップ、 $M$ は基準マッハ数です。

$\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) — 液体や弱圧縮性固体でよく使う形です。

p = (\gamma - 1)(\rho e - \rho q_\infty) - \gamma \pi_\infty

$\gamma$ は比熱比、 $e$ は単位質量当たり内部エネルギー、 $q_\infty$ は結合エネルギー定数、 $\pi_\infty$ は剛性圧力定数です。 $q_\infty = \pi_\infty = 0$ で理想気体に戻ります。音速は $c^2 = \gamma(p + \pi_\infty)/\rho$ です。

一般三次 EOS — Peng–Robinson と Van der Waals を一式にまとめた形です。

p = \frac{\rho R T}{1 - \rho b} - \frac{a(T)\,\rho^2}{(1 - \rho b r_1)(1 - \rho b r_2)}

$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 で 1 ステップを追う#

論文の 1D 線形音響だけを NumPy で回すと、陰的 1 ステップの仕事が数字で見えます。周期格子、スタガードセル、フーリエ解で陰的ステップを処理します。

import numpy as np
 
def acoustic_imex_step(rho, u, dt, dx, M):
    """1D 線形音響系の AP-IMEX 1 ステップ。
 
    rho: セル中心の密度摂動、長さ 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)
    # 巡回 Helmholtz: (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-sqrt(2), r2=-1+sqrt(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}$ を残します。二つの数字の間は四桁です。

陰的求解を陽的 1 ステップ(u_new = rhs に置換)に変えると、最初の反復で発散します。音響 CFL が 50 なのに、陽的安定条件は 1 未満を要求するからです。

再現中に踏んだ落とし穴

論文 5.1 節の等エントロピー渦の表を描くと、4 次 IMEX が $M = 10^{-4}$ 付近で次数低下を見せます。硬い問題で高次 RK がよく起こす現象です(Hairer–Wanner, 1996)。著者らは速度多項式の次数を $r+1$ に一段上げて避けます。費用はかかりますが、実用になります。

もう一つは、DG の四角形格子で低マッハ精度が落ちることです。Jung–Perrier(2024)が示したとおり、三角形は持ちこたえますが四角形は圧力に $\mathcal O(M)$ のノイズが入ります。論文はこの限界を認め、後続(エントロピー安定 DG、compatible FE)へ回しています。実務で四角形を使うなら、低マッハ修正(Thornber, 2008; Rieper, 2011)を別に足す必要があります。

OpenFOAM の rhoPimpleFoam 系はすでに圧力–運動量分離(SIMPLE/PISO)を使いますが、理想気体以外の EOS にはコード修正が要ります。論文の三次 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 音響で読んだ数字だけです。

[論文レビュー] 音速が発散しても時間ステップは生き残る — Orlando–Bonaventura(2024) 非理想気体向け AP-IMEX