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

. The HLLEM approximate Riemann solver can capture discontinuities sharply, maintain positive deﬁniteness, and satisfy the entropy condition automatically. These attractive properties make the HLLEM scheme widely used in the simulations of many compressible ﬂuid problems. However, in the simulations of low Mach incompressible ﬂow, 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 identiﬁed. This allows us to develop two modiﬁed methods to solve this low Mach number problem. The ﬁrst 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 diﬀerence 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 ﬂuid problems accurately and robustly.


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. erefore, it is urgently necessary to improve the performance of Godunovtype 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. ey found the results calculated by the numerical schemes contained pressure fluctuations of order Ma, while the actual physical pressure varied with Ma 2 . ey 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. e 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. is 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][10][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. eir 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]. ese all-speed methods not only exhibits high robustness for capturing shock waves but also can be applied to various low Mach number flows. e 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.
is 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. is 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. is 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.
is 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 twophase 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.
ere 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 Godunovtype schemes, the HLLEM scheme also suffers from the accuracy problem in the simulations of low Mach incompressible flows. is 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. rough 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 socalled 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. e 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.
rough 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.
is 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. en, 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.

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. e dimensionless compressible Euler equations can be written as where the equations have been scaled with the following reference quantities: where ρ is the density, u � (u, v) denotes the flow velocity, p represents the pressure, and a is the speed of sound. e 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): 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. e 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 Ma 2 , and therefore, the pressure can be written as which means the leading pressure and first-order pressure are constant in space.

Property 2.
e divergence constraint of zero order of flow velocity u (0) is zero, and thus Property 3. e second-order pressure p (2) is the balance force in the incompressible flow field and satisfies a Poisson equation of the following form:

Detailed Asymptotic Analysis of the HLLEM Scheme
3.1. e Expression of the HLLEM Scheme. Einfeldt [23] proposed the HLLEM scheme in 1988 which is now a widely used approximate Riemann solver. e interface flux of this scheme can be expressed as with where the vectors F and U represent the inviscid flux and conservative variable, respectively. ey can be defined as where e is the specific total energy, p represents the static pressure, ρ is the density, and u and v are velocity vector components. U � un x + vn y denotes the contravariant velocity, and n � (n x , n y ) T represents an outward unit vector perpendicular to the interface. S L and S R 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 a is the speed of sound.

Mathematical Problems in Engineering
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.

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 Ma 1 . ese results are in contrast with those of the compressible Euler equations. eir 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. e nomenclatures defined in this paper can be found in [5].
index vector for neighboring cells, l � (i, j ± 1), (i ± 1, j). A i : area of the reference cell. il: index for edge between cell i and cell j. δ il : length of the interface il. n il � (n x , n y ) T il : unit outer normal vector from cell i to j. t il � (−n y , n x ) T il : unit transverse vector from cell i to j.
Our analysis is conducted on a uniform Cartesian grid to reduce the complexity. So, we set δ il � Δx � Δy � δ and With the Mach number approaching zero, equation (11) and (13) can be reduced to the following forms: e 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 (ρu and ρv) transport, and the energy density(ρe) transport can be expressed as follows: e 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. e continuous equations in dimensionless form can be expressed as e dimensionless equations of horizontal and vertical momentum can be obtained by steps similar to the above: e nondimensional energy conservation equation can be written as All physical quantities are assumed to satisfy the following asymptotic three-term expansion: l∈](i) Order Ma 0 :
To prove that the checkerboard modes can be damped, we express the pressure modes as Assuming the checkerboard does not decrease in time, A is the maximum of the four values of A, B, C, D.
en, equation (30) can be expressed as e 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. erefore, the discrete system satisfies the uniformity of p (0) , i.e., p (0) � constant.

6
Mathematical Problems in Engineering Equations (39) and (40) are both not equal to zero, which implies the first-order pressure p (1) is not constant. erefore, the solutions of discrete Euler equations (19)-(22) support pressure fluctuations of order Ma 1 , 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 (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.

Two Low Mach Fixes for HLLEM Scheme
rough 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. ese two responsible terms are 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 p (1) 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.

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

e LM-HLLEM Scheme.
e correction term can be written as where θ is a scaling function and can be calculated from the following equation: 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 e local Mach number Ma can be calculated from the following equation, which is introduced by ornber et al. [38]: Mathematical Problems in Engineering

e 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. ere 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 p (1) in equations (31) and (32) As the Mach number approaches zero, the weight function θ becomes zero, so we can obtain l∈](i) According to equation (47), the first-order pressure p (1) can satisfy the following relations:

e LM-HLLEM2 Scheme and Its Asymptotic
Analysis. e 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.
e LM-HLLEM2 Scheme. Taking the first-order pressure p (1) in the horizontal direction as an example, i.e., equation (39), the responsible term is Δ il u (0) + Δ il V (0) (n y ) il where the first term Δ il u (0) is derived from the difference of ρu. e second term Δ il V (0) (n y ) il results from the third right eigenvector of the Jacobian matrix R 3 . erefore, we propose the following modifications for conservative variables U LM and the right eigenvector R LM 3 : where θ is defined in equation (43). e modified numerical flux can be written as 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 . e 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. erefore, 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. 8 Mathematical Problems in Engineering

4.2.2.
e 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 p (1) in equations (31) and (32) As the Mach number approaches zero, the scaling function θ becomes zero, so we can obtain the same results as described in equation (48), namely,

Asymptotic Analysis of Discrete
System. e 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 p (1) 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. us, the discrete system satisfies the uniformity of p (1) , i.e., p (1) � constant, which indicates that the discrete solutions of these two schemes both satisfy the pressure fluctuations of order Ma 2 : (53) e numerical results in Section 5.4 support this analytical result.
According to the following equation, the energy density (ρe) (0) in equation (37) can be replaced by the pressure p (0) . rough this replacement, we can obtain the evaluation equations of p (0) : Roe's average values can be simplified according to ρ (0) il � ρ (0) and a (0) il � a (0) . And taking into account Δ il p (0) � 0, equation (37) can be transformed into Following the simplification in [16], namely, assuming that there is no net influx over the boundary due to u (0) in the ghost cells and there is no gradient in p (1) across the boundary, equation (55) can be further transformed into By applying the fact that Δ il p (1) � 0, 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 a (1) il and ρ (1) il 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)

Mathematical Problems in Engineering
For the LM-HLLEM2 scheme: (60) According to the above analysis, the discrete divergence constraint is shown to be satisfied and p (0) is constant. Considering these facts, equations (59) and (60) can be transformed into a Poisson equation of the following form: 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. e asymptotic analysis of these two new schemes demonstrates that the discrete systems can satisfy Properties 1-3.
is means the differences between continuous systems and discrete systems are eliminated. erefore, the LM-HLLEM scheme and the LM-HLLEM2 scheme are both capable of reaching physically correct solutions at low Mach number.

Numerical Test Cases
e 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. e 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. e 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.

Modified Sod Shock Tube.
is 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. is 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 ρ L , u L , p L � (1.0, 0.75, 1.0), 0.125, 0.0, 0.1).
Additionally, the position of discontinuity is set at x � 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." Figure 2 shows the numerical results of four schemes and the exact results at t � 0.2 s. 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. e 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.

Double Rarefaction.
e double rarefaction test case used here contains a trivial contact wave of zero speed and two rarefaction waves. is test case is used to determine whether the proposed LM-HLLEM and LM-HLLEM2 schemes are positively conservative. e position discontinuity is set at x � 0.5, and the initial conditions can be found in [40], namely, e boundary conditions are "Transmissive." e spatial domain is discretized with 400 uniform points in the interval [0, 1]. e first-order solutions computed with three schemes at t � 0.15 s and their exact results are shown in Figure 3. e 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.

Gresho Vortex.
e Gresho vortex is a kind of rotating flow that uses pressure gradient to balance the centrifugal force. is flow mode was proposed by Gresho in 1990 [41].
is case is used here to verify the accuracy and dissipation property of the LM-HLLEM and LM-HLLEM2 schemes in incompressible flows simulations. e 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. e simulation region is set as Ω � [0, 1] × [0, 1], which is evenly discretized by 100 × 100 quadrilateral grids. e initial position of the vortex is (x 0 , y 0 ) � (0.5, 0.5), and the initial radius of the vortex is 0.4. e simulated Mach number is set to 0.01. e 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 p N (x, y) to represent the order of pressure fluctuations under the incompressible limit: e Gresho vortex is a steady-state solution of the incompressible Euler equation, so the dimensionless pressure field p N (x, y) 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 Ma 0 � 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.   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.

Low Mach Number Inviscid Flows around the NACA0012
Airfoil.
e 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. e definition of the dimensionless pressure field is also expressed in equation (64). e 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 Ma 1 , 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 Ma 2 exactly, which is consistent with the physical situations.

Low Mach Number Turbulent Flows around the NACA0012 Airfoil.
e 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.
e Spalart-Allmaras oneequation model [42] is employed to simulate the turbulence. e unit Reynolds number is Re � 6 × 10 6 . e 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]. e 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.

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. e setting of flow conditions can be found in [16]. e 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 × 10 6 . e simulation region is discretized by 369 × 165 quadrilateral grids, and the boundary condition of the outer boundary is "Farfield." e 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.
ese results indicate that the LM-HLLEM scheme and LM-HLLEM2 scheme can solve the transonic turbulence flows.
is case is further used to test the calculation efficiency of these two low Mach correction schemes. e 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. is efficiency test was conducted on one single Intel i5-4200M CPU core. e results are presented in Table 1, where "Time/Time HLLEM " 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.

Conclusions
A detailed asymptotic analysis on the HLLEM scheme is conducted in this paper to study its low Mach number behavior. rough 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. e asymptotic analysis of these schemes indicates that the differences between continuous systems  and discrete systems are eliminated. erefore, 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. ese 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
e data supporting the research results can be found in the article.

Disclosure
Hang Yu and Zhengyu Tian are co-first authors.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.