各相にそれぞれの方程式を与える — Baer–Nunziato 7方程式モデル

1986年、Baer と Nunziato は、固体推進薬の中の炎がどう爆轟へ移るかを解こうとしました。火薬は細かい粒の山でした。隙間へ熱いガスが染み、粒を焼きます。同じ点に固体とガスが、違う圧力、違う速度でいました。圧力一つ、速度一つを仮定する混合モデルでは、その瞬間が入りません。そこで二つの相に保存則を別々に与えました。

緩和率 μ と λ を 0 に置くと二つの圧力曲線は出会わず、値を上げると一つの点へ引かれます。

圧力曲線と速度曲線が別々に付くかを見ます。片方のスライダーだけ上げても、もう一方の曲線が元の位置にいるかを確認します。

火薬の粒から生まれた七つの方程式

拡散界面(二つの相の境界を有限の厚さで見る見方)には枝があります。最も単純な4方程式は、圧力も速度も温度も等しいと見ます。Baer–Nunziato は反対側です。平衡を仮定しません。各相は自分の密度、速度、圧力を持ちます。一つの相が占める割合、体積分率 $\alpha$ の輸送式が一行追加されます。

この一般性が双曲性を守ります。非平衡を許すとモデルは無条件に双曲型です。圧力を一つに束ねた6方程式は、条件によっては双曲性が壊れ、音速の二乗が負になる解も出します。BN はその罠を避けます。代償があります。式が七つになり、標準の保存形に書けない項が入ります。

七つの方程式を広げると

二つの相を $k=1,2$ と書きます。まず体積分率です。

左辺は、体積分率が界面速度 $u_I$ で運ばれることです。右辺は圧力緩和です。$\mu$ が大きいほど平衡が速いです。残りの六つは各相の質量、運動量、エネルギーです。

$\rho_k$ は相密度、$u_k$ は相速度、$p_k$ は相圧力です。$E_k=\rho_k e_k+\tfrac12\rho_k u_k^2$ は相の全エネルギー密度です。$k^*$ は相手の相です。$\lambda$ は速度の緩和率です。各相は自分の状態方程式で閉じます。よく使うのは剛性気体です。

$\gamma_k$ は比熱比です。$p_{\infty,k}$ は分子間の反発を表す圧力定数です。水のようにほとんど縮まない液体では、$p_\infty$ が数 GPa と大きいです。

非保存項と界面変数 — 本当に難しいところ

運動量とエネルギーの右辺に $\partial\alpha_k/\partial x$ があります。この項は発散形に束ねられません。$\partial_t U+\partial_x F(U)=0$ とは書けません。Godunov 型の有限体積法は、こういう非保存項の前で揺れます。$\alpha$ が不連続で跳ぶとき、どの値を掛けるかが曖昧です。離散化を誤ると、圧力と速度が一様でも界面を横切って偽の振動が飛びます。

試金石は Abgrall の基準です。圧力と速度が最初に一様なら、体積分率がどう敷かれていても、その一様さは保たれなければなりません。同じ圧力、同じ速度で流れるなら、界面は運ばれるだけで、新しい波を作ってはいけません。一様なはずなのに圧力が21%跳んだ記録は、非保存項の掛け方を誤るとその基準がすぐ壊れる例です。

界面の値 $u_I$、$p_I$ もモデルが選んだ値です。よく使う Saurel–Abgrall の選択は、質量加重の速度と体積加重の圧力です。

この選択はモデルを対称にします。すべてのセルで同じ式、同じ数値手法が使えます。

緩和:平衡へ向かう道

右辺の $\mu(p_1-p_2)$ と $\lambda(u_{k^*}-u_k)$ が緩和演算子です。$\mu$ と $\lambda$ を大きくすると、二つの相はほとんど即座に同じ圧力、同じ速度へ引かれます。緩和率を無限大へ送る極限が、より単純な平衡モデルです。

上の画面で緩和率を 0 に下げると、二つの曲線は出会いません。7方程式が許す非平衡です。$\mu$ を大きくするほど圧力曲線が急に一つの点へ集まります。圧力と速度は別々に緩和します。

階層:7 → 6 → 5 → 4#

緩和を瞬間と仮定し、平衡を強制すると、式の一つが代数制約に変わって消えます。速度を束ねると6方程式です。圧力まで束ねると Kapila の5方程式です。温度まで束ねると4方程式です。BN はこの階層の頂点です。

5方程式が界面捕捉に最もよく使われます。混合音速は体積分率に対して非単調な Wood 曲線なので、体積分率のソースが硬いです。BN はその硬さを緩和ソースとして切り離します。数値的にはより柔軟です。

Python — 圧力緩和演算子を積分する#

緩和の骨だけを外します。二つの相の部分質量と内部エネルギーは保存します。体積分率 $\alpha$ だけが $d\alpha/dt=\mu(p_1-p_2)$ で動きます。剛性気体の圧力が $\alpha$ とともにどう平衡へ着くかを積分します。

import numpy as np
 
# 剛性気体(stiffened gas)二相の「圧力緩和」演算子。
# 各相の部分質量 m_k と内部エネルギー E_k は保存され、
# 体積分率 alpha だけが緩和 ODE  da/dt = mu (p1 - p2) で変化する。
 
def phase_pressures(alpha, mass, energy, gamma, p_inf):
    frac = np.array([alpha, 1.0 - alpha])
    rho  = mass / frac          # 相密度 = 部分質量 / 体積分率
    e    = energy / mass        # 比内部エネルギー
    return (gamma - 1.0) * rho * e - gamma * p_inf
 
def relax_pressure(alpha0, mass, energy, gamma, p_inf,
                   mu=4.0e-8, dt=1.0e-3, nmax=200000, tol=1.0):
    alpha, hist = alpha0, []
    for n in range(nmax):
        p = phase_pressures(alpha, mass, energy, gamma, p_inf)
        hist.append((n * dt, alpha, p[0], p[1]))
        gap = p[0] - p[1]
        if abs(gap) < tol:                       # 圧力差が消えたら終了
            break
        alpha += dt * mu * gap                    # da/dt = mu (p1 - p2)
        alpha = min(max(alpha, 1e-4), 1.0 - 1e-4)
    return alpha, np.array(hist)
 
# 水のような相(1, 剛性気体) + 空気のような相(2)
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])      # Pa
mass   = np.array([480.0, 0.55])     # 部分質量  (kg/m^3)
energy = np.array([3.93e8, 1.10e6])  # 部分内部エネルギー (J/m^3)
 
alpha_eq, hist = relax_pressure(0.50, mass, energy, gamma, p_inf)
p0  = phase_pressures(0.50, mass, energy, gamma, p_inf)
peq = phase_pressures(alpha_eq, mass, energy, gamma, p_inf)
print(f"初期: alpha=0.500  p1={p0[0]/1e5:8.3f} bar  p2={p0[1]/1e5:7.3f} bar")
print(f"平衡: alpha={alpha_eq:.3f}  p1={peq[0]/1e5:7.3f} bar  p2={peq[1]/1e5:7.3f} bar")
print(f"{len(hist)} ステップ,  最終圧力差 {abs(peq[0]-peq[1]):.3e} Pa")

実行結果はこうです。

初期: alpha=0.500  p1= 324.000 bar  p2=  8.800 bar
平衡: alpha=0.506  p1=  8.906 bar  p2=  8.906 bar
75 ステップ,  最終圧力差 9.012e-01 Pa

モデルを押すと、何を平衡に束ねたか、音速曲線がどこで非単調に折れるかが変わります。5式へ下ろした瞬間、体積分率ソースが硬くなる場所がその画面です。

出力の324 barは、二つの大きな項の差です。$\alpha=0.5$ で液体密度は480を0.5で割った960です。比内部エネルギーは $3.93\times10^{8}/480=8.1875\times10^{5}$ です。$(\gamma-1)\rho e$ は $2.6724\times10^{9}$ Pa です。$\gamma p_\infty$ は $4.4\times6\times10^{8}=2.64\times10^{9}$ Pa です。引くと $3.24\times10^{7}$ Pa、324 bar です。各項の1.2%だけが残ります。

傾きもその式から出ます。$\alpha=0.5$ で液体圧力の $\alpha$ 微分は約 $-5.34\times10^{9}$ Pa です。体積分率 0.001 は液体圧力で約53 bar です。印字された移動は0.500から0.506、差は0.006です。気体圧力は $4.4\times10^{5}/(1-\alpha)$ なので、0.5で8.8 bar、0.506付近で8.9 bar です。平衡は、気体が立っている場所へ液体が下りる側です。

最初の一歩の大きさは、コードに書いた係数です。$\mu=4\times10^{-8}$、時間刻みは $10^{-3}$、初期差は324 barから8.8 barを引いた $3.152\times10^{7}$ Pa です。体積分率の増分は、その三つを掛けた約 $1.26\times10^{-3}$ です。印字された総移動0.006の約5分の1です。75ステップは、そのあと差が縮み増分が細かく切られた結果です。許容誤差1 Pa は最後になってやっと意味を持ちます。

体積分率は0.500から0.506へ0.6%しか動いていないのに、液体圧力は324 barから8.9 barへ崩れました。剛性気体で置いた水はほとんど縮まないので、わずかな体積変化が巨大な圧力変化を呼びます。緩和平衡がガス圧力の近くに座る理由も同じです。硬い液体はほとんど膨張しなくても圧力を合わせられます。

最初の緩和ステップが体積率を0.001動かす失敗#

最初の一歩が体積分率をすでに $1.26\times10^{-3}$ 動かします。液体圧力はその移動だけで数十 bar 落ちます。非保存項の $p_I$ は、この崩壊を界面勾配に掛けます。緩和とフラックスを同じ段に置くと、偽の波の大きさがこの最初の一歩で決まります。

移流と圧力をフラックスの中で完全に割ると、その積はどこへ行くのか。Baer–Nunziato フラックスを移流と圧力に割った TV 分離は、その分離が非保存項まで覆うかをまだ開けたままです。324 bar が一歩で崩れるとき、分けた圧力フラックスは界面を静かに置くでしょうか。