[論文レビュー] 低マッハで音を殺さない方法 — AP-IMEX波動スキーム

低マッハ（low Mach）の流れを圧縮性コードに入れると、音速が速いほど絵が鮮明になりそうに思えます。逆です。マッハ数が下がるほど、標準の風上（upwind）は健全な波を格子から消します。Arun・Das Gupta・Samantaray（2019）は線形波動方程式系でその崩壊を切り分け、IMEX-RK でマッハ数に依らない安定な漸近保存（AP）スキームを組みました。音響演算子を陰的に解く理由は、下の数字に出ています。

$\varepsilon$ のスライダーを左へ動かすと、赤い風上曲線は傾き $-1$ で立ち上がり、シアンの AP 曲線は平らなままです。

$\varepsilon=10^{-3}$ では二つの曲線の比が数百倍まで開きます。マッハ数を下げて大きくなるのは波ではなく、風上の数値粘性です。

ε が 0 に近づくと格子が波を食べる#

スケールを整えた等エントロピー Euler 方程式には、小さな数が一つ入っています。マッハ数 $\varepsilon$ です。圧力項は運動量方程式に $\nabla p / \varepsilon^2$ として入ります。

$\rho$ は密度、$\mathbf{u}$ は速度、$p=\rho^\gamma$ は圧力です。$\varepsilon$ は基準速度を基準音速で割った値です。

$\varepsilon \to 0$ で解は非圧縮極限へ向かいます。特異摂動（singular perturbation）です。著者らは解析を短くするため、線形波動系を使います。

$\varrho$ は線形化した密度摂動、$a$ はその状態の音速です。実効音速は $c=a/\varepsilon$ です。マッハ数が小さくなると、この音響速度が爆発します。

well-prepared 空間 — 生き残るべき解#

極限 $\varepsilon\to 0$ で解は好きな方向へ漏れません。準備の整った（well-prepared）空間に集まります。密度は定数、速度の発散は 0 です。

この空間が音響演算子 $L(U)=(a\nabla\cdot\mathbf{u},\,a\nabla\varrho)$ の核（kernel）です。離散解も同じ場所へ着地すべきです。Dellacherie（2010）は、この不変性を破るスキームが低マッハで精度を失うことを示しました。核を保てば漸近精度（AA）です。安定制限が $\varepsilon$ に依らなければ漸近保存（AP）です。

風上粘性は音速に比例して大きくなる

犯人は風上の数値粘性です。波動系を特性変数で対角化すると、速度 $\pm c$ のスカラー移流が二つになります。各モードに Rusanov 風上を掛けると、1ステップの増幅が出ます。

$\nu$ は音響 CFL、$\theta$ はセルあたりの位相、$k$ は波数です。

1ステップの損失は小さいです。積み上がる方が問題です。物理時間 $T$ まで行くには音響ステップが $M = T/\Delta t = T c /(\nu\Delta x)$ 回必要です。$\theta$ が小さいとき 1ステップの減衰は $\tfrac{1}{2}\nu(1-\nu)\theta^2$ で、$M$ 回重ねると対数減衰はこうなります。

マッハ数を一桁下げると、実効粘性も一桁大きくなります。上の図で赤い曲線の傾きが $-1$ なのは、この比例です。

演算子分割 — 対流は陽的に、音響は陰的に

方程式を時間スケールで割ると、演算子は二つです。

$H$ は時間スケール $O(1)$ の対流です。$L/\varepsilon$ は時間スケール $O(\varepsilon)$ の音響です。硬さ（stiffness）はすべて $L/\varepsilon$ にあります。その部分だけ陰的に解きます。対流は陽的に残します。この混合が IMEX（Implicit-Explicit）です。同じ分離をマッハ 0 まで持っていったとき時間刻みが持つかは、IMEX TVD の再現が続けて扱います。

音響を陰的にすれば中心差分でも安定です。風上粘性は外れます。$\Delta t$ を決めるのは速い音響ではなく遅い対流です。このモデルの中では、$\varepsilon$ が小さくても刻みを縮める理由がありません。

IMEX-RK が硬さを受け止める#

時間積分は IMEX ルンゲ・クッタ（IMEX-RK）です。陽的表 $(\tilde A,\tilde b)$ が $H$ を、陰的表 $(A,b)$ が $L/\varepsilon$ を受け持ちます。各ステージで音響は、密度に対する楕円型方程式へ凝縮されます。圧力ポアソンに似た形です。

$\beta$ は IMEX 表の対角係数、$\nabla^2$ はラプラシアンです。この楕円問題の可逆性は、変分の鞍点（saddle point）で保証されます。離散では循環行列（circulant matrix）が同じ結論を与えます。一意解、$\varepsilon$ に依らない一様安定、well-prepared 空間の不変性。三つです。

Python — 1モードの振幅減衰を再現する#

下のコードはフーリエモード一つの風上増幅を計算し、物理時間 $T$ の後に残る振幅をマッハ数ごとに出します。陰的スキームはモードを保存するので、振幅は 1 です。

import numpy as np
 
a  = 1.0        # 線形化音速
k  = 2 * np.pi  # [0,1] 上の単一フーリエモード
T  = 1.0        # 固定物理終了時刻
nu = 0.5        # 音響CFL数
 
def upwind_amp_factor(eps, N):
    """波動系Rusanovにおける1特性モードの増幅 |g|。"""
    c  = a / eps              # 音響速度 c = a/ε
    dx = 1.0 / N
    th = k * dx               # セルあたり位相 θ = kΔx
    re = 1.0 - nu * (1.0 - np.cos(th))
    im = nu * np.sin(th)
    return np.hypot(re, im)   # |g| <= 1
 
def retained_amplitude(eps, N):
    """物理時間 T 後のモード振幅（風上 vs 陰的）。"""
    c  = a / eps
    dx = 1.0 / N
    dt = nu * dx / c          # 音響CFLステップ
    M  = T / dt               # 陽的ステップ数
    explicit = upwind_amp_factor(eps, N) ** M
    implicit = 1.0            # 中心-陰的はモードを保存
    return M, explicit, implicit
 
if __name__ == "__main__":
    N = 64
    print(f"{'eps':>8} {'steps M':>10} {'explicit':>12} {'AP':>6}")
    for eps in [0.4, 0.2, 0.1, 0.05, 0.02, 0.01]:
        M, ex, ap = retained_amplitude(eps, N)
        print(f"{eps:8.3f} {M:10.0f} {ex:12.3e} {ap:6.2f}")

出力はこうなります。

     eps    steps M     explicit     AP
   0.400        320    6.804e-01   1.00
   0.200        640    4.629e-01   1.00
   0.100       1280    2.143e-01   1.00
   0.050       2560    4.591e-02   1.00
   0.020       6400    4.430e-04   1.00
   0.010      12800    1.980e-07   1.00

Mach $\varepsilon$ を下げると、赤い振幅は物理時間 $T=1$ に届く前に床へ着きます。CFL $\nu$ を 1 側へ寄せると $\nu(1-\nu)$ が小さくなり、減衰は一瞬和らぎます。$\varepsilon$ をさらに下げると、その利得は $1/\varepsilon$ の累積に埋もれます。$N$ を大きくしても $\theta^2\propto N^{-2}$ と $\Delta x\propto N^{-1}$ は一部しか打ち消し合いません。

コードの数を式へ入れると、ステップ数は閉じた形になります。$a=1$、$T=1$、$\nu=0.5$、$N=64$ なら $\Delta x=1/64$ で $c=1/\varepsilon$ です。$\Delta t=\nu\Delta x/c$ なので $M=T/\Delta t=128/\varepsilon$ です。表の 320、640、1280、2560、6400、12800 は、$\varepsilon=0.4,0.2,0.1,0.05,0.02,0.01$ をこの式に入れた値です。$\varepsilon$ を半分にすると $M$ は 2 倍です。10 分の 1 にすると $M$ も 10 倍です。

対数減衰 $\tfrac12 T c(1-\nu)k^2\Delta x$ に同じ数を入れると、指数は $0.154/\varepsilon$ です。$k=2\pi$ だから係数がそうなります。$\varepsilon=0.4$ で指数は 0.386、$\mathrm{e}^{-0.386}=0.680$ です。出力の $6.804\times 10^{-1}$ と小数三位まで一致します。$\varepsilon=0.1$ では $\mathrm{e}^{-1.542}=0.214$ で、出力の $2.143\times 10^{-1}$ と合います。$\varepsilon=0.01$ で指数は 15.4 です。$\mathrm{e}^{-15.4}$ は $2.0\times 10^{-7}$ の近くで、出力は $1.980\times 10^{-7}$ です。$\theta=k\Delta x=2\pi/64\approx 0.098$ なので、小さな角度の仮定は表の一番下まで正確ではありません。表の端が近似より少し余計に死ぬ理由です。$\nu=0.5$ は $\nu(1-\nu)=0.25$ で、この積が最大になる CFL です。だからこの列は、1ステップの減衰がいちばん強い条件の記録です。

漸近保存（AP）と漸近精度（AA）#

二つは別の問いです。AP は、離散スキームの $\varepsilon\to 0$ 極限が極限方程式の正しい離散化かを問います。安定制限は $\varepsilon$ に依ってはいけません。AA はもう一歩です。格子は有限で、$\varepsilon$ は小さいが 0 ではないとき、解が well-prepared 空間の近くに留まるかです。風上は両方に失敗します。IMEX の陰解法は両方を通ります。分かれる場所は空間差分の次数ではなく、時間離散化です。

ν=0.5 の最悪減衰でも AP 列は 1.00 です#

この表の AP 列は六行とも 1.00 です。陰解法が風上粘性を外すので、指数 $0.154/\varepsilon$ が振幅に掛からないからです。$\nu=0.5$ は 1ステップの減衰因子が最大なのに、$M=128/\varepsilon$ 回積んだあとも 1.00 です。この線形波動系では、$L/\varepsilon$ だけを陰的にする選択が合っています。

その選択が低マッハの方程式すべてへそのまま移るわけではありません。音響を陰的に入れても、時間刻みが再び大きく縮む系があります。相対速度に二つ目の $1/M^2$ が残る場所は、エンタルピーに残った硬さに書いてあります。