The Riemann Problem: The Heart of CFD#

In the Finite Volume Method (FVM), accuracy and stability hang on the numerical flux at a cell face. The key is a Riemann problem: two constant states, one jump.

Push the pressure-ratio slider from near 1 toward 1000:1 and both outer waves switch from rarefaction to shock.

The contact in the middle stays put as the ratio grows. Both outer waves are shocks together only once the star pressure $p^*$ sits above the pressure on each side.

What is the Riemann Problem?#

The Riemann problem is an initial value problem where the initial conditions consist of two constant states separated by a single discontinuity:

$$\mathbf{U}(x, 0) = \begin{cases} \mathbf{U}_L & \text{if } x < 0 \\ \mathbf{U}_R & \text{if } x > 0 \end{cases}$$

For the 1D Euler equations, the analytical solution to this problem consists of three waves:

1. Left-moving wave (rarefaction or shock)
2. Contact discontinuity
3. Right-moving wave (rarefaction or shock)

Godunov's Idea (1959)#

Godunov's key insight is simple yet powerful:

“

To find the flux at a cell boundary, treat the average values of the two adjacent cells as the left and right states and solve the Riemann problem.

The update formula for the Finite Volume Method:

$$\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x}\left(\hat{\mathbf{F}}_{i+1/2} - \hat{\mathbf{F}}_{i-1/2}\right)$$

where $\hat{\mathbf{F}}_{i+1/2}$ is the numerical flux determined from the solution to the Riemann problem.

Exact Riemann Solver#

An exact Riemann solver finds the pressure in the star region, $p^*$, using the Newton-Raphson iterative method. Rankine-Hugoniot relations or isentropic relations are applied depending on the type of wave (shock or rarefaction).

Shock wave relation ($p^* > p_K$, $K = L$ or $R$):

$$f_K(p^*) = (p^* - p_K)\left[\frac{A_K}{p^* + B_K}\right]^{1/2}$$

where $A_K = \frac{2}{(\gamma+1)\rho_K}$ and $B_K = \frac{\gamma-1}{\gamma+1}p_K$.

Rarefaction wave relation ($p^* \leq p_K$):

$$f_K(p^*) = \frac{2c_K}{\gamma - 1}\left[\left(\frac{p^*}{p_K}\right)^{\frac{\gamma-1}{2\gamma}} - 1\right]$$

$p^*$ satisfies the following condition:

$$f_L(p^*) + f_R(p^*) + (u_R - u_L) = 0$$

The split is one comparison: $p^*$ against $p_K$. If $p^*$ is above $p_K$, that side is a shock. $f_K$ is the pressure jump times $[A_K/(p^*+B_K)]^{1/2}$. $A_K=2/((\gamma+1)\rho_K)$, so a denser fluid changes speed less for the same pressure jump. $B_K=(\gamma-1)p_K/(\gamma+1)$ fattens the denominator a little past $p^*$. If $p^*$ is at or below $p_K$, that side is a rarefaction. The factor $2c_K/(\gamma-1)$ grows with sound speed, and it grows again as $\gamma$ nears 1. The pressure-ratio exponent is $(\gamma-1)/(2\gamma)$. Lower $\gamma$, smaller exponent, and $f_K$ cares less about the pressure ratio.

The two $f$ values plus $u_R-u_L$ have to sum to zero. If $u_R$ is below $u_L$, fluid piles into the middle. Then $f_L+f_R$ is positive and $p^*$ is pushed above both side pressures. Both outer waves are shocks. The 1000:1 end of the slider is that limit. The low-side $p_K$ sits under the star pressure, so that $f_K$ stays on the shock branch. The high side can remain a rarefaction, depending on the ratio. Both sides are shocks only when the velocity difference closing on the contact is large enough. A cell treated as its own shock tube is the calculation in Godunov's shock tube inside one cell.

Approximate Riemann Solver: Why is it Needed?#

Exact solvers are expensive because they require iterative calculations. Most CFD codes use approximate Riemann solvers.

Roe Solver (1981)#

Roe linearizes the non-linear Riemann problem into a linear Riemann problem for the matrix $\hat{\mathbf{A}}$:

$$\hat{\mathbf{F}}_{i+1/2} = \frac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \frac{1}{2}\sum_{k=1}^{3}|\hat{\lambda}_k|\hat{\alpha}_k\hat{\mathbf{r}}_k$$

The Roe average is used to construct $\hat{\mathbf{A}}$:

$$\hat{u} = \frac{\sqrt{\rho_L}\,u_L + \sqrt{\rho_R}\,u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}, \quad \hat{H} = \frac{\sqrt{\rho_L}\,H_L + \sqrt{\rho_R}\,H_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}$$

Pros: Captures contact discontinuities accurately. Cons: Can violate entropy conditions (requires an entropy fix).

HLLC Solver (Toro, 1994)#

HLLC improves the HLL solver by restoring the contact discontinuity. It estimates three wave speeds: $S_L$, $S^*$, and $S_R$:

$$\hat{\mathbf{F}}_{i+1/2} = \begin{cases} \mathbf{F}_L & \text{if } S_L > 0 \\ \mathbf{F}^*_L & \text{if } S_L \leq 0 < S^* \\ \mathbf{F}^*_R & \text{if } S^* \leq 0 \leq S_R \\ \mathbf{F}_R & \text{if } S_R < 0 \end{cases}$$

Star region flux:

$$\mathbf{F}^*_K = \mathbf{F}_K + S_K(\mathbf{U}^*_K - \mathbf{U}_K)$$

HLLC is easy to implement and robust, making it the most widely used in compressible multiphase flow codes.

AUSM+ (Liou, 1996)#

The AUSM family splits the flux into convective and pressure terms:

$$\hat{\mathbf{F}}_{i+1/2} = \dot{m}_{1/2}\boldsymbol{\Psi}_{1/2} + \hat{p}_{1/2}\mathbf{D}$$

where $\dot{m}_{1/2}$ is the mass flow rate at the cell boundary and $\boldsymbol{\Psi}$ is the convected quantity. It is stable even in low-Mach flows, making it easy to extend to all-speed schemes.

Extension to Multiphase Flows#

In multiphase flows, since the EOS on both sides of the interface are different, Riemann solvers must be generalized to multi-material versions.

Common approaches:

• Ghost Fluid Method: Constructs ghost cells for the fluid on the opposite side near the interface and solves independent single-material Riemann problems for each phase.
• HLLC for multi-material: Determines $p^*$ and $u^*$ by applying each EOS on both sides of the star region.

The next article compares interface capturing techniques (VOF, Level Set, Diffuse Interface).

Splitting advection from pressure once the flow is multiphase continues in TV splitting of the Baer–Nunziato flux.