9. Compressible and High-Speed Flow Computation

Introduction and Context

Compressible-flow computation is required whenever variations in density materially influence the motion, pressure field, temperature, or energy of a fluid. In low-speed liquid flows and many low-Mach-number gas flows, density may be treated as constant without serious error. In high-speed aerodynamics, however, pressure disturbances travel at a finite speed and density changes are coupled to momentum and energy transport. Shock waves, expansion fans, stagnation-temperature changes, entropy production, and chemical nonequilibrium may then determine the flow physics and the engineering performance of a system.

The governing equations are the compressible Navier–Stokes equations, or their inviscid form, the Euler equations. Their numerical solution is more demanding than that of incompressible flow because the equations are hyperbolic in supersonic regions, mixed elliptic–hyperbolic in transonic regions, and strongly nonlinear near shocks. A successful computational method must therefore represent wave propagation, impose boundaries according to the direction of information travel, preserve conservation, and prevent nonphysical oscillations.

This chapter develops the computational foundations of subsonic, transonic, supersonic, and hypersonic flow simulation. It emphasizes conservative formulations, characteristic boundary conditions, Riemann solvers, flux-splitting methods, shock-capturing schemes, entropy generation, and high-temperature effects. The discussion assumes familiarity with finite-volume discretization, numerical stability, and the basic conservation laws introduced in earlier CFD study.

Learning Objectives

After studying this chapter, the reader should be able to:

1. Explain how Mach number and wave propagation determine the character of compressible flow.

2. Formulate the Euler equations in conservative variables and interpret their fluxes.

3. Distinguish subsonic, transonic, supersonic, and hypersonic computational regimes.

4. Apply characteristic reasoning to inflow, outflow, wall, and far-field boundary conditions.

5. Compare approximate Riemann solvers, flux-vector splitting, and flux-difference splitting.

6. Select shock-capturing methods while controlling numerical diffusion and oscillations.

7. Interpret shocks, expansion waves, entropy generation, and stagnation-property changes.

8. Recognize when thermally perfect or chemically reacting gas models are necessary.

9. Diagnose common failures in high-speed CFD and make defensible solver choices.

Governing Equations and the Meaning of Compressibility

For a calorically perfect ideal gas, the conservative state vector in two dimensions is

U = [rho, rho u, rho v, rho E]^T,

where rho is density, u and v are velocity components, and E is total specific energy. The inviscid Euler equations are

∂U/∂t + ∂F/∂x + ∂G/∂y = 0,

with

F = [rho u, rho u²+p, rho uv, u(rho E+p)]^T,

G = [rho v, rho uv, rho v²+p, v(rho E+p)]^T.

The equation of state is p = (gamma−1)rho e, where e is internal energy per unit mass and E = e + (u²+v²)/2. For air treated as a calorically perfect gas, gamma is often taken as approximately 1.4 over moderate temperature ranges. At elevated temperatures this approximation becomes inadequate because molecular vibrational modes, dissociation, ionization, and variable specific heats alter the thermodynamic closure.

The Mach number, M = V/a, compares flow speed V with the local speed of sound a = (gamma p/rho)^(1/2). It is not merely a label for speed. It determines the number and direction of characteristic waves and therefore influences both the physics and the numerical boundary treatment.

In a one-dimensional inviscid flow, the characteristic speeds are u−a, u, and u+a. The acoustic characteristics carry pressure and velocity disturbances, while the convective characteristic transports entropy and, in multidimensional settings, vorticity-related information. At a subsonic boundary, some characteristics enter and some leave the domain. At a supersonic inflow, all characteristics enter; at a supersonic outflow, all leave. This fact provides the rigorous basis for counting boundary conditions.

Flow Regimes and Computational Consequences

Subsonic flow has M<1. Information can propagate both upstream and downstream, so an outlet condition can influence the interior solution. In a duct, pressure, velocity, and density are strongly coupled through acoustic waves. Pressure-based boundary specifications must therefore be used carefully, since over-specification can reflect waves and contaminate the solution.

Transonic flow contains both subsonic and supersonic regions, often separated by a shock. The governing character changes within the domain, and a small alteration in geometry, back pressure, or numerical dissipation can move a shock substantially. Transonic computation is especially sensitive to grid clustering, limiter behavior, convergence acceleration, and the treatment of the energy equation.

Supersonic flow has M>1. In a uniform region, disturbances generally travel downstream within the flow. Compression waves may coalesce into shocks, while expansion waves spread into Prandtl–Meyer fans. Boundary conditions become more directional: an inviscid supersonic outlet generally requires no imposed thermodynamic data, whereas a supersonic inlet requires the complete incoming state.

Hypersonic flow is conventionally associated with M greater than about 5, although the relevant physics depends on altitude, length scale, and temperature. Strong bow shocks, severe aerodynamic heating, vibrational excitation, dissociation, ionization, catalytic wall effects, and radiative phenomena may occur. A perfect-gas Euler solution can still provide useful preliminary information about shock location and pressure, but it may be unreliable for heat transfer, surface chemistry, and thermal loads.

Finite-Volume Formulation and Conservation

The finite-volume method integrates the conservation law over each control volume. For a cell i, the semi-discrete equation is

dUi/dt = −(1/Vi) sum over faces [A_f F_f*],

where F_f* is a numerical flux and A_f is the oriented face area. Conservation follows because the flux leaving one cell enters its neighbor with the opposite sign. This property is essential for shock computation: a discontinuity must satisfy the integral conservation laws even when its internal thickness is represented by only a few cells.

A numerical flux is normally constructed from left and right reconstructed states at a face. These states are inserted into a local Riemann problem, which asks how two constant states separated by an interface evolve. Exact Riemann solutions are expensive and may require iterative thermodynamic calculations, so practical CFD uses approximate solvers.

The numerical method must balance three requirements. It must be conservative, stable in the presence of discontinuities, and sufficiently accurate in smooth regions. Excessive dissipation smears contact surfaces and weak shocks; insufficient dissipation produces oscillations or negative pressure and density. Higher-order reconstruction, slope limiters, positivity preservation, and suitable time integration are therefore used together rather than as isolated improvements.

Riemann Solvers and Numerical Fluxes

The local Riemann problem for the Euler equations contains nonlinear waves: a left-going acoustic wave, a contact wave, and a right-going acoustic wave. Approximate solvers differ in how they estimate wave speeds and intermediate states.

The HLL solver represents the solution using only the fastest left- and right-moving waves. It is robust and relatively simple, but it does not explicitly represent the contact wave and can be excessively diffusive near shear layers and material interfaces. HLLC restores the contact wave and usually gives sharper resolution of contact discontinuities and boundary layers in inviscid calculations.

Roe’s approximate solver linearizes the flux Jacobian about a Roe-averaged state. It resolves contact waves sharply and is efficient for many engineering problems. Its weaknesses include possible nonphysical states near strong rarefactions, entropy-violating expansion shocks, and difficulties when density or pressure approaches zero. Entropy fixes and positivity safeguards are commonly added.

The exact or approximate acoustic solvers used in the AUSM family separate convective and pressure transport. Such methods can provide low numerical diffusion at contact surfaces and perform well over a broad Mach-number range, although their formulation and parameter choices require care near low-speed limits and strong shocks.

Flux-vector splitting methods divide the physical flux according to the signs of characteristic speeds. Steger–Warming and van Leer splitting are classical examples. They are naturally upwind and robust, but may introduce excessive diffusion, particularly at low Mach numbers or across contact discontinuities. Flux-difference splitting instead compares reconstructed states and approximates the wave contributions directly, as in Roe-type methods.

No single solver is universally best. HLLC is often a good general-purpose choice for shock-containing ideal-gas flows; Roe may be attractive when sharp wave resolution is important and safeguards are available; more dissipative HLL-type approaches may be preferable when robustness is the primary concern.

Shock Capturing, Reconstruction, and Limiters

A shock-capturing method does not explicitly track a shock as a moving mesh surface. Instead, it represents the discontinuity as a narrow numerical transition whose position and integrated jump satisfy the conservation laws. The shock thickness depends on grid resolution, reconstruction order, limiter behavior, and flux dissipation.

First-order upwind schemes are robust but smear shocks and contact surfaces. Second-order methods reconstruct piecewise-linear states and can achieve much better accuracy in smooth regions. Near discontinuities, however, an unrestricted high-order reconstruction creates Gibbs-type oscillations. Total-variation-diminishing limiters, such as minmod, van Leer, and Superbee, constrain slopes to prevent new extrema. More compressive limiters preserve discontinuities sharply but may reduce smooth-region accuracy or create sensitivity on irregular grids.

MUSCL reconstruction is widely used in finite-volume compressible CFD. Primitive, conservative, or characteristic variables may be reconstructed. Characteristic reconstruction is often advantageous near shocks because it limits wave amplitudes rather than component variables independently. Modern high-order schemes include WENO and discontinuous Galerkin methods, which are useful for complex wave fields but require careful positivity, troubled-cell, and time-step treatments.

The Courant–Friedrichs–Lewy condition restricts the time step. For an explicit method, a representative condition is

Delta t <= CFL min over cells [Delta x/(|u|+a)],

with the multidimensional form based on the sum of directional wave speeds. Strong shocks, small cut cells, and high acoustic speeds can make the allowable time step very small. Implicit schemes, local time stepping, dual-time stepping, and multigrid methods improve efficiency, but convergence acceleration must not be confused with physical time accuracy.

Characteristic Boundary Conditions

Boundary conditions should be imposed on incoming characteristic information rather than by prescribing every primitive variable independently. For a one-dimensional subsonic inlet, one acoustic characteristic may leave the domain and the remaining characteristics enter; consequently, only the appropriate number of independent conditions should be specified. A common practical arrangement specifies total pressure and total temperature together with flow direction, while the outgoing characteristic is obtained from the interior solution.

At a subsonic outlet, the pressure is often prescribed and the remaining information is extrapolated from the interior. If the outlet is too close to a shock, recirculation region, or strong wave system, simple extrapolation can reflect disturbances. The computational domain should be extended or a nonreflecting characteristic boundary condition should be used.

At a supersonic inlet, all characteristics enter and the full thermodynamic and velocity state must be specified. At a supersonic outlet, all characteristics leave and no physical state variable should normally be imposed. Imposing an outlet pressure in this situation can generate an artificial reflected wave.

For an inviscid slip wall, normal velocity is zero and the normal momentum flux is balanced by wall pressure. For a viscous no-slip wall, velocity is zero relative to the wall, while thermal conditions may be adiabatic, isothermal, catalytic, or coupled to conjugate heat transfer. Far-field boundaries should reproduce the correct characteristic behavior and minimize reflections; simple fixed primitive values are acceptable only when the boundary is sufficiently distant and disturbances are weak.

Shocks, Expansions, and Entropy Generation

Across a stationary normal shock, mass, momentum, and total enthalpy are conserved, but total pressure decreases. Static pressure, temperature, and density increase, while Mach number falls from supersonic to subsonic. The loss of total pressure is a measure of irreversible entropy generation. Numerically, a shock-capturing method must reproduce the Rankine–Hugoniot jump conditions in the grid-integrated sense.

Oblique shocks occur when supersonic flow is turned toward itself, as over a compression corner or a wedge. Their strength depends on upstream Mach number, turning angle, and shock angle. If the required turning is too large, the attached solution may disappear and a detached bow shock forms. Expansion turning is isentropic in ideal inviscid theory and produces a centered Prandtl–Meyer fan. The flow accelerates, pressure and temperature decrease, and total pressure remains constant.

A numerical shock contains artificial entropy production associated with dissipation. Some entropy production is necessary for stability, but excessive production causes inaccurate drag, stagnation-pressure loss, and shock position. Entropy residuals, total-pressure contours, and grid-refinement studies are useful diagnostics. A solution that appears smooth may still be inaccurate if the shock is over-diffused.

High-Temperature and Nonequilibrium Effects

The calorically perfect gas model assumes constant specific heats and fixed chemical composition. A thermally perfect gas allows temperature-dependent specific heats while retaining fixed composition. This distinction becomes important when post-shock temperatures are high enough to excite molecular rotation and vibration.

At still higher temperatures, oxygen and nitrogen dissociate, and subsequent ionization produces charged species. A chemically reacting model must then solve species conservation equations, use reaction-rate mechanisms, and close the system with mixture thermodynamics and transport properties. Vibrational relaxation may occur on a time scale comparable to the flow residence time, producing thermal nonequilibrium. Separate translational–rotational and vibrational–electronic temperatures may be required.

These models introduce stiffness because chemical and relaxation time scales can be much shorter than the fluid-dynamic time scale. Implicit integration, operator splitting, adaptive chemistry, and positivity-preserving species treatment are common remedies. The model should be selected according to the engineering question: perfect-gas Euler calculations may estimate external pressure distributions, whereas predicting heat flux or ablation requires appropriate transport, radiation, surface chemistry, and wall models.

Practical Workflow for High-Speed CFD

A defensible simulation begins with the physical regime and the quantity of interest. The analyst should identify expected Mach number, Reynolds number, temperature range, shock structure, wall condition, and whether chemistry is active. A suitable equation set and thermodynamic model should then be selected.

The mesh must resolve geometric features, boundary layers where viscous effects matter, and expected shock paths. Grid alignment can reduce numerical diffusion, but excessive alignment with an assumed shock may hide errors if the shock moves. Initial calculations should often use a robust first-order or more dissipative solver, followed by higher-order reconstruction after the flow field is established.

Verification requires checking residual trends, conservation errors, positivity, grid sensitivity, and comparison with analytical or benchmark solutions. Validation requires experimental data or trusted reference calculations. For transonic and hypersonic problems, integrated force and heat-transfer predictions may converge more slowly than residuals. Engineering decisions should therefore be based on monitored outputs, not residual reduction alone.

Illustrative Case Study: Transonic Flow over a Supercritical Aerofoil

Background

Consider an illustrative computational investigation of a supercritical aerofoil at a moderate angle of attack and a freestream Mach number slightly below unity. The objective is to predict lift, drag, and the location of the upper-surface shock.

Situation and Problem

A preliminary calculation using a first-order upwind scheme predicts a broad shock and excessive wave drag. A second calculation using an unrestricted second-order scheme produces oscillations immediately downstream of the shock and a nonphysical pressure pattern.

Analysis

The first result is dominated by numerical diffusion, which spreads the shock and weakens the pressure jump. The second result violates the monotonicity required near a discontinuity. Because the flow is transonic, the correct solution contains both subsonic and supersonic regions, and the shock location is coupled to the downstream pressure field. The far-field boundary must therefore be sufficiently distant, and the outlet condition must allow outgoing information without artificial reflection.

Approach and Solution

The Euler equations are solved using a conservative finite-volume method with MUSCL reconstruction, a bounded limiter, and an HLLC-type flux. Freestream total conditions and flow direction are imposed through characteristic far-field treatment. The wall is treated as inviscid and adiabatic. The mesh is refined near the aerofoil and along the expected shock path, and calculations are repeated on successively finer grids. A more dissipative HLL calculation is used as a robustness comparison.

Outcome

The bounded second-order solution produces a compact shock without spurious oscillations. The predicted shock position changes modestly under refinement, while the first-order result shows a larger change in drag and pressure distribution. The HLLC solution resolves the pressure recovery and contact-like features more sharply, but the HLL comparison provides a useful conservative estimate of numerical uncertainty.

Learning Points

The case demonstrates that solver selection cannot be separated from reconstruction and boundary treatment. A formally higher-order method is not automatically more accurate near discontinuities. Grid studies, integrated-output monitoring, and comparison between solvers are essential for judging whether a transonic prediction is physically credible.

Limitations, Risks, and Professional Practice

Common errors include prescribing all primitive variables at every boundary, imposing outlet pressure in a fully supersonic outflow, using a perfect-gas model for high-temperature chemistry, and judging accuracy from residuals alone. Other risks include negative density or pressure from aggressive reconstruction, carbuncle-type shock instabilities on aligned grids, excessive low-Mach dissipation, and unrecognized reflection from truncated far fields.

Compressible CFD also has professional and safety implications. Predictions of aerodynamic heating, propulsion performance, blast loading, or re-entry environments may influence design decisions with serious consequences. Assumptions, uncertainty ranges, mesh quality, solver settings, and validation evidence should be documented. Results should not be presented with unjustified precision, and proprietary or experimental data must be handled according to applicable confidentiality and export-control requirements.

Current Developments and Best Practices

Current research emphasizes high-order shock-capturing schemes, entropy-stable and kinetic-energy-preserving discretizations, positivity-preserving algorithms, adaptive mesh refinement, and improved all-Mach-number methods that remain accurate from low-speed to hypersonic conditions. Machine-learning tools are being explored for turbulence closure, reduced-order modeling, mesh adaptation, and surrogate prediction, but they require carefully defined training domains and independent physical validation.

Best practice is to select the simplest model that answers the engineering question, use conservative discretization, impose characteristic-consistent boundaries, begin with a robust solver, and increase order only after stability is established. Verification and validation should be treated as part of the simulation, not as an optional final step.

Chapter Summary

Compressible-flow computation is governed by the interaction of conservation laws, finite-speed wave propagation, and thermodynamic closure. Mach number determines the characteristic structure and therefore the appropriate boundary conditions. Conservative finite-volume methods combined with approximate Riemann solvers, flux splitting, reconstruction, and limiting provide the principal tools for computing shocks and expansions. Subsonic, transonic, supersonic, and hypersonic regimes differ in both physical behavior and numerical sensitivity. Shocks conserve mass, momentum, and total enthalpy while generating entropy and reducing total pressure; expansion fans are ideally isentropic. At high temperatures, variable specific heats, vibrational relaxation, dissociation, ionization, and surface chemistry may require advanced nonequilibrium models. Reliable practice depends on model selection, positivity, boundary treatment, grid refinement, verification, validation, and transparent reporting of uncertainty.

Key Terms

Mach number: Ratio of flow speed to local speed of sound.

Speed of sound: Small-disturbance propagation speed determined by the thermodynamic state.

Conservative variables: Density, momentum components, and total energy used in conservation-law formulations.

Euler equations: Inviscid compressible conservation equations.

Characteristic: Direction and speed along which information propagates.

Riemann problem: Initial-value problem involving two constant states separated by an interface.

Riemann solver: Numerical approximation to the solution of a local Riemann problem.

Flux-vector splitting: Division of physical flux into components associated with positive and negative wave speeds.

Flux-difference splitting: Upwind evaluation based on state differences and characteristic waves.

Shock wave: Compressional discontinuity across which entropy increases.

Expansion fan: Continuous, approximately isentropic region produced by supersonic turning away from itself.

Total pressure: Stagnation pressure, which decreases across an irreversible shock.

Entropy generation: Irreversible thermodynamic production associated with losses.

Shock capturing: Representation of a discontinuity within a conservative numerical transition zone.

Limiter: Reconstruction device that suppresses spurious extrema near discontinuities.

CFL condition: Time-step restriction imposed by numerical wave propagation.

Calorically perfect gas: Gas model with constant specific heats and fixed composition.

Thermally perfect gas: Gas model with temperature-dependent specific heats and fixed composition.

Chemical nonequilibrium: State in which reaction rates and composition do not instantaneously follow temperature and pressure.

Positivity preservation: Numerical property maintaining physically admissible density, pressure, and species concentrations.

Review Questions

Short Answer Questions

1. Why does Mach number influence the number of boundary conditions required in a compressible-flow calculation?

2. State two differences between a shock wave and a Prandtl–Meyer expansion fan.

3. What is the principal advantage of HLLC over the basic HLL solver?

4. Why are slope limiters used in high-resolution shock-capturing methods?

5. When does a calorically perfect gas model become questionable?

Long Answer Questions

1. Derive and interpret the conservative variables and fluxes in the two-dimensional Euler equations.

2. Explain characteristic boundary conditions for subsonic inflow, subsonic outflow, supersonic inflow, and supersonic outflow.

3. Compare Roe, HLL, HLLC, and flux-vector-splitting approaches in terms of resolution, robustness, and limitations.

4. Discuss the numerical and physical sources of entropy generation across a computed shock.

5. Describe a complete verification and validation strategy for a transonic aerofoil simulation.

Analytical / Application Questions

1. A supersonic duct calculation imposes static pressure at the outlet even though all characteristics leave the domain. Predict the likely numerical consequence and propose a better treatment.

2. A hypersonic vehicle calculation predicts surface pressure accurately but gives highly uncertain heat flux. Analyze which physical models may be missing and identify additional validation data required.

3. Two transonic simulations give different shock locations: one uses first-order HLL and the other uses second-order HLLC with limiting. Develop a diagnostic plan to determine whether the difference is physical, grid-related, or numerical.

22. Multiple Choice Questions with Answers

Question 1

In a fully supersonic outflow, which boundary treatment is generally appropriate?

A. Prescribe static pressure and temperature

B. Prescribe all primitive variables

C. Allow all characteristics to leave without imposing physical data

D. Impose total pressure only

Correct Answer: C. Allow all characteristics to leave without imposing physical data

Question 2

What is the principal weakness of the basic HLL Riemann solver?

A. It cannot conserve mass

B. It does not explicitly resolve the contact wave

C. It requires chemical reaction rates

D. It is restricted to incompressible flow

Correct Answer: B. It does not explicitly resolve the contact wave

Question 3

Across an ideal stationary normal shock, which quantity decreases?

A. Static pressure

B. Density

C. Total pressure

D. Static temperature

Correct Answer: C. Total pressure

Question 4

The main purpose of a slope limiter is to:

A. Increase the speed of sound

B. Suppress nonphysical oscillations near discontinuities

C. Eliminate the energy equation

D. Force every flow to be subsonic

Correct Answer: B. Suppress nonphysical oscillations near discontinuities

Question 5

A subsonic outlet commonly requires which physical information to be specified?

A. The complete state vector

B. One incoming characteristic, often represented by pressure

C. No boundary information under any condition

D. All outgoing wave amplitudes

Correct Answer: B. One incoming characteristic, often represented by pressure

Question 6

Which method usually restores the contact wave omitted by HLL?

A. HLLC

B. First-order central differencing

C. Minmod reconstruction alone

D. The CFL condition

Correct Answer: A. HLLC

Question 7

Why may a perfect-gas model be inadequate for hypersonic heat-transfer prediction?

A. Mach number has no relevance to heat transfer

B. High temperatures can cause variable specific heats, dissociation, and ionization

C. The Euler equations always predict zero pressure

D. Hypersonic flow contains no shocks

Correct Answer: B. High temperatures can cause variable specific heats, dissociation, and ionization

Question 8

In an ideal Prandtl–Meyer expansion, total pressure is expected to:

A. Remain constant

B. Decrease sharply because of entropy generation

C. Become zero

D. Equal static pressure everywhere

Correct Answer: A. Remain constant

Question 9

What is the key conservation advantage of a finite-volume method at a shock?

A. It removes all numerical dissipation

B. Fluxes cancel between neighboring cells, preserving integral conservation

C. It avoids the need for a thermodynamic equation of state

D. It guarantees exact shock thickness

Correct Answer: B. Fluxes cancel between neighboring cells, preserving integral conservation

Question 10

If a second-order reconstruction creates oscillations behind a shock, the most appropriate immediate remedy is to:

A. Remove conservation from the discretization

B. Use a bounded limiter or a more dissipative robust flux

C. Prescribe supersonic outlet pressure

D. Increase the angle of attack

Correct Answer: B. Use a bounded limiter or a more dissipative robust flux