Letting Sound and Matter Flow Separately — Acoustic-Convective Splitting (Lagrange-Projection)#

Sound still runs at 340 m/s next to someone walking at 5 km/h. A compressible solver stuffs both speeds into one step. The ratio is the nuisance. In a slow flow the acoustic waves outrun the matter by hundreds. A direct Godunov step is chained to the fast waves. The slow convection you wanted gets smeared.

Drop the Mach slider and both the split of the pressure pulse and the stall of the entropy blob change together.

At 0.02 the amber pulse peels into two waves that race outward. The teal blob barely moves. The three triangles sit at $u\!-\!c$, $u$, and $u\!+\!c$. Lower the Mach number and only the middle marker stops. The outer gap stays $2c$.

ten Eikelder et al. (2017) pull the two apart. Acoustic part, convective part, each with its own scheme, one after the other. The notes below put that Lagrange-Projection split on the single-phase Euler equations and check it on an acoustic pulse in a slow stream.

“

Paper: M.F.P. ten Eikelder, F. Daude, B. Koren, A.S. Tijsseling, "An acoustic-convective splitting-based approach for the Kapila two-phase flow model", Journal of Computational Physics 331 (2017) 188–208. DOI: 10.1016/j.jcp.2016.11.031

The Jacobian Splits in Two#

Write the 1D Euler equations in primitive variables $\mathbf{W}=(\rho,u,p)^T$ and they take a quasi-linear form.

Here $\rho$ is density, $u$ is velocity, $p$ is pressure. The coefficient matrix $\mathbf{B}$ splits cleanly into two pieces.

$\mathbf{A}$ holds every pressure term — the acoustic part. $\mathbf{C}=u\mathbf{I}$ simply carries matter along — the convective part. The split really amounts to peeling the $u\partial_x$ term out of the Lagrangian derivative $\mathrm{D}/\mathrm{D}t = \partial_t + u\partial_x$.

The Additive Split of the Eigenvalues#

Why this split is powerful becomes clear once you look at the eigenvalues. The wave speeds of the full system are $u-c,\ u,\ u+c$. They add up exactly from an acoustic and a convective contribution.

$c=\sqrt{\gamma p/\rho}$ is the speed of sound. The acoustic speeds are $\pm c$, independent of the flow velocity. The convective speed is $u$ across the board. In the low-Mach limit $M=u/c\to 0$, the acoustic band keeps its width $2c$ while the convective speed shrinks to zero. That gap is the stiffness. The triangle spacing in the figure above is this same addition. As $M\to 0$, only $\lambda^c$ vanishes. $\lambda^a=\pm c$ stays.

Acoustic Step — HLLC in Lagrangian Coordinates#

The acoustic part is the expansion and compression that pressure drives. The paper moves it into mass (Lagrangian) coordinates and solves it with an HLLC-type Riemann solver. The starred state at face $j+1/2$ comes down to two formulas.

$a_{j+1/2}=\max(\rho_j c_j,\ \rho_{j+1}c_{j+1})$ is the acoustic impedance (density times sound speed). Those starred velocities and pressures alone drive the update of the Eulerian variables. Each cell's expansion rate compresses into a single coefficient.

$R_j$ tells you how much the acoustic step stretched or squeezed the cell volume. Density follows immediately as $\rho^{n+1-}_j = \rho^n_j / R_j$.

The background in the run is $c_0=1$, $\rho_0=1$, $\gamma=1.4$, so $p_0=\rho_0 c_0^2/\gamma=1/1.4$. The pressure pulse has amplitude $10^{-3}$. The density anomaly has amplitude $5\times 10^{-2}$. Their centers sit at $x=0.5$ and $x=0.25$, so they do not overlap at $t=0$. The density peak therefore starts on the entropy blob, near $1.05$. If $R_j$ is above 1 the cell expands and density drops. Below 1, it compresses. Because $\Delta t/\Delta x$ multiplies the jump in face $u^*$, a CFL of 0.8 still keeps $R_j$ close to 1. With 400 cells on a domain of length 1, $\Delta x=1/400$.

Convective Step — Upwind Projection#

Now the starred velocity $u^*$ from the acoustic step carries the matter. Update the conserved quantities $\varphi\in\{\rho,\rho u,\rho E\}$ with a plain upwind scheme.

$\varphi^{\text{up}}_{j+1/2}$ is $\varphi_j$ when $u^*\!\ge\!0$, else $\varphi_{j+1}$. The $R_j$ from the previous step comes back here, and mass, momentum, and energy are conserved exactly. Walk through the two steps in order and one time step is done.

Positive $u^*$ takes the face value from the left cell. Negative $u^*$ takes it from the right. $R_j\,\varphi^{n+1-}$ undoes the volume change first. The face flux moves the matter after that. Swap the order and the sums of mass, $\rho u$, and $\rho E$ drift apart. Periodic faces do not leak. With the $R_j$ feedback in place, the momentum sum keeps only roundoff. How wide a plain Roe flux spreads at low Mach is measured separately in the three artificial-diffusion scalings of Hope–Collins.

Python — Acoustic Pulse in a Slow-Moving Flow#

Here is the whole pipeline in numpy: periodic boundaries, an ideal gas, three conserved quantities per cell. The initial condition lays a small pressure pulse and a density anomaly on top of a background flow $u_0=Mc_0$.

import numpy as np
 
GAMMA = 1.4  # ideal-gas heat capacity ratio
 
def primitives(rho, mom, Ene):
    """Conserved -> primitive variables (velocity, pressure, sound speed)."""
    u = mom / rho
    e = Ene / rho - 0.5 * u * u            # specific internal energy
    p = (GAMMA - 1.0) * rho * e
    c = np.sqrt(GAMMA * p / rho)
    return u, p, c
 
def acoustic_faces(rho, u, p, c):
    """HLLC acoustic state u*, p* at face j+1/2 (paper eq. 34)."""
    rp, up, pp, cp = (np.roll(a, -1) for a in (rho, u, p, c))
    a = np.maximum(rho * c, rp * cp)       # acoustic impedance a = max(rho*c) (eq. 32)
    ustar = 0.5 * (u + up) + (p - pp) / (2 * a)
    pstar = 0.5 * (p + pp) + 0.5 * a * (u - up)
    return ustar, pstar
 
def lagrange_projection_step(rho, mom, Ene, dx, dt):
    """Acoustic step -> convective step (paper eq. 38, 40)."""
    u, p, c = primitives(rho, mom, Ene)
    uf, pf = acoustic_faces(rho, u, p, c)          # j+1/2
    uf_m, pf_m = np.roll(uf, 1), np.roll(pf, 1)    # j-1/2
    lam = dt / dx
    # 1) acoustic step: expansion/compression driven by pressure
    R = 1.0 + lam * (uf - uf_m)                     # eq. (39)
    rho1 = rho / R
    mom1 = (mom - lam * (pf - pf_m)) / R
    Ene1 = (Ene - lam * (pf * uf - pf_m * uf_m)) / R
    # 2) convective step: carry matter at u* (upwind projection)
    def project(phi1):
        phi_f = np.where(uf >= 0, phi1, np.roll(phi1, -1))
        phi_f_m = np.roll(phi_f, 1)
        return R * phi1 - lam * (uf * phi_f - uf_m * phi_f_m)
    return project(rho1), project(mom1), project(Ene1)
 
def run_acoustic_pulse(mach, n=400, cfl=0.8, tmax=0.25):
    x = (np.arange(n) + 0.5) / n
    dx = 1.0 / n
    c0, rho0 = 1.0, 1.0
    p0 = rho0 * c0**2 / GAMMA
    u0 = mach * c0
    dp = 1e-3 * np.exp(-((x - 0.5) / 0.03)**2)     # acoustic pressure pulse
    ds = 5e-2 * np.exp(-((x - 0.25) / 0.03)**2)    # entropy (density) anomaly -> contact wave
    rho = rho0 + dp / c0**2 + ds
    u = np.full(n, u0)
    p = p0 + dp
    mom = rho * u
    Ene = p / (GAMMA - 1) + 0.5 * rho * u * u
    t, m0 = 0.0, mom.sum()
    while t < tmax:
        _, _, c = primitives(rho, mom, Ene)
        dt = min(cfl * dx / np.max(np.abs(mom / rho) + c), tmax - t)
        rho, mom, Ene = lagrange_projection_step(rho, mom, Ene, dx, dt)
        t += dt
    return rho, mom, m0, mom.sum()
 
for M in (0.02, 0.2):
    rho, mom, m0, m1 = run_acoustic_pulse(M)
    print(f"M={M}: momentum error={abs(m1-m0)/abs(m0):.1e}, rho_max={rho.max():.4f}")
# M=0.02: momentum error=2.2e-16, rho_max=1.0497
# M=0.2:  momentum error=0.0e+00, rho_max=1.0452

Momentum is conserved to machine precision. Pressure stays bounded and density stays positive. Splitting does not break conservation because the convective step feeds the $R_j$ back in.

The printed $\rho_{\max}$ is 1.0497 at $M=0.02$ and 1.0452 at $M=0.2$. The initial crest near 1.05 barely erodes. The pressure pulse adds only $dp/c_0^2=10^{-3}$ to density, so the peak belongs to the contact. Momentum errors of $2.2\times 10^{-16}$ and 0 say the split itself is not what breaks conservation. The step is still limited by $|u|+c$, so raising $M$ from 0.02 to 0.2 barely changes the step count. With $c_0=1$, $|u|+c$ is 1.02 versus 1.2. Convection can speed up by ten and acoustics still own most of the step.

What Actually Diverges at Low Mach#

The split earns its keep at low Mach. A direct scheme's stable step stays chained to $|u|+c$, even when the physics of interest lives on $u$.

Move the Mach slider and the amber band stays locked at width $2c$ while the teal speed $u$ shrinks.

The acoustic width does not depend on $M$. Only $u$ converges to zero. At $c=1$ and $M=0.05$, $u=0.05$ and $u+c=1.05$, so $(u+c)/u=21$. The split can advance convection on a different step from acoustics. Slow physics stops getting dragged by the fast waves. The same split also dodges the extra numerical diffusion a direct Godunov scheme picks up at low Mach.

A Critical Look#

Three snags. The split here is first order in time. Second order means a Strang wrap — acoustic, convective, acoustic — and another pass on the bill. The paper's real stage is the Kapila five-equation two-phase model. The non-conservative term $K\partial_x u$ in the volume-fraction equation, and the positivity proof, are the actual hurdles. This single-phase cut skips them. The impedance $a=\max(\rho c)$ is robust and also diffusive. On a strong shock that diffusion can smear the contact. The paper therefore times accuracy and cost against the direct scheme on their own.

In a code base this split is not exotic. OpenFOAM's pressure-based rhoPimpleFoam and the all-Mach family already treat acoustics and convection implicitly and explicitly from the same root. Lagrange-Projection is the version written in Riemann-solver language.

Drop R_j and conservation dies first#

Leave the acoustic $R_j$ out of the convective projection and the updates for mass, momentum, and energy stop agreeing. The run above reaches $2.2\times 10^{-16}$ only with that feedback in place. A Strang wrap can raise the order. The sequence grows to acoustic–convective–acoustic, and one extra pass comes with it.

Step up to the Kapila five-equation model and $K\partial_x u$ in the volume fraction is still there. The single-phase cut never touches that non-conservative term. Setting $a=\max(\rho c)$ survives shocks and thickens contacts. Another way to lift acoustics out of the time integrator is Peluchon's acoustic-transport IMEX. At $M=0.05$, $(u+c)/u=21$. Whether that factor of 21 sits in the acoustic solver or stays on the convective CFL is the implementation fork.