方程式を減らしたら、むしろ壊れなくなった — 4方程式モデルと界面平衡条件

一様な圧力が残るのはΓそのものを風上で運ぶときで、保存したρΓをρで割ると折れます。

方程式を一つ消したら、コードはむしろ壊れなくなりました。圧縮性多相流(異なる流体が混ざる圧縮性の流れ)では、ふつう逆に習います。相(phase、物質の異なる形態)ごとに圧力・速度・温度を別々に解く7方程式が最も安全で、式を減らすほど危うい、と。Collis(2025)はその順番をひっくり返します。最も短い4方程式の方が持ちこたえます。鍵は界面平衡です。

Naiveのトグルを入れて right gas のγを1.66か1.09へ押すと、界面の赤い圧力線が上がり、緑の線は平らなままです。

right gas γ = 1.66CFL = 0.60

Uniform velocity, uniform pressure, only γ jumps across the shaded interface. The naive scheme grows a pressure spike at the interface; the consistent one keeps P flat to machine precision.

赤い山はγのジャンプとともに大きくなります。900ステップでも緑が床に張り付いているかが、この条件が生きているかどうかです。

平衡のはしご — 7・6・5・4方程式#

拡散界面法(diffuse interface、境界面を数セルににじませて解く方法)には、にじんだセルの中で二相の状態を決める規則が要ります。その規則が「何を平衡とみなすか」です。

7方程式(Baer–Nunziato): 相ごとに速度・圧力・温度がすべて違います。最も一般的で、最も高いです。
6方程式: 速度は共有します。圧力は別に置き、緩和(relaxation)で合わせます。
5方程式: 速度と圧力を共有します。温度だけ別です。
4方程式: 速度・圧力・温度をすべて共有します。最も短いです。

物理的な界面の厚さはナノメートルです。数値界面はそれを格子数セルに膨らませたものです。その狭い帯の中では、熱・力学平衡は事実上すぐに終わります。界面付近の物理には4方程式が最も近いです。

T^{(1)} = T^{(2)}, \quad P^{(1)} = P^{(2)}, \quad \mathbf{u}^{(1)} = \mathbf{u}^{(2)}

上付き $(1),(2)$ は二相です。4方程式はこの三つの等式をセルの中で全部強制します。

方程式が少ないほどなぜ頑健か

5方程式以上は、たいてい体積分率(volume fraction、セル内である相が占める体積の割合)や比熱比に付く $1/(\gamma-1)$ のような余分な輸送式を追加します。PDEの水準では重複です。離散化すると互いにずれます。同じ量を二つの経路で計算して、値が割れます。

4方程式はその余分を外します。質量・運動量・エネルギーだけを保存します。圧力は混合状態方程式で閉じます。式が少ないので、ずれる場所も少ないです。費用も下がります。

問題は一つです。余分な式が持っていた余裕が消えると、保存型の離散化が界面で圧力をずらします。

界面平衡条件が破れる瞬間

界面平衡条件(IEC、interface equilibrium condition)は単純です。圧力・温度・速度が最初に一様なら、界面が通過してもその一様さは破れてはいけません。

破ったスキームは、界面の周りに無かった圧力振動を作ります。振動は時間とともに育ちます。密度比が大きいと、負の密度や負の内部エネルギーが出て、コードが止まります。

なぜ保存型がこれを破るのか。理想気体の混合で $\Gamma \equiv 1/(\gamma-1)$ と置くと、全エネルギーは次の形です。

E = P\,\Gamma + \tfrac{1}{2}\,\rho\,u^2

$E$ は全エネルギー、 $P$ は圧力、 $\rho$ は密度、 $u$ は速度です。 $\Gamma$ は流体ごとに違い、界面でジャンプします。

一様な $u,P$ のもとで一次風上を一セル進めます。 $\nu = u\,\Delta t/\Delta x$ はCFL数です。保存型のエネルギー更新は $+P$ の項が打ち消され、次のように残ります。

E_i^{n+1} = E_i - \nu\,(E_i - E_{i-1})

全エネルギーの定義を代入すると、右辺が割れます。

E_i^{n+1} = P_0\big[\,\Gamma_i - \nu(\Gamma_i - \Gamma_{i-1})\,\big] + \tfrac{1}{2}\,u_0^2\,\rho_i^{n+1}

新しい圧力が $P_0$ のままであるには、 $\Gamma$ があの括弧どおりに更新されなければなりません。非保存の風上です。保存型は $\rho\Gamma$ をまとめて運んだあと、再び $\rho$ で割ります。二つの結果は界面で割れます。一様な圧力なのに非保存経路がずれると圧力が揺れる、その場所がAbgrallが残した基準です。その差が圧力振動です。

先ほどγを動かしたとき赤が伸びたのが、この差です。IEC-consistentに切り替えると、緑は900ステップでも平らです。

Python — 振動の正体と処方#

同じ計算です。一様な $u,P$ で、 $\gamma$ だけがジャンプする界面を一次風上で進めます。違うのは二つのスキームだけです。

import numpy as np
 
def advect_interface(scheme, gamma_R=1.66, nu=0.6, steps=300, N=200):
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    ramp = np.clip((x - 0.28) / 0.06, 0.0, 1.0)   # 数セルににじんだ界面
    gamma = 1.4 + (gamma_R - 1.4) * ramp          # 空気(1.4) -> 選んだ気体
    rho = 1.0 + (0.2 - 1.0) * ramp                # 密度もジャンプする
    u0, P0 = 1.0, 1.0
    Gam = 1.0 / (gamma - 1.0)                      # 混合物の剛性 1/(gamma-1)
    mom = rho * u0
    E = P0 * Gam + 0.5 * rho * u0**2               # E = P*Gamma + 0.5*rho*u^2
 
    for _ in range(steps):
        rhoG = rho * Gam                           # 更新前の rho*Gamma
        r_up, m_up, e_up = np.roll(rho, 1), np.roll(mom, 1), np.roll(E, 1)
        rho = rho - nu * (rho - r_up)              # 質量:保存型風上
        mom = mom - nu * (mom - m_up)              # 運動量:保存型
        E   = E   - nu * (E   - e_up)              # エネルギー:保存型
        if scheme == "naive":
            rhoG = rhoG - nu * (rhoG - np.roll(rhoG, 1))  # rho*Gamma を保存輸送
            Gam = rhoG / rho                              # 再び rho で割る
        else:  # consistent
            Gam = Gam - nu * (Gam - np.roll(Gam, 1))      # Gamma を直接、非保存風上
 
    P = (E - 0.5 * mom**2 / rho) / Gam             # 圧力を復元
    return float(np.max(np.abs(P - P0)))
 
for s in ("naive", "consistent"):
    print(f"{s:11s}  max|P-P0| = {advect_interface(s):.3e}")

出力は二桁ではなく、十三桁まで割れます。

naive        max|P-P0| = 3.83e-02
consistent   max|P-P0| = 4.44e-15

naiveは界面で圧力が約4%ずれます。consistentは浮動小数の丸めまで平らです。差を作ったのは一行です。 $\Gamma$ を rho で割って取り戻すか、 $\Gamma$ そのものを運ぶかです。

$\Gamma$ を数字にすると、ジャンプの背丈が見えます。 $\Gamma = 1/(\gamma-1)$ なので、空気 $\gamma=1.4$ は $\Gamma=2.5$ です。ヘリウム $\gamma=1.66$ は $\Gamma=1/0.66 \approx 1.515$ です。SF6 $\gamma=1.09$ は $\Gamma=1/0.09 \approx 11.11$ です。空気からヘリウムへの差は約0.985です。空気からSF6へは約8.61です。SF6のジャンプはヘリウムの約8.74倍です。γを1.09へ押すと赤が伸びる倍率です。

鋭い一セルを、コードの $\nu=0.6$ で一度だけ進めた場合です。右セルは $\rho=0.2$ 、 $\Gamma\approx 1.515$ で、風上の隣は空気 $\rho=1$ 、 $\Gamma=2.5$ です。 $\rho\Gamma$ を保存型で運んでから $\rho$ で割ると $\Gamma_{\mathrm{naive}}\approx 2.384$ です。 $\Gamma$ を直接進めると $\Gamma_{\mathrm{cons}}\approx 2.106$ です。上のエネルギー式を圧力に戻すと $P/P_0 = \Gamma_{\mathrm{cons}}/\Gamma_{\mathrm{naive}} \approx 0.883$ なので、一ステップの $|P-P_0|\approx 0.117$ です。同じ一セルをSF6にすると $|P-P_0|\approx 0.692$ です。出力の $3.83\times 10^{-2}$ は、この鋭いセルではありません。ランプ幅 $0.06$ は $\Delta x=1/200$ で12セルであり、300ステップのあいだ界面は $\nu \cdot 300 \cdot \Delta x = 0.9$ だけ動きます。にじんだ界面を長く進めると、最大誤差は一ステップのジャンプより小さくなり、 $0.0383$ に留まります。鋭い一ステップの $0.117$ とこの $0.0383$ の比は約 $3.05$ です。12セルに割れても、誤差は一桁パーセントに残ります。consistentの $4.44\times 10^{-15}$ は、経路の差が0なので丸めだけが残った値です。

論文が実際に行ったこと

上のトイは理想気体の $\gamma$ を一つ押しただけです。Collis(2025)の舞台はもっと広いです。

NASG状態方程式(Noble–Abel Stiffened Gas、液体まで扱う状態式)で水・空気・ヘリウム・SF6をまとめて閉じます。混合圧力はAmagatの法則の二次式で解けます。
一次風上の代わりに ENO系(WENO/TENO)の高次再構成を使います。その再構成を4方程式の熱・力学平衡と噛み合わせて、IECを機械精度で満たします。
強い衝撃が大きい密度比の界面を叩く場合のために、正値保存リミッター(positivity-preserving limiter)を4方程式へ広げます。密度と音速の二乗が負へ漏れないように止めます。
余分な体積分率・比熱比の方程式なしに、保存性もスカラーフィルタも落とさずに、これを達成しています。

限界もあります。液体は単一成分、気体は理想気体成分に限っています。成分が増えると、圧力・温度の平衡を反復ソルバで解く必要があります。まだプレプリント(査読前)です。

Γをρで割って取り戻すと、平らな圧力が一ステップで折れる

鋭い一セルでは、ヘリウムがすでに $|P-P_0|\approx 0.117$ です。SF6は $0.692$ です。300ステップの出力 $0.0383$ は、そのジャンプを12セルのランプが分け持った結果です。 $\rho\Gamma$ を保存型で運んでから $\rho$ で割る一行を残すだけで、一様だった $P_0$ は界面ですぐに折れます。式を減らした代わりに、この一行を非保存へ替える選択が、4方程式の中身です。

5方程式は速度と圧力を共有したまま、余分な式を一つ残します。その余分が界面の音速を崩す側なら、IECのために式を消したことが病も消したのか、混合物の状態式の中へ隠しただけなのかが残ります。水と空気を半分ずつ混ぜたときに音速が崩れる5式が、4方程式の閉包の中にまだいるのかが次の分かれ目です。