Compression Schemes -- Riemann solvers, upwinding, shock capturing

Why incompressible schemes fail for compressible flow

The incompressible Navier-Stokes equations assume constant density and use pressure-based SIMPLE or PISO coupling. Compressible flow changes the mathematics fundamentally. Density is no longer constant. The equation of state p = rho R T couples pressure and temperature. The system supports acoustic waves (pressure disturbances propagating at the speed of sound) and discontinuities (shocks, contact surfaces, expansion fans).

The Mach number regime determines the appropriate numerical approach:

At Mach greater than 0.3, density changes become significant and the pressure-based algorithm loses accuracy by solving a pressure-correction equation for incompressible flow. Density-based solvers instead solve the coupled mass, momentum, and energy equations simultaneously using Riemann solvers to compute face fluxes from left and right states.

The Euler equations

The Euler equations govern inviscid flow when viscous terms are neglected:

$$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = 0$$

The vector U contains the conserved variables rho, rho u, and rho E. The flux F is the convective flux. In 1D, F equals rho u, rho u squared plus p, and u times rho E plus p.

The Euler equations form a first-order hyperbolic system with three characteristic waves: left-running acoustic at speed u minus a, contact at speed u, and right-running acoustic at speed u plus a. These waves carry information at finite speed, fundamentally different from incompressible Navier-Stokes where pressure disturbances propagate at infinite speed through an elliptic pressure equation.

The Riemann problem

The Riemann problem is the initial value problem for the Euler equations, starting with piecewise-constant data separated by a single discontinuity. The solution structure is universal for the Euler equations regardless of left and right states: two acoustic waves emanate from the origin, one contact discontinuity travels at velocity u-star, and two constant regions separate three waves.

Riemann solvers

Approximate Riemann solvers replace the exact wave structure with a simplified model that computes fluxes in constant time.

Roe's solver (1981)

Roe's solver linearizes the flux Jacobian across the discontinuity using Roe-averaged variables. The wavespeeds are u-tilde plus or minus a-tilde. The flux combines eigenvalues and eigenvectors of the Roe-averaged Jacobian. Advantages: exact contact-resolution, cheap computation, second-order accurate. Disadvantages: can produce non-physical entropy-violating solutions requiring an entropy fix. OpenFOAM rhoCentralFoam uses Roe as default. Fluent and Star-CCM+ use Roe with entropy fix.

AUSM family (Advection Upstream Splitting Method)

Van Leer's AUSM family avoids Riemann eigenstructure entirely, instead splitting the scalar convective velocity into positive and negative parts. AUSM+ properties: contact-sharper than Roe with less dissipation, robust at all Mach numbers including low-Mach when preconditioned, no eigenvalue computation (faster per-flux evaluation), default in Fluent for density-based formulation, not in base OpenFOAM (available via third-party).

HLLC (Harten-Lax-van Leer with Contact)

HLLC is a compromise between the cheap HLL solver and exact Riemann solver. HLL has no contact wave. HLLC adds the missing contact wave, resolving it exactly while being cheaper and more robust than Roe. It is available in Fluent, Star-CCM+, and OpenFOAM rhoCentralFoam.

Comparison: which solver to choose?

For most compressible flow applications, the choice is between Roe and AUSM. Roe defaults in OpenFOAM rhoCentralFoam. AUSM defaults in Fluent. For shocks with strong detonation or combustion, use AUSM+up or SLAU2. For transonic external aerodynamics, Roe or HLLC suffices. The practical difference on integrated forces is typically 2 to 5% for well-resolved shocks. Near the shock location, wall-shear and heat-transfer differences reach 10 to 20%.

Upwinding and flux limiting

Upwinding in compressible flow means solving the Riemann problem at each face instead of just taking the upstream value. At low Mach numbers (M less than 0.3), acoustic wavespeed greatly exceeds flow speed, making the upwind direction ambiguous since acoustic waves travel in both directions. The pressure-based algorithm is more robust in this regime, treating pressure as elliptic with infinite wave speed. At transonic and supersonic speeds, upwinding becomes unambiguous and the density-based Riemann solver is more accurate.

The Sod shock tube as a test case

The Sod shock tube is the canonical validation case for compressible flow solvers, testing every essential feature: right-running supersonic shock, contact discontinuity traveling at flow velocity, and left-running centered expansion fan. Pre-shock conditions are rho equals 0.125, p equals 0.1, u equals 0. Post-shock conditions are rho equals 1.0, p equals 1.0, u equals 0. An analytical solution exists for all three zones. Interstitial pressure p-star and velocity u-star are found from the Riemann solver solution.

OpenFOAM rhoCentralFoam handles this class of problem. It is density-based with shock-capturing. The scheme selection in fvSchemes:

ddtSchemes { default  Euler; }
divSchemes { default  none;
             div(phi,U)      Gauss rhoCentral;
             div(phi,k)      Gauss vanLeer;
             div(phi,omega)  Gauss vanLeer;
             div(phir,p)     Gauss rhoCentral;
}

Cross-references

• FVM Overview -- upwinding for incompressible convection
• TVD Limiters -- shock-capturing limiters alongside Riemann solvers
• Numerical Dispersion -- shock-smearing from numerical diffusion
• OpenFOAM Shock Fluid -- rhoCentralFoam, Riemann solvers, Sod shock tube
• OpenFOAM Compressible Fluid -- rhoSimpleFoam, buoyant solvers
• Toro, E. F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer.