Giving Each Phase Its Own Equations — The Baer–Nunziato Seven-Equation Model

In 1986 Baer and Nunziato wanted the moment a flame in a solid propellant turns into a detonation. The powder was a pile of fine grains. Hot gas seeped through the gaps and burned them. At one point, solid and gas sat at different pressures and different speeds. A mixture model with one pressure and one velocity could not hold that instant. So each phase got its own conservation laws.

Setting the relaxation rates μ and λ to zero keeps the two pressure curves from meeting, and raising them pulls both curves into one point.

Watch whether pressure and velocity lock on separately. Raise only one slider and check that the other curve stays put.

Seven Equations Born from Powder Grains#

Diffuse-interface models, which give the boundary between two phases a finite thickness, come in branches. The simplest four-equation model assumes equal pressure, velocity, and temperature. Baer–Nunziato stands at the other end. No equilibrium is assumed. Each phase keeps its own density, velocity, and pressure. One more equation transports the volume fraction $\alpha$, the share of volume one phase occupies.

That generality protects hyperbolicity. Allowing non-equilibrium makes the model unconditionally hyperbolic. A six-equation model that forces a single pressure can lose hyperbolicity and produce a negative square of sound speed. BN avoids that trap. The price is seven equations, plus terms that will not fit the standard conservation form.

Unfolding the Seven Equations#

Write the phases as $k=1,2$. Volume fraction comes first.

The left side carries volume fraction at the interface velocity $u_I$. The right side is pressure relaxation. Larger $\mu$ means faster equilibrium. The other six equations are mass, momentum, and energy in each phase.

$\rho_k$ is the phasic density, $u_k$ the phasic velocity, $p_k$ the phasic pressure. $E_k=\rho_k e_k+\tfrac12\rho_k u_k^2$ is the phasic total-energy density. $k^*$ is the other phase. $\lambda$ is the velocity-relaxation rate. Each phase closes with its own equation of state. Often that is a stiffened gas.

$\gamma_k$ is the ratio of specific heats. $p_{\infty,k}$ is a pressure constant for repulsion between molecules. A liquid that barely compresses, like water, has $p_\infty$ of several GPa.

Non-Conservative Terms and Interface Variables — The Truly Hard Part#

The momentum and energy equations carry $\partial\alpha_k/\partial x$ on the right. That term will not fold into a divergence. It cannot be written as $\partial_t U+\partial_x F(U)=0$. Godunov-type finite volume stumbles in front of such a non-conservative term. When $\alpha$ jumps across a discontinuity, it is unclear which value to multiply. A bad discretization throws fake oscillations across the interface even when pressure and velocity start uniform.

The test is Abgrall's criterion. If pressure and velocity start uniform, they must stay uniform no matter how volume fraction is laid out. Two phases streaming at one pressure and one velocity should only transport the interface, not make a new wave. The record of a 21 percent pressure spike on a flow that should have stayed uniform is what happens when the non-conservative term is multiplied wrong.

The interface values $u_I$ and $p_I$ are modeling choices. The usual Saurel–Abgrall pick is a mass-weighted velocity and a volume-weighted pressure.

That choice makes the model symmetric. Every cell can use the same equations and the same numerics.

Relaxation: The Road to Equilibrium#

The sources $\mu(p_1-p_2)$ and $\lambda(u_{k^*}-u_k)$ are the relaxation operators. Large $\mu$ and $\lambda$ drag the phases to one pressure and one velocity almost at once. Sending the rates to infinity is the simpler equilibrium model.

On the figure above, a rate of zero means the two curves never meet. That is the non-equilibrium the seven equations allow. Larger $\mu$ pulls the pressure curves into one point more steeply. Pressure and velocity relax apart from each other.

The Hierarchy: 7 → 6 → 5 → 4#

Assume relaxation is instantaneous, force equilibrium, and one equation disappears into an algebraic constraint. Tie the velocities and the model has six equations. Tie the pressures as well and it is Kapila's five-equation model. Tie the temperatures and it is four equations. BN sits at the top of that hierarchy.

The five-equation model is the one most used for interface capturing. Its mixture sound speed is non-monotone in volume fraction, the Wood curve, so the volume-fraction source is stiff. BN peels that stiffness off into a relaxation source. Numerically it is the more flexible of the two.

Python — Integrating the Pressure Relaxation Operator#

Strip relaxation down to the bone. Partial masses and internal energies of both phases stay fixed. Only the volume fraction $\alpha$ moves, by $d\alpha/dt=\mu(p_1-p_2)$. Integrate how a stiffened-gas pressure reaches equilibrium as $\alpha$ changes.

import numpy as np
 
# Pressure-relaxation operator for two stiffened-gas phases.
# Each phase's partial mass m_k and internal energy E_k are conserved;
# only the volume fraction alpha evolves by  da/dt = mu (p1 - p2).
 
def phase_pressures(alpha, mass, energy, gamma, p_inf):
    frac = np.array([alpha, 1.0 - alpha])
    rho  = mass / frac          # phasic density = partial mass / volume fraction
    e    = energy / mass        # specific internal energy
    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:                       # stop once the gap vanishes
            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)
 
# Water-like phase (1, stiffened gas) + air-like phase (2)
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])      # Pa
mass   = np.array([480.0, 0.55])     # partial mass  (kg/m^3)
energy = np.array([3.93e8, 1.10e6])  # partial internal energy (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"initial: alpha=0.500  p1={p0[0]/1e5:8.3f} bar  p2={p0[1]/1e5:7.3f} bar")
print(f"equil.:  alpha={alpha_eq:.3f}  p1={peq[0]/1e5:7.3f} bar  p2={peq[1]/1e5:7.3f} bar")
print(f"{len(hist)} steps,  final pressure gap {abs(peq[0]-peq[1]):.3e} Pa")

The run prints this.

initial: alpha=0.500  p1= 324.000 bar  p2=  8.800 bar
equil.:  alpha=0.506  p1=  8.906 bar  p2=  8.906 bar
75 steps,  final pressure gap 9.012e-01 Pa

Pressing a model shows which equilibrium was forced and where the sound-speed curve bends out of monotone. The stiff volume-fraction source shows up the moment the hierarchy drops to five equations.

The printed 324 bar is two huge terms canceling. At $\alpha=0.5$ liquid density is 960. Specific internal energy is $3.93\times10^{8}/480=8.1875\times10^{5}$. Then $(\gamma-1)\rho e$ is $2.6724\times10^{9}$ Pa and $\gamma p_\infty$ is $2.64\times10^{9}$ Pa. The gap is $3.24\times10^{7}$ Pa, 324 bar, the leftover 1.2 percent.

At $\alpha=0.5$ liquid pressure changes by about $-5.34\times10^{9}$ Pa per unit $\alpha$. A shift of 0.001 is about 53 bar. The printed move is 0.500 to 0.506. Gas pressure $4.4\times10^{5}/(1-\alpha)$ is 8.8 bar at 0.5 and about 8.9 bar near 0.506. Equilibrium is the liquid dropping to the gas.

The first step uses $\mu=4\times10^{-8}$, time step $10^{-3}$, and gap $3.152\times10^{7}$ Pa. The product is about $1.26\times10^{-3}$, near a fifth of the printed 0.006. The 75 steps are the shrinking tail. The 1 Pa tolerance matters only at the end.

Volume fraction moved only 0.6 percent, from 0.500 to 0.506, while liquid pressure fell from 324 bar to 8.9 bar. Water modeled as a stiffened gas barely compresses, so a tiny volume change makes a huge pressure change. Equilibrium sits near the gas pressure for the same reason. The stiff liquid can match pressure without expanding much.

The first relaxation step that already spends 0.001 in volume fraction#

The first step already moves volume fraction by $1.26\times10^{-3}$. Liquid pressure drops by tens of bar on that move alone. The interfacial pressure in the non-conservative term multiplies this collapse by the interface gradient. Put relaxation and the flux in the same stage, and the size of a fake wave is set on this first step.

Where does that product go if advection and pressure are split cleanly apart in the flux? The TV split that separates advection from pressure in the Baer–Nunziato flux still leaves open whether the split covers the non-conservative term. When 324 bar collapses in one step, does the separated pressure flux leave the interface quiet?