Implementing the 5-Equation Model#

Post-shock water sits at pressure $10^{9}$. The water ahead of it sits at $10^{5}$. An air bubble lies between them, and the volume fraction $\alpha_1$ jumps from nearly 0 to nearly 1. A 1D 5-equation diffuse interface model carries that jump with five equations. The two fluids share one velocity and one pressure. Only the volume fraction has its own transport.

Drop the artificial-diffusion slider to 0 and the volume-fraction front breaks across a single cell.

Raise the slider and the front spreads over several cells. The $\alpha$ in those cells goes straight into the denominator of the mixture pressure.

Summary of Governing Equations#

In 1D the conserved variables and the flux are these.

The fifth flux component is 0. $\alpha_1$ is non-conservative, so it is updated on its own.

Five equations is the model after both phases have been reduced to one velocity and one pressure. Putting an equation set back on each phase is the Baer–Nunziato seven-equation model.

Equation of State: Stiffened Gas EOS#

Both fluids use a stiffened gas.

Mixture pressure comes from this rule.

$p_k$ is the pressure of fluid $k$, $e_k$ its internal energy, and $\gamma_k$, $p_{\infty,k}$ the EOS constants. The denominator is $\sum \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))

pressure is the mixture rule rewritten on conserved variables. Inside sound_speed, $c_k^2=\gamma_k(p+p_{\infty,k})/\rho_k$ is the pure-fluid sound speed.

HLLC Flux Calculation#

HLLC covers the four conservation laws. The contact speed $S^\star$ is returned with the flux and later advects $\alpha$.

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)#

Time stepping is forward Euler. The step is $\mathrm{CFL}\,\Delta x/\max(|u|+c)$, with CFL 0.5.

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

The $\alpha_1$ update is not a conservative flux. The sign of the face speed uface picks the upwind side.

Initial Conditions: Shock-Gas Bubble Interaction#

Ahead of position 0.25 the water has already passed through the shock: $\rho=1100$, $u=50$, $p=10^{9}$. Between 0.4 and 0.6 is the air bubble, $\rho=1$, $u=0$, $p=10^{5}$. Everywhere else is ambient water, $\rho=1000$, $u=0$, $p=10^{5}$.

# 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

compute_total_energy is not in this listing. The EOS already says how to recover energy from pressure: $\rho e=(p+\gamma p_\infty)/(\gamma-1)$.

Considerations#

1. Limiting the Range of $\alpha$#

If $\alpha_1$ is exactly 0 or 1, the conserved product $\alpha_1\rho_1$ cannot be divided by $\alpha_1$. The code clips $\alpha_1$ between $10^{-8}$ and $1-10^{-8}$.

2. Mixture Sound Speed: Wood's Formula#

A mixed cell can be slower than either pure fluid. That is what the formula says, and the CFL reads the same speed.

The pure states in the initial condition already fix the number. Bubble air has $\rho=1$, $p=10^{5}$, $\gamma=1.4$, $p_\infty=0$. Then $c=\sqrt{1.4\times 10^{5}}\approx 374\,\mathrm{m/s}$. Ambient water has $\rho=1000$, $p=10^{5}$, $\gamma=4.4$, $p_\infty=6\times 10^{8}$, so $c\approx 1625\,\mathrm{m/s}$. Post-shock water has $\rho=1100$, $p=10^{9}$, so $c\approx 2530\,\mathrm{m/s}$. The group $p+p_\infty$ is about 2.67 times the ambient water value, while density is only $1100/1000=1.10$ times larger. Kinetic energy at $u=50$ is $0.5\times 1100\times 50^{2}=1.375\times 10^{6}$. That is about 0.13% of the internal-energy density $\rho e\approx 1.07\times 10^{9}$ in the same cell. The shock's energy sits in $p_\infty$, not in the 50 m/s.

Mix bubble air and ambient water at $\alpha_1=0.5$ and a common pressure $10^{5}$. In Wood's formula, $1/(\rho c^2)$ is almost entirely the air term. $\rho_1 c_1^2=1.4\times 10^{5}$, while $\rho_2 c_2^2\approx 2.64\times 10^{9}$. Mixture density is $0.5\times 1+0.5\times 1000=500.5$, and $c_\mathrm{mix}\approx 24\,\mathrm{m/s}$. Not 374, and not 1625.

The pressure denominator is sensitive in the same place. Pure water has $\sum\alpha/(\gamma-1)=1/3.4\approx 0.294$. Leak $\alpha_1=0.01$ into that water and the denominator becomes $0.01/0.4+0.99/3.4\approx 0.316$, about 7.5% larger. Pressure is divided by that denominator, so an $\alpha$ error in a smeared interface shows up as a pressure error.

The CFL step is $\Delta t=0.5\,\Delta x/\max(|u|+c)$, and $\Delta x$ is the spacing of linspace(0, 1, 1000), which is $1/999$. The mixture cell does not set it. Behind the shock, $|u|+c\approx 50+2530=2580$, so $\Delta t\approx 1.9\times 10^{-7}$. A cell at 24 m/s would allow a step about 110 times larger. A low mixture speed does not tighten the time step. What it does is leave interface acoustics almost frozen inside that same $\Delta t$.

3. Extension to Higher-Order Accuracy#

The time integrator above is first order, and the spatial update uses the cell value as it stands. Second order wants MUSCL plus a slope limiter. Above that, WENO. Time accuracy goes up with Runge–Kutta. None of that removes the non-conservative $\alpha$ update or Wood's speed.

Interpreting the Results#

Once the shock enters the bubble, a transmitted shock leaves through the downstream side and a reflected rarefaction runs back to the left. The bubble collapses and its internal pressure rises. A jet on the downstream wall needs 2D or 3D. Underwater blasts and shock-wave lithotripsy are built on this picture.

Upwinding α alone lets contact pressure drift first#

The fifth flux is 0, and $\alpha$ is upwinded by the sign of the face speed. The 7.5% swing in the pressure denominator from an $\alpha_1=0.01$ leak into water means one bad interface cell can move the equilibrium pressure. This first-order update does not enforce that equilibrium on its own. The discretization of the non-conservative term is what decides whether pressure stays put.

The test for that equilibrium is a contact where nothing else should happen. Where the pressure leaks on that contact is the question in the Abgrall criterion and the non-conservative term.