HLLC Review

Reference book: [Riemann Solvers and Numerical Methods for Fluid Dynamics] by Eleuterio F. Toro

​ Fundamentally, the HLC method belongs to the Godunov-type numerical schemes, primarily solving the

hyperbolic equations which are also known as Euler Equation:

$$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = 0$$
​ When solving the above model equation, the Godunov flux $F_{i+1/2}$ is crucial. The HLLC method, like other Godunov-type methods, aims to computate this flux. HLLC method treats it as a Riemann problem.

​ The HLLC method assumes that three wave types emerge after a Riemann problem including the right wave , left wave and contact wave. They propagate at wave speed $S_R, S_L, S*$ respectively.

​ Note that it is not an exact wave structure.

​ Under this assumption, an approximate solution is given by :

$$\tilde{U}(x) = \begin{cases} U_L, & \frac{x}{t} \leq S_L \\ U_{*L}, & S_L \leq \frac{x}{t} \leq S_* \\ U_{*R}, & S_* \leq \frac{x}{t} \leq S_R \\ U_{R}, & S_R \leq \frac{x}{t} \end{cases}$$
​ The corresponding HLLC flux is :
$$F_{i+1/2}(x) = \begin{cases} U_L, & \frac{x}{t} \leq S_L \\ U_{*L}, & S_L \leq \frac{x}{t} \leq S_* \\ U_{*R}, & S_* \leq \frac{x}{t} \leq S_R \\ U_{R}, & S_R \leq \frac{x}{t} \end{cases}$$
​ By applying Rankine-Hugoniot conditions(conservation law), there are:
$$F_{*L} = F_{L} + S_{L}(U_{*L} - U_L) \\ F_{*R} = F_{R} + S_{R}(U_{*R} - U_R)$$
​ Note that there are four unknown quantities, that is $F_{*L}$, $U_{*L}$, $F_{*R}$, $U_{*R}$. $F_{L}$ and $F_R$ are actually $F_i$ and $F_{i+1}$.

​ For the left-state equation:

$$S_LU_{*L} - F_{*L} = S_LU_{L} - F_{L}$$
​ Focus on the components:
$$S_L \rho_{*L} - \rho_{*L}u_{*L} = S_L\rho_L - \rho_Lu_L \\ S_L \rho _{*L} u_{*L} - (\rho_{*L} u_{*L}^2 + p_{*L}) = S_L \rho _{L} u_{L} - (\rho_{L} u_{L}^2 + p_{L}) \\ S_{L} E_{*L} - u_{*L}(E_{*L} + p_{*L}) = S_L E_L - u_L(E_L + p_L)$$
​ From the first two equations, $p_{*L}$ can be derived as:
$$p_{*L} = \rho_{L}(S_L - u_L)(u_{*L} - u_L) + p_L$$
​ Similarly, $p_{*R}$ can be obtained as:
$$p_{*R} = \rho_{R}(S_R - u_R)(u_{*R} - u_R) + p_R$$
​ Consider that the contact wave introduces critical conditions:
$$p_{*L} = p_{*R} \\ u_{*L} = u_{*R} = S*$$
​ This implies the pressure is equal and the normal velocity matches the contact wave speed on both sides.

​ With these extra conditions, $S_{*}$ can be solved:

$$S_{*} = \frac{p_R - p_L + \rho_Lu_L(S_L - u_L) - \rho_Ru_R(S_R - u_R)}{\rho_L(S_L - u_L) - \rho_R(S_R - u_R)}$$
​ Substituting $S_*$ into equations (6), $U_{*L}$ components are:
$$\rho_L(\frac{S_L - u_L}{S_L - S_*}) \\ \rho_L(\frac{S_L - u_L}{S_L - S_*})S_* \\ \rho_L(\frac{S_L - u_L}{S_L - S_*})[\frac{E_L}{\rho_L}+(S_* - u_L)(S_*+\frac{p_L}{\rho_L(S_L-u_L)})]$$
​ $U_{*R}$ could be obtained easily by changing the subscript L to R;

​ And obviously, the flux $F_{*L}$ and $F_{*R}$ can be derived from the equations (4)

​ Before using the above results, $S_L$ and $S_R$ must be provided first.

​ So here comes an important problem: how to estimate the wave speed?

Wave speed estimate

Davis suggested:

$$S_L = u_L-a_L, S_R = u_R+a_R \\ S_L = \min(u_L-a_L,u_R-a_R), S_R = \max(u_L-a_L,u_R-a_R)$$

Davis and Einfeldt proposed to use Roe average eigenvalues:

$$S_L = \tilde{u} - \tilde{a}, S_R = \tilde{u} + \tilde{a} \\ \tilde{u} = \frac{\sqrt{\rho_L}u_L+\sqrt{\rho_R}u_R}{\sqrt{\rho_L}+\sqrt{\rho_R}},\tilde{a} = \sqrt{[(\gamma-1)(\tilde{H} - \frac{1}{2}{\tilde{u}}^2)]} \\ \tilde{H} = \frac{\sqrt{\rho_L}H_L+\sqrt{\rho_R}H_R}{\sqrt{\rho_L}+\sqrt{\rho_R}}$$
Einfeldt also proposed HLLE:
$$S_L = \bar{u} - \bar{d}, S_R = \bar{u} + \bar{d} \\ \bar{d}^2 = \frac{\sqrt{\rho_L}a_L^2+\sqrt{\rho_R}a_R^2}{\sqrt{\rho_L}+\sqrt{\rho_R}} + \eta_2(u_R - u_L)^2 \\ \eta_2 = \frac{1}{2}\frac{\sqrt{\rho_L}\sqrt{\rho_R}}{(\sqrt{\rho_L}+\sqrt{\rho_L})^2}$$
Davis suggested that about Rusanov flux:
$$F_{i+1/2} = \frac{1}{2}(F_L+F_R) - \frac{1}{2}S^+(U_R - U_L)\\ S^+ = \max(|u_L - a_L|, |u_R - a_R|, |u_L+a_L|, |u_R+a_R|)$$
Toro gave a pressure-based method:
$$S_L = u_L - a_Lq_L,S_R = u_R+a_Rq_R \\ q_K = \begin{cases} 1 ,&if & p^*\leq p_K \\ [1+\frac{\gamma + 1}{2\gamma}(\frac{p^*}{p_K} - 1)]^{1/2} ,& if & p^* > p_K \end{cases}$$

Personal annotations

Consider the one-dimension equations with sectional area:

$$\frac{\partial{\rho A}}{\partial t} + \frac{\partial{\rho u A}}{\partial x} = 0 \\ \frac{\partial{\rho u A}}{\partial t} + \frac{\partial{(\rho u^2+p) A}}{\partial x} = 0 \\ \frac{\partial{\rho A e_t}}{\partial t} + \frac{\partial{(\rho A e_t+pA) u}}{\partial x} = 0 \\ $$
The corresponding $U$ and $F$ :
$$U = [\rho A ,\rho u A ,EA] \\ F = [\rho u A, (\rho u^2+p)A, u(E+P)A]$$
To keep the conservation of discrete format, reformulate the HLLC scheme for the one-dimensional equations.

According to the formula (6) :

$$S_L \rho_{*L}A_{*L} - \rho_{*L}u_{*L}A_{*L} = S_L\rho_LA_L - \rho_Lu_LA_L \\ S_L \rho _{*L} u_{*L}A_{*L} - (\rho_{*L} u_{*L}^2 + p_{*L})A_{*L} = S_L \rho _{L} u_{L}A_L - (\rho_{L} u_{L}^2 + p_{L})A_L \\ S_{L} E_{*L}A_{*L} - u_{*L}(E_{*L} + p_{*L})A_{*L} = S_L E_L A_L - u_L(E_L + p_L)A_L$$

So, from the first two equations:

$$p_{*L}A_{*L} = [\rho_{L}(S_L - u_L)(u_{*L} - u_L) + p_L]{A_L}$$
Correspondingly, right state is:
$$p_{*R}A_{*R} = [\rho_{R}(S_R - u_R)(u_{*R} - u_R) + p_R]{A_R}$$
Introduce the condition of contact wave(Pay attention to the difference compared to the constant sectional area):
$$p_{*L}A_{*L} = p_{*R}A_{*R} \\ u_{*L} = u_{*R} = S*$$
Then the wave speed of star region is:
$$S_{*} = \frac{p_RA_R - p_LA_L + \rho_LA_Lu_L(S_L - u_L) - \rho_RA_Ru_R(S_R - u_R)}{\rho_LA_L(S_L - u_L) - \rho_RA_R(S_R - u_R)}$$
Till now, it can be demonstrated that the difference between constant area and variable area lies in the need to incorporate area into those intensive properties such as pressure and density.

With $S_{*}$, $U_{*L}$ is:

$$\rho_LA_L(\frac{S_L - u_L}{S_L - S_*}) \\ \rho_LA_L(\frac{S_L - u_L}{S_L - S_*})S_* \\ \rho_LA_L(\frac{S_L - u_L}{S_L - S_*})[\frac{E_L}{\rho_L}+(S_* - u_L)(S_*+\frac{p_L}{\rho_L(S_L-u_L)})]$$
$F_{*L}$ still satisfies that：
$$F_{*L} = F_{L} + S_{L}(U_{*L} - U_L)$$
Right star region could be obtained in the same way.