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

In 1986, Mark Baer and John Nunziato set out to explain how a flame inside a solid propellant turns into a detonation. The explosive was a pile of tiny grains. Hot gas seeps through the gaps and burns the grains. At a single point, solid and gas coexisted at different pressures and different velocities. The mixture models of the day, which assume one pressure and one velocity, could not capture that moment. So the two of them built a model that gives each phase its own full set of conservation laws. This post walks through the structure of that seven-equation model — why seven equations are needed, why the non-conservative terms are a headache, and how relaxation gives birth to simpler models. At the end we integrate the pressure relaxation operator directly and watch the two phases march to equilibrium.

Seven Equations Born from Powder Grains#

Diffuse-interface methods — which smear the boundary between two phases over a finite thickness — come in branches. The simplest four-equation form assumes the two phases share pressure, velocity, and temperature alike. Baer–Nunziato (BN) sits at the opposite end. It assumes no equilibrium at all. Each phase carries its own density, velocity, and pressure. On top of that comes one more transport equation for the volume fraction $\alpha$ — the fraction of volume that one phase occupies.

Why does this generality matter? Allowing non-equilibrium makes the model unconditionally hyperbolic. The six-equation model, which forces a single pressure, loses hyperbolicity under certain conditions and produces unphysical solutions like a negative speed of sound squared. BN avoids that trap. It pays for it, though. The count rises to seven, and terms creep in that cannot be written in standard conservative form.

Unfolding the Seven Equations#

Label the two phases $k = 1, 2$. First, the volume-fraction transport equation.

The left side carries the volume fraction along at the interface velocity $u_I$; the right side is pressure relaxation, with $\mu$ the relaxation rate (larger means faster equilibrium). The remaining six are the mass, momentum, and energy balances of each phase.

Here $\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$ the total energy density of the phase, and $k^*$ the partner phase. $\lambda$ is the velocity relaxation rate. Each phase is closed by its own equation of state. Whether for propellant or cavitation, the stiffened gas law is a common choice.

$\gamma_k$ is the heat-capacity ratio and $p_{\infty,k}$ a pressure constant representing intermolecular repulsion. A nearly incompressible liquid such as water has a large $p_\infty$ of several GPa.

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

The right sides of the momentum and energy equations carry $\partial \alpha_k / \partial x$. This term does not fold into a divergence. That means the system cannot be written as a conservation law $\partial_t U + \partial_x F(U) = 0$. Standard Godunov-type finite-volume schemes stumble on such non-conservative terms, because it is ambiguous which value of $\alpha$ to multiply by as it jumps across a discontinuity. Discretize it wrong and spurious oscillations erupt across the interface even in a uniform pressure and velocity field.

The touchstone that filters out this trap is Abgrall's criterion. If pressure and velocity start uniform, that uniformity must be preserved no matter how the volume fraction is distributed. When two phases flow at the same pressure and velocity, the interface should merely be transported, not spawn new waves. The discretization of the non-conservative terms is designed to satisfy exactly this condition.

The interface variables $u_I$ and $p_I$ are also products of modeling. The most common Saurel–Abgrall choice is a mass-weighted velocity and a volume-weighted pressure.

This choice makes the model symmetric, so the same equations and the same numerical method apply in every cell.

Relaxation: The Road to Equilibrium#

The two right-side sources, $\mu(p_1-p_2)$ and $\lambda(u_{k^*}-u_k)$, are the relaxation operators. Crank $\mu$ and $\lambda$ up and the two phases get pulled to a common pressure and velocity almost instantly. That is the charm of BN: the limit of infinite relaxation rate is exactly a simpler equilibrium model.

Try the simulation below. The two phases start at different pressures and velocities, and the relaxation rates $\mu$ and $\lambda$ drag them toward equilibrium.

Drop the relaxation rates to zero and the two curves never meet. That is the non-equilibrium state the seven-equation model permits. The larger $\mu$, the more steeply the pressure curves collapse to a point. Notice too that pressure and velocity relax independently of each other.

The Hierarchy: 7 → 6 → 5 → 4#

Assume the relaxation is instantaneous — force the equilibrium — and one equation collapses into an algebraic constraint and disappears. Lock velocity to equilibrium and you get six equations; add pressure and you get Kapila's five; add temperature and you get four. BN sits at the top of this hierarchy.

Below, click each model to compare what is assumed in equilibrium and what is lost.

The five-equation Kapila model is the most widely used for interface capturing, but its mixture sound speed is non-monotonic in the volume fraction (the Wood curve), which makes the volume-fraction source stiff. BN is numerically more flexible in that it peels that stiffness off into a relaxation source and handles it separately.

Python — Integrating the Pressure Relaxation Operator#

Let us isolate the heart of relaxation. The partial masses and internal energies of the two phases are conserved, and only the volume fraction $\alpha$ moves, by $d\alpha/dt = \mu(p_1 - p_2)$. We integrate how the stiffened-gas pressures vary with $\alpha$ and reach equilibrium.

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

The striking part is that the volume fraction moved only 0.6%, from 0.500 to 0.506, yet the liquid pressure collapsed from 324 bar to 8.9 bar. Water modeled as a stiffened gas is nearly incompressible, so a tiny change in volume drives an enormous change in pressure. That is also why the relaxed equilibrium settles near the gas pressure: the stiff liquid can match pressures while barely expanding.

Worth Remembering#

• Seven equations = non-equilibrium. Giving each phase its own pressure, velocity, and energy buys you unconditional hyperbolicity.
• The non-conservative terms are the core difficulty. $p_I \partial_x \alpha$ does not fold into conservative form. Abgrall's criterion (preserving uniform pressure and velocity) is the touchstone for the discretization.
• Relaxation builds the hierarchy. The $\mu, \lambda \to \infty$ limit yields the 6 → 5 → 4-equation models. BN sits on top, giving the most general picture.