Table of Contents
- When Does Compressibility Matter? Ma > 0.3 and Prandtl-Glauert
- Governing Equations: Euler and Full Compressible N-S
- Mach Number Regimes: Subsonic to Hypersonic
- Shock Waves: Rankine-Hugoniot and Normal/Oblique Shocks
- Nozzle Flows: Choked Flow and Supersonic Expansion
- Numerical Challenges: Shock Capturing and Upwind Schemes
- OpenFOAM and Fluent: Density-Based Solvers
- Engineering Applications
When Does Compressibility Matter? Ma > 0.3 and the Prandtl-Glauert Rule
A flow must be treated as compressible whenever density variations are significant — typically when the Mach number exceeds 0.3. At this threshold, the maximum density change from the isentropic relation is about 4.5%, which is significant enough to affect aerodynamic forces, pressure distributions, and thermal states in many engineering applications.
The Prandtl-Glauert correction gives a quick estimate of compressibility effects on aerodynamic coefficients in the subsonic regime. For a thin airfoil at low angle of attack, the lift coefficient in compressible flow relates to the incompressible value as:
This factor diverges as Ma → 1.0 (the "Prandtl-Glauert singularity"), which is unphysical — the formula breaks down well before Ma = 1 because it is a linear theory valid only for small perturbations and purely subsonic flow. The more accurate Kármán-Tsien and Laitone corrections extend the range slightly, but for Ma > 0.7 on practical geometries, full compressible CFD is required.
So at Ma = 0.5 — far below supersonic — I already need a compressible solver? That really surprised me.
Yes, and here's why it matters more than the 0.3 threshold alone suggests. The free-stream Mach number might be 0.5, but locally on the suction side of an airfoil or inside a converging duct, the flow can accelerate to Ma = 0.8 or even higher — and that's where compressibility really kicks in. The pressure coefficient \(C_p\) is much more sensitive to density changes at Ma = 0.8 than at Ma = 0.5. Commercial transports cruise at Ma = 0.82–0.85, and even though the free stream is subsonic, there are supersonic patches on the upper wing surface that terminate in shocks. These shocks control where boundary-layer separation occurs, which directly affects drag. An incompressible solver would simply miss all of this — no shocks, wrong pressure distribution, wrong drag. The practical rule: if any part of your flow might reach Ma > 0.5 locally, use a compressible solver to be safe.
What about a high-pressure gas pipeline — the velocity is slow, but pressure is very high. Does that need a compressible model?
Interesting case. Natural gas pipelines typically flow at 5–15 m/s; at the speed of sound in natural gas (~400–450 m/s), that's Ma ≈ 0.01–0.04 — dynamically incompressible. However, over long pipelines, the pressure drops from 70 bar at the compressor to 35 bar at delivery. The density at those two points differs by nearly a factor of 2 because of the ideal gas law \(\rho = p/RT\). Even though each cross-section looks locally incompressible, the bulk density varies dramatically along the pipe. This is called low-Mach-number compressible flow — compressibility driven by thermodynamic state changes rather than dynamic pressure effects. You need a compressible model to get correct density, velocity, and mass flow rate distributions. It requires different numerical treatment than shock-dominated high-speed flows, and dedicated pipeline simulation codes handle it efficiently with implicit pressure-correction methods.
Governing Equations: Euler and Full Compressible Navier-Stokes
For compressible flows, density \(\rho\) is a primary variable, and the equation of state closes the system. The governing equations written in conservative form are:
The Euler Equations (Inviscid Compressible Flow)
The Euler equations neglect viscosity and are appropriate for high-Reynolds-number flows away from boundary layers:
Here \(E = e + |\mathbf{u}|^2/2\) is total specific energy, \(e\) is specific internal energy. For a calorically perfect ideal gas: \(p = \rho R T\) and \(p = (\gamma-1)\rho e\) where \(\gamma = c_p/c_v\).
The Full Compressible Navier-Stokes Equations
Adding viscous stresses and heat conduction gives the full system:
Where \(\boldsymbol{\tau} = \mu\left[(\nabla\mathbf{u}+\nabla\mathbf{u}^T) - \frac{2}{3}(\nabla\cdot\mathbf{u})\mathbf{I}\right]\) is the viscous stress tensor and \(\mathbf{q} = -k\nabla T\) is heat flux. Unlike incompressible flow, all five equations (continuity, 3 momentum, energy) are tightly coupled through density, pressure, and temperature simultaneously.
Why do we write the equations in "conservative form"? I see this phrase everywhere in CFD papers but textbooks don't always explain why it matters.
Conservative form means the equation is written as the time derivative of a conserved quantity (mass density, momentum density, energy density) plus the divergence of its flux. When you discretize this on a finite volume mesh, the flux leaving one cell face enters the adjacent cell exactly — mass, momentum, and energy are conserved to machine precision at the discrete level. This matters crucially near shocks. A shock must satisfy the Rankine-Hugoniot jump conditions, which are nothing more than the statement that mass, momentum, and energy are conserved across the discontinuity. If you use a non-conservative form — for instance, writing the momentum equation in terms of velocity directly — conservation is violated at the discrete level, and your numerical shock may travel at the wrong speed, have the wrong strength, or smear into a non-physical gradient. The Lax-Wendroff theorem guarantees that a conservative consistent scheme converges to the correct weak solution, which is what you need when discontinuities are present. That's why every serious compressible flow code writes the equations in conservative form.
Mach Number Regimes: Subsonic to Hypersonic
Compressible flows are classified into four regimes, each with distinct physics and modeling requirements:
Pressure disturbances propagate upstream and influence the entire domain. No wave drag. Aerodynamic characteristics change gradually with Mach number. Examples: commercial aircraft approach/climb, fans, compressor inlet stages. Prandtl-Glauert correction applicable for 0.3 < Ma < 0.7.
Mixed subsonic/supersonic regions coexist. Normal shocks terminate supersonic patches on airfoil suction surfaces. Wave drag rises sharply near Ma = 1. Computationally the most challenging regime. Most commercial transports fly here at cruise.
Entirely supersonic. Oblique shocks and Prandtl-Meyer expansion fans dominate. Bow shocks form ahead of blunt bodies. Wave drag is the main drag component. Examples: supersonic fighters, cruise missiles, rocket nozzle exit planes.
Shock layer temperatures cause vibrational excitation, dissociation, and ionization of air — ideal gas law fails. Radiative heating is significant above Ma ≈ 12. Requires two-temperature thermodynamic models. Examples: atmospheric reentry vehicles, scramjet intakes, ICBM warheads.
Why is transonic flow especially hard to compute, compared to fully supersonic? I'd expect supersonic to be harder because of the shocks.
Counterintuitive, but transonic is harder for several compounding reasons. First, the mixed-type problem: the governing equations are elliptic in subsonic regions (information propagates in all directions) and hyperbolic in supersonic patches (information propagates only downstream along characteristics). No single numerical scheme is naturally suited to both types simultaneously. Second, shock-induced boundary-layer separation: the normal shock at the foot of a supersonic bubble on an airfoil interacts with the boundary layer, and whether that interaction causes separation depends sensitively on Reynolds number, turbulence model, and local mesh resolution. RANS models notoriously disagree on transonic shock/boundary-layer interaction, and even well-validated models can be off by 20–30% on separation extent. Third, transonic buffet: near the design Mach number, small changes in angle of attack cause the shock to oscillate dynamically — a self-sustaining oscillation that shakes the aircraft structure. Predicting buffet onset requires unsteady URANS or DES simulation, not steady RANS. Purely supersonic flows, by contrast, have well-defined shock positions, no upstream influence, and are computationally more predictable despite the sharp features.
Shock Waves: Rankine-Hugoniot Conditions
A shock wave is a thin region (a few molecular mean-free-paths thick, effectively a discontinuity at CFD scales) where flow properties change abruptly. Shocks form whenever a supersonic flow must decelerate or turn, and the flow cannot communicate this requirement upstream through pressure waves.
Normal Shock: Rankine-Hugoniot Jump Conditions
Conservation of mass, momentum, and energy across a stationary normal shock gives exact algebraic relations between upstream (subscript 1) and downstream (subscript 2) states:
The post-shock Mach number is always less than 1 — a normal shock always converts supersonic flow to subsonic. Total pressure \(p_0\) decreases across the shock due to entropy generation. A normal shock at Ma = 2 causes about 28% total pressure loss; at Ma = 5, the loss reaches 96%. This is why efficient supersonic inlets use a series of oblique shocks to decelerate flow progressively with minimal loss.
Oblique Shocks and the θ-β-Ma Relation
When supersonic flow encounters a wedge or compression corner, it deflects through an oblique shock at angle \(\beta\) to the upstream flow, turning through deflection angle \(\theta\). The governing relation is:
For given Ma₁ and \(\theta\), two solutions exist: the weak shock (smaller \(\beta\), post-shock flow remains supersonic — the physically preferred solution) and the strong shock (larger \(\beta\), post-shock subsonic). If \(\theta\) exceeds the maximum deflection angle for the given Mach number, no attached oblique shock is possible and a detached bow shock forms.
Shock waves convert kinetic energy into heat — does that mean a supersonic vehicle heats up from the shock itself, or from friction in the boundary layer?
Both contribute, but in different ways. The shock raises the static temperature of the gas dramatically — at orbital reentry speed (Ma ≈ 25), the stagnation temperature behind the bow shock reaches about 8,000 K, far beyond what any metal can survive. But the heat that reaches the vehicle surface is transferred by convection and radiation through the hot gas layer — the shock itself stands off from the surface at some distance. The rate of heat transfer to the surface is determined by the temperature gradient in the boundary layer and the thermal conductivity and radiation intensity of the shock layer gas. On a blunt reentry vehicle, the bow shock stands well off the nose, and the thick, hot shock layer acts as a natural insulator (this is why blunt shapes are preferred over sharp noses for reentry — counterintuitive but true). On slender bodies at lower speeds like a supersonic fighter at Ma = 2, aerodynamic heating is modest (skin temperatures around 100–200°C at the leading edge) and comes primarily from viscous dissipation in the turbulent boundary layer, not shock heating.
Nozzle Flows: Choked Flow and Supersonic Expansion
The de Laval nozzle (converging-diverging nozzle) demonstrates how flow accelerates from subsonic to supersonic through a throat. The area-Mach relationship for isentropic flow is:
Where \(A^*\) is the throat area at Ma = 1. This equation has two solutions for each area ratio — one subsonic and one supersonic — and which applies depends on the downstream back pressure.
Choked flow occurs when the throat reaches Ma = 1 and mass flow rate reaches its maximum possible value. Reducing back pressure below the critical value cannot increase mass flow — information cannot travel upstream past a sonic throat. The choked mass flow rate is:
I've seen choked flow mentioned for safety relief valves in pressure vessels. Is that the same thing? Why does it matter for safety design?
Exactly the same phenomenon — and it's critical for pressure relief valve sizing. When gas escapes through a relief valve orifice and the upstream-to-downstream pressure ratio exceeds the critical value (about 1.9 for air with γ = 1.4), the flow at the throat goes sonic and locks at its maximum mass flow rate. At that point, the discharge rate depends only on upstream stagnation conditions and orifice area — not on the downstream (atmospheric) pressure. This is useful for design because the maximum relief capacity is well-defined and predictable. The danger is sizing the valve using incompressible flow formulas when the pressure ratio clearly exceeds the critical value — you'll calculate a discharge rate up to 30% higher than the actual choked value, meaning the valve is undersized. In a high-pressure gas vessel upset scenario, an undersized relief valve can't depressurize fast enough, leading to failure of the vessel. This is why process codes like API 520/521 explicitly require choked flow calculations when the pressure ratio exceeds the critical value.
In a rocket engine nozzle, why does the flow have to go through the throat at exactly Mach 1? Can't it just stay subsonic all the way through a converging nozzle?
The area-Mach equation tells you that for any isentropic flow to pass through a geometric throat (minimum area point), it must pass through Ma = 1 at that throat — or it must remain subsonic throughout. If you have purely subsonic flow in a converging-diverging nozzle with the right back pressure, the flow accelerates to a maximum velocity at the throat, then decelerates in the diverging section — no Ma = 1, no supersonic region. This is what happens at low pressure ratios across the nozzle. For a rocket, the combustion pressure is so much higher than ambient that the pressure ratio easily exceeds the critical value, driving the throat to Ma = 1, and the diverging section then further accelerates the flow supersonically. The nozzle shape in the diverging section is designed (using the method of characteristics) to produce the desired exit Mach number with minimum entropy generation — a smooth fan of Prandtl-Meyer expansion waves rather than shock diamonds.
Numerical Challenges: Shock Capturing and Upwind Schemes
The central numerical challenge in compressible CFD is capturing shocks accurately without spurious oscillations while maintaining accuracy in smooth flow regions.
Upwind Schemes
In supersonic flow, information travels only downstream along characteristics. Upwind schemes respect this by biasing the finite difference stencil in the direction of information propagation. First-order upwinding is highly dissipative — shocks are smeared over several cells. Higher-order extensions (MUSCL — Monotone Upstream-centred Scheme for Conservation Laws) achieve second or third order accuracy in smooth regions while applying limiters near shocks.
Roe's Approximate Riemann Solver
The flux at each cell face is computed by solving an approximate Riemann problem — the wave propagation problem across the left/right state discontinuity. Roe's scheme (1981) linearizes the Euler equations around a density-weighted "Roe average" state, giving a flux that identifies all wave speeds correctly (shock, contact discontinuity, rarefaction fan). It is among the most widely used flux schemes in compressible aerodynamics. An "entropy fix" is required to prevent non-physical expansion shocks at sonic points where the scheme otherwise produces numerically admissible but physically wrong solutions.
AUSM (Advection Upstream Splitting Method)
AUSM and its variants (AUSM+, AUSM+-up) split the flux into separate convective and pressure contributions. AUSM schemes are generally more robust than Roe's scheme at low Mach numbers — Roe's scheme adds excessive dissipation as Ma → 0 — making AUSM+-up particularly attractive for codes that must handle both subsonic and supersonic regions in the same computation.
TVD Schemes and Flux Limiters
Total Variation Diminishing (TVD) schemes guarantee that total variation of the numerical solution does not increase in time, preventing spurious oscillations from growing. Flux limiters locally reduce the scheme to first-order near discontinuities while maintaining higher-order accuracy in smooth regions. Common limiters include van Leer (smooth, second-order, well-balanced), minmod (most diffusive, most robust), and superbee (least diffusive but can clip smooth extrema).
Why do central differencing schemes create oscillations near shocks? I understand upwinding intuitively but not why central schemes fail at shocks specifically.
A shock is essentially an infinitely steep gradient — a discontinuity. Central differencing approximates derivatives by averaging left and right cell values symmetrically, without knowing which direction information comes from. Near a steep gradient, this symmetric stencil encounters a large jump, and the resulting polynomial interpolation tries to represent that jump smoothly — this is the Gibbs phenomenon, identical to what happens when you approximate a step function with a finite-term Fourier series. The oscillations (Gibbs wiggles) appear on both sides of the shock and can easily be large enough to drive density or pressure negative, instantly crashing the simulation. Upwind schemes "look upstream" — they use information only from the direction the signal actually comes from physically — so they cannot generate artificial oscillations from the downstream side. The cost is numerical dissipation, which smears the shock over 2–4 cells. The challenge in designing high-order compressible schemes is combining central accuracy in smooth regions with upwind stability near shocks — that's exactly what TVD, ENO (Essentially Non-Oscillatory), and WENO (Weighted ENO) schemes achieve through locally adaptive blending.
In OpenFOAM's rhoCentralFoam, how does the Kurganov-Tadmor scheme actually work? Is it an upwind scheme or central?
It's a central-upwind scheme — a blended approach. The Kurganov-Tadmor (KT) method computes cell-face fluxes by taking a weighted average of left and right flux estimates, with the weights determined by the local maximum wave speed (essentially the maximum eigenvalue of the flux Jacobian, which includes both the flow velocity and the speed of sound). Near shocks where the wave speed difference is large, the scheme adds proportionally more dissipation — it becomes more upwind-like. In smooth subsonic regions where the eigenvalues are moderate, the scheme is closer to central and preserves second-order accuracy. In OpenFOAM you control this through the fvSchemes dictionary's interpolationSchemes settings; the reconstruct(rho), reconstruct(U), etc. entries use MUSCL reconstruction with limiters, and you can choose vanLeer, minmod, or SuperBee depending on how much you want to prioritize sharp shock resolution vs. smooth-region accuracy. For most practical supersonic cases, vanLeer is a good default — robust but not excessively diffusive.
OpenFOAM and Fluent: Density-Based Solvers
OpenFOAM Compressible Solvers
| Solver | Algorithm | Best For | Notes |
|---|---|---|---|
rhoCentralFoam | Kurganov-Tadmor (density-based explicit) | High-speed flows with shocks, Ma > 1 | Excellent shock resolution; small explicit Δt (CFL < 1); best for shock tubes, supersonic nozzles |
rhoPimpleFoam | PIMPLE (pressure-based) | Subsonic/transonic compressible transient | Larger Δt than rhoCentralFoam; less suited to strong shocks; good for 0.3 < Ma < 1.2 |
rhoSimpleFoam | SIMPLE (pressure-based) | Steady compressible subsonic-transonic | Steady-state; turbomachinery at moderate Mach numbers |
sonicFoam | PISO (pressure-based) | Transonic/supersonic transient | Legacy solver; handles sonic transitions; rhoPimpleFoam generally preferred |
ANSYS Fluent: Density-Based Solver
Fluent's density-based solver is the correct choice for compressible flows at Ma > 0.3. Key settings:
- Solver type: Density-Based — solves continuity, momentum, and energy simultaneously in a coupled system
- Implicit time-stepping — allows CFL > 1 and is the standard choice for steady-state compressible simulations
- Flux scheme: Roe FDS (Flux Difference Splitting) or AUSM — both are upwind schemes appropriate for compressible flows
- High-order reconstruction: 3rd-order MUSCL for smoother shock profiles and better accuracy in smooth regions
- Equation of state: Ideal gas (default for air at moderate temperatures); real gas models (Peng-Robinson, NIST Refprop) available for supercritical fluids and very high temperatures
Engineering Applications of Compressible Flow Analysis
Jet Engines and Gas Turbines
Modern high-bypass turbofan engines have fan tip speeds near Ma = 1.3–1.5 and core compressor blades at transonic conditions. CFD predicts shock losses on blade suction surfaces, shock/boundary-layer interaction leading to separation, and total pressure loss through the compressor stages. The engine inlet must decelerate flight Mach number (~0.85) to compressor face Mach number (~0.5) with minimum total pressure recovery — a transonic diffuser design problem where shock strength management is critical.
Supersonic Aircraft and Missiles
Fighter aircraft at Ma = 1.5–2.5 require CFD for wave drag prediction, intake shock system design, and assessment of kinetic heating on leading edges and control surfaces. Supersonic cruise missiles use area-rule body shaping to minimize wave drag. Sonic boom prediction requires propagating the near-field pressure signature from CFD simulation through an atmospheric model to ground level — an active area of design for overland supersonic business jet certification.
Rocket Nozzles and Spacecraft Propulsion
Rocket nozzle contours are optimized to produce uniform supersonic exit flow with the target Mach number and minimal entropy generation. CFD evaluates over- and under-expanded operation (when ambient pressure differs from exit pressure at off-design altitudes), shock diamond patterns in the plume, and side-load forces during transient startup. Reentry vehicle aerodynamics involves hypersonic flow with real-gas chemistry — dissociation and ionization significantly alter heat loads compared to ideal-gas predictions.
Industrial Gas Systems
Safety relief valves, pressure regulators, control valves, and compressor check valves all involve compressible flow at pressure ratios that may exceed the critical value. CFD determines accurate discharge coefficients for valve certification, identifies regions of high velocity and erosion risk, and assesses noise generation from shock cells in high-velocity jets. Industrial piping system designers increasingly use CFD for valve sizing validation beyond the limitations of standard compressible flow discharge coefficient tables.