Skip to content
cfd-lab:~/ja/posts/five-equation-model● online
$_
Younglin's CFD Blog
毎日の記録 · 計算流体力学
NOTE #006DAY SAT 논문리뷰DATE 2026.03.07READ 2 min readWORDS 904#Multiphase#Diffuse-Interface#Python#Implementation

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

5-equation diffuse interface model を Python で直接実装し、1D 衝撃波と気泡の相互作用問題をシミュレーションします。

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

理論だけではなかなか実感が湧かないものです。 この記事では、1D 5-equation diffuse interface model を Python で直接実装し、衝撃波が気体（ガス）の気泡を貫通する問題をシミュレーションします。

支配方程式の整理#

1D 5-equation model の保存変数とフラックスは以下の通りです。

U=[α1ρ1α2ρ2ρuEα1],F=[α1ρ1uα2ρ2uρu2+p(E+p)u0]\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}

5番目の方程式（α1\alpha_1）は非保存型であるため、別に処理します。

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

状態方程式：Stiffened Gas EOS#

両方の流体で stiffened gas EOS を使用します。

pk=(γk−1)ρkek−γkp∞,kp_k = (\gamma_k - 1)\rho_k e_k - \gamma_k p_{\infty,k}

混合則（mixture rule）を用いて全体の圧力を求めます。

p=ρe−∑kαkγkp∞,kγk−1∑kαkγk−1p = \frac{\rho e - \sum_k \alpha_k \frac{\gamma_k p_{\infty,k}}{\gamma_k - 1}}{\sum_k \frac{\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))

HLLC フラックスの計算#

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次精度）#

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

初期条件：衝撃波と気泡#

# 初期条件：水中の空気の気泡に衝撃波が当たる問題
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

注意事項#

1. α\alpha の範囲制限#

α1\alpha_1 が正確に 0 や 1 になると、EOS の計算でゼロ除算（division by zero）が発生します。 常に ϵ<α1<1−ϵ\epsilon < \alpha_1 < 1 - \epsilon (ϵ∼10−8\epsilon \sim 10^{-8}) でクリッピング（clipping）を行ってください。

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

混合領域における音速は、直感に反して非常に低くなることがあります。 これは物理的に正しい現象であり、CFL 条件に影響を与えます。

1ρc2=α1ρ1c12+α2ρ2c22\frac{1}{\rho c^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}

α1=0.5\alpha_1 = 0.5 付近では、cmixc_{mix} は両方の流体の音速よりも 遥かに小さく なります。

3. 高次精度への拡張#

上記のコードは 1次精度です。実際の研究では以下が用いられます。

  • MUSCL 再構築 ＋ 勾配制限関数（slope limiter）で 2次精度
  • WENO 再構築で 5次精度以上
  • ルンゲ・クッタ法 による時間積分で時間精度の向上

結果の解釈#

シミュレーションを実行すると、衝撃波が気泡を圧縮しながら以下の現象が起こります。

  1. 透過衝撃波（transmitted shock） が気泡の後方へ進行
  2. 反射膨張波（reflected rarefaction） が左側へ移動
  3. 気泡が急激に 圧縮 され、内部圧力が上昇
  4. 後方の壁面で ジェット（jet） が形成される（2D/3D の場合）

この現象は、水中爆발や衝撃波結石破砕術（SWL）などにおいて極めて重要な物理です。

次のステップ#

次の記事では、数値解析の基本概念を ゲームのように 学べるインタラク티브なコンテンツを用意しました。数値的な安定性や CFL 条件などを直接体験することができます。

人工拡散を 0 にすると界面が非物理的に sharp になる — THINC 補正が必要な理由が体感できる。

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

役に立ったらシェアしてください。

GitHub でソースを見る ↗誤字・改善提案 ↗

関連記事

界面捕捉法の比較：VOF vs Level Set vs Diffuse Interface圧縮性多相流における界面追跡の主要な3つの手法について、原理、長所・短所、および適用例を比較します。
Mar 6, 2026·3 min read
20℃の水も沸く — キャビテーションの物理減圧で起こる沸騰、キャビテーション数、気泡崩壊の衝撃波。
Apr 30, 2026·5 min read
SPHで流体を解く ― カーネルから人工粘性まで格子なしで粒子に流体方程式を解かせるSPHの原理、コード、落とし穴
Apr 29, 2026·5 min read