Abstract
This paper presents a new flux splitting scheme for the Euler equations. The proposed scheme termed TV-HLL is obtained by following the Toro-Vazquez splitting (Toro and Vázquez-Cendón, 2012) and using the HLL algorithm with modified wave speeds for the pressure system. Here, the intercell velocity for the advection system is taken as the arithmetic mean. The resulting scheme is more accurate when compared to the Toro-Vazquez schemes and also enjoys the property of recognition of contact discontinuities and shear waves. Accuracy, efficiency, and other essential features of the proposed scheme are evaluated by analyzing shock propagation behaviours for both the steady and unsteady compressible flows. The accuracy of the scheme is shown in 1D test cases designed by Toro.
1. Introduction
Today, upwind schemes undoubtedly have become the main spatial discretization techniques used for solving the Euler/Navier-Stokes equations. Indeed, an interesting feature of these schemes is that the discretization of the equations on a mesh is performed according to the direction of propagation of information on that mesh. They are categorized [1] into two major family schemes, namely, flux vector splitting (FVS) and flux difference splitting (FDS).
Concerning the flux difference splitting schemes (FDS), they are based on the difference in the decomposition of fluxes, constructed on an approximated solution of the local Riemann problem between two adjacent states. There are several FDS schemes formulated in the literature among which are the Roe and HLL numerical schemes. The HLL scheme developed by Harten-Lax-Van Leer is a direct approximation of the numerical flux to compute Godunov flux [2]. It has been shown that the HLL scheme is very efficient and robust, has an entropy satisfaction property, resolves isolated shock exactly, and preserves positivity [2–4]. Unfortunately, the main demerit of the HLL scheme is that it cannot resolve contact discontinuity exactly. Another well-known FDS scheme is Roe’s scheme [5] largely used because of its accuracy, quality, and mathematical clarity. Unfortunately, Roe’s scheme admits rarefaction shocks that do not satisfy the entropy condition and sometimes gives rise to spurious solutions such as carbuncle phenomena and odd-even decoupling [6]. To improve the quality of solutions, Quirk proposed a strategy to use combined fluxes so that a dissipative approach can be used in the shock regions [7].
In contrast, the flux vector splitting (FVS) has proven to be a simple and useful technique for arriving at upwind differencing and is preeminently suited for use in implicit schemes. They are considered as the first level of upwind schemes. Their approach achieves upwinding by splitting the flux vector into two components that are the positive and the negative parts according to the sign of the eigenvalues of the coefficient matrix. The identification of the upwind directions is done with less effort than in the Godunov-type methods. The upwinding is also performed for multidimensional flows. Unfortunately, the simplicity of these splittings comes at the price of reduced accuracy due to numerical diffusion. Many FVS schemes have also been formulated in the literature such as Steger-Warming and Van Leer schemes. Concerning the Steger-Warming [8] and Anderson et al. [9] schemes, it has been shown that they cannot resolve intermediate characteristic fields.
In order to combine the advantages presented by the above family schemes (FDS and FVS), Liou and Steffen [10] developed a hybrid flux splitting method that is able to recognize contact waves. This scheme, termed AUSM (advection upstream splitting method), has received much attention by the computational community. It consists of two steps: the first step is to recognize that the inviscid flux vector consists of two physical parts, namely, convective and pressure terms, and the next step is to discretize the two components separately. This scheme is able to capture shock waves and contact discontinuities. However, its drawbacks also surface. To overcome the flaws of the AUSM scheme, Liou and Wada [11] proposed the AUSMDV scheme. It eliminates numerical overshoots behind shock waves shown in AUSM though not completely. Unfortunately, it allows carbuncle phenomena due to numerical instability. The above disadvantage, presented by AUSMDV scheme, motivated Liou [12] to propose the AUSM+ scheme that features the following properties: exact resolution of 1D contact and shock discontinuities, positivity preserving of scalar quantity such as the density, free of “carbuncle phenomenon,” algorithm simplicity, and easy extension to treat other hyperbolic systems. Unfortunately, it shows numerical overshoots behind strong shock waves and numerical oscillations near regions of small convection velocity or small pressure gradient, such as around a wedge or near a wall. Hence, Kim et al. [13] proposed the AUSMPW scheme that features the following: no carbuncle phenomena, elimination of overshoots, high resolution of discontinuities, reduced numerical dissipation, constancy of total enthalpy, efficiency, and robustness. In addition, Kim et al. [14] proposed the AUSMPW+ to increase the accuracy and computational efficiency of AUSMPW in capturing an oblique shock without compromising robustness. The extension of AUSM schemes family to all speed regimes has been introduced by Liou [15] and the resulting scheme is called AUSM+-up. Similarly, Edwards [16] developed a low-diffusion flux-splitting for Navier-Stokes calculations. Later, the interest for the multiphase flows is developed by Edwards et al. [17] and Ihm and Kim [18]. It is well known that numerical fluxes for conservative methods for solving hyperbolic equations are expected to be not just robust but also able to resolve intermediate characteristic fields accurately. This requirement is not met by centred fluxes, nor by incomplete Riemann solvers or classical flux vector splitting methods. The inability to resolve intermediate characteristic fields badly affects the correct resolution of contact waves, material interfaces, shear waves, vortices, and ignition fronts. Moreover, when extending these schemes to solve viscous problems, this shortcoming manifests itself when computing shear layers. It is in this vein that Toro and Vázquez-Cendón proposed the TV and TV-AWS schemes [19] that are simple, robust, and accurate when compared with existing methods and also enjoy a very desirable property: recognition of contact and shear waves in general and exact preservation of isolated stationary contacts.
In this paper, an improved version of the TV schemes termed TV-HLL is presented. The proposed scheme is more accurate than the TV schemes and also recognizes contact discontinuities and shear waves.
The paper is organized as follows. Section 2 gives the preliminaries related to this paper. In Section 3, the TV and HLL schemes are presented. In Section 4, the TV-HLL scheme is developed. Section 5 is devoted to some 1D numerical tests to assess both the robustness and accuracy of the proposed scheme.
2. Preliminaries
Let us consider the one-dimensional system of conservation laws for ideal-gas flows given by where and are, respectively, vectors of conservative quantities, and fluxes given by where is the density, is the velocity component in the -direction, and is the static pressure. is the total specific energy defined by , is the internal energy, and the ideal-gas equation of state is assumed to be with being the ratio of specific heats.
The finite volume method applied on (1) leads to where , , and . and are, respectively, the time step and the grid size. is the numerical flux vector defined by two neighbouring cells of left and right cells:
In the finite-volume formulation (3), the differences among all the numerical schemes lie essentially in the definition of the numerical flux evaluated at the cell interface. The performance of the numerical flux will vary, depending on what set of properties it must meet.
3. The TV and HLL Schemes
3.1. The TV Scheme [19]
The first step in designing the TV schemes [19] is to split the flux vector (2) as follows: with the corresponding advection and pressure fluxes defined as
For an ideal gas, the pressure flux (6) becomes
The system (1) is replaced by two subsystems (the advection and the pressure systems, respectively; see [19] for more details) that are
The objective is to compute the numerical flux as
In order to compute the numerical fluxes and , Toro and Vázquez-Cendón [19] proposed a careful study of (8) and (9) separately.
3.1.1. The Advection System [19]
The system (8) can be written in a quasilinear form as with and being the Jacobian matrix given by
The analysis of matrix shows that the system is weakly hyperbolic as there is no complete set of linearly independent eigenvectors. Indeed, its eigenvectors are and the corresponding eigenvalues are
3.1.2. The Pressure System [19]
The system (9) can be written in a quasilinear form as with being the Jacobian matrix given by
The analysis of matrix shows that it has three real eigenvalues given by with the corresponding eigenvectors given by
3.1.3. Numerical Fluxes
Concerning the advection system, Toro and Vázquez-Cendón [19] first proposed (TV scheme) to evaluate the numerical flux as where is the intercell velocity given by with and .
Secondly, they proposed [19] (TV-AWS scheme) to express the advection numerical flux by with where
Concerning the pressure system, they proposed [19] to evaluate the pressure numerical flux as with given by (20) and given by
Therefore, the numerical flux is given by (10) with defined by (19) or (21) and by (24).
3.2. The HLL Scheme [3]
Harten et al. [3] proposed to consider the Riemann problem with two waves where only nonlinear waves are used in the computation (see Figure 1).
Considering (1), the numerical flux is given by where and are estimates of wave speeds [20], given by with and , respectively, being the left and right sound speeds defined by
4. The TV-HLL Scheme
The basic idea in designing the proposed scheme, TV-HLL, is to use the HLL algorithm to evaluate the numerical pressure flux with modified wave speeds and the arithmetic mean to define the intercell velocity for the advection numerical flux. Concerning the pressure system, the central idea is to assume a wave configuration for the solution that consists of two waves separating three constants’ states. Hence, the new numerical pressure flux is given by where is given by (7).
The most well-known approach for estimating bounds for the minimum and maximum signal velocities present in the solution of the Riemann problem is to provide, directly, and . Our suggestion, following the Davis idea for the wave speeds estimates [20] and the expression of the eigenvalues (17) of the Jacobian matrix for the pressure system, is with and given by (27), respectively, and and defined as Here, the quantities and are given by
The quantities appearing in (31) represent the maximum and the minimum eigenvalues of the Jacobian matrix of the pressure system. and are the sound speeds given by (28).
Concerning the advection numerical flux, we propose
For the definition of the intercell velocity, many choices are possible. One of the choices is the one suggested in [19], but our suggestion is to use a simple choice which takes into account the two adjacent velocities. It is given by
Therefore, the final numerical flux proposed in this paper is given by (10) with defined by (33) and by (29).
5. Numerical Results and Discussion
The performance of the TV-HLL scheme is illustrated in the 1D Euler equations for ideal gas with .
5.1. Results for 1D Problem Tests
Six tests have been solved, which are generally tested by Toro [21]. The tests are selected to investigate the performance of TV-HLL solver. In all chosen tests, data consists of two constant states, and , separated by a discontinuity at a position . The states are given in Table 1. The exact Riemann and numerical solutions are found in the spatial domain using 100 cells. The Courant number coefficient is taken as 0.9. For each test problem, an initial location of discontinuity, , and the output time are selected; these are stated in the caption of each figure. Boundary conditions are transmissive. For these 1D problem tests, the numerical results are compared with those of five other schemes (exact Riemann, TV, TV-AWS, TV-HLL, and Roe).
5.1.1. Results for Test 1
Test 1, which is known as a modified version of Sod’s test, has solutions consisting of a right shock wave, a right travelling contact wave, and a left expansion wave with a sonic point inside. This test is devised to assess the entropy satisfaction property of numerical methods. Figures 2 and 3 show solution profile for the density, velocity, pressure, and internal energy at the output time units. The TV, TV-AWS, TV-HLL, Roe, and exact Riemann solvers are used here.
(a)
(b)
(a)
(b)
Discussion for Test 1. As expected, Roe’s solver exhibits an expansion shock near the sonic point because an entropy fix is not used. It can be noticed that the TV-HLL scheme resolves the sonic point as well as the TV and TV-AWS schemes. In the results from the TV, TV-AWS, and TV-HLL schemes, the resolution of the shock and the right travelling contact is comparable with that of exact Riemann scheme. Concerning the left rarefaction, the proposed scheme (TV-HLL) seems to be the most accurate.
5.1.2. Results for Test 2
Test 2 has solutions consisting of two symmetric expansion waves and a trivial contact wave of zero speed; the star region between the nonlinear waves is close to vacuum, which makes this problem a suitable test for assessing the performance of numerical methods for low-density flows. The results of this test are shown in Figures 4 and 5.
(a)
(b)
(a)
(b)
Discussion for Test 2. It can be noticed that Roe’s scheme fails in this test. The TV-HLL scheme as well as TV scheme performs better near the minimum for density. The pressure and velocity distributions from TV, TV-AWS, and TV-HLL are virtually identical. From the internal energy plot (see Figure 5), it can be observed that the TV-HLL scheme is less dissipative than TV and TV-AWS schemes. The TV-HLL scheme gives more accurate results than those obtained with TV and TV-AWS schemes; in the vicinity of the trivial contact, where both density and pressure are close to zero, the results are somewhat erroneous; see internal energy plots. Most of the approximated Riemann solvers present a high dissipation for this test case.
5.1.3. Results for Test 3
Test 3 is designed to assess the robustness and accuracy of numerical methods; its solution consists of a strong running shock wave of shock Mach number 198, a contact surface, and a left rarefaction wave. Figures 6 and 7 show the results for five numerical schemes.
(a)
(b)
(a)
(b)
Discussion for Test 3. Compared to the exact Riemann solution, density levels remain relatively low behind the shock wave, as it can be easily seen in the density plots. From the velocity distribution, it can be seen that the TV-HLL captures better the left rarefaction wave when compared with the other schemes.
5.1.4. Results for Test 4
Test 4, as test 3, is also designed to assess the robustness of numerical methods; its solution consists of three strong discontinuities travelling to the right. It is an interaction of two strong impinging shock waves. The left shock wave moves to the right very slowly, which adds another difficulty to numerical methods.
Discussion for Test 4. The results from the TV-HLL scheme, shown in Figures 8 and 9, are overall comparable with those obtained with other schemes. The four schemes behave similarly well, with a slight improvement in results obtained with TV-HLL scheme. The proposed scheme resolves the slowly moving left shock slightly more sharply than the other schemes.
(a)
(b)
(a)
(b)
5.1.5. Results for Test 5
Test 5 is designed to assess the accuracy of numerical methods to capture stationary contact wave. It concerns the exact recognition of isolated contacts.
Discussion for Test 5. In Figures 10 and 11, all the schemes (TV, TV-AWS, TV-HLL, Roe, and exact Riemann) pass the test successfully. The ability of the TV-HLL scheme to resolve contact discontinuity is more straightforwardly demonstrated here. It is important to notice that many FVS schemes and the HLL scheme do not resolve the contact discontinuity.
(a)
(b)
(a)
(b)
5.1.6. Results for Test 6
Test 6 is designed to assess the accuracy of numerical methods to capture moving contact wave.
Discussion for Test 6. The result is shown in Figure 12 where it can be observed how all the used schemes lead to identical result.
5.1.7. Other 1D Shock-Tube Problems
Two two-material shock-tube problems of two gases modeled by ideal gas EOSs with different and with very different shock strengths are considered. In Case I-A, which is taken from Fedkiw et al. [22], Hu and Khoo [23], and Hu et al. [24], the shock is moderately strong. The initial condition is and the final time is units. In Case I-B, which is taken from Hu and Khoo [23] and Abgrall and Karni [25], the shock is very strong. The initial condition is and the final time is units.
Discussion for Case I-A. From Figures 13 and 14, it can be easily seen that the TV-HLL scheme is more close to the exact Riemann solution when compared to the other schemes (see the pressure and velocity distributions). The TV-AWS scheme is more dissipative than the other schemes (see the internal energy plots).
(a)
(b)
(a)
(b)
Discussion for Case I-B. From Figures 15 and 16, it can be observed that the solutions presented by all the schemes are overall comparable. There is a slight improvement with the TV-HLL scheme (see the velocity and pressure distributions).
(a)
(b)
(a)
(b)
6. Conclusions
The new scheme termed TV-HLL provided more accurate and efficient results. The idea behind the new scheme consisted firstly of a new splitting of the flux vector [19], followed by the HLL algorithm [3] applied on the pressure system and the use of modified wave speeds estimates. It has been evaluated on several well-known test cases and found to be more robust and accurate when compared to the TV, TV-AWS, and Roe schemes.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The first author acknowledges Professor Eleuterio Toro for his greatest contribution in the computational physics field.