[論文レビュー] 圧力を二つに分けて解き、一つに戻す — Saurel(2009) 圧縮性多相流の緩和法

5方程式拡散界面モデルの罠を、6方程式の圧力非平衡 + stiff緩和で迂回したSaurel・Petitpas・Berry(2009)のレビュー。

アイダホ国立研究所での核燃料安全シミュレーションが止まりました。液体ナトリウムに挟まれた小さなガス気泡セルで、混合音速が突如25 m/sまで落ち込んだのです。圧力平衡を仮定した5方程式多相流モデルの限界でした。Richard Saurelらは2009年、一歩引いた答えを示します — 二相の圧力をあえて別々に解いておき、各タイムステップの最後に強制的に揃えるのです。本稿では、その6方程式非平衡 + 緩和アルゴリズムがなぜ機能するのか、界面圧力 $\hat{p}_I$ のインピーダンス加重平均はどこから来るのか、Wood と frozen の二つの音速曲線が何を別々に見せるのかを追います。

一行のポイント#

平衡5方程式は音速で破綻します。6方程式に一段下がって解き、stiff緩和で再び揃えてやります。インピーダンス平均 $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1+Z_2)$ という一行がすべてを支えます。

Wood音速という落とし穴#

Kapila(2001)の5方程式モデルは圧力平衡を前提とします。1セル内で二相が同じ圧力を持つという仮定です。一見きれいですが、混合音速が奇妙にふるまいます。混合音速 $c_{eq}$ は Wood(1930) の関係

\frac{1}{\rho c_{eq}^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}

に従います。ここで $\alpha_k$ は相 $k$ の体積分率、 $\rho_k$ は密度、 $c_k$ は純粋音速です。問題はこの曲線が単調ではないこと。水(1500 m/s)と空気(340 m/s)を含むセルの平衡音速は、両端どちらより遥かに低い約24 m/sまで沈み込みます。

下のシミュレーションで自分の手で操作してみましょう。スライダーで $\rho_k$ や $c_k$ を変えて、二曲線の差をご覧ください。

ρ₁ = 1000 kg/m³c₁ = 1500 m/sρ₂ = 10.0 kg/m³c₂ = 340 m/scursor α₁ = 0.500log scale on y-axis

Wood (equilibrium) sound speed dips far below either pure-fluid value near intermediate α; frozen speed stays inside the [c₂, c₁] band.

数値拡散領域では仮想的な混合状態が短く現れ、音波はその中で急に遅くなります。波の到達時刻がずれ、衝撃波の追跡が崩れます。さらに $\alpha$ が0や1に近いところで体積分率の正値性を保つのも厄介です。

迂回路 — 圧力をあえて違うままにする#

Saurelは一歩引きます。圧力平衡の前提を外し、二相が別々の圧力 $p_1$ , $p_2$ を持てるようにします。7方程式 Baer–Nunziato モデルから速度緩和の極限だけを先に取った結果が6方程式モデルです。

\partial_t \alpha_1 + u\,\partial_x \alpha_1 = \mu(p_1 - p_2)

\partial_t (\alpha_k \rho_k) + \partial_x(\alpha_k \rho_k u) = 0

\partial_t (\rho u) + \partial_x(\rho u^2 + \alpha_1 p_1 + \alpha_2 p_2) = 0

\partial_t (\alpha_k \rho_k e_k) + \partial_x(\alpha_k \rho_k e_k u) + \alpha_k p_k\,\partial_x u = \pm p_I \mu (p_1 - p_2)

ここで $\mu$ は圧力緩和率(stiff)、 $p_I$ は界面圧力、 $e_k$ は相ごとの内部エネルギーです。このモデルの音速は

c_f^2 = Y_1 c_1^2 + Y_2 c_2^2

と単調になります。 $Y_k = \alpha_k \rho_k / \rho$ は質量分率です。Woodが作った落とし穴が消えます。

界面圧力 $\hat{p}_I$ — インピーダンス加重平均#

内部エネルギー方程式の右辺に現れる $p_I$ は、7方程式モデルの非対称項を速度平衡の極限で整理した結果で、

p_I = \frac{Z_2 p_1 + Z_1 p_2}{Z_1 + Z_2}, \quad Z_k = \rho_k c_k

という綺麗な形に落ちます。 $Z_k$ は音響インピーダンスです。固い相ほど大きな重みを持ちます。直感的には、二相が界面で仕事をやり取りするとき、固い側ほど動かず、柔らかい側ほど動く比率です。この一行がエネルギー保存とエントロピー不等式を同時に保証します。

緩和ステップ — 単一変数の方程式一本#

数値計算では $\mu \to \infty$ の極限を直接時間積分しません。代わりに分割法で双曲部を先に解き、緩和項のみのODEを別途積分します。緩和ステップでは $u$ , $\rho E$ , $\alpha_k \rho_k$ が保存され、各相の比体積 $v_k = 1/\rho_k$ のみが変化します。Stiffened gas EOS

p_k = (\gamma_k - 1)\rho_k e_k - \gamma_k p_{\infty,k}

を用いると、積分は閉形式に落ちます。

v_k(p) = v_k^0\,\frac{p^0 + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}{p + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}

飽和条件 $\sum_k (\alpha\rho)_k v_k(p) = 1$ を解けば平衡圧 $p$ が一発で求まります。二分法一回で終わり。これがアルゴリズムの心臓部です。

再初期化 — 保存エネルギーで再アンカー#

緩和で得た $p$ , $\alpha_k$ は各相EOSとは整合しますが、保存された混合エネルギーとはずれることがあります。非保存の内部エネルギー方程式が衝撃波で誤差を貯め込むためです。Saurelはもう一手間加えます。保存された混合全エネルギーから混合EOS

p(\rho, e, \alpha_1, \alpha_2) = \frac{\rho e - \big(\frac{\alpha_1 \gamma_1 p_{\infty,1}}{\gamma_1 - 1} + \frac{\alpha_2 \gamma_2 p_{\infty,2}}{\gamma_2 - 1}\big)}{\frac{\alpha_1}{\gamma_1 - 1} + \frac{\alpha_2}{\gamma_2 - 1}}

で圧力を計算し直し、その圧力で各相の $e_k$ を再初期化します。出発は非保存ですが、衝撃波付近では本物の保存系のように振る舞います。コストはセルあたりEOS一回分のみ。

コード — 1セルで二圧力が出会う瞬間#

緩和ステップだけNumPyで30行ほど追ってみましょう。Stiffened gasと二分法でサチュレーション式を解くコアループです。

import numpy as np
 
def stiffened_gas_volume(p, p0, v0, gamma, p_inf, p_hat_I):
    num = p0 + gamma * p_inf + (gamma - 1.0) * p_hat_I
    den = p   + gamma * p_inf + (gamma - 1.0) * p_hat_I
    return v0 * num / den
 
def saturation_residual(p, alpha_rho, p0, v0, gamma, p_inf, p_hat_I):
    total = 0.0
    for k in range(2):
        v_k = stiffened_gas_volume(p, p0[k], v0[k], gamma[k], p_inf[k], p_hat_I)
        total += alpha_rho[k] * v_k
    return total - 1.0
 
def stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z, tol=1e-9):
    p_hat_I = (Z[1] * p0[0] + Z[0] * p0[1]) / (Z[0] + Z[1])
    lo, hi = 1.0, 1.0e10
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        r = saturation_residual(mid, alpha_rho, p0, v0, gamma, p_inf, p_hat_I)
        if r > 0:
            lo = mid       # 圧力が低すぎる -> 体積の和が1超
        else:
            hi = mid
        if hi - lo < tol * mid:
            break
    return 0.5 * (lo + hi), p_hat_I
 
# 水(alpha=0.7) + 空気(alpha=0.3), 圧力が異なる
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])
rho0   = np.array([1000.0, 1.0])
alpha0 = np.array([0.7, 0.3])
p0     = np.array([5.0e5, 1.0e5])
Z      = np.array([1500.0 * 1000.0, 340.0 * 1.0])
 
alpha_rho = alpha0 * rho0
v0        = 1.0 / rho0
 
p_eq, p_hat_I = stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z)
print(f"インピーダンス平均 ^pI = {p_hat_I:.3e} Pa")
print(f"緩和後の平衡圧 p     = {p_eq:.3e} Pa")

水のインピーダンス( $Z_1 = 1.5\times 10^6$ )は空気( $Z_2 = 340$ )の約4400倍です。よって $\hat{p}_I$ は実質的に水の圧力とほぼ等しくなります。緩和後の平衡圧は、二相がわずかに圧縮・膨張し合った結果です。

観察 — 二圧力が出会う時間#

1セルだけでは動的感覚がつかめません。1Dグリッド上の簡略化ODEモデルを動かしてみましょう。各セルの $p_1, p_2$ がインピーダンス平均 $\hat{p}_I$ に向かって指数的に収束します — $\mu \cdot \Delta t$ が大きいほど速く。

Z₁ = 1.50 MPa·s/mZ₂ = 3.40 kPa·s/mμ·Δt = 0.050p₁(left) = 5.00 barp₂(left) = 1.00 barp₁(right) = 1.00 barp₂(right) = 5.00 bar

Stiff relaxation drives both p₁ and p₂ toward the impedance-weighted average p̂ᴵ = (Z₂p₁ + Z₁p₂)/(Z₁+Z₂). Increase μ·Δt for faster pull-in.

$Z_1$ を大きくすると $\hat{p}_I$ が相1の圧力側へより強く引き寄せられます。インピーダンスの比率がそのまま重みです。 $\mu \cdot \Delta t$ を1付近にすれば、1ステップで実質平衡に到達します — これがstiff緩和の意味です。Saurel論文ではこの手続きを毎タイムステップ末尾に無条件で適用し、5方程式モデルのstiff極限を回復します。

限界 — 7方程式の影#

6方程式は万能ではありません。速度平衡は依然仮定です。多孔質媒体や粉末爆薬のように二相速度が大きく違う問題では、7方程式 Baer–Nunziato モデルに戻る必要があります。衝撃波が極めて強いとき、非保存の内部エネルギー誤差も無視できず、Petitpasの人工熱交換補正が入ります(論文5節)。さらに $\mu \to \infty$ は界面応答が無限に速いという仮定であり、表面張力や物質移動が入る場合は別モデルが要ります。

それでも6方程式 + 緩和は、現在の圧縮性多相シミュレーションで最も広く使われる土台です。AP-IMEX時間積分(2026-05-03の記事)やMUSCL-THINC-BVD界面再構成(2026-05-02の記事)は、いずれもこの6方程式の上で動いています。

持ち帰る三つの事実#

5方程式の罠は音速にある。 Wood平衡音速の非単調性が界面で音波を誤配します。
6方程式 + 緩和はstiff極限で5方程式を復元する。 非保存系を恐れず、一度解き、一度揃え直すこと。
インピーダンス平均 $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1+Z_2)$ がすべてを支える。 エネルギー保存・エントロピー不等式・単調音速 — 一行で全部です。