時間刻みを決めるのはその24 m/sではなく、衝撃波背後の水の約2530 m/sです。

5-Equation Model を直接実装する#

衝撃波の背後の水は圧力 $10^{9}$ 、その手前の水は $10^{5}$ です。あいだに空気の気泡があり、体積分率 $\alpha_1$ はほぼ 0 とほぼ 1 を行き来します。1D の 5-equation diffuse interface model は、この跳びを式五つで運びます。二つの流体は速度を一つ、圧力を一つ分け合い、体積分率だけが別に流れます。

人工拡散のスライダーを 0 まで下げると、体積分率の界面は一セルで切れます。

속도 u 1.00인공 확산 0.020초기 폭 0.20

5-방정식 모델의 핵심: 부피분율 α의 이류. upwind 만으로는 계면이 점차 번지므로 실제 코드는 THINC / interface sharpening을 추가한다.

スライダーを上げると界面は数セルに広がります。そのセルの $\alpha$ が、混合圧力の分母にそのまま入ります。

支配方程式の整理

1D の保存変数とフラックスは次のとおりです。

\mathbf{U} = \begin{bmatrix} \alpha_1 \rho_1 \\ \alpha_2 \rho_2 \\ \rho u \\ E \\ \alpha_1 \end{bmatrix}, \quad \mathbf{F} = \begin{bmatrix} \alpha_1 \rho_1 u \\ \alpha_2 \rho_2 u \\ \rho u^2 + p \\ (E+p)u \\ 0 \end{bmatrix}

五行目のフラックスは 0 です。 $\alpha_1$ は非保存型なので、別に進めます。

\frac{\partial \alpha_1}{\partial t} + u \frac{\partial \alpha_1}{\partial x} = 0

速度一つ、圧力一つへ縮めたのがこの五つの式です。各相に方程式を戻す側は Baer–Nunziato の 7 方程式です。

状態方程式：Stiffened Gas EOS#

両方の流体が stiffened gas です。

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

混合圧力はこの規則で求めます。

p = \frac{\rho e - \sum_k \alpha_k \frac{\gamma_k p_{\infty,k}}{\gamma_k - 1}}{\sum_k \frac{\alpha_k}{\gamma_k - 1}}

$p_k$ は流体 $k$ の圧力、 $e_k$ はその内部エネルギー、 $\gamma_k$ と $p_{\infty,k}$ は EOS の定数です。分母は $\sum \alpha_k/(\gamma_k-1)$ です。

Python による実装#

初期設定と EOS#

import numpy as np
import matplotlib.pyplot as plt
 
# ドメイン設定
N = 1000
x = np.linspace(0, 1, N)
dx = x[1] - x[0]
 
# EOS パラメータ
# 流体 1: 空気 (理想気体)
gamma1, pinf1 = 1.4, 0.0
# 流体 2: 水 (stiffened gas)
gamma2, pinf2 = 4.4, 6e8
 
def pressure(alpha1, arho1, arho2, E, u):
    """保存変数から混合圧力を計算します。"""
    alpha2 = 1.0 - alpha1
    rho = arho1 + arho2
    ke = 0.5 * rho * u**2
    rhoe = E - ke
 
    num = rhoe - alpha1*gamma1*pinf1/(gamma1-1) - alpha2*gamma2*pinf2/(gamma2-1)
    den = alpha1/(gamma1-1) + alpha2/(gamma2-1)
    return num / den
 
def sound_speed(alpha1, arho1, arho2, p):
    """混合音速を計算します (Wood の公式)。"""
    alpha2 = 1.0 - alpha1
    rho = arho1 + arho2
    rho1 = arho1 / np.maximum(alpha1, 1e-12)
    rho2 = arho2 / np.maximum(alpha2, 1e-12)
 
    c1sq = gamma1 * (p + pinf1) / rho1
    c2sq = gamma2 * (p + pinf2) / rho2
 
    # Wood の公式
    inv_rhoc2 = alpha1/(rho1*c1sq) + alpha2/(rho2*c2sq)
    return np.sqrt(1.0 / (rho * inv_rhoc2))

pressure は上の混合則を保存変数で書き直したものです。sound_speed の $c_k^2=\gamma_k(p+p_{\infty,k})/\rho_k$ が、純粋流体の音速です。

HLLC フラックスの計算#

保存式四つだけが HLLC です。 $\alpha$ を運ぶ接触速度 $S^\star$ も一緒に返します。

def hllc_flux(UL, UR, alpha1L, alpha1R):
    """4つの保存方程式に対する HLLC 数値フラックスを計算します。"""
    arho1L, arho2L, momL, EL = UL
    arho1R, arho2R, momR, ER = UR
 
    rhoL = arho1L + arho2L
    rhoR = arho1R + arho2R
    uL = momL / rhoL
    uR = momR / rhoR
    pL = pressure(alpha1L, arho1L, arho2L, EL, uL)
    pR = pressure(alpha1R, arho1R, arho2R, ER, uR)
    cL = sound_speed(alpha1L, arho1L, arho2L, pL)
    cR = sound_speed(alpha1R, arho1R, arho2R, pR)
 
    # 波速の推定
    SL = min(uL - cL, uR - cR)
    SR = max(uL + cL, uR + cR)
    Sstar = (pR - pL + rhoL*uL*(SL-uL) - rhoR*uR*(SR-uR)) \
          / (rhoL*(SL-uL) - rhoR*(SR-uR))
 
    FL = np.array([arho1L*uL, arho2L*uL, momL*uL + pL, (EL+pL)*uL])
    FR = np.array([arho1R*uR, arho2R*uR, momR*uR + pR, (ER+pR)*uR])
 
    if SL >= 0:
        return FL, Sstar
    elif Sstar >= 0:
        coeff = rhoL * (SL - uL) / (SL - Sstar)
        UstarL = coeff * np.array([
            arho1L/rhoL,
            arho2L/rhoL,
            Sstar,
            EL/rhoL + (Sstar - uL)*(Sstar + pL/(rhoL*(SL-uL)))
        ])
        return FL + SL*(UstarL - UL), Sstar
    elif SR > 0:
        coeff = rhoR * (SR - uR) / (SR - Sstar)
        UstarR = coeff * np.array([
            arho1R/rhoR,
            arho2R/rhoR,
            Sstar,
            ER/rhoR + (Sstar - uR)*(Sstar + pR/(rhoR*(SR-uR)))
        ])
        return FR + SR*(UstarR - UR), Sstar
    else:
        return FR, Sstar

時間積分（1次精度）#

時間前進は前進オイラーです。 $\Delta t$ は $\mathrm{CFL}\,\Delta x/\max(|u|+c)$ で、CFL の既定値は 0.5 です。

def solve(U, alpha1, dx, t_end, cfl=0.5):
    """前進オイラー法による時間積分を行います。"""
    t = 0.0
    while t < t_end:
        rho = U[0] + U[1]
        u = U[2] / rho
        p = pressure(alpha1, U[0], U[1], U[3], u)
        c = sound_speed(alpha1, U[0], U[1], p)
 
        dt = cfl * dx / np.max(np.abs(u) + c)
        if t + dt > t_end:
            dt = t_end - t
 
        # セル界面でのフラックスを計算
        flux = np.zeros((4, N+1))
        uface = np.zeros(N+1)
 
        for i in range(N+1):
            iL = max(i-1, 0)
            iR = min(i, N-1)
            UL_local = U[:4, iL]
            UR_local = U[:4, iR]
            flux[:, i], uface[i] = hllc_flux(
                UL_local, UR_local, alpha1[iL], alpha1[iR]
            )
 
        # 保存変数の更新
        for k in range(4):
            U[k, 1:-1] -= dt/dx * (flux[k, 2:-1] - flux[k, 1:-2])
 
        # alpha1 の更新 (非保存型, アップウィンド)
        for i in range(1, N-1):
            if uface[i] >= 0:
                alpha1[i] -= dt/dx * uface[i] * (alpha1[i] - alpha1[i-1])
            else:
                alpha1[i] -= dt/dx * uface[i] * (alpha1[i+1] - alpha1[i])
 
        alpha1 = np.clip(alpha1, 1e-8, 1-1e-8)
        t += dt
 
    return U, alpha1

$\alpha_1$ の更新は保存フラックスではありません。面速度 uface の符号で風上を選びます。

初期条件：衝撃波と気泡

位置 0.25 より手前は、衝撃波を通過した水です。 $\rho=1100$ 、 $u=50$ 、 $p=10^{9}$ 。0.4 と 0.6 のあいだは空気の気泡で、 $\rho=1$ 、 $u=0$ 、 $p=10^{5}$ 。残りは周囲の水、 $\rho=1000$ 、 $u=0$ 、 $p=10^{5}$ です。

# 初期条件：水中の空気の気泡に衝撃波が当たる問題
U = np.zeros((4, N))
alpha1 = np.zeros(N)
 
for i in range(N):
    if x[i] < 0.25:
        # 衝撃波背後の水 (高圧)
        rho2 = 1100.0; u_val = 50.0; p_val = 1e9
        alpha1[i] = 1e-8
        arho1 = alpha1[i] * 1.0  # 無視できる程度の空気
        arho2 = (1 - alpha1[i]) * rho2
    elif 0.4 < x[i] < 0.6:
        # 空気の気泡
        rho1 = 1.0; u_val = 0.0; p_val = 1e5
        alpha1[i] = 1 - 1e-8
        arho1 = alpha1[i] * rho1
        arho2 = (1 - alpha1[i]) * 1000.0
    else:
        # 周囲の水
        rho2 = 1000.0; u_val = 0.0; p_val = 1e5
        alpha1[i] = 1e-8
        arho1 = alpha1[i] * 1.0
        arho2 = (1 - alpha1[i]) * rho2
 
    rho = arho1 + arho2
    E = compute_total_energy(alpha1[i], arho1, arho2, p_val, u_val)
 
    U[0, i] = arho1
    U[1, i] = arho2
    U[2, i] = rho * u_val
    U[3, i] = E

compute_total_energy はこの断片にありません。圧力からエネルギーへ戻す式は EOS 側にあります。 $\rho e=(p+\gamma p_\infty)/(\gamma-1)$ です。

注意事項

1. $\alpha$ の範囲制限#

$\alpha_1$ がちょうど 0 または 1 だと、保存量 $\alpha_1\rho_1$ を $\alpha_1$ で割れません。コードは $\alpha_1$ を $10^{-8}$ と $1-10^{-8}$ のあいだに切ります。

2. 混合音速：Wood の公式#

混合セルの音速は、どちらの純粋流体よりも低くなり得ます。直感とは逆ですが、式はそう言っています。CFL もこの音速を見ます。

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

初期条件の純粋状態をこの式に入れると、数が出ます。気泡の空気は $\rho=1$ 、 $p=10^{5}$ 、 $\gamma=1.4$ 、 $p_\infty=0$ です。 $c^2=\gamma(p+p_\infty)/\rho$ なら $c=\sqrt{1.4\times 10^{5}}\approx 374\,\mathrm{m/s}$ です。周囲の水は $\rho=1000$ 、 $p=10^{5}$ 、 $\gamma=4.4$ 、 $p_\infty=6\times 10^{8}$ なので $c\approx 1625\,\mathrm{m/s}$ です。衝撃波背後の水は $\rho=1100$ 、 $p=10^{9}$ で $c\approx 2530\,\mathrm{m/s}$ です。 $p+p_\infty$ は周囲の水の約 2.67 倍なのに、密度は $1100/1000=1.10$ 倍だけです。 $u=50$ の運動エネルギーは $0.5\times 1100\times 50^{2}=1.375\times 10^{6}$ です。同じセルの内部エネルギー密度 $\rho e\approx 1.07\times 10^{9}$ の約 0.13% です。この衝撃のエネルギーは速度より $p_\infty$ 側にあります。

$\alpha_1=0.5$ で気泡の空気と周囲の水を、圧力 $10^{5}$ のまま混ぜると、Wood の式の $1/(\rho c^2)$ はほとんど空気の項です。 $\rho_1 c_1^2=1.4\times 10^{5}$ 、 $\rho_2 c_2^2\approx 2.64\times 10^{9}$ です。混合密度は $0.5\times 1+0.5\times 1000=500.5$ で、 $c_\mathrm{mix}\approx 24\,\mathrm{m/s}$ です。374 でも 1625 でもありません。

圧力の分母も同じ場所で敏感です。純粋な水の $\sum \alpha/(\gamma-1)$ は $1/3.4\approx 0.294$ です。水に $\alpha_1=0.01$ だけ漏れても、分母は $0.01/0.4+0.99/3.4\approx 0.316$ となり、約 7.5% 増えます。圧力はこの分母で割るので、広がった界面の $\alpha$ の誤差は圧力の誤差になります。

CFL の $\Delta t=0.5\,\Delta x/\max(|u|+c)$ で、 $\Delta x$ は linspace(0, 1, 1000) の間隔 $1/999$ です。制限するのは混合セルではありません。衝撃波の背後は $|u|+c\approx 50+2530=2580$ なので $\Delta t\approx 1.9\times 10^{-7}$ です。音速 24 m/s のセルだけを見ると、間隔は約 110 倍大きくなります。混合音速が低くても、時間刻みはそのセルが締めません。同じ $\Delta t$ のあいだ、界面の音響がほとんど止まって見える方が影響です。

3. 高次精度への拡張#

上の時間積分は 1 次です。空間もセル値をそのまま使います。2 次が要るなら MUSCL と slope limiter、その上は WENO です。時間は Runge–Kutta で上げます。次数を上げても、 $\alpha$ の非保存更新と Wood の音速は残ります。

結果の解釈

衝撃が気泡に入ると、下流側へ透過衝撃波（transmitted shock）が出て、左へ反射膨張波（reflected rarefaction）が戻ります。気泡は急に押しつぶされ、内部の圧力が上がります。ジェット（jet）は 2D か 3D で、下流側の壁にできます。水中爆発と衝撃波結石破砕は、この絵を使います。

風上の α だけでは接触面の圧力が先に揺れます

五番目のフラックスは 0 で、 $\alpha$ は面速度の符号で風上になります。水に $\alpha_1=0.01$ が漏れただけで圧力の分母が 7.5% 動くという計算は、界面一セルの誤差が平衡圧力に触れるという意味です。この 1 次の更新は、その平衡を別に強制しません。非保存項の離散化が圧力を保つかどうかが、先に分かれます。

平衡を測る問題は、何も起きてはいけない接触です。そこで圧力がどこへ漏れるかは、Abgrall 基準と非保存項が問います。

5-Equation Model 実装ガイド：Python で作る多相流ソルバー

5-Equation Model を直接実装する#

支配方程式の整理

状態方程式：Stiffened Gas EOS#

Python による実装#

初期設定と EOS#

HLLC フラックスの計算#

時間積分（1次精度）#

初期条件：衝撃波と気泡

注意事項

1. $\alpha$ の範囲制限#

2. 混合音速：Wood の公式#

3. 高次精度への拡張#

結果の解釈

風上の α だけでは接触面の圧力が先に揺れます

関連記事

5-Equation Model 実装ガイド：Python で作る多相流ソルバー

5-Equation Model を直接実装する#

支配方程式の整理

状態方程式：Stiffened Gas EOS#

Python による実装#

初期設定と EOS#

HLLC フラックスの計算#

時間積分（1次精度）#

初期条件：衝撃波と気泡

注意事項

1. α\alphaα の範囲制限#

2. 混合音速：Wood の公式#

3. 高次精度への拡張#

結果の解釈

風上の α だけでは接触面の圧力が先に揺れます

関連記事

1. $\alpha$ の範囲制限#