Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 7451240, 21 pages
https://doi.org/10.1155/2020/7451240

Research Article

On Asymptotic Behavior of HLL-Type Schemes at Low Mach Numbers

Hang Yu, Zhengyu Tian, Fan Yang, and Hua Li

College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, Hunan, China

Correspondence should be addressed to Zhengyu Tian; yuhang18@nudt.edu.cn

Received 18 May 2020; Accepted 24 July 2020; Published 17 August 2020

Academic Editor: Luca Chiapponi

Copyright © 2020 Hang Yu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The HLLEM approximate Riemann solver can capture discontinuities sharply, maintain positive definiteness, and satisfy the entropy condition automatically. These attractive properties make the HLLEM scheme widely used in the simulations of many compressible fluid problems. However, in the simulations of low Mach incompressible flow, the accuracy of HLLEM solver cannot be guaranteed. In the current study, a detailed discrete asymptotic analysis is conducted on the HLLEM scheme and the responsible terms for the loss of accuracy are identified. This allows us to develop two modified methods to solve this low Mach number problem. The first method is to add a low Mach number correction term on the basis of the original HLLEM scheme. The second is to simply rescale the responsible terms with a Mach number-based function. The asymptotic analysis of these two low Mach correction methods shows that the difference between the continuous system and the discrete system disappears, which means the resulting LM-HLLEM and LM-HLLEM2 schemes are both capable of obtaining physically correct solutions in low Mach limit. The results obtained from various test cases demonstrate that both these two HLLEM-type schemes can simulate incompressible and compressible fluid problems accurately and robustly.

1. Introduction

Due to the robustness and clear physical interpretation, Godunov-type Riemann solvers have been widely applied in practical engineering problem calculation and theoretical research. But most Godunov-type approximate Riemann solvers cannot ensure that the simulation results are always physically correct for incompressible or weakly compressible flows [1]. In practice, most of the problems not only include compressible flow regions but also incompressible flow areas where the fluid flows at very low speeds. Schemes developed to capture shock waves usually encounter accuracy problems in the simulations of incompressible flows. Therefore, it is urgently necessary to improve the performance of Godunov-type schemes in low Mach incompressible flow simulations.

Much work has been done by a large number of scholars to try to solve this accuracy problem of Godunov-type Riemann solvers under the low Mach limit. Guillard et al. [2, 3] conducted an asymptotic analysis of the Godunov-type schemes at low Mach numbers. They found the results calculated by the numerical schemes contained pressure fluctuations of order , while the actual physical pressure varied with . They applied a preconditioned technique to rescale the Roe scheme [4] and thus modified the numerical dissipation to solve this problem. Rieper [5] used a Mach number-based function to rescale the velocity jump perpendicular to the cell interface and developed a low Mach modified Roe scheme called the LMRoe scheme. Oßwald et al. [6] further improved the LMRoe scheme which exhibits a reduction in numerical dissipation when the Mach number is low. The preconditioning method proposed by Turkel [7] also improves the accuracy and efficiency in the calculations of low Mach flow cases [8]. However, a large number of test cases showed that when low Mach flow and high Mach flow both exist in the simulated region, the cutoff problem will inevitably occur with the preconditioning method. This problem reduces the accuracy of this method in the simulations of the transition flow region and the low/high-speed mixed flow regions. Li et al. [9–11] proposed an All-Speed-Roe scheme to solve this cutoff problem. In addition, an asymptotic analysis was conducted by them to demonstrate that the asymptotic behavior of this new scheme is the same as that of the compressible governing equations at the low Mach limit. Chen et al. [12] split the jump of the left and right states into the density diffusion part and the velocity diffusion part to improve the accuracy and efficiency of upwind schemes. By multiplying a scaling function on the velocity diffusion part, the upwind schemes can be extended to the calculations of low Mach number incompressible flows. Dellacherie et al. [13, 14] have also done a series of research in order to solve this low Mach limit problem and put forward a corresponding correction method. By applying their low Mach correction, a new low Mach Godunov scheme for solving linear wave equations was developed, and they further generalized it to nonlinear Euler cases. Their further research showed that this method can be used in other numerical schemes and can also be used in the verification of other low Mach number schemes, like the AUSM-type methods [15] and the LMRoe scheme [5]. Based on the work of Dellacherie et al. [13, 14], Xie et al. [16, 17] managed to develop an all-speed HLLC-type Riemann solver and an all-speed Roe-type Riemann solver, which are not only endowed with a high resistance against shock anomalies but also provide accurate solutions in low Mach number systems. Based on the simple and robust AUSM family schemes [15, 18], Shima and Kitamura developed all-speed flux functions, namely, SLAU [19], SLAU2 [20], and SD-SLAU [21]. These all-speed methods not only exhibits high robustness for capturing shock waves but also can be applied to various low Mach number flows.

The HLL scheme was developed by Harten et al. [22] in 1983, which is a unified framework for constructing approximate Riemann solutions. Using this method, the approximate Riemann solutions at different resolutions for various discontinuities can be obtained easily. As indicated in [23], under certain wavespeed estimates, the HLL-type schemes satisfy the entropy condition automatically and can maintain positive definiteness. By means of a linear distribution method, Einfeldt [23] modified the intermediate states of the Riemann problem and obtained a HLL-type approximate Riemann solver, namely, the HLLEM scheme. This scheme is a complete Riemann solver, so it can distinguish both linear and nonlinear waves. In addition, the HLLEM scheme is able to resolve linear waves at the cost of minimal diffusion. This scheme can also maintain positive definiteness and satisfy the entropy condition automatically, like the HLLC scheme [24]. Due to these attractive features, many scholars have conducted studies on the HLLEM scheme. This scheme has been widely extended to various applications like magnetohydrodynamics [25], multiphase flows [26], and chemically reacting flows [27]. It is worth noting that Dumbser and Balsara [28] recently developed a simple and universal formulation of the HLLEM scheme which performs well in both general conservative and nonconservative hyperbolic systems, namely, the HLLI scheme. This new HLLI scheme inherits many desirable properties of the HLL Riemann solver and has wider applicability than the original HLLEM scheme. Its desirable features have been appreciated, and it was shown to be feasible to compute shallow water-type equations and two-phase flow models, as well as for gas dynamics with real equation of state, magnetohydrodynamics, and nonlinear elasticity [28]. Based on this HLLI scheme, the hybridized SLAU2-HLLI and AUSMPW+HLLI [29] have been proposed to solve the magnetohydrodynamics problems. There are also improved works that concern the shock instability problem of the HLLEM scheme, which can be referred in [30, 31]. However, like many other Godunov-type schemes, the HLLEM scheme also suffers from the accuracy problem in the simulations of low Mach incompressible flows. This low Mach number problem limits the performance of the HLLEM scheme to accurately simulate all-speed flows.

Some scholars have attempted to improve the low Mach performance of the HLLEM scheme. Park et al. [32] proposed a preconditioned HLLEM method to solve low Mach viscous flows, and they introduced additional limitations and modified global cutoff terms to enhance the robustness of the preconditioning system. Through a correction item added to the original HHLEM scheme, Yu et al. [33] improved the simulation accuracy of the low Mach problem. By adopting a controlling function to the momentum equations at low speeds, Qu et al. [34] developed the so-called HLLEMS-AS scheme which is capable of obtaining physically correct solutions in incompressible flows. However, most of the current analyses of HLLEM-type schemes for the low Mach problem are not so complete. The reason why the HLLEM scheme fails to be accurate at low Mach number flows needs to be further clarified. And a proper modification based on the correct analysis of HLLEM scheme at low Mach incompressible flows is still rare.

In this paper, a detailed asymptotic analysis of HLLEM scheme is conducted to study its low Mach number behavior. Through this analysis, the terms which are responsible for resulting in the loss of accuracy of HLLEM solver at low Mach limit are identified. Based on the results of this analysis, two low Mach number correction methods are put forward to rescale the responsible terms, resulting in the LM-HLLEM scheme and the LM-HLLEM2 scheme. Compared with the existing fixes for the HLLEM scheme in low Mach number situations like Park et al.’s preconditioned HLLEM method [32] and Qu et al.’s HLLEMS-AS scheme [34], our newly proposed schemes have much simpler formulas and only the actual responsible terms for the loss of accuracy are rescaled. In this sense, the proposed low Mach HLLEM-type schemes are effective for use as reliable tools for all-speed flow computations.

This paper is organized in the following way. In Section 2, the essential results of the asymptotic analysis on compressible Euler equations are briefly reviewed. In Section 3, the HLLEM scheme and its discrete asymptotic analysis are introduced. Then, two low Mach number fix methods along with their asymptotic analysis are discussed in Section 4. Following on in Section 5, we will examine the performance of these two HLLEM-type schemes with several numerical tests. Finally, the conclusion of this paper will be given in Section 6.

2. Low Mach Number Behavior of Compressible Euler Equations

Following Guillard and Viozat [2], the asymptotic analysis of compressible Euler equations is briefly reviewed in this section.

The dimensionless compressible Euler equations can be written as

where the equations have been scaled with the following reference quantities:

where is the density, denotes the flow velocity, represents the pressure, and is the speed of sound. The superscript “” represents the reference quantities of each variable, and is an arbitrary length scale. is the reference Mach number.

As in [5, 34], we expand the physical quantities , , and into powers of the reference Mach number to find the asymptotic solutions of system (1):

Inserting equation (3) into equation (1) and sorting the dimensionless equations according to the powers of , an enlarged system of equations is obtained. Only the main results are summarized below. The detailed asymptotic analysis of this set of equations can be found in [35, 36].

Property 1. In the low Mach number limit, the physical pressure scales with the square of the Mach number , and therefore, the pressure can be written as

which means the leading pressure and first-order pressure are constant in space.

Property 2. The divergence constraint of zero order of flow velocity is zero, and thus

Property 3. The second-order pressure is the balance force in the incompressible flow field and satisfies a Poisson equation of the following form:

3. Detailed Asymptotic Analysis of the HLLEM Scheme

3.1. The Expression of the HLLEM Scheme

Einfeldt [23] proposed the HLLEM scheme in 1988 which is now a widely used approximate Riemann solver. The interface flux of this scheme can be expressed as

with

where the vectors and represent the inviscid flux and conservative variable, respectively. They can be defined as

where is the specific total energy, represents the static pressure, is the density, and and are velocity vector components. denotes the contravariant velocity, and represents an outward unit vector perpendicular to the interface.

and denote the left and right wave speeds, respectively, and can be calculated from

where the symbol “” represents Roe’s averaged values at the interface and is the speed of sound.

where and denote the 2nd and 3rd right eigenvectors of Jacobian matrix, and are the wave strengths, and and represent the antidiffusion terms, which control the dissipation terms of linearly degenerate waves. One should note that the definition in equation (12) follows the improvement by Park and Kwon [37] instead of its original version in [23]. Such a definition enhances the resolution of the HLLEM scheme for contact discontinuities.

3.2. Discrete Asymptotic Analysis of the HLLEM Scheme

Guillard and Viozat [2] first introduced progressive analysis to analyze low-Mach performance in numerical schemes. By using asymptotic analysis, they proved that the solutions of discrete system contain pressure fluctuations of order . These results are in contrast with those of the compressible Euler equations. Their analysis explains why the approximations of the compressible fluid flow equations cannot converge to the solutions that approach the incompressible flows. In [5], Rieper offered a more detailed and coherent asymptotic analysis of the Roe scheme. Here, we follow the work of Rieper [5] to conduct an asymptotic analysis on the HLLEM scheme to study its low Mach behavior. The nomenclatures defined in this paper can be found in [5].

 : index vector for cell of reference, . : index set for neighboring cells, . : index vector for neighboring cells, . : area of the reference cell. : index for edge between cell and cell . : length of the interface . : unit outer normal vector from cell to . : unit transverse vector from cell to . : difference operator. : Roe’s averaged values of and . : velocity with Cartesian (global) coordinates. : normal component of . : transverse component of .

Our analysis is conducted on a uniform Cartesian grid to reduce the complexity. So, we set and .

With the Mach number approaching zero, equation (11) and (13) can be reduced to the following forms:

The expressions of the 2nd and 3rd right eigenvectors are transformed into

and the corresponding wave strengths can be expressed as

For two-dimensional cases, the HLLEM scheme can be written as

with the HLLEM flux function

Complete equations for density transport, momentum density ( and ) transport, and the energy density transport can be expressed as follows:

The same reference quantities in [5] are used here to scale the equations (19)–(22) and sort the nondimensional equations according to the order of .

The continuous equations in dimensionless form can be expressed as

The dimensionless equations of horizontal and vertical momentum can be obtained by steps similar to the above:

The nondimensional energy conservation equation can be written as

All physical quantities are assumed to satisfy the following asymptotic three-term expansion:

Then, insert the asymptotic expansions of all physical quantities into the semidiscrete equations (23)–(26) and sort the equations according to the power of .

Order :

Order :

Order :

Equations (28) and (29) imply that has a four-field checkerboard solution as depicted in Figure 1.

Figure 1: Checkboard structure.

To prove that the checkerboard modes can be damped, we express the pressure modes as , , , and . Assuming the checkerboard does not decrease in time, is the maximum of the four values of . Then, equation (30) can be expressed as

The coefficient of equation (38) is positive, and thus equation (38) cannot be equal to zero unless , which indicates that the checkerboard structure will disappear in large time limit. Therefore, the discrete system satisfies the uniformity of , i.e., .

Roe’s average values can be simplified according to and . And taking into account , the order momentum equations (31) and (32) can be reduced to

Equations (39) and (40) are both not equal to zero, which implies the first-order pressure is not constant. Therefore, the solutions of discrete Euler equations (19)–(22) support pressure fluctuations of order , namely,

Equation (41) shows that the low-Mach asymptotic behavior of the HLLEM scheme is inconsistent with that of the continuous system (Property 1). So, the discrete terms (horizontal direction) and (vertical direction) are responsible for the loss of accuracy at low Mach numbers. Since Property 1 cannot be satisfied, Properties 2 and 3 also cannot be satisfied naturally.

4. Two Low Mach Fixes for HLLEM Scheme

Through the discrete asymptotic analysis of the HLLEM scheme in the previous section, we have determined which responsible terms make the HLLEM scheme unable to converge to a physically correct solution when the Mach number approaches zero. These two responsible terms are (horizontal direction) and (vertical direction). In this section, two simple methods of correction are proposed to remove these responsible terms, resulting in two new HLLEM-type schemes. Asymptotic analysis of these two schemes shows that the modified first-order pressure can be kept constant, which is consistent with Property 1. Based on this result, it can be proved that Properties 2 and 3 can be satisfied.

4.1. The LM-HLLEM Scheme and Its Asymptotic Analysis

Inspired by Dellacherie et al. [13, 14] and following on from our previous work [33], the first low Mach number fix is to simply add a correction term on the original HLLEM scheme to rescale the responsible terms.

4.1.1. The LM-HLLEM Scheme

The correction term can be written as

where is a scaling function and can be calculated from the following equation:

where and are Roe’s averaged values of and , denoting the speed of sound and density at the grid interface. is the local Mach number.

According to the low Mach fix in equation (42), the LM-HLLEM scheme in two-dimensional space can be written as

The local Mach number can be calculated from the following equation, which is introduced by Thornber et al. [38]:

4.1.2. The Asymptotic Analysis of the LM-HLLEM Scheme

In this section, an asymptotic analysis will be conducted on the LM-HLLEM scheme to complete the construction of this scheme. There is no need to repeat each step of the asymptotic analysis, so we will merely concentrate on the parts where the low Mach number fix works.

By applying the low Mach correction, the first-order pressure in equations (31) and (32) becomes

As the Mach number approaches zero, the weight function becomes zero, so we can obtain

According to equation (47), the first-order pressure can satisfy the following relations:

4.2. The LM-HLLEM2 Scheme and Its Asymptotic Analysis

The second low Mach correction method is multiplying the responsible terms with the scaling function . Compared with the first low Mach number fix in previous section, this approach is simpler and more straightforward.

4.2.1. The LM-HLLEM2 Scheme

Taking the first-order pressure in the horizontal direction as an example, i.e., equation (39), the responsible term is where the first term is derived from the difference of . The second term results from the third right eigenvector of the Jacobian matrix . Therefore, we propose the following modifications for conservative variables and the right eigenvector :

where is defined in equation (43).

The modified numerical flux can be written as

In [34], Qu et al. applied a controlling function not only to the conservative variables and third right eigenvectors but also to the second right eigenvectors . The resulting HLLEMS-AS scheme is capable of obtaining physically correct solutions in incompressible flows. However, the asymptotic analysis of the HLLEM scheme in Section 3 showed that has no effect on its low Mach behavior. Therefore, there is no need to modify this term. Compared with Qu et al.’s HLLEMS-AS scheme, this LM-HLLEM2 scheme is more concise and preserves the original form of the HLLEM scheme as much as possible.

4.2.2. The Asymptotic Analysis of LM-HLLEM2 Scheme

Here, we also conduct an asymptotic analysis on the LM-HLLEM2 scheme.

By applying the second low Mach fix, the first-order pressure in equations (31) and (32) becomes

As the Mach number approaches zero, the scaling function becomes zero, so we can obtain the same results as described in equation (48), namely,

4.3. Asymptotic Analysis of Discrete System

The two low Mach number correction methods proposed in Sections 4.1 and 4.2 result in the same conclusion, namely, the first order of pressure can be remedied to zero. Based on this result, the remaining parts of the asymptotic analysis on the LM-HLLEM and LM-HLLEM2 schemes can be completed.

Equations (48) and (52) indicate that the first-order pressure of the cells adjacent to the reference cell is equal. Thus, the discrete system satisfies the uniformity of , i.e., , which indicates that the discrete solutions of these two schemes both satisfy the pressure fluctuations of order :

The numerical results in Section 5.4 support this analytical result.

According to the following equation, the energy density in equation (37) can be replaced by the pressure . Through this replacement, we can obtain the evaluation equations of :

Roe’s average values can be simplified according to and . And taking into account , equation (37) can be transformed into

Following the simplification in [16], namely, assuming that there is no net influx over the boundary due to in the ghost cells and there is no gradient in across the boundary, equation (55) can be further transformed into

By applying the fact that , equation (56) can be simplified into

and therefore, the discrete divergence constraint is satisfied:

where is the discrete divergence operator. So, Property 2 is satisfied in these two modified discrete systems.

Assuming that and are constant, the momentum equations in equations (35) and (36) can be simplified. Taking equation (35) as an example, this equation can be transformed into the following forms, for the LM-HLLEM scheme:

For the LM-HLLEM2 scheme:

According to the above analysis, the discrete divergence constraint is shown to be satisfied and is constant. Considering these facts, equations (59) and (60) can be transformed into a Poisson equation of the following form:

with denoting the Laplacian operator in discrete form. Equation (61) means these two discrete systems satisfy Property 3.

In this section, two low Mach number fix methods are introduced. As a result, two modified HLLEM-type schemes are developed, namely, the LM-HLLEM scheme and the LM-HLLEM2 scheme. The asymptotic analysis of these two new schemes demonstrates that the discrete systems can satisfy Properties 1–3. This means the differences between continuous systems and discrete systems are eliminated. Therefore, the LM-HLLEM scheme and the LM-HLLEM2 scheme are both capable of reaching physically correct solutions at low Mach number.

5. Numerical Test Cases

The effectiveness of these two low Mach fixes are verified by six numerical test cases in this section. Two one dimensional test cases are conducted to test the LM-HLLEM and LM-HLLEM2 scheme’s ability of sharply capturing shocks and contact discontinuities. The Gresho vortex test case and inviscid flows around NACA0012 airfoil test case are applied to verify whether these two schemes are capable of obtaining physically correct solutions for incompressible flows. The two turbulence test cases are conducted to examine whether the proposed two HLLEM-type schemes can be applied to the simulations of low Mach number and transonic turbulence flows or not. Moreover, the transonic turbulent case is further used to test the calculation efficiency of these two low Mach correction schemes.

5.1. Modified Sod Shock Tube

This case is derived from the classical Sod shock tube problem [39], which contains a left sonic rarefaction wave, a right travelling contact wave, and a right shock wave. This Sod shock tube test case is used to evaluate if these two newly developed low Mach schemes can capture different flow states as clearly as the original HLLEM scheme and whether these numerical schemes meet the entropy condition. We refer to [40] and set the initial conditions as

Additionally, the position of discontinuity is set at , the computational domain is evenly discretized with 400 points, and the boundary condition of the boundaries is set as “Nonreflecting.”

Figure 2 shows the numerical results of four schemes and the exact results at . As shown in the density and Mach number profiles, both the LM-HLLEM and LM-HLLEM2 schemes can distinguish the contact discontinuity, rarefaction wave, and shock wave clearly like in the original HLLEM scheme. In addition, there is no “rarefaction shock” near the sonic point, indicating that these two new low Mach schemes inherit the characteristic of the HLLEM scheme satisfying the entropy condition. The Roe scheme is used as a comparison because the entropy condition is not satisfied, so an obvious nonphysical rarefaction shock is generated near the sonic point.

(a)

(b)

Figure 2: The profiles of density (a) and Mach number (b) at .

5.2. Double Rarefaction

The double rarefaction test case used here contains a trivial contact wave of zero speed and two rarefaction waves. This test case is used to determine whether the proposed LM-HLLEM and LM-HLLEM2 schemes are positively conservative. The position discontinuity is set at , and the initial conditions can be found in [40], namely,

The boundary conditions are “Transmissive.” The spatial domain is discretized with 400 uniform points in the interval . The first-order solutions computed with three schemes at and their exact results are shown in Figure 3. The results obtained by the LM-HLLEM scheme and the LM-HLLEM2 scheme are basically consistent with those calculated by the HLLEM scheme and agree well with the exact results. In addition, the results indicate that the proposed two low Mach schemes are positively conservative. In contrast, the Roe scheme fails to compute this problem, which means that it cannot preserve the positivity of the solution.

(a)

(b)

Figure 3: The profiles of density (a) and Mach number (b) at .

5.3. Gresho Vortex

The Gresho vortex is a kind of rotating flow that uses pressure gradient to balance the centrifugal force. This flow mode was proposed by Gresho in 1990 [41]. This case is used here to verify the accuracy and dissipation property of the LM-HLLEM and LM-HLLEM2 schemes in incompressible flows simulations.

The simulation conditions including the definitions of initial conditions, initial background pressure, and initial pressure field can be found in [5], which will not be repeated in this paper. The simulation region is set as , which is evenly discretized by quadrilateral grids. The initial position of the vortex is , and the initial radius of the vortex is 0.4. The simulated Mach number is set to 0.01. The boundary conditions of the simulated area are “Farfield (for lower and upper boundaries)” and “Periodic (for right and left boundaries),” respectively.

A dimensionless pressure field is defined as to represent the order of pressure fluctuations under the incompressible limit:

The Gresho vortex is a steady-state solution of the incompressible Euler equation, so the dimensionless pressure field defined above should not change with time in the process of calculating the low Mach flow fields.

Figure 4 shows the initial state of the dimensionless pressure field and the state calculated by four schemes under the condition . After one period, both the LM-HLLEM scheme and the LM-HLLEM2 scheme successfully maintained the vortex structure which agrees well with the initial result, indicating the low dissipation properties of these two schemes in the calculation of low Mach problems. On the contrary, due to the excessive numerical dissipation when the Mach number is low, the vortex structures calculated by the original Roe scheme and HLLEM scheme are completely dissipated after only one cycle. Figure 5 depicts the normalized pressure distribution on the centerline of the calculation region. As can be seen, the dimensionless pressure distributions calculated by the LM-HLLEM scheme and the LM-HLLEM2 scheme are very close to the exact result, but the distributions obtained by the Roe scheme and HLLEM scheme have a large deviation from the exact result. Differences in the results of dimensionless pressure distributions mean that the newly developed schemes are more accurate in simulating low Mach flows.

(a)

(b)

(c)

(d)

(e)

Figure 4: The initial state of the dimensionless pressure field and the state after one period calculated by four schemes. (a) Initial state. (b) Result from Roe scheme. (c) Result from HLLEM scheme. (d) Result from LM-HLLEM scheme. (e) Result from LM-HLLEM2 scheme.

Figure 5: The dimensionless pressure distribution on the centerline of the calculation region after one cycle.

5.4. Low Mach Number Inviscid Flows around the NACA0012 Airfoil

The NACA0012 airfoil is widely used in the evaluation of low Mach schemes. In this section, three schemes are used to calculate the dimensionless pressure field of the NACA0012 airfoil when the inlet Mach numbers are 0.1, 0.01, and 0.001, respectively. The definition of the dimensionless pressure field is also expressed in equation (64). The angle of attack is zero, the simulation region is discretized by 240  120 Cartesian grids, and the residual is calculated from the L2-norm of density. Only when the residual drops by at least five orders of magnitude the calculation is considered to have converged. In this case, the air is considered to be inviscid.

Figures 6–8 show the dimensionless pressure field of the NACA0012 airfoil at inlet Mach Numbers of 0.001, 0.01 and 0.1, respectively. It can be found that the results calculated by the LM-HLLEM scheme and LM-HLLEM2 scheme are all close to the solutions of the incompressible flow, but the original HLLEM scheme fails to achieve the physically correct results in the calculations of the low Mach flows. As shown in Figure 9, the slope of the unmodified HLLEM scheme curve is close to 1, which means the pressure fluctuations of the HLLEM scheme scale with , which is why the original HLLEM scheme cannot guarantee physical solutions. In contrast, with the adoption of the low Mach number corrections, the pressure fluctuations of the LM-HLLEM and LM-HLLEM2 schemes can scale with exactly, which is consistent with the physical situations.

(a)

(b)

(c)

Figure 6: The dimensionless pressure fields of NACA0012 airfoil when . (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.

(a)

(b)

(c)

Figure 7: The dimensionless pressure fields of NACA0012 airfoil when . (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.

(a)

(b)

(c)

Figure 8: The dimensionless pressure fields of NACA0012 airfoil when . (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.

Figure 9: The relationship between the pressure fluctuation of three schemes and the inlet Mach number.

5.5. Low Mach Number Turbulent Flows around the NACA0012 Airfoil

The distribution of pressure coefficient on the upper surface of NACA0012 airfoil in turbulent flow will be calculated in this section to verify the accuracy of the LM-HLLEM and LM-HLLEM2 schemes in the simulations of low Mach turbulent flow. The Spalart–Allmaras one-equation model [42] is employed to simulate the turbulence. The unit Reynolds number is . The settings of the incoming Mach number and angle of attack, the definition of the residual, the judgment of convergence, and the discrete method of the simulation region are all the same as those used in the previous section and will not be repeated here. Figure 10 shows the comparison between the numerical calculation results and the experimental data [43]. The pressure distributions of the upper surface calculated by the LM-HLLEM scheme and the LM-HLLEM2 scheme all agree well with the experimental data, indicating that these two newly developed numerical schemes are both capable of simulating essentially incompressible turbulent flows.

(a)

(b)

(c)

Figure 10: Upwall pressure coefficient for NACA0012 airfoil. (a) Ma₀ = 0.0001. (b) Ma₀ = 0.01. (c) Ma₀ = 0.1.

5.6. Transonic Turbulent Flows around the RAE2822 Airfoil

We use this test case to assess whether the LM-HLLEM scheme and the LM-HLLEM2 scheme can accurately simulate transonic turbulent flows and thus make a smooth transition between the calculations of incompressible and compressible flows.

The setting of flow conditions can be found in [16]. The inlet Mach number is 0.729, and the angle of attack is . Turbulence is simulated by means of the two-equation model [44], and the unit Reynolds number is . The simulation region is discretized by 369  165 quadrilateral grids, and the boundary condition of the outer boundary is “Farfield.” The second order MUSCL reconstruction method with minmod limiter is employed.

What is shown in Figure 11 is the Mach number isolines. It can be seen from the figure that the results of the original HLLEM scheme and the newly developed two schemes are almost indistinguishable. In Figure 12, the surface pressure coefficients computed by the LM-HLLEM scheme and LM-HLLEM2 scheme both agree well with the experimental data. These results indicate that the LM-HLLEM scheme and LM-HLLEM2 scheme can solve the transonic turbulence flows.

(a)

(b)

(c)

Figure 11: The Mach number isolines for RAE2822 airfoil. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.

Figure 12: Pressure coefficient on the surface of RAE2822 airfoil .

This case is further used to test the calculation efficiency of these two low Mach correction schemes. The reason why this case was chosen to test the computational efficiency is because it is turbulent and the computation is relatively large compared to the others, so it can provide more accurate and referential results. This efficiency test was conducted on one single Intel i5-4200M CPU core. The results are presented in Table 1, where “” in this table refers to the ratio of the CPU time required to calculate one time step by each scheme to that required by the HLLEM scheme. Since both of these new schemes introduce additional items, their computational cost is higher than the original HLLEM scheme. As can be seen from Table 1, compared to the original HLLEM scheme, the LM-HLLEM scheme requires an additional 13% of CPU time to calculate one time step, while the LM-HLLEM2 scheme takes a longer time, increasing by 21%, thus indicating that the LM-HLLEM format is more efficient.

Table 1: Comparison of calculation time of three schemes.

Scheme
HLLEM
LM-HLLEM
LM-HLLEM2

1.0
1.13
1.21

6. Conclusions

A detailed asymptotic analysis on the HLLEM scheme is conducted in this paper to study its low Mach number behavior. Through this analysis, we have identified the responsible terms which result in the loss of accuracy of the HLLEM Riemann solver. Two low Mach number fix methods are proposed to rescale these responsible terms at low Mach number, resulting in two new HLLEM-type schemes, namely, the LM-HLLEM scheme and the LM-HLLEM2 scheme. The asymptotic analysis of these schemes indicates that the differences between continuous systems and discrete systems are eliminated. Therefore, the LM-HLLEM scheme and LM-HLLEM2 scheme are both capable of obtaining physically correct solutions at low speeds.

In conclusion, the LM-HLLEM scheme and LM-HLLEM2 scheme inherit many attractive properties of HLLEM scheme, and they are capable of obtaining physically correct solutions in incompressible and weakly compressible cases. Comparing these two modification methods, it can be found that the LM-HLLEM scheme is more efficient, and this advantage will be more prominently highlighted in the cases of heavy calculation burden. For the LM-HLLEM2 scheme, it retains the form of HLLEM scheme to the maximum, so it is quite straightforward to add these low Mach number corrections to existing codes. These two proposed low Mach HLLEM-type schemes are both attractive for their favorable properties and could be extended to compute flow problems ranging from low Mach incompressible to compressible flows.

Data Availability

The data supporting the research results can be found in the article.

Disclosure

Hang Yu and Zhengyu Tian are co-first authors.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (grant no. 11472004).

References

G. Volpe, “Performance of compressible flow codes at low mach numbers,” AIAA Journal, vol. 31, no. 1, pp. 49–56, 1993.

H. Guillard and C. Viozat, “On the behaviour of upwind schemes in the low mach number limit,” Computers & Fluids, vol. 28, no. 1, pp. 63–86, 1999.

H. Guillard and A. Murrone, “On the behavior of upwind schemes in the low mach number limit: II. Godunov type schemes,” Computers & Fluids, vol. 33, no. 4, pp. 655–675, 2004.

P. L. Roe, “Approximate Riemann solvers, parameter vectors, and difference schemes,” Journal of Computational Physics, vol. 43, no. 2, pp. 357–372, 1981.

F. Rieper, “A low-mach number fix for Roe’s approximate Riemann solver,” Journal of Computational Physics, vol. 230, no. 13, pp. 5263–5287, 2011.

K. Oßwald, A. Siegmund, P. Birken, V. Hannemann, and A. Meister, “L2Roe: a low dissipation version of Roe's approximate Riemann solver for low Mach numbers,” International Journal for Numerical Methods in Fluids, vol. 81, no. 2, pp. 71–86, 2016.

E. Turkel, “Preconditioned methods for solving the incompressible and low speed compressible equations,” Journal of Computational Physics, vol. 72, no. 2, pp. 277–298, 1987.

J. R. Edwards and M.-S. Liou, “Low-diffusion flux-splitting methods for flows at all speeds,” AIAA Journal, vol. 36, no. 9, pp. 1610–1617, 1998.

X.-s. Li and C.-w. Gu, “An all-speed roe-type scheme and its asymptotic analysis of low mach number behaviour,” Journal of Computational Physics, vol. 227, no. 10, pp. 5144–5159, 2008.

X.-s. Li, C.-w. Gu, and J.-z. Xu, “Development of roe-type scheme for all-speed flows based on preconditioning method,” Computers & Fluids, vol. 38, no. 4, pp. 810–817, 2009.

X.-s. Li and C.-w. Gu, “Mechanism of roe-type schemes for all-speed flows and its application,” Computers & Fluids, vol. 86, pp. 56–70, 2013.

S.-s. Chen, C. Yan, S. Lou, and B.-x. Lin, “An improved entropy-consistent euler flux in low mach number,” Journal of Computational Science, vol. 27, pp. 271–283, 2018.

S. Dellacherie, “Analysis of godunov type schemes applied to the compressible euler system at low mach number,” Journal of Computational Physics, vol. 229, no. 4, pp. 978–1016, 2010.

S. Dellacherie, J. Jung, P. Omnes, and P.-A. Raviart, “Construction of modified godunov-type schemes accurate at any mach number for the compressible euler system,” Mathematical Models and Methods in Applied Sciences, vol. 26, no. 13, pp. 2525–2615, 2016.

M.-S. Liou, “A sequel to ausm, part ii: AUSM+-up for all speeds,” Journal of Computational Physics, vol. 214, no. 1, pp. 137–170, 2006.

W. Xie, R. Zhang, J. Lai, and H. Li, “An accurate and robust HLLC-type Riemann solver for the compressible Euler system at various Mach numbers,” International Journal for Numerical Methods in Fluids, vol. 89, no. 10, pp. 430–463, 2019.

W. Xie, Y. Zhang, Q. Chang, and L. Hua, “Towards an accurate and robust roe-type scheme for all mach number flows,” Advances in Applied Mathematics and Mechanics, vol. 11, no. 1, pp. 132–167, 2019.

M.-S. Liou, “A sequel to AUSM: AUSM+,” Journal of Computational Physics, vol. 129, no. 2, pp. 364–382, 1996.

E. Shima and K. Kitamura, “Parameter-free simple low-dissipation ausm-family scheme for all speeds,” AIAA Journal, vol. 49, no. 8, pp. 1693–1709, 2011.

K. Kitamura and E. Shima, “Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes,” Journal of Computational Physics, vol. 245, pp. 62–83, 2013.

E. Shima and K. Kitamura, “Multidimensional numerical noise from captured shock wave and its cure,” AIAA Journal, vol. 51, no. 4, pp. 992–998, 2013.

A. Harten, P. D. Lax, and B. v. Leer, “On upstream differencing and godunov-type schemes for hyperbolic conservation laws,” SIAM Review, vol. 25, no. 1, pp. 35–61, 1983.

B. Einfeldt, “On godunov-type methods for gas dynamics,” SIAM Journal on Numerical Analysis, vol. 25, no. 2, pp. 294–318, 1988.

E. F. Toro, M. Spruce, and W. Speares, “Restoration of the contact surface in the hll-riemann solver,” Shock Waves, vol. 4, no. 1, pp. 25–34, 1994.

M. Wesenberg, “Efficient MHD riemann solvers for simulations on unstructured triangular grids,” Journal of Numerical Mathematics, vol. 10, no. 1, pp. 37–71, 2002.

N. T. Nguyen and M. Dumbser, “A new hllem-type riemann solver for compressible multi-phase flows with surface tension,” in Proceedings of the 2015 International Conference on Advanced Computing and Applications (ACOMP), pp. 112–117, IEEE, Ho Chi Minh, Vietnam, November 2015.

Y. Wada, “Approximate Riemann solvers: Improvement of the hllem scheme and its generalization for chemically reacting flows,” Tech. Rep., 1992, NASA STI/Recon Technical Report N 94, Washington, DC, USA.

M. Dumbser and D. S. Balsara, “A new efficient formulation of the hllem riemann solver for general conservative and non-conservative hyperbolic systems,” Journal of Computational Physics, vol. 304, pp. 275–319, 2016.

K. Kitamura and D. S. Balsara, “Hybridized SLAU2-HLLI and hybridized AUSMPW+−HLLI Riemann solvers for accurate, robust, and efficient magnetohydrodynamics (MHD) simulations, part I: one-dimensional MHD,” Shock Waves, vol. 29, no. 5, pp. 611–627, 2019.

W. Xie, W. Li, H. Li, Z. Tian, and S. Pan, “On numerical instabilities of godunov-type schemes for strong shocks,” Journal of Computational Physics, vol. 350, pp. 607–637, 2017.

S. Simon and J. Mandal, “Strategies to cure numerical shock instability in the Hllem Riemann solver,” International Journal for Numerical Methods in Fluids, vol. 89, no. 12, pp. 533–569, 2019.

S. H. Park, J. E. Lee, and J. H. Kwon, “Preconditioned Hllem method for flows at all mach numbers,” AIAA Journal, vol. 44, no. 11, pp. 2645–2653, 2006.

H. Yu, H. Li, J. Lai, and Y. Zhang, “A simple modification of HLLEM approximate Riemann solver applied to the compressible euler system at low mach number,” IOP Conference Series: Materials Science and Engineering, vol. 751, Article ID 012002, 2020.

F. Qu, J. Chen, D. Sun, J. Bai, and C. Yan, “A new all-speed flux scheme for the euler equations,” Computers & Mathematics with Applications, vol. 77, no. 4, pp. 1216–1231, 2019.

S. Klainerman and A. Majda, “Compressible and incompressible fluids,” Communications on Pure and Applied Mathematics, vol. 35, no. 5, pp. 629–651, 1982.

R. Klein, “Semi-implicit extension of a godunov-type scheme based on low mach number asymptotics I: one-dimensional flow,” Journal of Computational Physics, vol. 121, no. 2, pp. 213–237, 1995.

S. H. Park and J. H. Kwon, “On the dissipation mechanism of godunov-type schemes,” Journal of Computational Physics, vol. 188, no. 2, pp. 524–542, 2003.

B. Thornber, A. Mosedale, D. Drikakis, D. Youngs, and R. J. R. Williams, “An improved reconstruction method for compressible flows with low mach number features,” Journal of Computational Physics, vol. 227, no. 10, pp. 4873–4894, 2008.

G. A. Sod, “A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws,” Journal of Computational Physics, vol. 27, no. 1, pp. 1–31, 1978.

E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, Berlin, Germany, 1999.

P. M. Gresho and S. T. Chan, “On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. Part 2: implementation,” International Journal for Numerical Methods in Fluids, vol. 11, no. 5, pp. 621–659, 1990.

P. Spalart and S. Allmaras, “A one-equation turbulence model for aerodynamic flows,” in Proceedings of the 30th Aerospace Sciences Meeting and Exhibit, p. 439, Reno, NV, USA, 1992.

N. Gregory and C. O'. Reilly, Low-speed Aerodynamic Characteristics of Naca0012 Aerofoil Section, Including the Effects of Upper-Surface Roughness Simulating Hoar Frost, 1970, ARC R & M (3726), London, UK.

F. R. Menter, “Influence of freestream values on k-omega turbulence model predictions,” AIAA Journal, vol. 30, no. 6, pp. 1657–1659, 1992.