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

非平衡7方程式モデルの構造、非保存項、そして緩和が生むモデル階層

1986年、Mark Baer と John Nunziato は、固体推進薬の中で炎がどのように爆轟へ転じるかを解こうとしました。火薬は微小な粒の山でした。その隙間を高温のガスが染み込み、粒を燃やします。同じ一点に、固体とガスが異なる圧力、異なる速度で共存していました。単一の圧力・単一の速度を仮定する当時の混合モデルでは、この瞬間を捉えられません。そこで二人は、各相にそれぞれの保存則を丸ごと与えるモデルを立てました。本稿はその7方程式モデルの構造をたどり、なぜ方程式が七つも必要なのか、非保存項がなぜ厄介なのか、そして緩和(relaxation)がどのように単純なモデルを生むのかを見ていきます。最後に圧力緩和演算子を直接積分し、二相が平衡へ向かう過程を再現します。

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

拡散界面(diffuse interface、二相の境界を有限の厚みでにじませて捉える見方)法には枝分かれがあります。最も単純な4方程式は、二相が圧力・速度・温度まですべて等しいとみなします。Baer–Nunziato(BN)は正反対の端に立ちます。いかなる平衡も仮定しません。 各相は自分だけの密度・速度・圧力を持ちます。さらに、一方の相が占める体積の割合、すなわち体積分率(volume fraction) $\alpha$ の輸送方程式が一本加わります。

なぜこの一般性が重要なのでしょうか。非平衡を許すと、モデルは**無条件に双曲型(hyperbolic)**になります。単一圧力を強制した6方程式モデルは、ある条件下で双曲性を失い、音速の二乗が負になるような非物理的な解を生みます。BN はその落とし穴を避けます。ただし代償を払います。方程式が七つに増え、標準的な保存形では書けない項が紛れ込みます。

七つの方程式を広げると#

二相を $k = 1, 2$ と書きましょう。まず体積分率の輸送方程式です。

\frac{\partial \alpha_1}{\partial t} + u_I \frac{\partial \alpha_1}{\partial x} = \mu\,(p_1 - p_2)

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

\frac{\partial (\alpha_k \rho_k)}{\partial t} + \frac{\partial (\alpha_k \rho_k u_k)}{\partial x} = 0

\frac{\partial (\alpha_k \rho_k u_k)}{\partial t} + \frac{\partial (\alpha_k \rho_k u_k^2 + \alpha_k p_k)}{\partial x} = p_I \frac{\partial \alpha_k}{\partial x} + \lambda\,(u_{k^*} - u_k)

\frac{\partial (\alpha_k E_k)}{\partial t} + \frac{\partial \big(\alpha_k (E_k + p_k)\,u_k\big)}{\partial x} = p_I u_I \frac{\partial \alpha_k}{\partial x} + \lambda u_I (u_{k^*} - u_k)

ここで $\rho_k$ は相密度、 $u_k$ は相速度、 $p_k$ は相圧力、 $E_k = \rho_k e_k + \tfrac12 \rho_k u_k^2$ は相の全エネルギー密度、 $k^*$ は相手の相です。 $\lambda$ は速度緩和率です。各相は自分の状態方程式で閉じます。推進薬でもキャビテーションでも、よく使われるのは剛性気体(stiffened gas)です。

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

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

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

上の運動量・エネルギー方程式の右辺に $\partial \alpha_k / \partial x$ があります。この項は発散形にまとまりません。 つまり保存則 $\partial_t U + \partial_x F(U) = 0$ の形には書けないということです。標準的な Godunov 型有限体積法は、こうした非保存(non-conservative)項の前でつまずきます。不連続をまたいで $\alpha$ が跳ぶとき、どの値を掛けるべきかが曖昧だからです。離散化を誤ると、一様な圧力・速度場でも界面をまたいで偽の振動が生じます。

この落とし穴をふるい分ける試金石が Abgrall の判定基準です。圧力と速度が最初に一様なら、体積分率がどう分布しても、その一様性は保たれなければなりません。 二相が同じ圧力・速度で流れるなら、界面はただ輸送されるだけで、新しい波を生んではいけません。非保存項の離散化は、この条件を正確に満たすよう設計されます。

界面変数 $u_I$ 、 $p_I$ もモデル化の産物です。最もよく使われる Saurel–Abgrall の選択は、質量加重速度と体積加重圧力です。

u_I = \frac{\alpha_1 \rho_1 u_1 + \alpha_2 \rho_2 u_2}{\alpha_1 \rho_1 + \alpha_2 \rho_2}, \qquad p_I = \alpha_1 p_1 + \alpha_2 p_2

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

緩和:平衡へ向かう道#

右辺の二つのソース、 $\mu(p_1-p_2)$ と $\lambda(u_{k^*}-u_k)$ が緩和演算子です。 $\mu$ 、 $\lambda$ を大きくすると、二相はほぼ瞬時に同じ圧力・速度へ引き寄せられます。これが BN の妙味です。緩和率を無限大にする極限が、そのまま単純な平衡モデルになります。

下のシミュレーションで実際に操作してみましょう。二相は異なる圧力・速度から出発し、緩和率 $\mu$ ・ $\lambda$ が両者を平衡へ引き寄せます。

압력 완화율 μ: 18속도 완화율 λ: 14초기 압력차 Δp: 1.6초기 속도차 Δu: 0.8

종료 시 압력차 0.000 bar · 속도차 0.000 m/s

μ·λ를 0으로 내리면 두 상은 영영 평형에 도달하지 못한다 — 이것이 7-방정식의 “비평형” 상태다.

緩和率を 0 に下げると、二つの曲線は決して交わりません。これが7方程式モデルが許す非平衡状態です。 $\mu$ を大きくするほど、圧力曲線はより急に一点へ収束します。圧力と速度が互いに独立して緩和する点にも注目してください。

階層:7 → 6 → 5 → 4#

緩和を瞬時(instantaneous)と仮定する、つまり平衡を強制すると、方程式の一つが代数的制約に変わって消えます。速度を平衡に縛れば6方程式、さらに圧力も縛れば Kapila の5方程式、温度まで縛れば4方程式です。BN はこの階層の頂点に立ちます。

下で各モデルを押して、何を平衡と仮定し、何を失うのかを比べてみましょう。

→

Baer–Nunziato (full)

없음 — 압력·속도·온도 모두 독립

미지수: α, (ρ, u, p)₁, (ρ, u, p)₂

장점
가장 일반적. 무조건 쌍곡형(hyperbolic). 상마다 독립 EOS.

대가
비보존 항 + 두 속도장. 계면 변수 PI·uI 모델링이 까다롭다.

평형을 하나씩 더 부과할수록 방정식은 줄지만 — 일반성도 함께 줄어든다.

5方程式(Kapila)が界面捕捉に最も広く使われますが、混合音速が体積分率に対して非単調(Wood 曲線)であるため、体積分率ソースが硬い(stiff)という難点があります。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

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

覚えておきたいこと#

7方程式 = 非平衡。 各相にそれぞれの圧力・速度・エネルギーを与えた代償として、無条件の双曲性を得ます。
非保存項が核心の難点。 $p_I \partial_x \alpha$ は保存形にまとまりません。Abgrall の判定基準(一様な圧力・速度の保存)が離散化の試金石です。
緩和が階層を作る。 $\mu, \lambda \to \infty$ の極限が 6 → 5 → 4 方程式モデルです。BN はその頂点で、最も一般的な描像を与えます。