Fewer Equations, and It Stopped Blowing Up — The Four-Equation Model and the Interface Equilibrium Condition

Uniform pressure survives only if Γ itself is upwinded; dividing a conserved ρΓ by ρ is what bends it.

Delete one equation and the code stops crashing. Compressible multiphase flow is usually taught the other way: a seven-equation model, with its own pressure, velocity, and temperature per phase, is the safe one, and fewer equations means more risk. Collis (2025) flips that order. The shortest model, four equations, is the one that holds. The hinge is interface equilibrium.

Turn the Naive toggle on and push right-gas γ to 1.66 or 1.09: the red pressure line at the interface rises, and the green line stays flat.

right gas γ = 1.66CFL = 0.60

Uniform velocity, uniform pressure, only γ jumps across the shaded interface. The naive scheme grows a pressure spike at the interface; the consistent one keeps P flat to machine precision.

The red peak grows with the γ jump. Whether the green line is still on the floor at step 900 is whether the equilibrium condition is alive.

The Ladder of Equilibria — 7, 6, 5, 4 Equations#

A diffuse-interface method needs a rule for the two phases inside a smeared cell. That rule is what counts as equilibrium.

Seven equations (Baer–Nunziato): each phase has its own velocity, pressure, and temperature. Most general, most expensive.
Six equations: one velocity. Pressures stay separate and meet by relaxation.
Five equations: shared velocity and pressure. Temperature stays separate.
Four equations: velocity, pressure, and temperature all shared. Shortest.

A physical interface is nanometers thick. The numerical one is that sheet inflated across a few cells. Inside that band, thermal and mechanical equilibrium is effectively immediate. Near the interface, four equations are the closer physics.

T^{(1)} = T^{(2)}, \quad P^{(1)} = P^{(2)}, \quad \mathbf{u}^{(1)} = \mathbf{u}^{(2)}

Superscripts $(1),(2)$ are the two phases. The four-equation model enforces all three equalities inside the cell.

Why Fewer Equations Are More Robust#

Five equations and up usually add an extra transport equation, for volume fraction or for something like $1/(\gamma-1)$ . At the PDE level those equations repeat information. Once discretized they drift apart. The same quantity, two paths, two values.

Four equations drop the extra. Only mass, momentum, and energy are conserved. Pressure closes through a mixture equation of state. Fewer equations, fewer places to drift. Cheaper, too.

One problem remains. The slack those extra equations used to provide is gone, and a conservative discretization then pulls the pressure off at the interface.

The Moment the Interface Equilibrium Condition Breaks#

The interface equilibrium condition (IEC) is plain. If pressure, temperature, and velocity start uniform, the interface passing through must not break that uniformity.

A scheme that breaks it grows a pressure oscillation that was not there. The oscillation grows in time. At a large density ratio, negative density or negative internal energy shows up and the run stops.

Why a conservative scheme does this. For an ideal-gas mixture, set $\Gamma \equiv 1/(\gamma-1)$ . Total energy is

E = P\,\Gamma + \tfrac{1}{2}\,\rho\,u^2

$E$ is total energy, $P$ pressure, $\rho$ density, $u$ velocity. $\Gamma$ jumps at the interface because the fluids differ.

Push one cell with first-order upwind under uniform $u$ and $P$ . $\nu = u\,\Delta t/\Delta x$ is the CFL number. The conservative energy update cancels the $+P$ piece and leaves

E_i^{n+1} = E_i - \nu\,(E_i - E_{i-1})

Substitute the energy definition and the right-hand side splits.

E_i^{n+1} = P_0\big[\,\Gamma_i - \nu(\Gamma_i - \Gamma_{i-1})\,\big] + \tfrac{1}{2}\,u_0^2\,\rho_i^{n+1}

For the new pressure to stay $P_0$ , $\Gamma$ has to update exactly like that bracket. That is a non-conservative upwind. The conservative path advects $\rho\Gamma$ as a whole and divides by $\rho$ again. The two answers split at the interface. A uniform pressure that still shakes when the non-conservative path is inconsistent is the spot Abgrall's criterion marks. The split is the pressure oscillation.

The red line that grew when γ moved is this split. Switch to IEC-consistent and the green line is still flat at step 900.

Python — The Oscillation and Its Cure#

Same calculation. Uniform $u$ and $P$ , only $\gamma$ jumps, first-order upwind. Only the scheme differs.

import numpy as np
 
def advect_interface(scheme, gamma_R=1.66, nu=0.6, steps=300, N=200):
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    ramp = np.clip((x - 0.28) / 0.06, 0.0, 1.0)   # interface smeared over a few cells
    gamma = 1.4 + (gamma_R - 1.4) * ramp          # air (1.4) -> chosen gas
    rho = 1.0 + (0.2 - 1.0) * ramp                # density jumps too
    u0, P0 = 1.0, 1.0
    Gam = 1.0 / (gamma - 1.0)                      # mixture stiffness 1/(gamma-1)
    mom = rho * u0
    E = P0 * Gam + 0.5 * rho * u0**2               # E = P*Gamma + 0.5*rho*u^2
 
    for _ in range(steps):
        rhoG = rho * Gam                           # rho*Gamma before the update
        r_up, m_up, e_up = np.roll(rho, 1), np.roll(mom, 1), np.roll(E, 1)
        rho = rho - nu * (rho - r_up)              # mass: conservative upwind
        mom = mom - nu * (mom - m_up)              # momentum: conservative
        E   = E   - nu * (E   - e_up)              # energy: conservative
        if scheme == "naive":
            rhoG = rhoG - nu * (rhoG - np.roll(rhoG, 1))  # advect rho*Gamma conservatively
            Gam = rhoG / rho                              # then divide by rho again
        else:  # consistent
            Gam = Gam - nu * (Gam - np.roll(Gam, 1))      # advect Gamma directly (non-conservative)
 
    P = (E - 0.5 * mom**2 / rho) / Gam             # recover pressure
    return float(np.max(np.abs(P - P0)))
 
for s in ("naive", "consistent"):
    print(f"{s:11s}  max|P-P0| = {advect_interface(s):.3e}")

The printout splits by thirteen digits, not two.

naive        max|P-P0| = 3.83e-02
consistent   max|P-P0| = 4.44e-15

Naive pressure is off by about 4% at the interface. Consistent is flat down to floating-point roundoff. One line did it. Recover $\Gamma$ by dividing by rho, or advect $\Gamma$ itself.

Put numbers on $\Gamma = 1/(\gamma-1)$ . Air at $\gamma=1.4$ has $\Gamma=2.5$ . Helium at $\gamma=1.66$ has $\Gamma=1/0.66 \approx 1.515$ . SF6 at $\gamma=1.09$ has $\Gamma=1/0.09 \approx 11.11$ . Air to helium is a jump of about 0.985. Air to SF6 is about 8.61. The SF6 jump is about 8.74 times the helium jump. That factor is why the red line grows when γ goes to 1.09.

One step of the listing's $\nu=0.6$ on a sharp cell. The right cell is $\rho=0.2$ , $\Gamma\approx 1.515$ ; the upwind neighbor is air, $\rho=1$ , $\Gamma=2.5$ . Advect $\rho\Gamma$ conservatively and divide by $\rho$ , and $\Gamma_{\mathrm{naive}}\approx 2.384$ . Advect $\Gamma$ directly and $\Gamma_{\mathrm{cons}}\approx 2.106$ . Invert the energy relation above and $P/P_0 = \Gamma_{\mathrm{cons}}/\Gamma_{\mathrm{naive}} \approx 0.883$ , so one step gives $|P-P_0|\approx 0.117$ . The same sharp cell with SF6 gives $|P-P_0|\approx 0.692$ . The printed $3.83\times 10^{-2}$ is not that sharp cell. The ramp width $0.06$ is 12 cells at $\Delta x=1/200$ , and over 300 steps the interface travels $\nu \cdot 300 \cdot \Delta x = 0.9$ . A smeared interface run that long keeps the maximum under the one-step jump, at $0.0383$ . The sharp-step $0.117$ over that $0.0383$ is about $3.05$ . Spread across 12 cells, the error is still a few percent. The consistent $4.44\times 10^{-15}$ is roundoff left behind once the two paths agree.

What the Paper Actually Did#

The toy above only pushes one ideal-gas $\gamma$ . Collis (2025) is wider.

A NASG equation of state (Noble–Abel stiffened gas, an EOS that reaches liquids) closes water, air, helium, and SF6 together. Mixture pressure is a closed quadratic from Amagat's law.
High-order ENO reconstruction (WENO/TENO) replaces first-order upwind. The reconstruction is tied to the four-equation equilibrium assumption, so the IEC holds to machine precision.
A positivity-preserving limiter extended to four equations covers a strong shock hitting a large density-ratio interface. Density and the square of the sound speed are kept from going negative.
None of this adds a volume-fraction or specific-heat equation, and it does not drop conservation or a scalar filter.

Limits remain. The liquid is a single component and the gas is ideal-gas components. More components means an iterative solver for pressure and temperature equilibrium. The paper is still a preprint.

Recovering Γ by dividing by ρ bends a flat pressure in one step#

On one sharp cell, helium already has $|P-P_0|\approx 0.117$ . SF6 has $0.692$ . The 300-step printout $0.0383$ is that jump shared across a 12-cell ramp. Leave the one line that advects $\rho\Gamma$ and divides by $\rho$ , and a uniform $P_0$ bends at the interface immediately. Trading the extra equations for that non-conservative line is the whole of the four-equation choice.

Five equations share velocity and pressure and still keep one extra equation. If that extra piece is what collapses the interface sound speed, it is open whether deleting it to save the IEC also deletes the disease, or only hides it in the mixture EOS. Whether the five-equation wood speed that collapses when water and air are mixed half and half is still inside the four-equation closure is the next cut.