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

方程式を一つ消したら、コードがかえって壊れにくくなりました。圧縮性多相流(異なる流体が混ざる圧縮性の流れ)では、普通は逆に習います。相(phase・物質の異なる状態)ごとに圧力・速度・温度を別々に解く7方程式モデルが最も安全で、方程式を減らすほど不安定になる、と。Collis(2025)はその通念をひっくり返します。最も方程式の少ない4方程式モデルが、むしろ最も頑健だったのです。今日はその鍵となる界面平衡条件をたどります。なぜ「保存型」スキームが一様な圧力を揺らすのかを、Pythonで実際に確かめます。

平衡のはしご — 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$ のもとで1次風上を1ステップ進めます。 $\nu = u\,\Delta t/\Delta x$ はCFL数(1ステップで情報がセルを横切る割合)。保存型のエネルギー更新は $+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$ で割ります。二つの結果は界面で食い違う。その食い違いこそ圧力振動です。

下のシミュレーションで実際に操作してみましょう。

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.

Naiveをオンにすると、界面で赤い圧力線が跳ね上がります。right gas $\gamma$ を1.66(ヘリウム)や1.09(SF6)へ動かしてジャンプを大きくすると、振動も大きくなります。IEC-consistentに切り替えると、緑の線はステップ900まで平らなままです。

Python — 振動の正体と処方#

同じ計算をコードに移します。一様な $u,P$ で $\gamma$ だけがジャンプする界面を、1次風上で押します。異なるのは二つのスキームだけです。

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$ だけを扱いました。Collis(2025)の舞台ははるかに広い。

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

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

覚えておくこと

4方程式モデルは速度・圧力・温度を相の間ですべて共有します。方程式が最も少なく、物理に最も近い。
界面平衡条件(IEC):一様な圧力・速度は界面が通っても揺れてはならない。保存型スキームは混合物EOSのパラメータを誤って輸送し、この条件を破ります。
処方は、そのパラメータ( $1/(\gamma-1)$ )を平衡の仮定と整合的に輸送すること。Collis(2025)はこれを高次ENOへ引き上げ、余分な方程式なしに4方程式モデルを頑健にしました。