Compressible Flow: The Euler Equations and the Sod Shock Tube

Overview

So far the fluid tutorials have treated flow as incompressible — density is a constant. When speeds approach the speed of sound, that assumption breaks: the density itself becomes a conserved variable, and shocks, rarefactions, and contact discontinuities appear. This is compressible (gas) dynamics, governed by the Euler equations.

Rasmah implements a finite-volume solver for the Euler equations with the Rusanov (local Lax–Friedrichs) flux and explicit strong-stability-preserving time stepping, plus the exact Riemann solver needed to verify it on the canonical benchmark: the Sod shock tube (Sod 1978).

The Euler equations

In conservation form the compressible Euler system for a perfect gas is

\[\partial_t \begin{pmatrix}\rho\\ \rho u\\ \rho E\end{pmatrix} + \partial_x \begin{pmatrix}\rho u\\ \rho u^2 + p\\ u(\rho E + p)\end{pmatrix} = 0, \qquad p = (\gamma-1)\left(\rho E - \tfrac12 \rho u^2\right)\]

with ρ the density, u the velocity, E the total energy, p the pressure, and γ the ratio of specific heats. Rasmah’s EulerGas(γ) fixes the gas model (air, γ = 1.4).

The Sod shock tube

The Sod problem starts with a tube split by a membrane: a high-pressure, high-density state on the left and a low-pressure, low-density state on the right. When the membrane bursts, a shock runs right, a rarefaction fan runs left, and a contact discontinuity trails the shock. The exact solution is known in closed form, which makes it the standard test for any shock-capturing scheme.

using Rasmah
using LinearAlgebra

# Run the numerical solver: 1D Sod tube, t ∈ [0, 0.2].
xc, ρ, u, p = euler_sod_1d(200, 0.2)
println("grid points = ", length(xc), "   ρ ∈ [", round(minimum(ρ), digits=3), ", ", round(maximum(ρ), digits=3), "]")

grid points = 200   ρ ∈ [0.125, 1.0]

Verifying against the exact Riemann solution

The exact Sod solution is the reference sod_riemann(x, t, γ). The textbook post-shock density for the standard Sod initial data is the constant

\[\rho^* = 0.426319,\]

which the exact solution must reproduce on the plateau between the contact discontinuity and the shock. We check both claims:

# The exact solution's plateau density (excluding the undisturbed left state 1.0).
exact_ρ = [sod_riemann(x, 0.2, 1.4)[1] for x in range(0.5, 1.0, length=2000)]
plateau = maximum(filter(ρ -> ρ < 0.9, exact_ρ))
println("exact post-shock ρ* = ", round(plateau, digits=6))
println("textbook ρ*         = 0.426319")

exact post-shock ρ* = 0.426319
textbook ρ*         = 0.426319

The numerical solution is a first-order shock-capturing approximation: the rarefaction and contact regions match the exact solution closely, while the shock itself is smeared over a couple of cells (the price of the diffusive Rusanov flux, which is what keeps the scheme total-variation-stable):

err = maximum(abs(ρ[i] - sod_riemann(xc[i], 0.2, 1.4)[1]) for i in eachindex(xc))
println("max |ρ_num − ρ_exact| = ", round(err, digits=3), "   (shock smeared over a few cells)")

max |ρ_num − ρ_exact| = 0.088   (shock smeared over a few cells)

Refine the grid and the smeared error shrinks — the signature of a convergent shock-capturing scheme.

Next steps

• Navier–Stokes — the incompressible limit the Euler equations generalize.
• Stabilized flow — upwind/Petrov–Galerkin stabilization for convection-dominated problems, the incompressible analogue of the shock-capturing flux used here.