音波と物質を別々に流す — 音響・対流分離（Lagrange–Projection）#

音速は秒速340メートルです。時速5キロで歩く人の横でも、音は同じ速度です。圧縮性ソルバーはその二つを一つの時間ステップに押し込みます。厄介なのは比です。遅い流れでは音波が物質より数百倍速く走ります。直接のGodunov系は、その速い波に時間ステップを縛られます。遅い対流は数値拡散でぼけます。

Machのスライダーを下げると、圧力パルスが割れる速さと、エントロピーの塊が止まる位置が一緒に変わります。

0.02まで下げるとamberのパルスは左右に割れて走り、tealの塊はほとんど動きません。三角三つは $u\!-\!c$, $u$, $u\!+\!c$ です。遅くするほど真ん中だけが止まり、外側の間隔は $2c$ のままです。

ten Eikelderら（2017）は、この二つをまるごと引き離します。音響と対流を別のソルバーで交互に解きます。ここは単相（single-phase）EulerにそのLagrange–Projectionを載せ、「遅い流れの中の音響パルス」で数字を合わせます。

“

論文: M.F.P. ten Eikelder, F. Daude, B. Koren, A.S. Tijsseling, "An acoustic-convective splitting-based approach for the Kapila two-phase flow model", Journal of Computational Physics 331 (2017) 188–208. DOI: 10.1016/j.jcp.2016.11.031

ヤコビアンが二つに割れる

1次元Euler方程式を原始変数 $\mathbf{W}=(\rho,u,p)^T$ で書くと、準線形の形になります。

ここで $\rho$ は密度、$u$ は速度、$p$ は圧力です。係数行列 $\mathbf{B}$ はきれいに二つに割れます。

$\mathbf{A}$ は圧力項をすべて含む音響（acoustic）部分です。$\mathbf{C}=u\mathbf{I}$ は物質をそのまま運ぶ対流（convective）部分です。この分離は実のところ、ラグランジュ微分 $\mathrm{D}/\mathrm{D}t = \partial_t + u\partial_x$ から $u\partial_x$ 項を剥がす作業に他なりません。

固有値の加法的分解

この分離がなぜ強力かは、固有値を見ると分かります。全系の波速は $u-c,\ u,\ u+c$ です。これが音響と対流にちょうど足し合わさります。

$c=\sqrt{\gamma p/\rho}$ は音速です。音響の波速は $\pm c$ で、流速には依りません。対流の波速はすべて $u$ です。低マッハ極限 $M=u/c\to 0$ では、音響帯は幅 $2c$ のまま、対流速度だけがゼロへ縮みます。この差が剛性（stiffness）の正体です。上の画面の三角の間隔が、この加法そのものです。$M\to 0$ では $\lambda^c$ だけが消え、$\lambda^a=\pm c$ は残ります。

音響ステップ — ラグランジュ座標のHLLC#

音響部分は、圧力が生む膨張・圧縮です。論文はこれを質量（ラグランジュ）座標へ移し、HLLC型リーマンソルバーで解きます。面 $j+1/2$ のスター状態は二つの式で済みます。

$a_{j+1/2}=\max(\rho_j c_j,\ \rho_{j+1}c_{j+1})$ は音響インピーダンス（密度×音速）です。このスター速度・圧力だけでEuler変数の更新が行われます。各セルの膨張率は一つの係数に圧縮されます。

$R_j$ は、音響ステップがセル体積をどれだけ伸ばしたり縮めたりしたかを表します。密度はただちに $\rho^{n+1-}_j = \rho^n_j / R_j$ として得られます。

実行の背景は $c_0=1$, $\rho_0=1$, $\gamma=1.4$ なので、$p_0=\rho_0 c_0^2/\gamma=1/1.4$ です。圧力パルスの振幅は $10^{-3}$、密度異常の振幅は $5\times 10^{-2}$ です。中心は $x=0.5$ と $x=0.25$ で、初期には重なりません。密度の最大はエントロピーの塊側、およそ $1.05$ から始まります。$R_j$ が 1 より大きければそのセルは膨張し、密度は下がります。1 より小さければ圧縮です。$\Delta t/\Delta x$ が面の $u^*$ の差に掛かるので、CFL 0.8 でも $R_j$ は 1 の近くに留まります。格子 400、領域長さ 1 なら $\Delta x=1/400$ です。

対流ステップ — 風上投影

音響ステップが生んだスター速度 $u^*$ で、今度は物質を運びます。保存量 $\varphi\in\{\rho,\rho u,\rho E\}$ を単純な風上法で更新します。

$\varphi^{\text{up}}_{j+1/2}$ は $u^*\!\ge\!0$ なら $\varphi_j$、そうでなければ $\varphi_{j+1}$ です。前ステップの $R_j$ がここで戻ってきて、質量・運動量・エネルギーが厳密に保存されます。二つのステップを順に踏めば、一つの時間ステップが完了します。

$u^*$ が正なら面の値は左セルです。負なら右セルです。$R_j\,\varphi^{n+1-}$ が先に体積変化を戻し、そのあと面フラックスが物質を運びます。順序が逆だと質量と $\rho u$、$\rho E$ の総和がずれます。周期境界なので面フラックスは領域の外へ漏れません。$R_j$ のフィードバックが合っていれば、運動量の総和に残るのは丸め誤差だけです。低マッハで素のRoeがどれだけ広がるかは、Hope–Collinsの人工拡散の三つのスケーリングが別に測っています。

Python — 遅い流れの中の音響パルス#

パイプライン全体をnumpyで組みます。周期境界、理想気体、セルあたり三つの保存量です。初期条件は、背景流 $u_0=Mc_0$ の上に小さな圧力パルスと密度異常を載せます。

import numpy as np
 
GAMMA = 1.4  # 理想気体の比熱比
 
def primitives(rho, mom, Ene):
    """保存量 -> 原始変数（速度、圧力、音速）。"""
    u = mom / rho
    e = Ene / rho - 0.5 * u * u            # 比内部エネルギー
    p = (GAMMA - 1.0) * rho * e
    c = np.sqrt(GAMMA * p / rho)
    return u, p, c
 
def acoustic_faces(rho, u, p, c):
    """面 j+1/2 のHLLC音響状態 u*, p*（論文 式34）。"""
    rp, up, pp, cp = (np.roll(a, -1) for a in (rho, u, p, c))
    a = np.maximum(rho * c, rp * cp)       # 音響インピーダンス a = max(rho*c)（式32）
    ustar = 0.5 * (u + up) + (p - pp) / (2 * a)
    pstar = 0.5 * (p + pp) + 0.5 * a * (u - up)
    return ustar, pstar
 
def lagrange_projection_step(rho, mom, Ene, dx, dt):
    """音響ステップ -> 対流ステップ（論文 式38, 40）。"""
    u, p, c = primitives(rho, mom, Ene)
    uf, pf = acoustic_faces(rho, u, p, c)          # j+1/2
    uf_m, pf_m = np.roll(uf, 1), np.roll(pf, 1)    # j-1/2
    lam = dt / dx
    # 1) 音響ステップ: 圧力が生む膨張/圧縮
    R = 1.0 + lam * (uf - uf_m)                     # 式 (39)
    rho1 = rho / R
    mom1 = (mom - lam * (pf - pf_m)) / R
    Ene1 = (Ene - lam * (pf * uf - pf_m * uf_m)) / R
    # 2) 対流ステップ: 物質を u* で運ぶ（風上投影）
    def project(phi1):
        phi_f = np.where(uf >= 0, phi1, np.roll(phi1, -1))
        phi_f_m = np.roll(phi_f, 1)
        return R * phi1 - lam * (uf * phi_f - uf_m * phi_f_m)
    return project(rho1), project(mom1), project(Ene1)
 
def run_acoustic_pulse(mach, n=400, cfl=0.8, tmax=0.25):
    x = (np.arange(n) + 0.5) / n
    dx = 1.0 / n
    c0, rho0 = 1.0, 1.0
    p0 = rho0 * c0**2 / GAMMA
    u0 = mach * c0
    dp = 1e-3 * np.exp(-((x - 0.5) / 0.03)**2)     # 音響圧力パルス
    ds = 5e-2 * np.exp(-((x - 0.25) / 0.03)**2)    # エントロピー(密度)異常 -> 接触波
    rho = rho0 + dp / c0**2 + ds
    u = np.full(n, u0)
    p = p0 + dp
    mom = rho * u
    Ene = p / (GAMMA - 1) + 0.5 * rho * u * u
    t, m0 = 0.0, mom.sum()
    while t < tmax:
        _, _, c = primitives(rho, mom, Ene)
        dt = min(cfl * dx / np.max(np.abs(mom / rho) + c), tmax - t)
        rho, mom, Ene = lagrange_projection_step(rho, mom, Ene, dx, dt)
        t += dt
    return rho, mom, m0, mom.sum()
 
for M in (0.02, 0.2):
    rho, mom, m0, m1 = run_acoustic_pulse(M)
    print(f"M={M}: 運動量誤差={abs(m1-m0)/abs(m0):.1e}, rho_max={rho.max():.4f}")
# M=0.02: 運動量誤差=2.2e-16, rho_max=1.0497
# M=0.2:  運動量誤差=0.0e+00, rho_max=1.0452

運動量は機械精度で保存されます。圧力は有界に、密度は正のままです。分離しても保存が崩れないのは、対流ステップに $R_j$ を戻し入れているからです。

出力の $\rho_{\max}$ は $M=0.02$ で 1.0497、$M=0.2$ で 1.0452 です。初期の峰 1.05 からほとんど削れていません。圧力パルスが密度に足すのは $dp/c_0^2=10^{-3}$ だけなので、最大を作るのは接触波です。運動量誤差 $2.2\times 10^{-16}$ と 0 は、分離そのものが保存を壊していない数字です。時間ステップを $|u|+c$ で切っているので、$M$ を 0.02 から 0.2 に上げてもステップ数はほぼ同じです。$c_0=1$ なら $|u|+c$ は 1.02 と 1.2 です。対流が十倍速くなっても、時間ステップの大半は音響が持っていきます。

低マッハで分かれるもの

分離が生きる場所は低マッハです。直接法の安定な時間ステップは、常に $|u|+c$ に縛られます。着目する物理が $u$ の上にあってもです。

Machのスライダーを動かすと、amberの帯の幅は $2c$ に固定されたまま、tealの $u$ だけが縮みます。

音響帯の幅は $M$ に依らず常に $2c$ です。ゼロへ収束するのは対流速度 $u$ だけです。$c=1$、$M=0.05$ なら $u=0.05$、$u+c=1.05$ で、比 $(u+c)/u$ は 21 です。分離法は対流ステップを音響と別の時間ステップで進められます。遅い物理が速い波に引きずられません。低マッハで直接Godunovが被る過剰な数値拡散も、この分離で避けられます。

批判的に見ると

引っかかる点が三つあります。時間精度は1次です。2次を得るにはStrang分離で音響–対流–音響の順に包む必要があり、その分コストが増します。原論文の舞台はKapilaの5方程式二相流です。体積分率方程式の非保存項 $K\partial_x u$ と正値性の保証が本当の難所で、単相への縮約はそこを飛び越えます。インピーダンス $a=\max(\rho c)$ は頑健ですが散逸的です。強い衝撃波では、この散逸が接触面をぼかし得ます。論文も直接法に対する精度と効率を別に比べています。

実務では、この分岐は目新しくありません。OpenFOAMの圧力ベース rhoPimpleFoam やall-Mach系ソルバーが、音響と対流を陰的・陽的に分けて扱うのと同じ根です。Lagrange–Projectionは、その分離をリーマンソルバーの言葉で書いた版です。

R_jを戻さないと保存が先に死ぬ#

音響ステップの $R_j$ を対流投影へ戻さないと、質量・運動量・エネルギーの更新が互いにずれます。上の実行で誤差が $2.2\times 10^{-16}$ まで落ちたのは、そのフィードバックがあるときだけです。Strangで包めば次数は上がります。順序は音響–対流–音響に伸び、パスが一つ増えます。

Kapilaの5方程式へ進むと、体積分率の $K\partial_x u$ が残ります。単相に縮めるとその非保存項には触れません。インピーダンスを $a=\max(\rho c)$ にすると衝撃には耐えますが、接触面は厚くなります。時間積分から音響を外す別の切り方は、Peluchonの音響-輸送IMEXです。$M=0.05$ で $(u+c)/u=21$ なら、その21倍を音響ソルバーに預けるか、対流CFLに残すかが実装の分岐です。