The Riemann Problem: The Heart of CFD#

In the Finite Volume Method (FVM), the accuracy and stability of the solution depend heavily on how the numerical flux at the cell boundaries is calculated. The key to determining this numerical flux is the Riemann problem.

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$$

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.

In the next article, we will compare interface capturing techniques (VOF, Level Set, Diffuse Interface).

Crank the left/right pressure ratio toward 1000:1 to see both waves turn into shocks.