Shock Tube Problem (Riemann Solver)

The shock tube problem (represented by the 1D Riemann problem like the Sod problem) is the most fundamental benchmark for evaluating the accuracy and stability of numerical schemes for compressible flow. Because an exact analytical solution is available, it allows for quantitative evaluation of numerical diffusion, overshoot, and shock wave resolution by comparing with the numerical solution.

It is formulated as an initial value problem for the 1D Euler equations. There is a diaphragm (partition) in the center of the tube, with the left side in a high-pressure state $(\rho_L, u_L, p_L)$ and the right side in a low-pressure state $(\rho_R, u_R, p_R)$. When the diaphragm ruptures, three waves are generated.

1. Leftward traveling expansion wave (rarefaction fan)

2. Central contact discontinuity

3. Rightward traveling shock wave

The conservative form of the 1D Euler equations is as follows.

The total energy is $E = \frac{p}{\gamma - 1} + \frac{1}{2}\rho u^2$.

The relations derived from the conservation laws across the shock wave are the Rankine-Hugoniot relations. Letting the shock wave speed be $s$,

Here $[\cdot]$ represents the difference across the shock wave. Expressed in terms of pressure ratio,

where $M_s$ is the shock Mach number. Combining this with the continuity conditions for pressure and velocity at the contact discontinuity and the Riemann invariants within the expansion wave, all state variables in the four regions (left undisturbed, inside expansion fan, star region, right undisturbed) can be determined.

It's the region on either side of the contact discontinuity where pressure and velocity are equal ($p^* = p_L^* = p_R^*$, $u^* = u_L^* = u_R^*$) but density is discontinuous. The equation for this $p^*$ is nonlinear, so it is solved iteratively using the Newton-Raphson method.

Godunov's idea is to obtain the flux at each cell interface from the solution of a local Riemann problem. At the boundary between cell $i$ and $i+1$, solve the Riemann problem consisting of the left state $\mathbf{U}_i$ and right state $\mathbf{U}_{i+1}$, and determine the interface flux $\mathbf{F}_{i+1/2}$.

The exact Riemann solver requires iterative Newton method each time, which is computationally expensive. Therefore, many practical approximate Riemann solvers have been developed. Let's compare some representative ones.

HLL considers only two waves, so information about the contact discontinuity (the third wave) is lost. The 'C' in HLLC stands for Contact wave; by adding the third wave, it becomes capable of resolving the contact discontinuity. Toro's textbook has a detailed derivation.

Roe's idea is to linearize the nonlinear Euler equations at the interface. Construct the Jacobian matrix $\hat{A}$ using Roe-averaged states,

The Roe-averaged density and velocity are defined as follows.

It can sometimes generate non-physical solutions called expansion shocks. This violates the entropy condition and needs to be corrected with the Harten-Hyman entropy fix. Also, in low Mach number regions, excessive numerical diffusion can occur, leading to proposed fixes like the Low-Mach Roe correction (Rieper fix, etc.).

The standard approach is to achieve second-order accuracy using the MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) method. Extrapolate left and right states at cell interfaces via linear reconstruction.

Here $\phi(r)$ is the TVD limiter function, and $r$ is an indicator of solution smoothness. By choosing the limiter, it automatically switches to second-order accuracy in smooth regions and first-order accuracy near discontinuities.