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

Authors: Hang Yu; Zhengyu Tian; Fan Yang; Hua Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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 Ma, while the actual physical pressure varied with Ma2. 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 Section2, 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(1)∂∂t˜ρ˜+∇⋅ρ˜u˜=0,∂∂t˜ρ˜u˜+ρ˜u˜⋅∇u˜+1Ma2∇p˜=0,∂∂t˜p˜+u˜⋅∇p˜+ρ˜a˜2∇⋅u˜=0,where the equations have been scaled with the following reference quantities:(2)ρ˜=ρρ∗,u˜=uu∗,x˜=xl∗,t˜=tu∗l∗,p˜=pρ∗a∗2,where ρ is the density, u=u,v denotes the flow velocity, p represents the pressure, and a is the speed of sound. The superscript “∗” represents the reference quantities of each variable, and l∗ is an arbitrary length scale. Ma=u∗/a∗ is the reference Mach number.As in [5, 34], we expand the physical quantities ρ˜, p˜, and u˜ into powers of the reference Mach number to find the asymptotic solutions of system (1):(3)ρ˜=ρ˜0+Ma⋅ρ˜1+Ma2⋅ρ˜2+OMa2,p˜=p˜0+Ma⋅p˜1+Ma2⋅p˜2+OMa2,u˜=u˜0+Ma⋅u˜1+Ma2⋅u˜2+OMa2.Inserting equation (3) into equation (1) and sorting the dimensionless equations according to the powers of Ma, 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 numberMa2, and therefore, the pressure can be written as(4)px,t=p0t+Ma2p2x,t,which means the leading pressure and first-order pressure are constant in space.Property 2.
The divergence constraint of zero order of flow velocityu0 is zero, and thus(5)∇⋅u0=0.Property 3.
The second-order pressurep2 is the balance force in the incompressible flow field and satisfies a Poisson equation of the following form:(6)∇2p2=fx,u0,ρ0.


## 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(7)FHLLEM=SR+FUL−SL−FURSR+−SL−+SL−SR+SR+−SL−UR−UL−δ^2α^2R^2−δ^3α^3R^3,with(8)SL−≡min0,SL,SR+≡max0,SR,where the vectors F and U represent the inviscid flux and conservative variable, respectively. They can be defined as(9)U=ρρuρvρe,F=ρUρuU+pnxρvU+pnyρHU,where e is the specific total energy, p represents the static pressure, ρ is the density, and u and v are velocity vector components. U=unx+vny denotes the contravariant velocity, and n=nx,nyT represents an outward unit vector perpendicular to the interface.SL and SR denote the left and right wave speeds, respectively, and can be calculated from(10)SL=minUL−aL,U^−a^,SR=maxUR+aR,U^+a^,where the symbol “∧” represents Roe’s averaged values at the interface and a is the speed of sound.(11)R^2=1u^v^u^2+v^22,R^3=0−nynx−u^ny+v^nx,(12)α^2=Δρ−Δpa^2,α^3=ρ^−nyΔu+nxΔv.(13)δ^2=δ^3=a^a^+U^,where R^2 and R^3 denote the 2nd and 3rd right eigenvectors of Jacobian matrix, α^2 and α^3 are the wave strengths, and δ^2 and δ^3 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 Ma1. 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].i: index vector for cell of reference, i=i,j.νi: index set for neighboring cells, νi=i±1,j,i,j±1.l: index vector for neighboring cells, l=i,j±1,i±1,j.Ai: area of the reference cell.il: index for edge between cell i and cell j.δil: length of the interface il.nil=nx,nyilT: unit outer normal vector from cell i to j.til=−ny,nxilT: unit transverse vector from cell i to j.Δil=ϕi−ϕl: difference operator.ϕil: Roe’s averaged values of ϕi and ϕl.u=u,vT: velocity with Cartesian (global) coordinates.U=u⋅n: normal component of u.V=u⋅t: transverse component of u.Our analysis is conducted on a uniform Cartesian grid to reduce the complexity. So, we setδil=Δx=Δy=δ and Ai=ΔxΔy=δ2.With the Mach number approaching zero, equation (11) and (13) can be reduced to the following forms:(14)SL−≈−ail,SR+≈ail,δ^2=δ^3=1.The expressions of the 2nd and 3rd right eigenvectors are transformed into(15)R^2=1uilviluil2+vil22,R^3=0−nyilnxilVil,and the corresponding wave strengths can be expressed as(16)α^2=Δilρ−Δilpail2,α^3=ρilΔilV.For two-dimensional cases, the HLLEM scheme can be written as(17)ddtUi+1Ai∑l∈νiFUi,Ul,nilδil=0,with the HLLEM flux function(18)FUi,Ul,nil=FUi+FUl2⋅nil−ail2Ul−Ui−α^2R^2−α^3R^3.Complete equations for densityρ transport, momentum density (ρu and ρv) transport, and the energy densityρe transport can be expressed as follows:(19)Aiddtρi+12∑l∈νiρlulnilδil−12∑l∈νiΔilpail=0,(20)Aiddtρiui+12∑l∈νiρlululnil+plnxilδil−ail2∑l∈νiΔilρu−uilΔilρ−Δilpail2+ρilΔilVnyilδil=0,(21)Aiddtρivi+12∑l∈νiρlvlulnil+plnyilδil−ail2∑l∈νiΔilρv−vilΔilρ−Δilpail2−ρilΔilVnxilδil=0,(22)Aiddtρiei+12∑l∈νiρlel+plulnilδil−ail2∑l∈νiΔilρe−uil2+vil22Δilρ−Δilpail2−ρilVilΔilVδil=0.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 Ma.The continuous equations in dimensionless form can be expressed as(23)Ma012∑l∈νiΔilpail−Ma1Aiddtρi+12∑l∈νiρlulnilδil=0.The dimensionless equations of horizontal and vertical momentum can be obtained by steps similar to the above:(24)Ma−112∑l∈νiplnxilδil−Ma0ail2∑l∈νiΔilρu−uilΔilρ−Δilpail2+ρilΔilVnyilδil+Ma1Aiddtρiui+12∑l∈νiρlululnilδil=0,(25)Ma−112∑l∈νiplnyilδil−Ma0ail2∑l∈νiΔilρv−vilΔilρ−Δilpail2−ρilΔilVnxilδil+Ma1Aiddtρivi+12∑l∈νiρlvlulnilδil=0.The nondimensional energy conservation equation can be written as(26)Ma0ail2∑l∈νiΔilρe−Ma1Aiddtρiei+12∑l∈νiρlel+plulnilδil−Ma2ail2∑l∈νiuil2+vil22Δilρ−Δilpail2+ρilVilΔilVδil=0.All physical quantities are assumed to satisfy the following asymptotic three-term expansion:(27)ϕ=ϕ0+Ma1ϕ1+Ma2ϕ2+OMa2as Ma⟶0.Then, insert the asymptotic expansions of all physical quantities into the semidiscrete equations (23)–(26) and sort the equations according to the power of Ma.OrderMa−1:(28)∑l∈νipl0nxilδil=0,(29)∑l∈νipl0nyilδil=0.OrderMa0:(30)∑l∈νiΔilp0ail0=0,(31)∑l∈νipl1nxilδil−∑l∈νiail0Δilρ0u0−uil0Δilρ0−Δilp0ail02+ρil0ΔilV0nyilδil=0,(32)∑l∈νipl1nyilδil−∑l∈νiail0Δilρ0v0−vil0Δilρ0−Δilp0ail02−ρil0ΔilV0nxilδil=0,(33)∑l∈νiΔilρ0e0=0.OrderMa1:(34)Aiddtρi0+12∑l∈νiρl0ul0nilδil−12∑l∈νiΔilp1ail1=0.(35)Aiddtρiui0+12∑l∈νiρl0ul0ul0nilδil−12∑l∈νiail1Δilρ1u1−uil1Δilρ1−Δilp1ail12+ρil1ΔilV1nyilδil+12∑l∈νipl2nxilδil=0.(36)Aiddtρivi0+12∑l∈νiρl0vl0ul0nilδil−12∑l∈νiail1Δilρ1v1−vil1Δilρ1−Δilp1ail12−ρil1ΔilV1nxilδil+12∑l∈νipl2nyilδil=0,(37)Aiddtρiei0+12∑l∈νiρl0el0+pl0ul0nilδil−12∑l∈νiΔilρ1e1=0.Equations (28) and (29) imply that p0 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 asA=pi,j0, B=pi+1,j0, C=pi,j+10, and D=pi+1,j+10. Assuming the checkerboard does not decrease in time, A is the maximum of the four values of A,B,C,D. Then, equation (30) can be expressed as(38)2B+2C−4Aail0=0.The coefficient of equation (38) is positive, and thus equation (38) cannot be equal to zero unless A=B=C=D, which indicates that the checkerboard structure will disappear in large time limit. Therefore, the discrete system satisfies the uniformity of p0, i.e., p0=constant.Roe’s average values can be simplified according toρil0=ρ0 and ail0=a0. And taking into account Δilp0=0, the Ma0 order momentum equations (31) and (32) can be reduced to(39)∑l∈νipl1nxilδil=∑l∈νia0ρ0Δilu0+ΔilV0nyilδil,(40)∑l∈νipl1nyilδil=∑l∈νia0ρ0Δilv0−ΔilV0nxilδil.Equations (39) and (40) are both not equal to zero, which implies the first-order pressure p1 is not constant. Therefore, the solutions of discrete Euler equations (19)–(22) support pressure fluctuations of order Ma1, namely,(41)px,t=p0t+Ma1p1x,t.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 ∑l∈νiΔilu0+ΔilV0nyil (horizontal direction) and ∑l∈νiΔilv0−ΔilV0nxil (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.


## 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(7)FHLLEM=SR+FUL−SL−FURSR+−SL−+SL−SR+SR+−SL−UR−UL−δ^2α^2R^2−δ^3α^3R^3,with(8)SL−≡min0,SL,SR+≡max0,SR,where the vectors F and U represent the inviscid flux and conservative variable, respectively. They can be defined as(9)U=ρρuρvρe,F=ρUρuU+pnxρvU+pnyρHU,where e is the specific total energy, p represents the static pressure, ρ is the density, and u and v are velocity vector components. U=unx+vny denotes the contravariant velocity, and n=nx,nyT represents an outward unit vector perpendicular to the interface.SL and SR denote the left and right wave speeds, respectively, and can be calculated from(10)SL=minUL−aL,U^−a^,SR=maxUR+aR,U^+a^,where the symbol “∧” represents Roe’s averaged values at the interface and a is the speed of sound.(11)R^2=1u^v^u^2+v^22,R^3=0−nynx−u^ny+v^nx,(12)α^2=Δρ−Δpa^2,α^3=ρ^−nyΔu+nxΔv.(13)δ^2=δ^3=a^a^+U^,where R^2 and R^3 denote the 2nd and 3rd right eigenvectors of Jacobian matrix, α^2 and α^3 are the wave strengths, and δ^2 and δ^3 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 Ma1. 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].i: index vector for cell of reference, i=i,j.νi: index set for neighboring cells, νi=i±1,j,i,j±1.l: index vector for neighboring cells, l=i,j±1,i±1,j.Ai: area of the reference cell.il: index for edge between cell i and cell j.δil: length of the interface il.nil=nx,nyilT: unit outer normal vector from cell i to j.til=−ny,nxilT: unit transverse vector from cell i to j.Δil=ϕi−ϕl: difference operator.ϕil: Roe’s averaged values of ϕi and ϕl.u=u,vT: velocity with Cartesian (global) coordinates.U=u⋅n: normal component of u.V=u⋅t: transverse component of u.Our analysis is conducted on a uniform Cartesian grid to reduce the complexity. So, we setδil=Δx=Δy=δ and Ai=ΔxΔy=δ2.With the Mach number approaching zero, equation (11) and (13) can be reduced to the following forms:(14)SL−≈−ail,SR+≈ail,δ^2=δ^3=1.The expressions of the 2nd and 3rd right eigenvectors are transformed into(15)R^2=1uilviluil2+vil22,R^3=0−nyilnxilVil,and the corresponding wave strengths can be expressed as(16)α^2=Δilρ−Δilpail2,α^3=ρilΔilV.For two-dimensional cases, the HLLEM scheme can be written as(17)ddtUi+1Ai∑l∈νiFUi,Ul,nilδil=0,with the HLLEM flux function(18)FUi,Ul,nil=FUi+FUl2⋅nil−ail2Ul−Ui−α^2R^2−α^3R^3.Complete equations for densityρ transport, momentum density (ρu and ρv) transport, and the energy densityρe transport can be expressed as follows:(19)Aiddtρi+12∑l∈νiρlulnilδil−12∑l∈νiΔilpail=0,(20)Aiddtρiui+12∑l∈νiρlululnil+plnxilδil−ail2∑l∈νiΔilρu−uilΔilρ−Δilpail2+ρilΔilVnyilδil=0,(21)Aiddtρivi+12∑l∈νiρlvlulnil+plnyilδil−ail2∑l∈νiΔilρv−vilΔilρ−Δilpail2−ρilΔilVnxilδil=0,(22)Aiddtρiei+12∑l∈νiρlel+plulnilδil−ail2∑l∈νiΔilρe−uil2+vil22Δilρ−Δilpail2−ρilVilΔilVδil=0.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 Ma.The continuous equations in dimensionless form can be expressed as(23)Ma012∑l∈νiΔilpail−Ma1Aiddtρi+12∑l∈νiρlulnilδil=0.The dimensionless equations of horizontal and vertical momentum can be obtained by steps similar to the above:(24)Ma−112∑l∈νiplnxilδil−Ma0ail2∑l∈νiΔilρu−uilΔilρ−Δilpail2+ρilΔilVnyilδil+Ma1Aiddtρiui+12∑l∈νiρlululnilδil=0,(25)Ma−112∑l∈νiplnyilδil−Ma0ail2∑l∈νiΔilρv−vilΔilρ−Δilpail2−ρilΔilVnxilδil+Ma1Aiddtρivi+12∑l∈νiρlvlulnilδil=0.The nondimensional energy conservation equation can be written as(26)Ma0ail2∑l∈νiΔilρe−Ma1Aiddtρiei+12∑l∈νiρlel+plulnilδil−Ma2ail2∑l∈νiuil2+vil22Δilρ−Δilpail2+ρilVilΔilVδil=0.All physical quantities are assumed to satisfy the following asymptotic three-term expansion:(27)ϕ=ϕ0+Ma1ϕ1+Ma2ϕ2+OMa2as Ma⟶0.Then, insert the asymptotic expansions of all physical quantities into the semidiscrete equations (23)–(26) and sort the equations according to the power of Ma.OrderMa−1:(28)∑l∈νipl0nxilδil=0,(29)∑l∈νipl0nyilδil=0.OrderMa0:(30)∑l∈νiΔilp0ail0=0,(31)∑l∈νipl1nxilδil−∑l∈νiail0Δilρ0u0−uil0Δilρ0−Δilp0ail02+ρil0ΔilV0nyilδil=0,(32)∑l∈νipl1nyilδil−∑l∈νiail0Δilρ0v0−vil0Δilρ0−Δilp0ail02−ρil0ΔilV0nxilδil=0,(33)∑l∈νiΔilρ0e0=0.OrderMa1:(34)Aiddtρi0+12∑l∈νiρl0ul0nilδil−12∑l∈νiΔilp1ail1=0.(35)Aiddtρiui0+12∑l∈νiρl0ul0ul0nilδil−12∑l∈νiail1Δilρ1u1−uil1Δilρ1−Δilp1ail12+ρil1ΔilV1nyilδil+12∑l∈νipl2nxilδil=0.(36)Aiddtρivi0+12∑l∈νiρl0vl0ul0nilδil−12∑l∈νiail1Δilρ1v1−vil1Δilρ1−Δilp1ail12−ρil1ΔilV1nxilδil+12∑l∈νipl2nyilδil=0,(37)Aiddtρiei0+12∑l∈νiρl0el0+pl0ul0nilδil−12∑l∈νiΔilρ1e1=0.Equations (28) and (29) imply that p0 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 asA=pi,j0, B=pi+1,j0, C=pi,j+10, and D=pi+1,j+10. Assuming the checkerboard does not decrease in time, A is the maximum of the four values of A,B,C,D. Then, equation (30) can be expressed as(38)2B+2C−4Aail0=0.The coefficient of equation (38) is positive, and thus equation (38) cannot be equal to zero unless A=B=C=D, which indicates that the checkerboard structure will disappear in large time limit. Therefore, the discrete system satisfies the uniformity of p0, i.e., p0=constant.Roe’s average values can be simplified according toρil0=ρ0 and ail0=a0. And taking into account Δilp0=0, the Ma0 order momentum equations (31) and (32) can be reduced to(39)∑l∈νipl1nxilδil=∑l∈νia0ρ0Δilu0+ΔilV0nyilδil,(40)∑l∈νipl1nyilδil=∑l∈νia0ρ0Δilv0−ΔilV0nxilδil.Equations (39) and (40) are both not equal to zero, which implies the first-order pressure p1 is not constant. Therefore, the solutions of discrete Euler equations (19)–(22) support pressure fluctuations of order Ma1, namely,(41)px,t=p0t+Ma1p1x,t.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 ∑l∈νiΔilu0+ΔilV0nyil (horizontal direction) and ∑l∈νiΔilv0−ΔilV0nxil (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Δilu0+ΔilV0nyil (horizontal direction) and Δilv0−ΔilV0nxil (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 p1 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(42)FLM=θ−1ρ^a^20ΔUnxΔUny0,where θ is a scaling function and can be calculated from the following equation:(43)θMa=minMa,1,where a^ and ρ^ are Roe’s averaged values of a and ρ, denoting the speed of sound and density at the grid interface. Ma 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(44)FLM−HLLEM=FHLLEM+θ−1ρ^a^20ΔUnxΔUny0.The local Mach numberMa can be calculated from the following equation, which is introduced by Thornber et al. [38]:(45)Ma=maxui2+vi2ai,uj2+vj2aj.


#### 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 pressurep1 in equations (31) and (32) becomes(46)∑l∈νipl1nxilδil=∑l∈νiail0Δilρ0u0−uil0Δilρ0−Δilp0ail02+ρil0ΔilV0nyilδil+θ−1∑l∈νiail0ρil0ΔilU0nxilδil,∑l∈νipl1nyilδil=∑l∈νiail0Δilρ0v0−vil0Δilρ0−Δilp0ail02−ρil0ΔilV0nxilδil+θ−1∑l∈νiail0ρil0ΔilU0nyilδil.As the Mach number approaches zero, the weight functionθ becomes zero, so we can obtain(47)∑l∈νipl1nxilδil=0,∑l∈νipl1nyilδil=0.According to equation (47), the first-order pressure p1 can satisfy the following relations:(48)pi+1,j1−pi−1,j1=0,pi,j+11−pi,j−11=0.


### 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 pressurep1 in the horizontal direction as an example, i.e., equation (39), the responsible term is Δilu0+ΔilV0nyil where the first term Δilu0 is derived from the difference of ρu. The second term ΔilV0nyil results from the third right eigenvector of the Jacobian matrix R^3. Therefore, we propose the following modifications for conservative variables ULM and the right eigenvector R^3LM:(49)ULM=ρθ⋅ρuθ⋅ρvρe,R^3LM=0−θ⋅nyθ⋅nx−u^ny+v^nx,where θ is defined in equation (43).The modified numerical flux can be written as(50)FLM−HLLEM2=SR+FUL−SL−FURSR+−SL−+SL−SR+SR+−SL−URLM−ULLM−δ^2α^2R^2−δ^3α^3R^3LM.In [34], Qu et al. applied a controlling function not only to the conservative variables U and third right eigenvectors R^3 but also to the second right eigenvectors R^2. 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 R^2 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 pressurep1 in equations (31) and (32) becomes(51)∑l∈νipl1nxilδil=∑l∈νiail0θ⋅Δilρ0u0−uil0Δilρ0−Δilp0ail02+θ⋅ρil0ΔilV0nyilδil,∑l∈νipl1nyilδil=∑l∈νiail0θ⋅Δilρ0v0−vil0Δilρ0−Δilp0ail02−θ⋅ρil0ΔilV0nxilδil.As the Mach number approaches zero, the scaling functionθ becomes zero, so we can obtain the same results as described in equation (48), namely,(52)pi+1,j1−pi−1,j1=0,pi,j+11−pi,j−11=0.


### 4.3. Asymptotic Analysis of Discrete System

The two low Mach number correction methods proposed in Sections4.1 and 4.2 result in the same conclusion, namely, the first order of pressure p1 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 i,j is equal. Thus, the discrete system satisfies the uniformity of p1, i.e., p1=constant, which indicates that the discrete solutions of these two schemes both satisfy the pressure fluctuations of order Ma2:(53)px,t=p0t+Ma2p2x,t.The numerical results in Section5.4 support this analytical result.According to the following equation, the energy densityρe0 in equation (37) can be replaced by the pressure p0. Through this replacement, we can obtain the evaluation equations of p0:(54)ρiei0=1γ−1pi0.Roe’s average values can be simplified according toρil0=ρ0 and ail0=a0. And taking into account Δilp0=0, equation (37) can be transformed into(55)ddtp0+γp01Ai∑l∈νiul0nil2δil=12a01Ai∑l∈νiΔilp1δil.Following the simplification in [16], namely, assuming that there is no net influx over the boundary due to u0 in the ghost cells and there is no gradient in p1 across the boundary, equation (55) can be further transformed into(56)γp01Ai∑l∈νiul0nil2δil=12a01Ai∑l∈νiΔilp1δil.By applying the fact thatΔilp1=0, equation (56) can be simplified into(57)γp01Ai∑l∈νiul0nil2δil=0,and therefore, the discrete divergence constraint is satisfied:(58)γp0∇˜⋅u0=0,where ∇˜ is the discrete divergence operator. So, Property 2 is satisfied in these two modified discrete systems.Assuming thatail1 and ρil1 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:(59)Aiddtρiui0+12∑l∈νiρl0ul0ul0nilδil+12∑l∈νipl2nxilδil,−12∑l∈νiail1Δilρ1u1+ρil1ΔilV1nyil+θ−1∑l∈νiail1ρil1ΔilU1nxilδil=0.For the LM-HLLEM2 scheme:(60)Aiddtρiui0+12∑l∈νiρl0ul0ul0nilδil+12∑l∈νipl2nxilδil,−12∑l∈νiail1θ⋅Δilρ1u1+θ⋅ρil1ΔilV1nyilδil=0.According to the above analysis, the discrete divergence constraint is shown to be satisfied andp0 is constant. Considering these facts, equations (59) and (60) can be transformed into a Poisson equation of the following form:(61)∇˜2p2=fx,u0,ρ0,with ∇˜2 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 Properties1–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.


## 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(42)FLM=θ−1ρ^a^20ΔUnxΔUny0,where θ is a scaling function and can be calculated from the following equation:(43)θMa=minMa,1,where a^ and ρ^ are Roe’s averaged values of a and ρ, denoting the speed of sound and density at the grid interface. Ma 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(44)FLM−HLLEM=FHLLEM+θ−1ρ^a^20ΔUnxΔUny0.The local Mach numberMa can be calculated from the following equation, which is introduced by Thornber et al. [38]:(45)Ma=maxui2+vi2ai,uj2+vj2aj.


### 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 pressurep1 in equations (31) and (32) becomes(46)∑l∈νipl1nxilδil=∑l∈νiail0Δilρ0u0−uil0Δilρ0−Δilp0ail02+ρil0ΔilV0nyilδil+θ−1∑l∈νiail0ρil0ΔilU0nxilδil,∑l∈νipl1nyilδil=∑l∈νiail0Δilρ0v0−vil0Δilρ0−Δilp0ail02−ρil0ΔilV0nxilδil+θ−1∑l∈νiail0ρil0ΔilU0nyilδil.As the Mach number approaches zero, the weight functionθ becomes zero, so we can obtain(47)∑l∈νipl1nxilδil=0,∑l∈νipl1nyilδil=0.According to equation (47), the first-order pressure p1 can satisfy the following relations:(48)pi+1,j1−pi−1,j1=0,pi,j+11−pi,j−11=0.


## 4.1.1. The LM-HLLEM Scheme

The correction term can be written as(42)FLM=θ−1ρ^a^20ΔUnxΔUny0,where θ is a scaling function and can be calculated from the following equation:(43)θMa=minMa,1,where a^ and ρ^ are Roe’s averaged values of a and ρ, denoting the speed of sound and density at the grid interface. Ma 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(44)FLM−HLLEM=FHLLEM+θ−1ρ^a^20ΔUnxΔUny0.The local Mach numberMa can be calculated from the following equation, which is introduced by Thornber et al. [38]:(45)Ma=maxui2+vi2ai,uj2+vj2aj.


## 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 pressurep1 in equations (31) and (32) becomes(46)∑l∈νipl1nxilδil=∑l∈νiail0Δilρ0u0−uil0Δilρ0−Δilp0ail02+ρil0ΔilV0nyilδil+θ−1∑l∈νiail0ρil0ΔilU0nxilδil,∑l∈νipl1nyilδil=∑l∈νiail0Δilρ0v0−vil0Δilρ0−Δilp0ail02−ρil0ΔilV0nxilδil+θ−1∑l∈νiail0ρil0ΔilU0nyilδil.As the Mach number approaches zero, the weight functionθ becomes zero, so we can obtain(47)∑l∈νipl1nxilδil=0,∑l∈νipl1nyilδil=0.According to equation (47), the first-order pressure p1 can satisfy the following relations:(48)pi+1,j1−pi−1,j1=0,pi,j+11−pi,j−11=0.


## 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 pressurep1 in the horizontal direction as an example, i.e., equation (39), the responsible term is Δilu0+ΔilV0nyil where the first term Δilu0 is derived from the difference of ρu. The second term ΔilV0nyil results from the third right eigenvector of the Jacobian matrix R^3. Therefore, we propose the following modifications for conservative variables ULM and the right eigenvector R^3LM:(49)ULM=ρθ⋅ρuθ⋅ρvρe,R^3LM=0−θ⋅nyθ⋅nx−u^ny+v^nx,where θ is defined in equation (43).The modified numerical flux can be written as(50)FLM−HLLEM2=SR+FUL−SL−FURSR+−SL−+SL−SR+SR+−SL−URLM−ULLM−δ^2α^2R^2−δ^3α^3R^3LM.In [34], Qu et al. applied a controlling function not only to the conservative variables U and third right eigenvectors R^3 but also to the second right eigenvectors R^2. 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 R^2 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 pressurep1 in equations (31) and (32) becomes(51)∑l∈νipl1nxilδil=∑l∈νiail0θ⋅Δilρ0u0−uil0Δilρ0−Δilp0ail02+θ⋅ρil0ΔilV0nyilδil,∑l∈νipl1nyilδil=∑l∈νiail0θ⋅Δilρ0v0−vil0Δilρ0−Δilp0ail02−θ⋅ρil0ΔilV0nxilδil.As the Mach number approaches zero, the scaling functionθ becomes zero, so we can obtain the same results as described in equation (48), namely,(52)pi+1,j1−pi−1,j1=0,pi,j+11−pi,j−11=0.


## 4.2.1. The LM-HLLEM2 Scheme

Taking the first-order pressurep1 in the horizontal direction as an example, i.e., equation (39), the responsible term is Δilu0+ΔilV0nyil where the first term Δilu0 is derived from the difference of ρu. The second term ΔilV0nyil results from the third right eigenvector of the Jacobian matrix R^3. Therefore, we propose the following modifications for conservative variables ULM and the right eigenvector R^3LM:(49)ULM=ρθ⋅ρuθ⋅ρvρe,R^3LM=0−θ⋅nyθ⋅nx−u^ny+v^nx,where θ is defined in equation (43).The modified numerical flux can be written as(50)FLM−HLLEM2=SR+FUL−SL−FURSR+−SL−+SL−SR+SR+−SL−URLM−ULLM−δ^2α^2R^2−δ^3α^3R^3LM.In [34], Qu et al. applied a controlling function not only to the conservative variables U and third right eigenvectors R^3 but also to the second right eigenvectors R^2. 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 R^2 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 pressurep1 in equations (31) and (32) becomes(51)∑l∈νipl1nxilδil=∑l∈νiail0θ⋅Δilρ0u0−uil0Δilρ0−Δilp0ail02+θ⋅ρil0ΔilV0nyilδil,∑l∈νipl1nyilδil=∑l∈νiail0θ⋅Δilρ0v0−vil0Δilρ0−Δilp0ail02−θ⋅ρil0ΔilV0nxilδil.As the Mach number approaches zero, the scaling functionθ becomes zero, so we can obtain the same results as described in equation (48), namely,(52)pi+1,j1−pi−1,j1=0,pi,j+11−pi,j−11=0.


## 4.3. Asymptotic Analysis of Discrete System

The two low Mach number correction methods proposed in Sections4.1 and 4.2 result in the same conclusion, namely, the first order of pressure p1 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 i,j is equal. Thus, the discrete system satisfies the uniformity of p1, i.e., p1=constant, which indicates that the discrete solutions of these two schemes both satisfy the pressure fluctuations of order Ma2:(53)px,t=p0t+Ma2p2x,t.The numerical results in Section5.4 support this analytical result.According to the following equation, the energy densityρe0 in equation (37) can be replaced by the pressure p0. Through this replacement, we can obtain the evaluation equations of p0:(54)ρiei0=1γ−1pi0.Roe’s average values can be simplified according toρil0=ρ0 and ail0=a0. And taking into account Δilp0=0, equation (37) can be transformed into(55)ddtp0+γp01Ai∑l∈νiul0nil2δil=12a01Ai∑l∈νiΔilp1δil.Following the simplification in [16], namely, assuming that there is no net influx over the boundary due to u0 in the ghost cells and there is no gradient in p1 across the boundary, equation (55) can be further transformed into(56)γp01Ai∑l∈νiul0nil2δil=12a01Ai∑l∈νiΔilp1δil.By applying the fact thatΔilp1=0, equation (56) can be simplified into(57)γp01Ai∑l∈νiul0nil2δil=0,and therefore, the discrete divergence constraint is satisfied:(58)γp0∇˜⋅u0=0,where ∇˜ is the discrete divergence operator. So, Property 2 is satisfied in these two modified discrete systems.Assuming thatail1 and ρil1 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:(59)Aiddtρiui0+12∑l∈νiρl0ul0ul0nilδil+12∑l∈νipl2nxilδil,−12∑l∈νiail1Δilρ1u1+ρil1ΔilV1nyil+θ−1∑l∈νiail1ρil1ΔilU1nxilδil=0.For the LM-HLLEM2 scheme:(60)Aiddtρiui0+12∑l∈νiρl0ul0ul0nilδil+12∑l∈νipl2nxilδil,−12∑l∈νiail1θ⋅Δilρ1u1+θ⋅ρil1ΔilV1nyilδil=0.According to the above analysis, the discrete divergence constraint is shown to be satisfied andp0 is constant. Considering these facts, equations (59) and (60) can be transformed into a Poisson equation of the following form:(61)∇˜2p2=fx,u0,ρ0,with ∇˜2 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 Properties1–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(62)ρL,uL,pL=1.0,0.75,1.0,ρR,uR,pR=0.125,0.0,0.1.Additionally, the position of discontinuity is set atx=0.3, the computational domain x∈0,1 is evenly discretized with 400 points, and the boundary condition of the boundaries is set as “Nonreflecting.”Figure2 shows the numerical results of four schemes and the exact results at t=0.2s. 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.Figure 2
The profiles of density (a) and Mach number (b) att=0.2s.
(a)(b)


### 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 atx=0.5, and the initial conditions can be found in [40], namely,(63)ρL,uL,pL=1.0,−2.0,0.4,ρR,uR,pR=1.0,2.0,0.4.The boundary conditions are “Transmissive.” The spatial domain is discretized with 400 uniform points in the interval0,1. The first-order solutions computed with three schemes at t=0.15s 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.Figure 3
The profiles of density (a) and Mach number (b) att=0.15s.
(a)(b)


### 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 Ω=0,1×0,1, which is evenly discretized by 100×100 quadrilateral grids. The initial position of the vortex is x0,y0=0.5,0.5, 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 aspNx,y to represent the order of pressure fluctuations under the incompressible limit:(64)pNx,y=px,y−pminpmax−pmin.The Gresho vortex is a steady-state solution of the incompressible Euler equation, so the dimensionless pressure fieldpNx,y defined above should not change with time in the process of calculating the low Mach flow fields.Figure4 shows the initial state of the dimensionless pressure field and the state calculated by four schemes under the condition Ma0=0.01. 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.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.
(a)(b)(c)(d)(e)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.Figures6–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 Ma1, 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 Ma2 exactly, which is consistent with the physical situations.Figure 6
The dimensionless pressure fields of NACA0012 airfoil whenMa0=0.001. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)Figure 7
The dimensionless pressure fields of NACA0012 airfoil whenMa0=0.001. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)Figure 8
The dimensionless pressure fields of NACA0012 airfoil whenMa0=0.1. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)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 Re=6×106. 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.Figure 10
Upwall pressure coefficientCp for NACA0012 airfoil. (a) Ma0 = 0.0001. (b) Ma0 = 0.01. (c) Ma0 = 0.1.
(a)(b)(c)


### 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 2.31°. Turbulence is simulated by means of the k−ω two-equation model [44], and the unit Reynolds number is Re=6.5×106. 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 Figure11 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.Figure 11
The Mach number isolines for RAE2822 airfoil. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)Figure 12
Pressure coefficient on the surface of RAE2822 airfoilCp.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 Table1, where “Time/TimeHLLEM” 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.
SchemeHLLEMLM-HLLEMLM-HLLEM2Time/timeHLLEM1.01.131.21


## 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(62)ρL,uL,pL=1.0,0.75,1.0,ρR,uR,pR=0.125,0.0,0.1.Additionally, the position of discontinuity is set atx=0.3, the computational domain x∈0,1 is evenly discretized with 400 points, and the boundary condition of the boundaries is set as “Nonreflecting.”Figure2 shows the numerical results of four schemes and the exact results at t=0.2s. 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.Figure 2
The profiles of density (a) and Mach number (b) att=0.2s.
(a)(b)


## 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 atx=0.5, and the initial conditions can be found in [40], namely,(63)ρL,uL,pL=1.0,−2.0,0.4,ρR,uR,pR=1.0,2.0,0.4.The boundary conditions are “Transmissive.” The spatial domain is discretized with 400 uniform points in the interval0,1. The first-order solutions computed with three schemes at t=0.15s 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.Figure 3
The profiles of density (a) and Mach number (b) att=0.15s.
(a)(b)


## 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 Ω=0,1×0,1, which is evenly discretized by 100×100 quadrilateral grids. The initial position of the vortex is x0,y0=0.5,0.5, 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 aspNx,y to represent the order of pressure fluctuations under the incompressible limit:(64)pNx,y=px,y−pminpmax−pmin.The Gresho vortex is a steady-state solution of the incompressible Euler equation, so the dimensionless pressure fieldpNx,y defined above should not change with time in the process of calculating the low Mach flow fields.Figure4 shows the initial state of the dimensionless pressure field and the state calculated by four schemes under the condition Ma0=0.01. 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.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.
(a)(b)(c)(d)(e)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.Figures6–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 Ma1, 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 Ma2 exactly, which is consistent with the physical situations.Figure 6
The dimensionless pressure fields of NACA0012 airfoil whenMa0=0.001. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)Figure 7
The dimensionless pressure fields of NACA0012 airfoil whenMa0=0.001. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)Figure 8
The dimensionless pressure fields of NACA0012 airfoil whenMa0=0.1. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)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 Re=6×106. 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.Figure 10
Upwall pressure coefficientCp for NACA0012 airfoil. (a) Ma0 = 0.0001. (b) Ma0 = 0.01. (c) Ma0 = 0.1.
(a)(b)(c)


## 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 2.31°. Turbulence is simulated by means of the k−ω two-equation model [44], and the unit Reynolds number is Re=6.5×106. 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 Figure11 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.Figure 11
The Mach number isolines for RAE2822 airfoil. (a) HLLEM. (b) LM-HLLEM. (c) LM-HLLEM2.
(a)(b)(c)Figure 12
Pressure coefficient on the surface of RAE2822 airfoilCp.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 Table1, where “Time/TimeHLLEM” 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.
SchemeHLLEMLM-HLLEMLM-HLLEM2Time/timeHLLEM1.01.131.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.



---

*Source: 7451240-2020-08-17.xml*