[Paper Review] Half Water, Half Air, and Sound Crawls at 24 m/s — The K∇·u Term in the Murrone–Guillard Five-Equation Model

A one-line formula from a 1930 textbook#

In 1930, A. B. Wood of the British Admiralty Research Laboratory published a textbook on sound. Tucked inside is a one-line formula for the sound speed of bubbly water. It says the compressibility of a mixture is the volume-weighted sum of the phase compressibilities. Plug in a 50/50 air–water foam and you get about 24 m/s. That is far slower than water (1625 m/s) and far slower than air (342 m/s).

This post reads the Murrone–Guillard INRIA report RR-4778 (2003). The report studies a reduced five-equation model (a model obtained from the seven-equation model in the limit of instantaneous relaxation) for compressible two-phase flow. One question drives the post: what does the $K\nabla\cdot u$ term in the volume-fraction equation actually do? We answer it with a 1D Python run. Same tube, same pulse, with and without $K$, and we measure the speed of sound.

Reducing the equations by locking velocity and then relaxing pressure is laid out in the Saurel (2009) pressure-relaxation review. Today the measurement is the $K$ left in the volume-fraction equation after that relaxation is pushed all the way.

The α slider splits tube A's Wood pulse from tube B's fast pulse on the same screen.

Set the α slider to 0.5 and the pulse in tube A crawls while the pulse in tube B laps it dozens of times. Move α to about 1% and the blue curve has already dropped to 1/14 of the water sound speed.

The question the 2003 report reopened#

Wood's formula is old news. The real question is whether a numerical model reproduces it on its own. The seven-equation Baer–Nunziato model carries separate velocities and pressures for both phases. It is physically rich but expensive. If the phases lock to a common velocity and pressure quickly, you can drop variables.

Kapila et al. (2001) derived a five-equation model in this limit. Murrone and Guillard rederived it with an asymptotic expansion (expanding the solution in powers of a small parameter). They then worked out its hyperbolicity (real, finite wave speeds) and wave structure, and tested two numerical schemes. Cutting still further, until only interface equilibrium remains, is the choice taken by the four-equation model. Today the five equations stay, and only the $K$ line is switched on or off.

One figure in the report is the starting point here. It plots the air–water mixture sound speed against the air volume fraction. There are two curves. "5 Eq. Reduc" is the reduced model. "5 Eq. Trans" is a model that only transports the volume fraction. The two curves are far apart.

Five equations left after infinitely fast relaxation#

Send the velocity and pressure relaxation times of the seven-equation model to zero and this remains:

$\alpha_k$ is the volume fraction of phase $k$, $\rho_k$ its density, $u$ and $p$ the shared velocity and pressure, $E$ the mixture total energy. The first four equations are conservative. Only the last one is not, and its right-hand side $K\nabla\cdot u$ is the subject of this post.

Where K∇·u comes from#

In pressure equilibrium both phases see the same pressure change. Away from shocks, phase $k$ is compressed isentropically, so $D\ln\rho_k/Dt = (\rho_k c_k^2)^{-1} Dp/Dt$. Combine this with phase mass conservation, $D(\alpha_k\rho_k)/Dt = -\alpha_k\rho_k\nabla\cdot u$:

Add the two phase equations. Since $\alpha_1+\alpha_2=1$, the left side vanishes and a pressure equation drops out:

$c_W$ is Wood's sound speed. Substitute it back and $K$ is fixed:

$\rho_k c_k^2$ is the bulk modulus (resistance to compression) of phase $k$. Water is about 19,000 times stiffer than air. With air as phase 1, $K\approx\alpha_1$. When the mixture is squeezed, the air absorbs almost all of the volume change.

Wood's speed: compressibility by volume, inertia by density#

Wood's speed is low because two properties come from different phases. The soft air sets the compressibility. The heavy water sets the inertia. At 50% air, the compressibility is about half that of pure air. The density is about 400 times that of air. Since sound speed scales as $\sqrt{1/(\text{density}\times\text{compressibility})}$, it drops to roughly $1/\sqrt{200}$ of the air sound speed.

Set $K=0$ and the story changes. The volume fraction just rides with the flow. The mixture then behaves like a single stiffened gas (an equation of state that approximates stiff liquids like water). Its sound speed follows $\rho\xi c^2 = \sum_k \alpha_k\xi_k\rho_k c_k^2$ with $\xi_k = 1/(\gamma_k-1)$. Water's huge bulk modulus leaks in through an energy weighting. So sound stays fast even at 50% air.

Code: two models in the same tube#

Equations of state: air $\gamma=1.4,\ \pi=0$, water $\gamma=4.4,\ \pi=6\times10^8$ Pa. Fill a 1 m tube with a 50% air mixture. Place a Gaussian pressure pulse with 0.1% amplitude at the center. Run until the right-going pulse should have moved 0.3 m, then measure the speed from the peak position. Conservative variables use a Rusanov flux. The volume fraction uses central differences plus artificial diffusion of the same size.

import math
 
# stiffened gas: phase 1 = air, phase 2 = water
G = (1.4, 4.4)
PI = (0.0, 6.0e8)
P0 = 1.0e5
RHO0 = (1.2, 1000.0)
 
def wood_speed(a1):
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    comp = a1 / (G[0] * (P0 + PI[0])) + (1 - a1) / (G[1] * (P0 + PI[1]))
    return 1.0 / math.sqrt(rho * comp)
 
def isobaric_speed(a1):
    """alpha only advected (K = 0): the mixture is one stiffened gas"""
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    xi = a1 / (G[0] - 1) + (1 - a1) / (G[1] - 1)
    gpi = a1 * G[0] * PI[0] / (G[0] - 1) + (1 - a1) * G[1] * PI[1] / (G[1] - 1)
    gam = 1 + 1 / xi
    return math.sqrt(gam * (P0 + gpi / xi / gam) / rho)
 
def mix_pressure(a1, rhoe):
    a2 = 1 - a1
    num = rhoe - a1 * G[0] * PI[0] / (G[0] - 1) - a2 * G[1] * PI[1] / (G[1] - 1)
    return num / (a1 / (G[0] - 1) + a2 / (G[1] - 1))
 
def kapila_K(a1, p):
    k1, k2 = G[0] * (p + PI[0]), G[1] * (p + PI[1])
    a2 = 1 - a1
    return a1 * a2 * (k2 - k1) / (a1 * k2 + a2 * k1)
 
def pulse_speed(a_bg, use_K, n=400, travel=0.3, amp=1e-3):
    dx = 1.0 / n
    x = [(i + 0.5) * dx for i in range(n)]
    a1 = [a_bg] * n
    p = [P0 * (1 + amp * math.exp(-((xi - 0.5) / 0.03) ** 2)) for xi in x]
    m1 = [a_bg * RHO0[0]] * n
    m2 = [(1 - a_bg) * RHO0[1]] * n
    mom = [0.0] * n
    E = [a_bg * (pi + G[0] * PI[0]) / (G[0] - 1) + (1 - a_bg) * (pi + G[1] * PI[1]) / (G[1] - 1) for pi in p]
    c = wood_speed(a_bg) if use_K else isobaric_speed(a_bg)
    s = 1.05 * c                        # Rusanov speed = the model's own sound speed
    t_end = travel / c
    steps = math.ceil(t_end / (0.4 * dx / s))
    dt = t_end / steps
    g = lambda q, i: q[min(max(i, 0), n - 1)]   # zero-gradient ghost cells
    for _ in range(steps):
        u = [mom[i] / (m1[i] + m2[i]) for i in range(n)]
        pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] * u[i]) for i in range(n)]
        U = (m1, m2, mom, E)
        F = ([m1[i] * u[i] for i in range(n)], [m2[i] * u[i] for i in range(n)],
             [mom[i] * u[i] + pr[i] for i in range(n)], [(E[i] + pr[i]) * u[i] for i in range(n)])
        new = []
        for q, f in zip(U, F):
            fl = [0.5 * (g(f, i - 1) + g(f, i)) - 0.5 * s * (g(q, i) - g(q, i - 1)) for i in range(n + 1)]
            new.append([q[i] - dt / dx * (fl[i + 1] - fl[i]) for i in range(n)])
        na = []
        for i in range(n):
            K = kapila_K(a1[i], pr[i]) if use_K else 0.0
            dudx = (g(u, i + 1) - g(u, i - 1)) / (2 * dx)
            dadx = (g(a1, i + 1) - g(a1, i - 1)) / (2 * dx)
            lap = g(a1, i + 1) - 2 * a1[i] + g(a1, i - 1)
            na.append(a1[i] + dt * (-u[i] * dadx + K * dudx + s * lap / (2 * dx)))
        m1, m2, mom, E = new
        a1 = na
    pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] ** 2 / (m1[i] + m2[i])) for i in range(n)]
    ipk = max((i for i in range(n) if x[i] > 0.6), key=lambda i: pr[i])
    return c, (x[ipk] - 0.5) / t_end, steps
 
print("alpha_air  wood[m/s]  K=0[m/s]")
for a in [0.0, 1e-4, 1e-3, 0.01, 0.1, 0.5, 0.9, 0.99, 1.0]:
    print(f"{a:9.4f}  {wood_speed(a):9.1f}  {isobaric_speed(a):8.1f}")
amin = min((k / 1000 for k in range(1, 1000)), key=wood_speed)
print(f"Wood minimum: {wood_speed(amin):.2f} m/s at alpha_air = {amin:.3f}")
print(f"K at alpha_air=0.5: {kapila_K(0.5, P0):.6f}")
print("\n1D pulse, alpha_air = 0.5, N = 400")
for useK in (True, False):
    c, meas, st = pulse_speed(0.5, useK)
    tag = "with K" if useK else "K = 0 "
    print(f"{tag}  theory {c:7.2f}  measured {meas:7.2f}  ratio {meas / c:.3f}  steps {st}")

alpha_air  wood[m/s]  K=0[m/s]
   0.0000     1624.9    1624.9
   0.0001      956.6    1624.3
   0.0010      364.8    1618.9
   0.0100      118.6    1567.2
   0.1000       39.4    1228.3
   0.5000       23.6     745.3
   0.9000       39.2     581.7
   0.9900      112.4     540.9
   1.0000      341.6     341.6
Wood minimum: 23.65 m/s at alpha_air = 0.501
K at alpha_air=0.5: 0.499947
 
1D pulse, alpha_air = 0.5, N = 400
with K  theory   23.65  measured   23.55  ratio 0.996  steps 316
K = 0   theory  745.30  measured  742.19  ratio 0.996  steps 315

Run the code above right in your browser. The same Rusanov computation advances one step at a time in both tubes.

At the defaults (N = 400, α = 0.5), press Play and the measured column ends at 23.55 and 742.19. Drop N to 100 and both ratios fall together, yet the 31× gap between the models stays. Move α to time pulses at mixture ratios the output table never lists.

What 23.55 and 742.19 tell us#

Both runs share the grid, the flux, and the initial pulse. Only the $K$ line differs. Yet the pulse speeds split by a factor of 31.5: 23.55 m/s versus 742.19 m/s. Each matches its model's theoretical sound speed within 0.4%. That 0.4% is first-order Rusanov diffusion dragging the peak back.

So the scheme is not wrong. The two models solve different physics. The $K=0$ model transports this pulse accurately, just at the wrong speed. No amount of grid refinement brings 742 down to 23. Deleting $K$ and over-compressing water near an interface is the other face of the same term. Today's numbers are that term's footprint inside a uniform mixture.

Divide the two columns at $\alpha = 0.5$ and $745.3/23.6 \approx 31.6$. The pulse measurement is $742.19/23.55 \approx 31.5$, almost the same. Each measurement is 0.996 of theory. Subtract and the gaps are 0.10 m/s ($23.65-23.55$) and 3.11 m/s ($745.30-742.19$). The absolute miss is larger on the fast side. The ratio is not. First-order Rusanov drags the crest back by the same fraction. The step counts are 316 and 315. Each run sets its CFL from its own sound speed, so the slow wave and the fast wave take the same number of steps. What changes is the time, in the same ratio as the sound speeds.

$K$ is in the printout too. At $\alpha = 0.5$, $K = 0.499947$. If the two bulk moduli split without bound, this formula's $K$ becomes 0.5. The printed value is $5.3\times10^{-5}$ below 0.5, in line with water being about 19,000 times stiffer than air. At $\alpha = 0.01$ Wood's speed is 118.6 m/s, about 1/14 of pure water at 1624.9 m/s. At $\alpha = 0.0001$ it is already 956.6 m/s, 668 m/s below 1624.9. The $K=0$ column at that same spot is 1624.3 m/s. Bubbles cut the sound speed only when $K$ is present.

From 1% to 99% air: where the curves part#

The left column is sensitive to a pinch of air. At 0.01% the speed falls from 1625 m/s to 957 m/s, and at 1% it is 119 m/s. The minimum sticks near $\alpha\approx0.5$ because the bulk-modulus contrast is large.

The right column, at 99% air, is faster than pure air. One percent water has leaked its $\pi$ into the gas, and that is also where the report's "5 Eq. Trans" curve leaves the reduced model. A droplet sending a sound wave almost straight back is a separate story, in the acoustic-impedance post on sharp interfaces. Today's two curves ask, before any reflection off one interface, which sound speed a uniform mixture actually has.

What the report checks against experiment#

The report does not stop at sound-speed curves. It plots phase densities and internal energies from an air–water shock tube for both models side by side. It then plots the shock Hugoniot of an epoxy–spinel powder mixture, the relation between shock speed $U_s$ and particle speed $U_p$. The figure overlays the reduced five-equation model, the seven-equation model, and experimental points.

That comparison matters for a reason. $K\nabla\cdot u$ is non-conservative, so the model alone does not uniquely fix the shock jump conditions. The form of numerical viscosity can shift the post-shock state. Hence the need to check, figure by figure, how close the reduced model stays to the seven-equation model and to experiment. For weak waves, as in today's run, the issue disappears. For strong shocks it needs its own validation.

If the pulse is not near 23 m/s, the code dropped K#

Put a small pressure pulse in a uniform mixture at $\alpha=0.5$ and time it. Near 23 m/s, $K$ is alive. Near 745 m/s, the $K=0$ column in the table, the volume fraction is only being carried along. The sound speed inside the CFL limit belongs in the same check. A frozen speed in a code that should run at Wood's speed makes the time step tens of times too conservative, and feeding 23.6 m/s into a $K=0$ model breaks CFL. Where that stiffness goes, acoustics or transport, overlaps the split chosen in Peluchon's acoustic-transport IMEX.

At 0.1% air, Wood's speed is already 365 m/s, under a quarter of water at 1625 m/s. In early cavitation, or in a pump that ingests bubbles, that number sets when the pressure wave arrives. A cell that still prints 1600 m/s has wandered onto the side without $K$. Even at 0.01% the speed has fallen to 957 m/s, so keeping the pure-water value mis-times the arrival.