The Sod Gasdynamics Problem as a Tool for Benchmarking Face Flux Construction in the Finite Volume Method

1. Introduction

The term Riemann problem (Riemann, 1860; Castro-Orgaz and Hager, 2019) in basic gasdynamics refers to a setting where two adjacent regions of a compressible gas are initially separated by some barrier. Each region has a uniform set of conditions: density, velocity, and pressure. Suddenly, the barrier is removed and the energy imbalance between the two sides causes three waves such that some equilibrium average condition is encouraged. These waves can be a shock wave (normal shock), a contact discontinuity (contact surface), or an expansion fan (rarefaction fan). It is possible to have two shocks with no fans, two fans with no shocks, and one shock and one fan. In either case, a contact discontinuity would be separating the two other waves. There is also a possibility of reaching a vacuum state in zone surrounded by two contacts (one to the left and one to the right), which in turn separate two fans running away from each other (Kamm, 2015). However, the presence of vacuum is excluded in our work.

The simple Riemann problem is governed by the Euler equations (a system of time-dependent nonlinear hyperbolic partial differential equations), which admit discontinuities in the solution. For a computational fluid dynamics (CFD) solver built on the finite volume method (FVM) that is designed to handle compressible flows, the Riemann problem can be a good environment for testing the capability of the solver in forming fluxes at the finite volume face. The problem is focused on the convective/pressure flux in the one-dimensional space, thus eliminating complications of multi-dimensionality, inaccuracy of representing non-uniform boundary or initial conditions, and interference from other phenomena such as viscous effects and boundary layers, chemical reaction, mixing, multi-phase flow, heat transfer, thermal radiation, and source terms.

While the Riemann problem can be extended to more complicated environments (Hantke et al., 2020), the aim of this work is to solve a particular example of the basic Riemann problem known as the Sod problem with different numerical methods for building a flux vector at an cell interfaces, and analyze the results in comparison to the exact solution as well as to other numerical flux methods.

2. One-Dimensional Euler Equations

Euler formulation of a fluid motion neglects friction (viscosity) effects completely (Fletcher, 1991), leading to a reduced form of the Navier-Stokes equations. While Euler formulation may appear narrowly-restricted, there are many conditions where the effect of fluid viscosity can be safely neglected and this simplified formulation becomes accurate enough.

In vector form, the one-dimensional Euler equations governing inviscid fluids is (Hoffmann and Chiang, 1993)

$${\displaystyle \frac{\partial q}{\partial t}}+{\displaystyle \frac{\partial F}{\partial x}}=0$$

(1)

where $\left(q\right)$ is a vector of the conserved variables, and $\left(F\right)$ is a vector of flux (sum of convective flux in the x-direction and a pressure contribution).

$$q=\left\{\begin{array}{c}\rho \\ \rho u\\ \rho { e}_T\end{array}\right\}$$

(2)

In the above equation, $\left(\rho \right)$ is the fluid density, which is the mass per unit volume; $\left(u\right)$ is the fluid velocity, so $(\rho u)$ is the momentum per unit volume; $\left(e_T\right)$ is the total specific internal energy, so $(\rho e_T)$ is the total internal energy per unit volume. The density $\left(\rho \right)$, velocity $\left(u\right)$, and pressure $\left(p\right)$ are referred to as the primitive variables.

The total specific internal energy is the specific internal energy $\left(e\right)$ plus the kinetic energy per unit mass, thus

$$e_T=e+{\displaystyle \frac{1}{2}}{u}^2$$

(3)

The flux vector is

$$F=\left\{\begin{array}{c}\rho u\\ \rho u^2+p\\ \left(\rho { e}_T+p\right) u\end{array}\right\}=\left\{\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho u h_T\end{array}\right\}$$

(4)

where $\left(h_T\right)$ is the total specific enthalpy

$$h_T={ e}_T+{\displaystyle \frac{p}{\rho }}=e+{\displaystyle \frac{1}{2}}{u}^2+{\displaystyle \frac{p}{\rho }}$$

(5)

The specific enthalpy is

$$h=e+{\displaystyle \frac{p}{\rho }}=h_T-{\displaystyle \frac{1}{2}}{u}^2$$

(6)

For a thermally perfect gas, the pressure, density, and absolute temperature $\left(T\right)$ are related as

$$p=\rho R T$$

(7)

where $\left(R\right)$ is the specific gas constant. Also, for a thermally perfect gas, the specific heat capacity at constant volume $\left(C_v\right)$ and the specific heat capacity at constant pressure $\left(C_p\right)$both vary with the temperature only, and the following relation holds true:

$$C_p\left(T\right)=C_v\left(T\right)+R$$

(8)

The ratio of the specific heat capacities $\left(\gamma \right)$, also called the adiabatic index, is

$$\gamma ={\displaystyle \frac{C_p}{C_v}}$$

(9)

The following equations apply to a thermally perfect gas

$$C_v\left(T\right)={\displaystyle \frac{R}{\gamma \left(T\right)-1}}$$

(10)

$$C_p\left(T\right)={\displaystyle \frac{\gamma \left(T\right) R}{\gamma \left(T\right)-1}}$$

(11)

$$a=\sqrt{\gamma {\displaystyle \frac{ p}{\rho }} }=\sqrt{\gamma R T }$$

(12)

$${\displaystyle \frac{p}{\rho }}={\displaystyle \frac{a^2}{\gamma }}$$

(13)

$$h=e+{\displaystyle \frac{1}{\gamma }}{a}^2$$

(14)

$$h_T=e_T+{\displaystyle \frac{1}{\gamma }}{a}^2$$

(15)

$$\rho u^2+p=\rho a^2\left(M^2+{\displaystyle \frac{1}{\gamma }}\right)$$

(16)

where $\left(a\right)$ is the speed of sound in the gas (the sonic speed), and $(M=u/a)$ is the Mach number. If a thermally perfect gas has also the specialty that its specific enthalpy and its specific internal energy can be related as

$$h\left(T\right)=\gamma e\left(T\right)$$

(17)

with $\left(\gamma \right)$ being now a constant independent of the temperature, then the gas is described as calorically perfect and its specific heat capacities are also independent of the temperature (Fletcher, 2004). The calorically perfect gas assumption is usually reasonable for compressible flows at low speeds. This simplification permits analytical analysis for supersonic flows (NACA, 1953).

For a calorically perfect gas, the following formulas apply

$$p=\left(\gamma -1\right)\rho e$$

(18)

$$e={\displaystyle \frac{p}{\rho \left(\gamma -1\right)}}={\displaystyle \frac{1}{\gamma \left(\gamma -1\right)}}{a}^2$$

(19)

$$h={\displaystyle \frac{\gamma }{\gamma -1}} {\displaystyle \frac{ p}{\rho }}={\displaystyle \frac{1}{\gamma -1}}{a}^2$$

(20)

$$\gamma e_T={\displaystyle \frac{\gamma }{2}}{u}^2+{\displaystyle \frac{1}{\gamma -1}}{a}^2$$

(21)

$$h_T={\displaystyle \frac{1}{2}}{u}^2+{\displaystyle \frac{1}{\gamma -1}}{a}^2$$

(22)

$$\rho u h_T=\rho a^3\left({\displaystyle \frac{M^3}{2}}+{\displaystyle \frac{M}{\gamma -1}}\right)$$

(23)

This work considers a calorically perfect gas that is non-reacting and inviscid.

3. Sod’s Problem

A shock tube is a special case of the Riemann problem, in which the compressible gas on either side is initially stationary (Laney, 1998). Figure 1 gives a graphical illustration of the Riemann problem at the initial state. The Sod problem (Sod, 1978) is a shock tube problem with certain initial left and right densities and pressures, causing a right-going shock wave, a left-going expansion fan, and a right-going contact discontinuity in between. Table 1 contrasts the left and right conditions in the Sod problem at the initial state $(\mathbf{i}\mathbf{n}\mathbf{i}\mathbf{t}\mathbf{i}\mathbf{a}\mathbf{l}\mathbf{t}\mathbf{i}\mathbf{m}\mathbf{e}:\mathit{t}=0)$. The initial distribution of the primitive variables is demonstrated in Figure 2. The jump is located at a position $\mathit{x}=0.5$, which is midway between the left end of the problem domain $\left(\mathit{x}=0\right)$ and its right end $(\mathit{x}=1)$. Table 2 gives a summary of the fan-contact-shock wave pattern in the exact solution of the Sod problem.

At a time of $t=0.2$, the exact solution we obtained for the Sod problem is illustrated in Figure 3, expressed in terms of the primitive variables and the Mach number. The abrupt change in the pressure designates the location of a shock wave traveling to the right. The abrupt change in the density where no other variable showing any change designates the location of a contact discontinuity, which also travels to the right. The gradual change in all the variables designates the extent of an expansion fan, which travels to the left. The left front of the fan is called the head, and the right edge is called the tail.

4. Validation

The computation of the exact solution presented earlier in this work, as well as the numerical ones to be presented later, is done using computer codes written in GNU Octave, a scientific programming language (Eaton, 2020). As our exact solution is to be used as a reference for judging numerical ones, it is important to check its validity. This is supported here in two ways.

The first validity support of our implementation of the exact Riemann solver is through satisfying the Rankine-Hugoniot shock condition (Rankine 1870; Hugoniot 1887, 1889; Koval and Szabo, 2008) for the mass. The condition is used to compute a shock speed $\left(S\right)$ other than the one obtained from the exact Riemann solver (was given already in Table 2) as

$$S={\displaystyle \frac{\mathsf{\Delta }\rho u}{\mathsf{\Delta }\rho }}={\displaystyle \frac{\left(0.26557\right) \left(0.92745\right)-\left(0.125\right)\left(0\right) }{0.26557-0.125}}=1.75216$$

(24)

where the operator $\mathsf{\Delta }$ represents the difference across the shock wave. This speed matches the one given earlier in Table 2, which was found independently as

$$S=u_{RIGHT}+a_{RIGHT}\sqrt{1+{\displaystyle \frac{\gamma +1}{2 \gamma }}\mathrm{}\left({\displaystyle \frac{p_{CD}}{p_{RIGHT}}}-1\right) }$$

(25)

where $u_{RIGHT}$, $a_{RIGHT}$, and $p_{RIGHT}$ are the gas velocity, speed of sound, and pressure in the right undisturbed state; and $p_{CD}$ is the pressure around the contact discontinuity.

The second validity support of the provided exact solution of the Sod problem is through comparison with published exact and computational solutions at the same time $(t=0.2)$ taken from the Castro project validation data (Castro, 2020). Figure 5 shows these published results for the primitive variables. Visual inspection suggests they match with our solution given in Figure 3.

5. Face Flux Construction

The time integration of the conserved variables is done using Godunov scheme (Godunov, 1959), which is a conservative explicit finite volume formula. For a uniformly spaced one-dimensional array of computational cells having a cell width of $\left(\mathsf{\Delta }x\right)$, let the cell center positon index be $\left(i\right)$, then the left face flux vector is $\left(F_{i-0.5}\right)$ and the right face flux vector is $\left(F_{i+0.5}\right)$. Time integration of the conserved variables vector $\left(q\right)$ over one time step $\left(\mathsf{\Delta }t\right)$ is done as

$$q_i^{n+1}=q_i^n+{\displaystyle \frac{\mathsf{\Delta }t}{\mathsf{\Delta }x}}\mathrm{}\left(F_{i+0.5}^n-F_{i-0.5}^n\right)$$

(26)

The time index $\left(n\right)$ refers to the latest known variable values, while the index $(n+1)$ refers to the new values after time integration over one time step. The scheme is explicit, where the updating of a cell value is independent of updating other cells. This eliminates the need for solving a system of algebraic equations. Godunov scheme is of first order accuracy in time. The spatial accuracy depends on how the face flux vectors are constructed. It can be first order (lowest resolution) or second order (higher resolution).

Godunov scheme appears simple and straightforward, provided that the face flux vectors have been computed and are made available to the scheme.

The process of preparing face flux vectors at intercell faces is split into two steps:

• Step 1 (MUSCL): use neighbor cell-center values of the primitive variables $(\rho , u, p)$ to extrapolate and/or interpolate a ‘fictitious’ value at both sides of the facial interface. The left side of the face has the values denoted by: ${\rho }_L, u_L, p_L$. These are the locally left values. Similarly, the right values of the face are denoted by: ${\rho }_R, u_R, p_R$. These are the locally right values. It should be noted that although physically these locally left and locally right values are defined at the same spatial location, they are not necessarily equal. Once these fictitious face-left and face-right values are computed, a local fictitious Riemann problem is constructed. There is a jump in one or more flow variables across the face. Along with the interpolation/extrapolation, the preparation of the local Riemann problem at each face includes also the use of flux limiters (also called slope limiters) through a scheme called MUSCL, which stands for Monotonic Upstream-centered Scheme for Conservation Laws (van Leer, 1979). This scheme is a nonlinearized version of the so-called $\kappa$ 3-point scheme of parabolic interpolation/extrapolation (Waterson and Deconinck, 2007). Whereas the interpolation/extrapolation coefficients in the $\kappa$ scheme are solution-independent, the MUSCL scheme adds to the $\kappa$ scheme a solution-dependent limiter, which downgrades the interpolation/extrapolation level near sharp changes (like a shock wave) to use only the nearest cell-center value, while upgrades it to allow use of more cell-center values in regions of smooth variations. This adaptive switching allows the use of the high-resolution with its advantage of spatial accuracy without suffering from severe oscillations near sharp changes by ‘switching off’ the high-resolution capability and keeping only the first-order spatial accuracy. There are different flux limiters (Chu and Gao, 2013; Li et al., 2020). We here use the van Leer flux limiter (van Leer, 1974), which has intermediate characteristics among others in terms of the dissipation level (Bai et al., 2018).
• Step 2 (flux method): After a local fictitious Riemann problem is constructed, then comes the step of inferring a face flux vector based on the jump across the local left and local right flow variables. This step is the core of this work. We here consider 22 methods of face flux methods proposed in the literature, and compare their performance for the Sod problem as a test case.

To explain the MUSCL step, consider a generic variable $\left(v\right)$ that has a value in four adjacent cell centers as shown in Figure 4. The cell values are labeled in the direction of the increasing x-coordinate as $v_{MM}, v_M, v_P, \mathrm{a}\mathrm{n}\mathrm{d} v_{PP}$. The goal is to construct a limited interpolated/extrapolated fictitious left and right values at the face separating the cell $M$ and the cell $P$.

Let

$${\Delta }^M=v_M-v_{MM}$$

(27)

$${\Delta }^C=v_P-v_M$$

(28)

$${\Delta }^P=v_{PP}-v_P$$

(29)

Two consecutive gradient ratios are calculated, left one and right one, respectively, as

$$r_L={\displaystyle \frac{{\mathsf{\Delta }}^C}{{\mathsf{\Delta }}^M}}$$

(30)

$$r_R={\displaystyle \frac{{\mathsf{\Delta }}^C}{{\mathsf{\Delta }}^P}}$$

(31)

There is a possibility of a division by zero when calculating $r_L$ and $r_R$. However, this is avoided by applying the following condition

$$\mathrm{i}\mathrm{f} \left|{\mathsf{\Delta }}^M\right|\le ϵ, { r}_L=0$$

(32)

$$\mathrm{i}\mathrm{f} \left|{\mathsf{\Delta }}^P\right|\le ϵ, { r}_R=0$$

(33)

where $ϵ$ is a small threshold, which was decided by the computing program as $2.22\times {10}^{-16}$.

Then, two flux limiting values ${\varphi }_L$ and ${\varphi }_R$ are obtained for the left and right sides of the face, respectively, using the van Leer flux limiter function.

$${\varphi }_L={\displaystyle \frac{r_L+\left|r_L\right|}{1+\left|r_L\right|}}$$

(34)

$${\varphi }_R={\displaystyle \frac{r_R+\left|r_R\right|}{1+\left|r_R\right|}}$$

(35)

Now, the final stage is to compute the fictitious left and right face values $v_L$ and $v_R$, respectively, using the following version of the MUSCL scheme:

$$v_L=v_M+0.5 {\varphi }_L{\mathrm{}\mathsf{\Delta }}^M$$

(36)

$$v_R=v_P-0.5 {\varphi }_R {\mathsf{\Delta }}^P$$

(37)

In the case of unity flux limiters $\left({\varphi }_L={\varphi }_R=1\right)$, the above equations represent a pseudo linear-upwind interpolation (LUI), which is also called second order upwind (SOU). In that case, each of the fictitious face value is constructed by a standard linear extrapolation from the two nearest cell-center values on the same face side (either left or right). This pure extrapolation (no interpolation) feature corresponds to setting the parameter value $\kappa =-1$ in the linear version of the MUSCL scheme (the $\kappa$ scheme).

On the other hand, if the flux limiters are set to zero $\left({\varphi }_L={\varphi }_R=0\right)$, then a pseudo first order (classical upwinding) interpolation is obtained, where each fictitious face value is simply taken as the nearest cell-center value on the same face side (either left or right). The term ‘pseudo’ is used in the above text because in true upwinding schemes, taking the left or right cell-center values depends on the speed of the characteristic (Strikwerda, 2004), also called the wave speed or the information speed, relevant to the concerned variable at the face. Such checking of the direction of the information is not done here, and it is always assumed that for the left fictitious face value, the flow information (the wind) comes from the left, and for the right fictitious face value, it comes from the right.

6. Analyzed Flux Methods

Twenty two flux methods are tested in this work and they are listed here with a brief description on each. Giving detailed explanation or the mathematical formulation behind each method can make this section extremely long. Instead, the interested reader is cordially requested to visit the literature for deeper information about any of methods found of interest. In presenting the data, a short name (a code) is assigned to each method as a quick unique identifier. The methods are listed below by their codes. The sequence of introducing the methods here is the same one used later when presenting their graphical results and the detailed root-mean-square errors of their numerical results.

6.1. Riemann

This flux method is also called Godunov solver. Godunov is credited to formulate this exact (analytic) Riemann solver in the same reference cited before when addressing the Godunov time integration scheme. This flux method solves a Riemann problem exactly at each intercell face, with the local fictitious face-left and face-right conditions provided by the MUSCL scheme replacing the global initial physical left and right states of the gas in the Sod problem. This local Riemann solver identifies the state of the gas at the face interface, such as a zone between a left-going expansion fan and a right-going contact discontinuity, or a zone between a right-going shock wave and a right-going contact discontinuity, or a zone within an expansion fan whose tail and head go opposite to each other, which is called a transonic fan as it accelerates the gas from subsonic (Mach number < 1) to supersonic (Mach number > 1). Then, the face flux is constructed using the primitive values of the identified state. Thus, if the identified state has numerical values of the density, velocity, and pressure denoted by ${\rho }_0,u_0,\mathrm{a}\mathrm{n}\mathrm{d} p_0$, then the constructed face flux is simply

$F_{Face}=\left\{\begin{array}{c}{\rho }_0 u_0\\ {\rho }_0{u}_0^2+p_0\\ \left({\displaystyle \frac{\gamma p_0}{\left(\gamma -1\right)}}+{\displaystyle \frac{1}{2}}{{\rho }_{0}u}_0^2\right)u_0\end{array}\right\}$

6.2. Roe

This is the Roe approximate solver for face flux, which is also known as Roe upwind scheme (Roe, 1981). The nonlinear flux vector in Euler equations is replaced by a locally linearized approximate flux function, which utilizes matrix called a Roe-average Jacobian matrix. Roe flux method gives three equally-spaced waves. It cannot capture the spread of the expansion fan. Instead, it represents it as a zero-thickness wave. Denote the MUSCL returned face-left values of the density and velocity by ${\rho }_L$ and $u_L$, respectively; and the MUSCL returned face-right values of the density and velocity by ${\rho }_R$ and $u_R$, respectively. Then, Roe method introduces a customized square-root-density-weighted velocity defined as: $u_{Roe}=\left(u_L\sqrt{{\rho }_L}+u_R\sqrt{{\rho }_R}\right)/\left(\sqrt{{\rho }_L}+\sqrt{{\rho }_R}\right)$. A face-left specific total enthalpy is $h_{T,L}=\gamma p_L/\left(\gamma -1\right)+0.5 {\rho }_L{u}_L^2$, where $p_L$ is the face-left pressure as provided by the MUSCL scheme. Similarly, a face-right specific total enthalpy is $h_{T,R}=\gamma p_R/\left(\gamma -1\right)+0.5 {\rho }_R{u}_R^2$, where $p_R$ is the face-right pressure as provided by the MUSCL scheme. The Roe method introduces a square-root-density-weighted total specific enthalpy defined as: $h_{T,Roe}=\left(h_{T,L} \sqrt{{\rho }_L}+h_{T,R}\sqrt{{\rho }_R}\right)/\left(\sqrt{{\rho }_L}+\sqrt{{\rho }_R}\right)$. Finally, a Roe-specific speed of sound is defined as: $a_{Roe}=\sqrt{\left(\gamma -1\right)\left(h_{T,Roe}-0.5 u_{Roe}^2\right)}$. The three wave speeds in Roe method (also called Roe eigenvalues) are: $u_{Roe}, u_{Roe}-a_{Roe}, \mathrm{a}\mathrm{n}\mathrm{d} u_{Roe}+a_{Roe}.$The three Roe specific quantities ($u_{Roe}$, $h_{T,Roe}$, and $a_{Roe}$) are used to construct the face flux vector.

6.3. KNP

This is the central-upwind method of Kurganov, Noelle, and Petrova for the face flux (Kurganov at al., 2001). The term central-upwind here comes from the notice that the method is central in principle (thereby, offering high resolution with small dissipation), but still has an upwinding character by respecting the direction of the wave travel. The method does not attempt to approximate the exact flux as in the case of Roe method.

6.4. KT

This is another central-upwind method by Kurganov and Tadmor for forming the intercell face flux vector (Kurganov and Tadmor, 2000). It resembles the KNP method in structure, with the KT method (which is chronologically earlier) being simpler. For example, the KT method assigns fixed weights of 0.5 to the face-left and face-right MUSCL flux vectors, whereas the KNP method uses flow-based weights.

6.5. SW

This is the Steger and Warming flux-vector splitting (FVS) method for face flux (Steger and Steger, 1979 and 1981). The face flux vector is computed as a sum of two flux vector: a plus one and a minus one, as $F=F^{+}+F^{-}$. The subvector $F^{+}$is associated with the positive eigenvalues of the Jacobian matrix $\partial F/\partial q$, while the subvector $F^{-}$is associated with its negative eigenvalue es.

6.6. vanLeer

This is another flux vector splitting (FVS) method, made by van Leer (van Leer, 1982). It allows continuous slope of two formulation functions at the special values of Mach numbers -1, 0, and 1. A quadratic function of the fictitious face-left and face-right Mach numbers are introduced. This addresses a shortcoming in the in the chronologically earlier SW flux vector splitting method.

6.7. AUSM

This is the Advection Upstream Splitting Method (AUSM) for face flux construction (Liou and Steffen, 1993). The flux is divided into a pressure component and a convective component for better treatment.

6.8. AUSM⁺

This is a follow-up improvement (Liou, 1996) of the Advection Upstream Splitting Method (AUSM) for computing face fluxes by addressing observed shortcomings in some cases in which the original AUSM technique did not maintain its tradition of accuracy and efficiency.

6.9. AUSM⁺-up (for All Speeds)

It is an extension of AUSM⁺ to have the convergence rate nearly independent of the Mach number for the low-Mach subsonic limit up to a freestream Mach number of 0.5. This is achieved through proper scaling using asymptotic analysis (Liou, 2006). Also the convergence rate is nearly constant in transonic and supersonic flows.

6.10. AUFS

AUFS stands for Artificially Upstream Flux Vector Splitting method (Sun and Takayama, 2003). It employs two simultaneous flux-vector splitting mechanisms, one by dividing into pressure and non-pressure components, and another by using a method-specific Mach number (one that is based on the simple average of the face-left and face-right velocities, $u_{avg}=(u_L+u_R)/2$) as a splitting weight. The method also respects the sign of the average velocity of the face-left and face-right fictitious values. The method is relatively easy to implement.

6.11. HLL-Davis1

This is the Harten, Lax, van Leer (HLL) 2-wave model (Harten et al., 1983), with Davis first estimation of right and left wave speeds (Davis, 1988). Let ${\rho }_L$, $u_L$ , and $p_L$ be the face-left fictitious values of the density, velocity, and pressure, respectively; and ${\rho }_R$, $u_R$ , and $p_R$ be their face-right counterparts. Then, a face-left speed of sound is calculated as $a_L=\sqrt{\gamma p_L/{\rho }_L }$, and face-right speed of sound is calculated as $a_R=\sqrt{\gamma p_R/{\rho }_R }$. This version of the HLL family of methods sets the left wave speed as $S_L=u_L-a_L$, and the right wave speed as $S_R=u_R+a_R$.

6.12. HLL-Davis2

This is a second variant of the Harten, Lax, van Leer (HLL) 2-wave model, with Davis second estimation of right and left wave speeds. With the same meaning of $u_L,a_L, u_R,a_R$ as used in the HLL-Davis1 method earlier, this other version of the HLL family sets the left wave speed as $S_L=\min \left(u_L-a_L, u_R-a_R \right)$, and the right wave speed as $S_R=\max \left(u_L+a_L, u_R+a_R \right)$. The $\min \left( \right)$ function returns the smaller of the two arguments, and the $\max \left( \right)$ function returns the larger of them.

6.13. HLL-Roe

This is a third variant of the Harten, Lax, van Leer (HLL) 2-wave model, with two of the Roe eigenvalues used for the right and left wave speeds. Using same notations in the Roe method and the HLL-Davis1 method, the two wave speeds in the HLL-Roe method described here are: $S_L=u_{Roe}-a_{Roe}$ and $S_R=u_{Roe}+a_{Roe}$.

6.14. HLL-Einfelt

This is a fourth variant of the Harten, Lax, van Leer (HLL) 2-wave model, with Einfelt eigenvalues for the right and left wave speeds (Einfeldt, 1988). It is also called (HLLE) flux method. The left and right wave speeds here are: $S_L=u_{Roe}-u_{Einfelt}$ and $S_R=u_{Roe}+u_{Einfelt}$, respectively. The symbol $u_{Roe}$ refers to the Roe square-root-density-averaged velocity as described in the Roe method. The velocity $u_{Einfelt}$ is computed as: $u_{Einfelt}=\sqrt{\left(a_L^2\sqrt{{\rho }_L}+a_R^2\sqrt{{\rho }_R}\right)/\left(\sqrt{{\rho }_L}+\sqrt{{\rho }_R}\right) +0.5\sqrt{{\rho }_L}\sqrt{{\rho }_R}{ \left(u_R-u_L\right)}^2/{\left(\sqrt{{\rho }_L}+\sqrt{{\rho }_R}\right)}^2}$. A subscript $L$ refers to a face-left value, and a subscript $R$ refers to a face-right value. The face-left speed of sound $a_L$ and the face-right speed of sound $a_R$ are as discussed in the HLL-Davis1 method.

6.15. HLL-pBased

This is a fifth variant of the Harten, Lax, van Leer (HLL) 2-wave model, with pressure-based calculation for right and left wave speeds (Toro et al., 1994). The method introduces a customized threshold or criterion pressure $p^{*}$ used to determine the left and right wave speeds, where $p^{*}=0.5\left(p_L+p_R\right)+0.125 \left(u_R-u_L\right)\left(a_R+a_L\right)\left({\rho }_R+{\rho }_L\right)$. The left and right wave speeds here are: $S_L=u_L-f_L{a}_L$ and $S_R=u_R+f_R{a}_R$, respectively. The symbols in these formulas have the same meanings as those given for the HLL-Davis1 method. Two new symbols appear here ($f_L$ and $f_R$), which are factors decided based on $p^{*}$ and how it compares with the face-left pressure $p_L$ and the face-right pressure $p_R$, as follows: if $p^{*}\le p_L, f_L=1$, if $p^{*}>p_L, f_L=\sqrt{1+(p^{*}/p_L-1)(\gamma +1)/\left(2\gamma \right)}$, if $p^{*}\le p_R, f_R=1$, and if $p^{*}>p_R, f_R=\sqrt{1+(p^{*}/p_R-1)(\gamma +1)/\left(2\gamma \right)}$.

6.16. HLLC-Davis1

This is an extended 3-wave version (Toro et al., 1992) of the Harten, Lax, van Leer, Davis (HLL-Davis1) 2-wave model, where a third wave (contact discontinuity or contact surface) is included. The additional "C" letter in the method code stands for the added "contact" wave. This method for flux construction is comparable in accuracy and robustness to that of the exact Riemann flux method, while being simpler and computationally less demanding than it. Davis first estimation of right and left wave speeds is used.

6.17. HLLC-Davis2

This is an extended 3-wave version of the Harten, Lax, van Leer, Davis (HLL-Davis2) 2-wave model, where a third wave (contact discontinuity) is added.

6.18. HLLC-Roe

This is an extended 3-wave version of the Harten, Lax, van Leer, Roe (HLL-Roe) 2-wave flux model, where a third wave (contact discontinuity) is added.

6.19. HLLC-Einfelt

This is an extended 3-wave version of the Harten, Lax, van Leer, Einfelt (HLL-Einfelt or HLLE) 2-wave flux model, where a third wave (contact discontinuity) is added.

6.20. HLLC-pBased

This is an extended 3-wave version of the Harten, Lax, van Leer, pressure-based (HLL-pBased) 2-wave flux model, where a third wave (contact discontinuity) is added.

6.21. LF

The Lax-Friedrichs (LF) method, which is also called the Lax method (Lax, 1954) is a basic explicit three-point technique used for solving time-dependent partial differential equations. It is first-order accurate in space. This method gives the most dissipative flux among all basic stable explicit flux construction methods. The method does not require solving a Riemann problem (either exact or approximate) at each cell face. The face flux is constructed as $F_{face}=0.5 \left(F_L+F_R\right)-0.5\mathsf{\Delta }x/\mathsf{\Delta }t\left(q_R-q_L\right)$, where $F_L$ and $F_R$ are the flux vector after substituting with the face-left and face-right values of the primitive variables, respectively; $q_L$ and $q_R$ are the conserved variables vectors after substituting with the face-left and face-right values of the primitive variables, respectively; $\mathsf{\Delta }x$ is the cell width (the spatial step); and $\mathsf{\Delta }t$ is the time step.

6.22. Rusanov

This is a single-wave simplification of Roe 3-wave model (Rusanov, 1970). It is also called generalized Lax-Friedrichs flux method. Like the Lax-Friedrichs flux method, this method does not require solving a Riemann problem (either exact or approximate) at each cell face. However, unlike the Lax-Friedrichs flux method, it makes use of the local speed of sound. The flux calculation formula is nearly identical to the one given for the Lax-Friedrichs flux method, with the only difference is that the uniform virtual speed $\left(\mathsf{\Delta }x/\mathsf{\Delta }t\right)$ is replaced by a gas-dependent speed ${\lambda }_{max}=\max \left(\left|u_L\right|+a_L, \left|u_R\right|+a_R\right)$, where $\left|u_L\right|$ and $\left|u_R\right|$ are the absolute values of the face-left and face-right gas velocities, respectively; and $a_L$ and $a_R$ are the speeds of sounds obtained using the face-left and face-right densities and pressures, as mentioned for the HLL-Davis1 method.

7. Graphical Results

The Sod problem is solved numerically using each of the listed 22 flux methods in the previous section. No other change was done from one numerical simulation to another. Common simulation settings include a cell size of $\mathsf{\Delta }x=0.005$ (thus, the shock tube of unit length is divided into 200 equal cells), and a time step of $\mathsf{\Delta }t=0.001$, which was related to the cell size as

$$\mathsf{\Delta }t={\displaystyle \frac{{\mathrm{C}\mathrm{o}}_{max} \mathsf{\Delta }x}{S_{max}}}$$

(38)

where ${\mathrm{C}\mathrm{o}}_{max}$ is the desired upper limit of the Courant number (Courant–Friedrichs–Lewy, CFL, number) as a condition for spatio-temporal stability in explicit time integration of partial differential equations, requiring the Courant number (Co) to be below unity. As the Co decreases, stabilizing numerical dissipation increases (Chang and Wang, 2003). The parameter $S_{max}$ is an estimate of the maximum wave speed in the problem. We set $S_{max}=2$, which is a reasonable value for the highest possible sum of the local absolute gas velocity and the speed of sound in the problem, $\left|u\right|+a$. The target maximum Courant number was set to ${\mathrm{C}\mathrm{o}}_{max}=0.4$. This target was largely satisfied in all performed 22 simulations, as found by post-processing the simulation results.

In Figures 6–27, two plots are assigned per figure, and each figure represents results using one of the flux methods. The order of presenting the graphical results of the flux methods is the same one used when briefly describing them earlier. The first plot in each figure is the spatial distribution of the three primitive variables (density $\rho$, gas velocity $u$, and pressure $p$) at the final simulation time of 0.2 (thus, after 200 time steps). The other plot in each figure shows the spatial distribution of the local Mach number (which is initially zero everywhere) at the final simulation time of 0.2. The title of each plot includes the code (the short name) of the flux method used in the simulation. In addition, the maximum observed Co number (CFL number) is also listed in the plots titles.

The LX method in Figure 26 shows significant undesirable dissipation causing severe smearing of jumps in the solution. No other method showed such level of unacceptable performance.

AUSM (Figure 12) exhibits a small overshoot at the tail of the fan. AUSM⁺ (Figure 13) did not resolve this problem, but it disappears in the further-improved version of AUSM⁺-up (Figure 14).

The flux methods of Roe, KT, SW, AUFS, HLL-Davis1, HLL-Roe, HLL-Einfelt, HLL-pBased, HLLC-Davis1, HLLC-Roe, HLLC-Einfelt, HLLC-pBased, and Rusanov show small velocity oscillations (wiggles) to the left of the shock (i.e., after the gas is shocked). These wiggles are particularly noticeable with the SW method.

The flux methods of the exact Riemann solver (Figure 6), KNP (Figure 8), vanLeer (Figure 11), AUSM⁺-up (Figure 14), HLL-Davis2 (Figure 17), and HLLC-Davis2 (Figure 22) appear to perform the best in terms of being free from visual oscillations in the solution.

Within the expansion fan, near its head, the methods of AUSM (Figure 12), AUSM⁺ (Figure 13), and AUFS (Figure 15) do not capture well the spatial increase of the gas velocity. They produce a fan head that is shifted to the right of its exact location. On the other hand, the methods of KT (Figure 9), SW (Figure 10), and vanLeer (Figure 11) satisfactorily predict the propagation of the fan head.

Figure 6. Calculated final primitive variables and Mach number, using face flux of exact Riemann solver (also called Godunov flux).

Figure 6. Calculated final primitive variables and Mach number, using face flux of exact Riemann solver (also called Godunov flux).

Figure 7. Calculated final primitive variables and Mach number, using face flux of Roe.

Figure 7. Calculated final primitive variables and Mach number, using face flux of Roe.

Figure 8. Calculated final primitive variables and Mach number, using face flux of Kurganov, Noelle, and Petrova (KNP).

Figure 8. Calculated final primitive variables and Mach number, using face flux of Kurganov, Noelle, and Petrova (KNP).

Figure 9. Calculated final primitive variables and Mach number, using face flux of Kurganov and Tadmor (KT).

Figure 9. Calculated final primitive variables and Mach number, using face flux of Kurganov and Tadmor (KT).

Figure 10. Calculated final primitive variables and Mach number, using face flux of Steger and Warming (SW).

Figure 10. Calculated final primitive variables and Mach number, using face flux of Steger and Warming (SW).

Figure 11. Calculated final primitive variables and Mach number, using face flux of van Leer.

Figure 11. Calculated final primitive variables and Mach number, using face flux of van Leer.

Figure 12. Calculated final primitive variables and Mach number, using face flux of the Advection Upstream Splitting Method (AUSM).

Figure 12. Calculated final primitive variables and Mach number, using face flux of the Advection Upstream Splitting Method (AUSM).

Figure 13. Calculated final primitive variables and Mach number, using face flux of the first improved version of the Advection Upstream Splitting Method (AUSM⁺).

Figure 13. Calculated final primitive variables and Mach number, using face flux of the first improved version of the Advection Upstream Splitting Method (AUSM⁺).

Figure 14. Calculated final primitive variables and Mach number, using the face flux of the second improved version of the Advection Upstream Splitting Method (AUSM⁺-up, for all speeds).

Figure 14. Calculated final primitive variables and Mach number, using the face flux of the second improved version of the Advection Upstream Splitting Method (AUSM⁺-up, for all speeds).

Figure 15. Calculated final primitive variables and Mach number, using face flux of the Artificially Upstream Flux Vector Splitting (AUFS) method.

Figure 15. Calculated final primitive variables and Mach number, using face flux of the Artificially Upstream Flux Vector Splitting (AUFS) method.

Figure 16. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the first estimation of Davis for the right/left wave speeds.

Figure 16. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the first estimation of Davis for the right/left wave speeds.

Figure 17. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the second estimation of Davis for the right/left wave speeds.

Figure 17. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the second estimation of Davis for the right/left wave speeds.

Figure 18. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the eigenvalues of Reo for the right/left wave speeds.

Figure 18. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the eigenvalues of Reo for the right/left wave speeds.

Figure 19. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the eigenvalues of Einfelt for the right/left wave speeds.

Figure 19. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with the eigenvalues of Einfelt for the right/left wave speeds.

Figure 20. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with pressure-based calculation for the right/left wave speeds.

Figure 20. Calculated final primitive variables and Mach number, using face flux of Harten, Lax, van Leer (HLL), along with pressure-based calculation for the right/left wave speeds.

Figure 21. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the first estimation of Davis for the right/left wave speeds.

Figure 21. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the first estimation of Davis for the right/left wave speeds.

Figure 22. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the second estimation of Davis for the right/left wave speeds.

Figure 22. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the second estimation of Davis for the right/left wave speeds.

Figure 23. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the eigenvalues of Reo for the right/left wave speeds.

Figure 23. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the eigenvalues of Reo for the right/left wave speeds.

Figure 24. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the eigenvalues of Einfelt for the right/left wave speeds.

Figure 24. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with the eigenvalues of Einfelt for the right/left wave speeds.

Figure 25. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with pressure-based calculation for the right/left wave speeds.

Figure 25. Calculated final primitive variables and Mach number, using face flux of the improved version of Harten, Lax, van Leer (HLLC), along with pressure-based calculation for the right/left wave speeds.

Figure 26. Calculated final primitive variables and Mach number, using face flux of Lax and Friedrichs (LF).

Figure 26. Calculated final primitive variables and Mach number, using face flux of Lax and Friedrichs (LF).

Figure 27. Calculated final primitive variables and Mach number, using face flux of Rusanov.

Figure 27. Calculated final primitive variables and Mach number, using face flux of Rusanov.

8. Accuracy and Computational Efficiency

To quantity the closeness of each numerical solution to the exact one, the root-mean-square error (RMSE) for each primitive variable in addition to the Mach number were computed. Let $\left(v\right)$ denote any of these four quantities. At each cell-center, let the exact value be $\left(v_{ex}\right)$ while the obtained one by numerically solving Euler equations using a certain flux method be $\left(v_{num}\right)$. Then the RMSE for the generic variable $\left(v\right)$ is

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}\left(\mathrm{v}\right)=\sqrt{{\displaystyle \frac{\sum _{i=1}^n{\left(v_{num,i}-v_{ex,i}\right)}^2}{n}} }$$

(39)

where $\left(i\right)$ is an index for the cell-centers, and $\left(n\right)$ is the total number of cell-centers (which is 200 in this work).

The full list of these errors is given in Table 3. The flux methods are listed in the table in the same sequence of introducing them and presenting their graphical results earlier. The velocity seems to be the most challenging variable to predict, because it always has the highest RMSE value among the four quantities regardless of the flux method. Then comes the Mach number, which is based on the velocity.

The maximum magnitudes of the four quantities (density, gas velocity, pressure, and Mach number) in the exact solution are of the same order of magnitude. In addition, their individual RMSE values are also of the same order of magnitude. It may be then helpful to lump their RMSE values into a one measure of inaccuracy by simple addition, leading to a single aggregate RMSE for each flux method. This is done in Table 4. The methods are ordered ascendingly by that introduced aggregate RMSE. The method of AUFS comes first (aggregate RMSE = 0.05830). For all flux methods other than AUSM⁺-up, AUSM, SW and LF, the aggregate RMSE increases gradually with small increments from a method to the next better one up to AUSM⁺ (aggregate RMSE = 0.06651). This is followed by a mild increase with the two methods AUSM⁺-up (aggregate RMSE = 0.07285) and AUSM (aggregate RMSE = 0.07499). Then, as the end of the list, the aggregate RMSE jumps by a factor of about 4 for the SW and LF methods, thereby clearly delineating their solutions as improper.

While the accuracy of a flux method is very important when analyzing it, its computational demand is worthy of consideration. In large three-dimensional finite volume simulations, face flux construction is a process that happens millions of times to advance the solution by just one time step, and the solution may need thousands of time steps to complete. Even an extremely small time saving in the calculation of the flux vector at a face can thus have a useful impact when considering the total runtime of a realistic simulation. The elapsed clock time of the simulation excluding the processes of initializing variables, calculating an exact solution (as a reference to compare with), plotting results, or writing output data files was recorded for each flux method, and the results are presented in Table 5 in an ascending order. It is worthy of mentioning here that the relative time (one method relative to another) is more meaningful than the absolute time spent. The computing machine used in this work is a laptop Dell Inspiron (model N5010), with 6 GB of installed RAM memory and an Intel^® Core^TM i5 CPU, M 480 @ 2.67GHz. The SW method took the shortest time of 41.13235 s, while the Riemann exact solver took the longest time of 46.49466 s. The ratio of these two extrema values is 1.13037.

A percentage increase of the runtime (using the shortest recorded one as a reference) is included in Table 5. This relative measure of computational cost is calculated as

$$\mathrm{\%}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}={\displaystyle \frac{\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\left(\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{h}\mathrm{o}\mathrm{d}\right)-\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}\times 100\mathrm{\%}$$

(40)

The changes in the elapsed time are not large, with an extra time penalty not exceeding 13% when using the exact Riemann solver flux method (as the one found to be the most computationally demanding) compared to any other method. This can be attributed to the inherent simplicity of the problem.

9. Conclusions

Face flux construction is a core task in a finite volume solver. An end user of such a solver may be offered a large number of flux methods to select from for numerically handling convection (divergence) terms (OpenCFD, 2020). A good background in flux construction methods as well as an elementary analysis of their performance for a simple problem having a known analytical solution may help in making a good selection of a flux method in relevant simulations.

In this context, this work considered 22 flux methods and analyzed them while solving the Sod problem (a classical one-dimensional jump problem for an inviscid flow). Other than changing the flux method as the only discrete variable to examine its influence, the analysis is conducted at a set of unified conditions, where a fully-extrapolated version of the MUSCL scheme was used in which the classical van Leer flux limiter was employed. The one-dimensional grid was uniformly-spaced 200 cells, and the maximum CFL number was kept around 0.44. Qualitative and quantitative assessment of the accuracy of the resulting solution as well its estimated computational cost suggest that there is a number of methods that can offer similar levels of a good performance, including the methods of Kurganov-Noelle-Petrova, van Leer, Harten-Lax-van Leer-Davis2, and Harten-Lax-van Leer-Contact-Davis2. On the other hand, two flux methods may be avoided due to high oscillation or excessive dissipation.