Skip to content
cfd-lab:~/en/posts/five-equation-model● online
$_
Younglin's CFD Blog
daily journal · computational fluids
NOTE #006DAY SAT 논문리뷰DATE 2026.03.07READ 2 min readWORDS 351#Multiphase#Diffuse-Interface#Python#Implementation

5-Equation Model Implementation Guide: Building a Multiphase Flow Solver in Python

Directly implement a 5-equation diffuse interface model in Python and simulate a 1D shock-bubble interaction problem.

Implementing the 5-Equation Model#

Theory alone is often not enough to fully grasp a concept. In this post, we will directly implement a 1D 5-equation diffuse interface model in Python and simulate a problem where a shock wave passes through a gas bubble.

Summary of Governing Equations#

The conserved variables and fluxes for the 1D 5-equation model are as follows:

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}

The 5th equation (α1\alpha_1) is non-conservative and is handled separately:

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

Equation of State: Stiffened Gas EOS#

Both fluids use the stiffened gas EOS:

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

The mixture pressure is calculated using a 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 Implementation#

Initialization and EOS#

import numpy as np
import matplotlib.pyplot as plt
 
# Domain
N = 1000
x = np.linspace(0, 1, N)
dx = x[1] - x[0]
 
# EOS parameters
# Fluid 1: Air (ideal gas)
gamma1, pinf1 = 1.4, 0.0
# Fluid 2: Water (stiffened gas)
gamma2, pinf2 = 4.4, 6e8
 
def pressure(alpha1, arho1, arho2, E, u):
    """Compute mixture pressure from conserved variables."""
    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):
    """Compute mixture sound speed (Wood's formula)."""
    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's formula
    inv_rhoc2 = alpha1/(rho1*c1sq) + alpha2/(rho2*c2sq)
    return np.sqrt(1.0 / (rho * inv_rhoc2))

HLLC Flux Calculation#

def hllc_flux(UL, UR, alpha1L, alpha1R):
    """HLLC numerical flux for the 4 conservation equations."""
    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)
 
    # Wave speed estimates
    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

Time Integration (1st order)#

def solve(U, alpha1, dx, t_end, cfl=0.5):
    """Forward Euler time integration."""
    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
 
        # Compute fluxes at cell interfaces
        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]
            )
 
        # Update conserved variables
        for k in range(4):
            U[k, 1:-1] -= dt/dx * (flux[k, 2:-1] - flux[k, 1:-2])
 
        # Update alpha1 (non-conservative, upwind)
        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

Initial Conditions: Shock-Gas Bubble Interaction#

# Initial condition: shock hitting an air bubble in water
U = np.zeros((4, N))
alpha1 = np.zeros(N)
 
for i in range(N):
    if x[i] < 0.25:
        # Post-shock water (high pressure)
        rho2 = 1100.0; u_val = 50.0; p_val = 1e9
        alpha1[i] = 1e-8
        arho1 = alpha1[i] * 1.0  # negligible air
        arho2 = (1 - alpha1[i]) * rho2
    elif 0.4 < x[i] < 0.6:
        # Air bubble
        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:
        # Ambient water
        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

Considerations#

1. Limiting the Range of α\alpha#

If α1\alpha_1 becomes exactly 0 or 1, division by zero may occur in the EOS calculations. Always clip α1\alpha_1 to the range ϵ<α1<1−ϵ\epsilon < \alpha_1 < 1 - \epsilon (where ϵ∼10−8\epsilon \sim 10^{-8}).

2. Mixture Sound Speed: Wood's Formula#

The sound speed in a mixture can be surprisingly low, which is physically correct and affects the CFL condition.

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}

Near α1=0.5\alpha_1 = 0.5, cmixc_{mix} can be much smaller than the sound speed of either individual fluid.

3. Extension to Higher-Order Accuracy#

The code provided above is first-order accurate. In actual research:

  • Achieve 2nd-order accuracy with MUSCL reconstruction + slope limiter.
  • Use WENO reconstruction for 5th-order or higher.
  • Improve temporal accuracy with Runge-Kutta time integration.

Interpreting the Results#

Running the simulation shows the shock wave compressing the gas bubble, resulting in:

  1. A transmitted shock proceeding past the bubble.
  2. A reflected rarefaction wave moving to the left.
  3. Rapid compression of the bubble, causing internal pressure to rise.
  4. Formation of a jet at the rear wall (in 2D/3D).

This phenomenon is key in fields such as underwater explosions and shock wave lithotripsy.

Next Steps#

In the next post, I have prepared interactive content where you can learn basic numerical analysis concepts like a game. You will be able to experience numerical stability, the CFL condition, and more directly.

Set artificial diffusion to 0 and watch the interface go unphysically sharp — exactly why we need THINC-like sharpeners.

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

Share if you found it helpful.

View source on GitHub ↗Suggest an edit ↗

Related

Interface Capturing Methods: VOF vs Level Set vs Diffuse InterfaceComparing the principles, pros/cons, and applications of three representative interface capturing techniques in compressible multiphase flows.
Mar 6, 2026·4 min read
Water Boils at 20°C Too — The Physics of CavitationHow depressurization makes water boil, the cavitation number, and the 10,000-bar collapse.
Apr 30, 2026·5 min read
SPH Without Tears — Solving Fluids with Particles, from Kernel to Artificial ViscosityHow Smoothed Particle Hydrodynamics replaces grids with particles — equations, code, traps.
Apr 29, 2026·5 min read