OpenFOAM Shock-Capturing and Density-Based Solvers

Overview of Density-Based Shock-Capturing Solvers

When compressible flows contain discontinuities — shock waves, expansion fans, contact surfaces — pressure-based solvers become inefficient and often unstable. The root cause is that pressure-based methods decouple pressure from density through the equation of state and solve a Poisson-type pressure correction, which assumes smooth, elliptic pressure propagation. Shock waves violate the smoothness assumption entirely: density, pressure, and temperature change by finite amounts across an infinitesimally thin region (on the order of a few mean free paths, resolved as 1-3 cells in CFD).

Density-based solvers address this by solving the continuity, momentum, energy, and species equations as a strongly coupled system. The conserved variable vector $\mathbf{U}$ and the flux vector $\mathbf{F}$ are governed by the Euler (or Reynolds-Averaged) equations written in conservation form:

$$ \frac{\partial \mathbf{U}}{\partial t} + \nabla\cdot\mathbf{F}(\mathbf{U}) = \mathbf{S}(\mathbf{U}) $$

where for the Euler equations:

$$ \mathbf{U} = \begin{pmatrix} \rho \\ \rho\mathbf{U} \\ \rho E \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} \rho\mathbf{U} \\ \rho\mathbf{U}\mathbf{U} + p\mathbf{I} \\ (\rho E + p)\mathbf{U} \end{pmatrix} $$

The total energy $E = e + \frac{1}{2}\mathbf{U}\cdot\mathbf{U}$ couples the thermal and kinetic energy fields, and the Jacobian $\partial\mathbf{F}/\partial\mathbf{U}$ has eigenvalues $U_n$, $U_n+a$, and $U_n-a$, where $a$ is the speed of sound. These eigenvalues represent the convective and acoustic wave speeds — the fundamental characteristic waves in compressible flow.

**Compressible Flow :: The structure of compressible Euler and Navier-Stokes equations, characteristics, and Riemann problem theory.

**Compression Schemes :: High-resolution shock-capturing schemes, Riemann solvers, and flux limiters for discontinuous solutions.

**Numerical Schemes :: Discretization methods for compressible flow, from upwind schemes to implicit Runge-Kutta.

**Compressible Fluid Module :: Pressure-based compressible solvers (rhoSimpleFoam, rhoPimpleFoam) for subsonic/transonic flows without strong shocks.

rhoCentralFoam: Central-Upwind Scheme

rhoCentralFoam implements a node-centred finite-volume scheme based on the Kurganov-Tadmor (KT) central-upwind scheme (Kurganov & Tadmor, 2000). Unlike traditional upwind solvers (Roe, AUSM, Godunov variants), the KT scheme does not require the exact computation of left and right eigenvalues and eigenvectors. Instead, it uses:

1. An explicit central (Lax-Friedrichs) flux as the base. 2. Upwind/upwind-downwind corrections based on the local wave speeds. 3. A numerical speed of propagation $a^*$ computed as a local estimate:

$$ a^* = \max_i(|u_{i,n}| + a_i) $$

The scheme reconstructs cell-face values from neighboring cell centres using a slope limiter (van Albada, minmod, or van Leer available). The face flux for a conserved variable $\psi$ is:

$$ \mathbf{F}_{face} = \frac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \frac{a^*}{2}(\mathbf{U}_R - \mathbf{U}_L) + \frac{1}{a^* + |\mathbf{f}'|}(\mathbf{f}_R - \mathbf{f}_L) $$

where subscripts $L$ and $R$ denote left and right states after reconstruction, and $\mathbf{f}'$ represents the non-linear wave speed limits.

The advantage over exact Riemann solvers is simplicity and robustness — no eigenvalue decomposition, no Riemann iteration, no special handling of contact discontinuities. The trade-off is slightly more numerical dissipation, which can smear weak shocks by one additional cell width compared to Roe or HLLC.

Riemann Solvers and Roe's Flux

**Roe's flux decomposition** is the classic upwind approach for the Euler equations. Roe (1981) constructed a linearized Riemann solver based on a weighted average of the left and right states. The key idea: the Riemann problem at each cell face is approximated by a constant-coefficient linear system, where the coefficient matrix $\tilde{A}$ uses Roe-averaged quantities:

$$ \tilde{U} = \frac{U_L + U_R}{2}, \quad \tilde{a}^2 = (\gamma - 1)\left(\tilde{E} - \frac{1}{2}\tilde{U}\cdot\tilde{U}\right) $$

The flux difference is then decomposed into wave contributions:

$$ \mathbf{F}_R - \mathbf{F}_L = \sum_{k=1}^{5} |\tilde{\lambda}_k| \tilde{\alpha}_k \tilde{\mathbf{r}}_k $$

where $\tilde{\lambda}_k$ are the eigenvalues, $\tilde{\mathbf{r}}_k$ the right eigenvectors, and $\tilde{\alpha}_k$ the wave strength coefficients. This is exact (for the linearized system) and captures shocks sharply but requires careful handling of:

• **Entropy fix:** Near zero velocity, the eigenvalue $|U_n| = 0$ produces a singularity. An entropy fix adds a small perturbation to prevent unphysical expansion shocks.
• **Contact discontinuity:** Roe's solver can produce "carbuncle" type instabilities for strong shock-cell interactions. HLLC (Harten-Lax-van Leer-Contact) remedies this by adding a contact wave between the left and right shock waves.
• **HLLC flux:** HLLC reconstructs three waves (left fast, contact, right fast), recovering the correct contact discontinuity speed:

$$ u^*_{HLLC} = \frac{\rho_R u_R p_L - \rho_L u_L p_R + \rho_L \rho_R (u_L - u_R)a_{HLLC}}{\rho_R p_L - \rho_L p_R + \rho_L \rho_R (u_L - u_R)} $$

OpenFOAM's rhoCentralFoam does not implement Roe or HLLC directly; these are available in the density-based solver section of the finiteVolume solvers (rhoPimpleFoam with appropriate flux scheme). However, the mathematical foundations are shared.

Sod Shock Tube: Canonical Test Case

The Sod shock tube is the standard verification case for shock-capturing solvers. It consists of a long tube divided into two sections by a diaphragm:

When the diaphragm is removed at $t = 0$, four waves form: a right-propagating shock wave, a right-propagating contact discontinuity (shear), a left-propagating expansion fan, and the stationary reference frame. The exact solution can be computed analytically or numerically with high precision.

**Verification criteria for rhoCentralFoam:**

• Shock position at $t = 0.2$: within 0.5% of exact value.
• Density jump across contact: correctly captured (contact smeared over ~3-5 cells).
• Total energy conservation: within floating-point tolerance.

A well-configured simulation with second-order upwind reconstruction and appropriate time stepping achieves all three criteria. The solution is sensitive to the reconstruction scheme and the numerical dissipation coefficient; default settings in OpenFOAM typically yield acceptable results but may require refinement for high-resolution studies.

Upwinding vs Central Schemes: Understanding the Difference

Central schemes (Kurganov-Tadmor, Lax-Friedrichs) compute face fluxes symmetrically, using averages of left and right states plus a dissipation term proportional to the grid spacing and wave speed. Upwind schemes (Roe, HLLC, AUSM) decompose fluxes along characteristic waves and apply upwind-biased differences:

For most engineering applications involving shock-boundary layer interactions in external aerodynamics, the central-upwind scheme delivers sufficient accuracy. Upwind Riemann solvers are preferred when shock resolution is paramount: turbine cascade analysis, re-entry vehicle aerodynamics, and detonation waves.

Superimposed High-Speed Aerodynamics

Supersonic flow over wedges generates oblique shock waves, whose angle $\beta$ relative to the free-stream is determined by the $\theta$-$\beta$-$M$ relationship:

$$ \tan\theta = 2\cot\beta\,\frac{M^2\sin^2\beta - 1}{M^2(\gamma + \cos 2\beta) + 2} $$

where $\theta$ is the deflection angle and $M$ the upstream Mach number. This analytical relationship provides an excellent benchmark for validating shock-capturing solvers. The rhoCentralFoam solver on a structured O-grid around a wedge geometry typically resolves the shock angle within 1° for second-order reconstruction on a cell count of 100,000-500,000.

Low-Mach Preconditioning: Why Density-Based Solvers Struggle at Ma < 0.3

The fundamental problem with density-based solvers at low Mach numbers is that the eigenvalues of the flux Jacobian become highly unbalanced. The convective eigenvalue $U_n$ is proportional to $Ma$, while the acoustic eigenvalues $U_n \pm a$ are proportional to $1/\sqrt{\gamma}$ (O(1) in non-dimensional form). This creates a stiffness ratio of $O(1/Ma)$, causing:

1. Severe time step restriction for explicit schemes ($\Delta t \sim O(Ma)$). 2. Ill-conditioning of implicit linear systems. 3. Excessive numerical dissipation in the momentum equations relative to acoustic waves, which are not the physics of interest.

Low-Mach preconditioning modifies the time derivative terms to balance the eigenvalues, effectively rescaling the system for low-speed operation. However, this introduces complications: boundary conditions must be preconditioned consistently, and the approach adds algorithmic complexity that negates the simplicity of the density-based formulation. For $Ma < 0.3$, pressure-based solvers remain the clear choice — see Compressible Fluid Module.

Comparison with ANSYS Fluent

Fluent's coupled density-based solver offers more Riemann solver options, implicit time integration for larger time steps, and automatic low-Mach preconditioning. OpenFOAM's rhoCentralFoam trades some of this sophistication for algorithmic simplicity, transparent implementation, and the ability to modify any component (schemes, reconstruction, limiter) without black-box coupling.

References

• Kurganov, A. & Tadmor, E. (2000). New high-resolution central schemes for hyperbolic conservation laws. Journal of Computational Physics, 160(1), 241-282.
• Roe, P.L. (1981). Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 43(2), 357-372.
• Toro, E.F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics (3rd ed.). Springer. Chapters 5-15.
• Sod, G.A. (1978). A survey of several finite-difference methods for systems of nonlinear hyperbolic conservation laws. Journal of Computational Physics, 27(1), 1-31.
• OpenFOAM Foundation. (2023). OpenFOAM v2212/v14 User Guide, Section on rhoCentralFoam.
• Anderson, J.D. (1990). Computational Aerodynamics: An Introduction through Programming. AIAA Education Series. Chapter 9 on shock-capturing methods.