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

アイダホ国立研究所での核燃料安全シミュレーションが止まりました。液体ナトリウムに挟まれた小さなガス気泡セルで、混合音速が突如 25 m/s まで落ち込んだのです。圧力平衡を仮定した 5 方程式モデルの限界でした。2009 年、Saurel らは一歩引きます。二つの圧力は別々に解き、ステップの最後だけ揃えます。

$\rho_k$ と $c_k$ のスライダーを動かすと、Wood 曲線と frozen 曲線の開きが変わります。

水 1500 m/s、空気 340 m/s のあたりでは、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) の関係

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

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

コードに書いたセルを Wood の式へ入れると、本文の 24 m/s が戻ります。$\alpha_1 = 0.7$、$\rho_1 = 1000$、$c_1 = 1500$ なら $\rho_1 c_1^2 = 2.25\times 10^9$ です。空気側は $\alpha_2 = 0.3$、$\rho_2 = 1$、$c_2 = 340$ なので $\rho_2 c_2^2 = 1.156\times 10^5$ です。右辺は空気が $2.6\times 10^{-6}$ 付近、水は $3\times 10^{-10}$ でほとんど見えません。混合密度は $0.7\times 1000 + 0.3 = 700.3$ です。$1/(\rho c_{eq}^2)$ が空気の項に食われると $c_{eq}$ は約 24 m/s になります。同じセルの frozen 式 $c_f^2 = Y_1 c_1^2 + Y_2 c_2^2$ は逆です。質量のほとんどが水なので $c_f$ は 1500 m/s の脇に付きます。一セルで二つの音速が六十倍ほど開きます。この落ち込みは、水と空気を半々に混ぜた音速がさらに掘ります。

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

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

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

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

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

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

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

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

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

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

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

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

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

で圧力を計算し直し、その圧力で各相の $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 = (Z_2 p_1 + Z_1 p_2)/(Z_1 + Z_2)$ です。大きい重みは小さい側の圧力に付きます。$p_1 = 5\times 10^5$、$p_2 = 1\times 10^5$ を入れると、分子の $Z_1 p_2$ は $1.5\times 10^{11}$、$Z_2 p_1$ は $1.7\times 10^8$ で止まります。だから $\hat{p}_I$ は水の圧力の脇ではなく、空気圧 $1\times 10^5$ のすぐ隣、約 $1.001\times 10^5$ Pa に立ちます。固い相が界面を独占する絵ではありません。緩和後の共通圧力は、そのあと飽和条件の二分法が別に探します。

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

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

$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 時間積分や MUSCL-THINC-BVD 再構成 は、いずれもこの 6 方程式の上で動きます。

そのセルに5方程式を残すか#

Wood 音速が 24 m/s まで落ちたセルを 5 方程式に任せると、音波の到達が遅れ、衝撃の位置が揺れます。そのセルは 6 方程式へ下がり、二つの圧力を別々に持ち、ステップ末尾だけでインピーダンス平均に貼ります。

速度までずれると、この妥協は終わりです。多孔質や粉末爆薬は 7 方程式へ戻ります。式をさらに減らして界面を閉じるかは、4方程式の平衡が別に問います。