Jean-Christophe Boniface
email: jean-christophe.boniface@onera.fr

Rescaling of the Roe Scheme in Low Mach-Number Flow Regions II: Artificial Speed of Sound and Low Mach Number Fix

Keywords: Euler equations, Low Mach number, Compressible flow, Potential flow, Difference scheme, Finite-volume method, Roe scheme, Artificial speed of sound, Low Mach number fix, Incompressible time-scale, Acoustic time-scale, Wave equation, Checkerboard pressure-velocity decoupling, Stability, Implicit scheme, Characteristic time-step

We look at two simple modifications of the Roe scheme in the incompressible limit, based on different ideas: the Rossow's artificial speed of sound and the Rieper's low Mach number fix. Both schemes modify the eigenspaces of the dissipation matrix. The analysis emphasizes the properties of the dissipation matrix for the Von Neumann stability, the asymptotic behaviour and the solution accuracy in the incompressible limit. Numerical results in the very low-speed limit are discussed for robustness, consistency and accuracy issues of the numerical procedure. Possible occurrence of checkerboard pressure modes, when using a collocated arrangement for velocity components and pressure in the finite-volume scheme, and spurious acoustic modes, is also illustrated for both schemes.

Introduction

In the last twenty years, a number of numerical schemes based on the Roe's approximate Riemann solver [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes[END_REF] have been developed in order to compute consistent solutions of the Euler equations with the incompressible limit, while preserving the shock-capturing capability of the original Roe scheme. This is further motivated nowadays, where industry routinely uses CFD in the design process and there are many situations in which industrial flows are characterized by coexisting compressible and nearly incompressible flows. Local low-speed preconditioners became popular in years 90s with the pioneering works of Turkel [START_REF] Turkel | Preconditioned Methods For Solving the Incompressible And Low Speed Compressible Equations[END_REF], Choi and Merkel [START_REF] Choi | The Application of Preconditioning To Viscous Flows[END_REF], Weiss and Smith [START_REF] Weiss | Preconditioning Applied to Variable and Constant Density Flows[END_REF], Van Leer et al. [START_REF] Van Leer | Characteristic Time-Stepping or Local Preconditioning of the Euler Equations[END_REF] which were further analyzed and generalized by Turkel et al. [START_REF] Turkel | [END_REF][START_REF] Turkel | Leer Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes[END_REF]. This approach has optimal properties for the convergence of the steady-state problem, improving the conditioning of the flux Jacobian matrix by slowing down the acoustic wave speeds towards the the local velocity, and the accuracy, with a consistent approximation of low Mach number flows. However, the extension to unsteady flows is not trivial and without special care, the time accuracy may be lost [START_REF] Turkel | Vatsa Choice of Variables and Preconditioning for Time Dependent Problems AIAA paper[END_REF][START_REF] Turkel | Vatsa Local Preconditioners for Steady State and Dual Time Stepping ESAIM[END_REF].

Since then, the use of local preconditioners for time-dependent problems did not get a large acceptance by the CFD community. Additionally, low speed preconditioners modify the flux Jacobian matrix, and all boundary conditions based on the characteristic variables or the Riemann invariants must be completely reformulated [START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows AIAA paper[END_REF][START_REF] Turkel | Assessment of Preconditioning Methods For Multidimensional[END_REF]12]. For large aerodynamics codes developed for industrial internal and external flow configurations, a very large number of boundary conditions using the characteristic theory for compressible flows should be reformulated, which is practically not possible.

In the years 2000, a new density-based approach for compressible flows was considered, aiming at drastically simplifying the extension of compressible schemes to handle incompressible flows. This was initiated with the case of the Roe-Turkel scheme studied by Guillard and Viozat in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF], and followed with the All-Speed Roe-type scheme developed by Li and Gu only modifying the eigenvalues [START_REF] Li | An All-Speed Roe-Type Scheme and its Asymptotic Analysis of Low Mach Number Behaviour[END_REF][START_REF] Li | Xu Development of Roe-Type Scheme for All-Speed Flows Based on Preconditioning Method[END_REF], the artificial speed of sound according to Rossow [START_REF] Rossow | A Flux-Splitting Scheme for Compressible and Incompressible Flows[END_REF][START_REF] Rossow | A blended pressure/density based method for the computation of incompressible and compressible flows[END_REF][START_REF] Rossow | Efficient Computation of Compressible and Incompressible Flows[END_REF] and more recently the Rieper's low Mach number fix [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF], which of these approaches only modify the matrix-valued dissipation of the Roe scheme. With this new way of formulating the low speed preconditioning, the numerical scheme recovers the formulation for compressible flows, with a standard extension to unsteady flows, using usual time-marching algorithms and it is no longer required to reformulate characteristic-like boundary conditions. All these Roe-type schemes actually provide a consistent rescaling of the matrix dissipation with the conservation laws in the low speed limit, as pointed out by Turkel et al. [START_REF] Turkel | Leer Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes[END_REF], and ensure the proper scaling of the pressure field, following the asymptotic discrete analysis developed by Guillard and Viozat [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF]. This approach has likely motivated the extension to low Mach number flow of Godunov schemes [START_REF] Dellacherie | Analysis of Godunov Type Schemes Applied to The Compressible Euler System at Low Mach Number[END_REF][START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF], AUSM-type schemes [START_REF] Edwards | Low-Diffusion Flux-Splitting Methods for Flows at All Speeds[END_REF][START_REF] Mary | Deville An Algorithm for Unsteady Viscous Flows at all Speeds[END_REF][START_REF] Mary | Large Eddy Simulation of Flow Around an Airfoil Near Stall[END_REF] or more recently to approximate Riemann solver using HHL-type schemes [START_REF] Rieper | On the dissipation mechanism of upwind-schemes in the low Mach umber regime: A comparison between Roe and HLL[END_REF][START_REF] Pelanti | Low Mach number preconditioning techniques for Roe-type and HLLC-type methods for a two-phase compressible flow model[END_REF]. However, in modifying the matrix-valued dissipation of the Roe scheme, there are certain features of the numerical procedure that it would also be worth striving for.

First, the proper formulation of the Von Neuman stability condition for the explicit scheme is an essential feature of the numerical procedure, when no preconditioner is applied to the time-derivative. In the formulation of their Lax-Friedrichs-Turkel scheme, Birken and Meister [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF] have demonstrated that, in the low speed limit, the explicit scheme is stable if the local time-step satisfies a stringent stability condition, given by the spectral radius of the matrix dissipation. This is also the case of the Roe-Turkel scheme reformulated in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF] and with the Rossow's artificial speed of sound, as demonstrated in the paper. In the case of the low Mach number fix for instance, we will see that the stability condition is close to the usual Von Neumann condition for compressible flows, although not well posed for the classical Fourier analysis.

Second, some authors have argued that an accurate extension of the Roe scheme to the incompressible limit should satisfy the discrete divergence free constraint of the leading order velocity in the asymptotic limit. This is especially the case of the All-Speed Roe scheme developed in [START_REF] Li | An All-Speed Roe-Type Scheme and its Asymptotic Analysis of Low Mach Number Behaviour[END_REF] and the low Mach number fix [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF]. Some other approximate Riemann solvers were also designed to satisfy the discrete divergence free constraint, see for instance the HLLC scheme developed in [START_REF] Xie | An accurate and robust HLLC-type Riemann solver for the compressible Euler system at various Mach numbers[END_REF]. Note that the divergence constraint is not enforced by the Roe-Turkel scheme [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF] and the artificial speed of sound, as shown in the paper. As recently advocated by Guillard et al. [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF], "they are sound theoretical reasons to prefer a scheme that does not enforce the discrete divergence free constraint". Assuming isentropic flow, the functional analysis indicates that two limits of the continuous Euler equations exist, with on one hand the elliptic incompressible system satisfying the divergence free constraint for the slow incompressible time scale, and on the other hand, an hyperbolic acoustic system for the fast acoustic time scale. Then, the solution of the Euler equations in the incompressible limit is characterized by acoustic-incompressible interactions, in which the incompressible solution is only one component of the solution. This asymptotic behaviour with two time scales was first described by Klein for the compressible Euler equations applied to combustion science [START_REF] Klein | Semi-Implicit Extension of Godunov-Type Scheme Based on Low Mach Number Asymptotics I: One-Dimensional Flow[END_REF]. In [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF], it is shown that the All-Speed Roe scheme and the low Mach number fix have large inconsistencies in a vertex-centred triangular mesh. This stems from the fact that, from the incompressible theory, using a collocated arrangement for velocity and pressure on unstaggered grids without specific stabilization, undampted spurious pressure modes cannot be ruled out with a cell-centred finite-volume scheme. It is also shown in the paper that the low Mach number fix is intrinsically prone to pressure-velocity decoupling, regardless of the mesh used, structured or unstructured and this a consequence of enforcing the divergence free constraint of the velocity.

Third, when removing the preconditioning of the flux Jacobian matrix, the rescaled Roe schemes are subjected to a critical lack of robustness and accuracy near stagnation points, when the Euler equations are formulated for the primitive or the conservative variables. The mathematical framework used by Darmofal et. al [START_REF] Darmofal | The Importance of Eigenvectors for Local Preconditioners of the Euler Equations[END_REF] has demonstrated that the Roe scheme has a transient error growth at stagnation points due to the non-normality of the Jacobian matrix. This issue can only be circumvented by designing specific low speed preconditioners, which are not considered in the formulation of the rescaled Roe schemes mentioned above. This is still an open problem, discussed in the paper, especially noticeable in the very incompressible limit.

Thus, it was found interesting to thoroughly compare the Rossow's artificial speed of sound and the Rieper's low Mach number fix. Both schemes introduce very simple modifications of the first-order Roe scheme, based on different ideas, for the computation of low Mach number flows. Both automatically return the original Roe scheme at the sonic point according to the local Mach number, and therefore can also be used for the computation compressible flows with little effort. However, their resulting matrix-valued dissipation leads to very different behaviours of the discrete solutions in the incompressible limit. This contribution is organized as follows. In section 2, a general framework used for the analysis of roe-Type schemes in the low Mach-number limit is presented. Next, in section 3, the Rossow's artificial speed of sound and the Rieper's fix are introduced and the corresponding matrix valued dissipation is derived for the entropic symmetrizing variables. Following the results of section 3, the eigenvalues and spectral radius of the respective matrix dissipation are derived in section 4, and the Von Neumann condition for stability is discussed for both approaches. In section 5, properties of both schemes are compared in the low-speed limit, considering the Weiss-Smith decomposition and the normalized Euler equations for the continuous problem, in which the matrix-valued dissipation is explicitly accounted for in the asymptotic analysis. A complementary scheme truncation error analysis is carried out in section 6. Then, numerical experiments presented in section 7 illustrate the above analysis. Some comments and guidelines are given as concluding remarks in section 8.

General Framework for the Analysis of Roe-Type Schemes in the Low Mach-Number Limit

The following set of independent variables are introduced:

-W = [ρ, ρu, ρv, ρE] T (conservative variables used for shock-capturing),

-Q = [ρ, ρu, ρv, ρH] T (flux density vector used for shock-capturing),

-U = [ρ, u, v, p] T (primitive variables used for the asymptotic analysis in the incompressible limit),

d W0 = [dΦ, du, dv, dS ] T (symmetrizing variables used for the stability analysis in the incompressible limit),

where ρ is the density, V = [u, v] T are the velocity components, p is the pressure, E is the total energy per unit volume, H = E + p ρ is the total enthalpy per unit mass, with the differential variables dΦ = d p ρc proportional to the pressure and dS = d p-c 2 dρ ρc proportional to the entropy, c being the speed of sound.

We suppose satisfied the ideal gaz law. So c 2 = γ p ρ and H = c 2 (γ -1)

+ |V| 2 2
, γ being the ratio of the specific heats.

We consider the Euler equations in integral form on each computational cell Ω of the structured mesh, with boundary ∂Ω with unit outward normal n

= [n x , n y ] T d dt Ω W dV + ∂Ω F(W).n dS = 0 ( 1
)
where F is the flux vector. System (1) is discretized using the finite-volume method. This leads to the semi-discrete form

V d dt W + R(W, n) = 0, with R(W, n) = cell boundary F(W).n S
where V denotes the volume of the computational cell and S is the area of the cell interface. We will consider next modifications in the low-speed limit of the baseline first-order Roe scheme [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes[END_REF], with the numerical flux h discretizing the physical flux F across the cell interfaces

h = 1 2 [F(W R ) + F(W L )] - 1 2 | A roe | (W R -W L ) (2)
where W L and W R are the left and right reconstructed solution vectors in adjacent cells to the cell interface and A roe is the Roe matrix. Three different formulations of the dissipation vector

d = |A roe | (W R -W L )
will be considered for the analysis as described in the following.

First formulation of the dissipation vector

The Rossow's artificial speed of sound and the Rieper's low Mach number fix are both formulated expanding the dissipation vector in terms of jumps of the conservative variables ∆W = W R -W L projected in the right eigenvector basis of the matrix dissipation [START_REF] Van Leer | Progress in multi-dimensional upwind differencing[END_REF]. For the Roe's flux-splitting scheme, the resulting individual disturbances can be further expressed in term of jumps of the primitive variables ∆U = U R -U L as follows

d = ∆F 0 + ∆F + + ∆F - (3)
with:

∆F 0 = |λ 0 |               ∆ρ - ∆p c 2               1 V |V| 2 2               + ρ           0 ∆V -∆q n n V.∆V -q n ∆q n                         , ∆F -= |λ -| ∆p -ρc∆q n 2c 2           1 V -cn H -cq n           , ∆F + = |λ + | ∆p + ρc∆q n 2c 2           1 V + cn H + cq n          
, and q n = un x + vn y is the normal velocity across the cell interface, ∆q n = ∆un x + ∆vn y is the jump of the normal velocity and λ 0 = q n , λ + = q n + c, λ -= q nc are the wave speeds, eigenvalues of matrix A roe . Note that ∆F 0 is explicitly formulated for the jump of the tangential velocity q t with ∆V -∆q n n = ∆q t in the momentum equation and V.∆Vq n ∆q n = q t .∆q t in the total energy equation.

Second formulation of the dissipation vector

The second formulation of the dissipation vector corresponds to the decomposition according to Weiss-Smith [START_REF] Weiss | Preconditioning Applied to Variable and Constant Density Flows[END_REF]. This decomposition has been especially used by Li and Gu [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF] for the analysis of Roe-type schemes in the low-speed limit

d = |q n |           ∆ρ ∆(ρV) ∆(ρE)           + δq n           ρ ρV ρH           + δp           0 n q n           , (4)
where the first term is the basic scalar upwind dissipation, with for the original Roe scheme

δq n = (c -|q n |) ∆p ρc 2 + q n c ∆q n and δp = q n c ∆p + (c -|q n |)ρ∆q n
The second and third coefficients δq n and δp are related to the convected velocity and the pressure at the cell interface, respectively. In the low-speed limit, Li and Gu have shown that the two last terms of decomposition (4) are also related to the pressure checkerboard issue and the physical consistency of the dissipation vector. These coefficients are modified by the artificial speed of sound and the low Mach number fix.

Third formulation of the dissipation vector

As pointed out by Birken and Meister in [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF], the standard CFL condition used for the computation of compressible flows may be no longer valid in the incompressible limit. It was demonstrated that the stability condition is strongly concerned with the eigenvalues of the modified dissipation matrix. This is especially the case of modified Roe-type schemes, since in most of the cases, both the eigenvectors and eigenvalues of the dissipation matrix are modified. So the reformulation of the CFL condition becomes an essential feature of the numerical procedure. In the low speed limit, the local time step may be restricted with a severe stability constraint and the use of an implicit scheme is highly recommended. Turkel [START_REF] Turkel | Vatsa Local Preconditioners for Steady State and Dual Time Stepping ESAIM[END_REF]12] suggested to use the symmetrizing variables d W0 , which greatly simplifies the analysis and a much simpler formulation of the dissipation matrix can be easily derived. A necessary condition for stability also requires that the asymptotic behaviour of the entries of the dissipation matrix matches that of the terms of the Euler equations [START_REF] Wong | The Solution of the Euler Equations at Low Mach Number Using a Stabilized Finite Element Algorithm[END_REF], and this choice of variables also simplifies the scaling analysis of the dissipation matrix.

Then, for the stability analysis, the dissipation vector will be formulated for jumps ∆ W0 in d W0 variables, within the following change of variables

d0 = ∂ W0 ∂W d = ∂ W0 ∂W | A roe | ∂W ∂ W0 ∆ W0 = | Ã0 |∆ W0 , (5)
with

∂ W0 ∂W and ∂W ∂ W0
given in Appendix A.

The Jacobian matrix | Ã0 | in d W0 variables has the same eigenvalues as | A roe | in conservative variables by similarity transformation. Assuming the flow subsonic, the Jacobian matrix and the corresponding eigenspaces can be derived explicitly. We will see in the next sections that both the Rossow's artificial speed of sound and the Rieper's low Mach number fix modify the eigenspaces of the Jacobian matrix and the dissipation vector will be expressed with a modified dissipation matrix d0 = D0 ∆ W0 .

Additionally, the limit to incompressible equations has been investigated in [START_REF] Turkel | Leer Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes[END_REF]12] in symmetrizing d W0 variables. Within an asymptotic analysis for a constant mesh size, it was shown that the proper scaling of the matrix-valued dissipation should ensure a maximum allowed dissipation in the incompressible limit. In modifying the Jacobian matrix | Ã0 |, we will see that the Rossow's artificial speed of sound and the Rieper's low Mach number fix yield less dissipative schemes in the incompressible limit than the Roe-Turkel scheme.

Part of the asymptotic analysis also requires to formulate the dissipation vector in primitive variables U, when considering the normalized Euler equations, following the framework analysis introduced in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF]. The change of variables d W0 into dU = [dρ, du, dv, d p] T is straightforward.

Rescaling of the Roe Scheme According to Rossow and Rieper

In all the following, we assume the flow subsonic across the cell interface, i.e. M n ≤ 1, where M n = |q n | c is the directional Mach number. Then |λ -| = cq n and |λ + | = c + q n , and the dissipation vector formulated with (3) for the jumps of the primitive variables can be formulated explicitly.

Rossow's artificial speed of sound

In successive works [START_REF] Rossow | A Flux-Splitting Scheme for Compressible and Incompressible Flows[END_REF][START_REF] Rossow | A blended pressure/density based method for the computation of incompressible and compressible flows[END_REF][START_REF] Rossow | Efficient Computation of Compressible and Incompressible Flows[END_REF], Rossow introduces the interface Mach number M 0 = sgn( q n c ) min(M n , 1) and formulates the dissipation vector as an expansion with terms factored by M 0 and (1 -|M 0 |). So M 0 is bounded by ±1 in the supersonic range. For subsonic flows, the interface Mach number is then simply defined as M 0 = q n c and therefore |M 0 | = M n the directional Mach number.

Let d = [∆F ρ , ∆F ρV , ∆F ρE ] T be the components of the dissipation vector corresponding to the mass, momentum and energy conservation laws. The dissipation for the original Roe's flux-difference scheme is then formulated as follows

∆F ρ = |q n |∆ρ + ρM 0 ∆q n + 1 c (1 -|M 0 |)∆p ∆F ρV = |q n |V∆ρ + ρ|q n |∆V + ρVM 0 + ρcn(1 -|M 0 |) ∆q n + nM 0 + 1 c V(1 -|M 0 |) ∆p ∆F ρE = |q n | |V| 2 2 ∆ρ + ρ|q n |V.∆V + ρHM 0 + ρq n c(1 -|M 0 |) ∆q n + 1 (γ-1) |q n | + q n M 0 + H c (1 -|M 0 |) ∆p (6)
The rescaling of the Roe scheme in the incompressible limit was achieved by introducing an artificial speed of sound as first described in [START_REF] Rossow | A Flux-Splitting Scheme for Compressible and Incompressible Flows[END_REF]. This approach was further developed by Rossow in [START_REF] Rossow | A blended pressure/density based method for the computation of incompressible and compressible flows[END_REF] for a blended pressure/density based scheme and in [START_REF] Rossow | Efficient Computation of Compressible and Incompressible Flows[END_REF] considering a pure density-based approach. We shall consider here this later implementation. The rescaling is obtained by replacing the speed of sound c in expansions [START_REF] Turkel | [END_REF] with an artificial speed of sound c , excepted in the formulation of the total enthalpy, which must be formulated with the speed of sound according to the ideal gaz law. Thus, the corresponding dissipation for the rescaled Roe scheme with artificial speed of sound reads

∆F ρ = |q n |∆ρ + ρM 0 ∆q n + 1 c (1 -|M 0 |)∆p ∆F ρV = |q n |V∆ρ + ρ|q n |∆V + ρVM 0 + ρc n(1 -|M 0 |) ∆q n + nM 0 + 1 c V(1 -|M 0 |) ∆p ∆F ρE = |q n | |V| 2 2 ∆ρ + ρ|q n |V.∆V + ρHM 0 + ρq n c (1 -|M 0 |) ∆q n + 1 (γ-1) |q n | + q n M 0 + H c (1 -|M 0 |) ∆p (7)
So, all terms factored by M 0 dominating in the compressible regime are not modified by the correction. The artificial speed of sound is defined from a local preconditioner

c = c α 2 M 2 n + β 2 , with α = 1 2 (1 -β 2 ) and β 2 = min(max(M 2 , M 2 re f ), 1), ( 8
)
where M is the Mach number, and M re f is related to a reference Mach number. This formulation actually borrows the modified speed of sound arising from the preconditioning of the Jacobian matrix P A, with all eigenvalues slowed down towards the local flow velocity as the Mach number goes to zero (see for instance Turkel et al for a general form of preconditioners [START_REF] Turkel | Leer Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes[END_REF]). Then we have

c -→ V re f as M -→ 0, (9)
where V re f is some reference speed. This asymptotic behaviour has large effects on the pressure differences with all terms with coefficients 1 c (1

-|M 0 |) increased with a factor O( 1 M ) in the incompressible limit 1 c (1 -|M 0 |) O(1) while 1 c (1 -|M 0 |) O(M).
So the dominance of the pressure difference terms is strongly amplified in the matrix dissipation. The artificial speed of sound [START_REF] Turkel | Vatsa Choice of Variables and Preconditioning for Time Dependent Problems AIAA paper[END_REF] actually corresponds to an artificial pressure p defined through the following relationship

p = ρc 2 = γp 1 4 (1 -M 2 ) 2 M 2 n + M 2 (10)
assuming β 2 = M 2 . Let introduce in the two-dimensional case the local velocity flow angle θ with respect to the unit outward normal n to the cell interface. Then q n = |V|cos(θ), and so, for the directional Mach number, M n = M|cos(θ)|.

Upon substitution in expression [START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows AIAA paper[END_REF], this can be reformulated as

p = ρ |V| 2 2 1 2 (1 -M 2 ) 2 cos(θ) 2 + 2 ρ |V| 2 2 1 2 cos(θ) 2 + 2 when M re f < M << 1.
Then, we see that the artificial pressure at the the cell interfaces can be interpreted in the low speed limit as a smooth variation of the dynamic pressure according to the local velocity flow angle. Expressions [START_REF] Turkel | Leer Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes[END_REF] are then further expanded as indicated in section 2.2 to be recast in jumps of d W0 variables. By applying the pre-multiplication with ∂ W0 ∂W for the change of variables, the dissipation vector can be expressed in matrix dissipation form d0 = D0 ∆ W0 , with the modified matrix-valued dissipation in d W0 variables

D0 =                                      r(c -|q n |) + |q n | n x q n n y q n 0 n x q n n 2 x 1 r (c -|q n |) + |q n | n x n y 1 r (c -|q n |) 0 n y q n n x n y 1 r (c -|q n |) n 2 y 1 r (c -|q n |) + |q n | 0 as M -→ 0.
So we clearly see that the Rossow's artificial speed of sound introduces less consistent dissipation in the incompressible limit.

Rieper's low Mach number fix

A low Mach number fix was originally proposed by Rieper in [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF]. This idea stems from the observation that jumps of the normal velocity component at the cell interfaces are responsible of non-physical entropy production in the low speed limit. Rieper then suggested a smooth correction for vanishing jumps ∆q n in the wave strengths, occurring explicitly in the first formulation (3) of the dissipation vector

∆p ± ρc∆q n 2c 2 -→ ∆p ± ρcz∆q n 2c 2 ,
without modifying the jumps for the tangential velocity. The proposed fix z is simply an averaged Mach number across each cell interface, bounded by 1

z = min(1, M) with M = |q n | + |q t | c cell interface , q t = un y -vn x . (13)
So

z -→ 0 as M -→ 0, (14)
regardless of the cell face orientation. In our implementation, the Roe's average was used as averaging for the fix.

The same idea was considered in [START_REF] Oßwald | a Low Dissipation Version of Roe's Approximate Riemann Solver for Low Mach Numbers[END_REF] with the L 2 Roe scheme, by also applying the same fix to the jumps for the tangential velocity in order to further decrease the numerical dissipation for a specific test-case with decreasing isotropic turbulence. This later correction is not investigated in this paper. Note that in reference [START_REF] Oßwald | a Low Dissipation Version of Roe's Approximate Riemann Solver for Low Mach Numbers[END_REF], we also find z = min(1, max(M L , M R )), where M L/R are the Mach numbers computed with the right and left states to the cell interface. The Rieper's fix also modifies both the eigenvalues and eigenvectors of the matrix dissipation. The modified dissipation vector can be expressed for the jumps of the primitive variables as follows

∆F ρ = |q n |∆ρ + ρzM 0 ∆q n + 1 c (1 -|M 0 |)∆p ∆F ρV = |q n |V∆ρ + ρ|q n |∆V + ρVzM 0 + ρcn(z -|M 0 |) ∆q n + nM 0 + 1 c V(1 -|M 0 |) ∆p ∆F ρE = |q n | |V| 2 2 ∆ρ + ρ|q n |V.∆V + ρHzM 0 + ρq n c(z -|M 0 |) ∆q n + as M -→ 0.
Therefore, it can be seen that the Rieper's low Mach number fix yields a much less dissipative scheme, especially in the scaling of the pressure equation, with diminishing dissipation coefficients with an order O(M) compared to the Rossow's artificial speed of sound and up to O(M 2 ) compared to the Roe-Turkel scheme.

Stability Condition for the Explicit Scheme

The proper formulation of the stability condition for the explicit scheme is an essential feature of the numerical procedure. The stability condition is strongly concerned with the eigenvalues of the matrix-valued dissipation, as the rescaling has a major effect on the solution accuracy in the incompressible limit, ensuring a consistent balancing with the flux divergence. A different behaviour for the acoustic eigenvalues was found for the two investigated schemes, as described in the following, with consequences on the Von Neumann condition for stability of the explicit scheme in the low Mach-number range.

Rossow's artificial speed of sound

We find that matrix [START_REF] Turkel | Assessment of Preconditioning Methods For Multidimensional[END_REF] has the following eigenvalues

             µ 0 = |q n |
with multiplicity 2 for the 2D problem

µ ± = c 2r (r 2 + 1)(1 -M n ) + 2rM n ± (r 2 -1) 2 (1 -M n ) 2 + 4r 2 M 2 n (r = c c ) (17)
with the following properties:

1. The acoustic eigenvalues are positive in the subsonic range with µ + > 0 and

µ + µ -= c 2 r 2 r 2 (1 -M n ) + r(1 + r 2 )M n (1 -M n ) ≥ 0 when M n ≤ 1; 2.
Thus the symmetric matrix [START_REF] Turkel | Assessment of Preconditioning Methods For Multidimensional[END_REF] is positive semi-definite in the subsonic range; 3. The spectral radius is given by max(µ 0 , µ -, µ + ) = µ + (see Appendix B); 4. We find that µ + ≥ |q n | + c = ρ( A) the spectral radius of the original Jacobian matrix, ∀M n ≤ 1 (see Appendix B); 5. At the sonic line, c = c => r = 1 and µ + = ρ( A).

The right and left eigenvector matrices are explicitly derived in Appendix A, for the diagonalization in the symmetrizing variables d W0 variables

D0 = R0 Λ L0 with Λ = diag(µ + , µ -, |q n |, |q n |).
In the incompressible limit, M n -→ 0 and the acoustic eigenvalues behave as

µ ± c 2 (r 2 + 1) ± |r 2 -1| . Since r ≥ 1 µ + c 2 c = O( 1 M 2 ), µ -c V re f = O(1). ( 18
)
For the other extreme case, at the sonic line M n -→ 1 and

µ + = 2c, µ -= 0.
These above results show that Rossow's artificial speed of sound is closely related to the rescaling of the Roe scheme in the incompressible limit reformulated for the Roe-Turkel scheme in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF], with identical properties. This is likely related to the fact that the artificial speed of sound is defined from the same local preconditioner and that in modifying the pressure jumps, the fastest acoustic speed travels at almost infinite speed, with large consequences on the numerical stability. In particular, we see that the stringent stability condition that applies to the Roe-Turkel scheme has to be applied also to the artificial speed of sound in the asymptotic limit, with acoustics waves associated to µ + travelling at infinite speed of the order O( 1 M 2 ) while acoustic waves associated to µ -are slowed down to a reference flow velocity.

An important result for the stability in the low speed limit has been given by Birken and Meister, under an essential condition that the fastest wave speed of the dissipation matrix is of order O( 1 M 2 ) as M → 0 (see [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF], Lemma 4.1 and Theorem 4.2). The proof uses a subordinate matrix norm for the amplification matrix, given in one space-dimension with mesh spacing δx by

G(φ, σ) = Id + (cos(ξ) -1)σ D -i sin(ξ)σ A, ( 19
)
where Id is the identity matrix, D the dissipation matrix, A the flux Jacobian matrix, ξ the wave number (normalized by 1/δx), σ = ∆t δx and i 2 = -1. An essential mechanism of this proof relies in the fact that ρ( D) >> ρ( A) as M → 0, where ρ( D) and ρ( A) are the spectral radius of D and A, respectively. As ρ(A) = O( 1 M ), this condition is satisfied by the Rossow's artificial speed of sound with ρ( D) = µ + and result [START_REF] Rossow | A blended pressure/density based method for the computation of incompressible and compressible flows[END_REF]. Therefore the necessary Von Neumann condition is the same as for the Roe-Turkel scheme reformulated in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF], with for a fixed mesh

∆t ≤ h µ + O(hM 2 ) as M -→ 0, (20)
where h represents some characteristic cell distance and µ + is given by [START_REF] Rossow | A Flux-Splitting Scheme for Compressible and Incompressible Flows[END_REF]. With the above property 5., the Von Neumann condition [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF] returns the standard stability condition at the sonic line

∆t ≤ h |q n | + c as M -→ 1. ( 21
)
This necessary stability condition is well posed and allows the use of very large CFL numbers in the incompressible limit, provided that an implicit scheme is developed accordingly. Note that condition [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF] is much more restrictive than the stability condition used by Rossow, formulated by replacing the speed of sound by the artificial speed of sound in [START_REF] Dellacherie | Analysis of Godunov Type Schemes Applied to The Compressible Euler System at Low Mach Number[END_REF] with

∆t ≤ h |q n | + c O(h) as M -→ 0.
An implicit LU-SGS scheme was formulated for the steady-state problem from the backward Euler scheme considered in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF] for the implicit Roe-Turkel scheme, using a scalar or a matrix time-step. For the finite-volume scheme, the scalar local time-step is formulated for the fastest wave speed according to the CFL condition [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF] with

∆t cell = CFL h cell boundary µ + (22)
while the characteristic time-step matrix ∆t c is defined in the computational cell following requirements indicated in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF] (

∆t c ) cell =         cell boundary R∆t -1 L         -1 (23)
where R and L are the right and left eigenvector matrix of the flux Jacobian matrix with A = R diag(λ k )L and ∆t is the diagonal matrix ∆t = diag(∆t + , ∆t -, ∆t 0 , ∆t 0 )

with the characteristic time-steps

∆t + = CFL h µ * + , ∆t -= CFL h µ * - , ∆t 0 = CFL h µ * 0 . (24)
The star symbol indicates that the eigenvalues must be bounded away from zero and that µ ± must be modified to transition smoothly to the original |λ k | at the sonic line. We may simply define µ * ± as linear functions of the directional Mach number

µ * ± = (1 -M n )µ ± + M n |λ ± |.
The scalar time-step ( 22) and matrix R∆t -1 L for the matrix time-step ( 23) are evaluated at the cell interfaces using the Roe's average. The implicit scheme is stable for very large CFL numbers in the incompressible limit, with typically CFL O( 1M 2 ) in order to counterbalance the constraint of very small time steps ∆t O(M 2 ).

Rieper's low Mach number fix

In the case of the Rieper's fix, the corresponding matrix dissipation ( 16) has eigenvalues

             µ 0 = |q n | with multiplicity 2 for the 2D problem µ ± = c 2 1 + z ± (1 + z) 2 -4z(1 -M 2 n ) (25)
with the following properties:

1. The acoustic eigenvalues are positive in the subsonic range with µ + > 0 and

µ + µ -= c 2 z(1 -M 2 n ) ≥ 0 when M n ≤ 1; 2.
The spectral radius is given by max(µ 0 , µ -, µ

+ ) = µ + (see Appendix C); 3. However, µ + ≤ |q n | + c = ρ( A), ∀M n ≤ 1 (see Appendix C); 4.
At the sonic line, z = 1 and µ + = ρ( A); 5. The following inequalities hold:

µ -≤ |λ -| and µ -≤ |λ + |, ∀M n ≤ 1 (see Appendix C).
The right and left eigenvector matrices are explicitly derived in Appendix A, for the diagonalization in the symmetrizing variables d W0 variables. In the incompressible limit, with M n → 0, we see that

µ ± c 2 [(1 + z) ± |1 -z|] . With z ≤ 1 µ + c = O( 1 M ), µ -zc = |q n | + |q t | = O(1). ( 26
)
For the other extreme case, at the sonic line M n → 1 and

µ + = 2c, µ -= 0.
We note that:

-The Rieper's fix has no lower bound for acoustic eigenvalue µ -, with two vanishing eigenvalues µ 0 and µ -at a stagnation point. On the other hand, µ + has the same asymptotic behaviour as the wave speeds λ ± of the original Roe scheme, with the pressure jumps not altered by the fix;

-The fastest wave speed µ + is lower than ρ( A), with a behaviour of the same order of magnitude O( 1 M ) in the low-speed limit. Thus, the previous stability result according to Birken and Meister does not apply in theory to the rescaling of the Roe scheme using the Rieper's fix.

On one hand, it is expected that largest time steps can be used as µ + is one order smaller that the fastest wave speed found for the rescaling of the Roe scheme using the artificial speed of sound. However, as discussed next, a rigorous necessary condition for stability cannot be easily derived. On the other hand, the vanishing acoustic eigenvalue µ - may prevent from developing an efficient implicit scheme.

The main difficulty for the stability analysis of modified Roe-type schemes lies in the fact that the dissipation matrix is no longer an explicit function of the flux Jacobian matrix. This can be readily seen considering the onedimensional system. Assuming that matrices D and A have the same diagonalization basis (which is not true in the subsonic regime), with eigenvalues µ k and λ k respectively, then the eigenvalues of the amplification factor [START_REF] Rossow | Efficient Computation of Compressible and Incompressible Flows[END_REF] are given by

λ k G (φ, σ) = 1 + (cos(ξ) -1)σµ k -i sin(ξ)σλ k , and |λ k G (φ, σ)| 2 = (X k -1) 2 + Y 2 k , with X k = (1 -cos(ξ))σµ k , Y k = sin(ξ)σλ k .
After expanding and rearranging terms, we find

|λ k G (φ, σ)| 2 = 1 -2(1 -cos(ξ)) (1 -σµ k )σµ k + cos 2 ( ξ 2 )[(σµ k ) 2 -(σλ k ) 2 ] .
So a necessary and sufficient condition for the Von Neumann stability is

(1 -σµ k )σµ k + cos 2 ( ξ 2 )[(σµ k ) 2 -(σλ k ) 2 ] ≥ 0, ∀k, ∀ξ ∈]0, 2π[. ( 27
)
In the case of the original Roe scheme, µ k = |λ k | and we get the usual stability condition for the explicit scheme 1 -σ|λ k | ≥ 0, ∀k. Now assume that µ k >> |λ k |. Then a necessary and sufficient condition for the Von Neumann stability must be formulated with

1 -σµ k ≥ 0 or else ∆t ≤ δx µ k .
This assumption is actually not true for all k. This is only true for the spectral radius in the case of the Roe-Turkel scheme in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF], the preconditioned Lax-Friedrich scheme [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF] and the artificial speed of sound . This explains why using a matrix norm in [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF], a necessary condition could be formulated based on the spectral radius, with for all matrix norms ||M|| ≥ ρ(M). For the rescaling of Roe scheme with the Rieper's fix, we have in the incompressible limit µ + |λ + | c and O(1) µ -< |λ -| c with property 5. But in general, µ + may be larger or smaller than |λ + |, depending on the flow conditions at the cell interface, and furthermore, condition ( 27) is never satisfied for the very small or vanishing µ -. So a possible Von Neumann condition for stability, neither rigorously necessary nor sufficient, can be formulated considering that ( 27) is satisfied for k ∈ {0, +} in the incompressible limit with

∆t ≤ h µ + O(hM) as M -→ 0 ( 28
)
for a fixed grid spacing. Note that this condition returns the usual Von Neumann condition for the Roe scheme at the sonic line with the above property 4. With property 3., the usual stability condition for compressible flows yields a somewhat more restrictive CFL condition as

∆t ≤ h |q n | + c ≤ h µ + ∀M, with h |q n | + c O(hM) when M -→ 0. ( 29
)
Of course, the above discussion for the scalar case is only a conjecture, as the starting point is an assumption not valid in subsonic regime. However, this illustrates the difficulty in finding out a rigorous stability condition for the rescaling of the Roe scheme using the Rieper's fix. Evidence to support the practical validity of the standard CFL condition [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF] has been shown by many authors and has been also experienced by the author for the baseline explicit first-order scheme. However, considering some common stretched mesh or grid convergence studies brings into question this stability condition not only in the very incompressible limit for Mach numbers, but also in the transonic regime.

It turned out very difficult to develop a robust implicit scheme for the Rieper's fix. The main reason is the unbounded acoustic wave speed µ -, together with very small linear wave speeds µ 0 = |q n | in the incompressible limit. This results in a singular dissipation matrix and the LU-SGS implicit scheme for the steady-state problem is unstable whatever the formulation of the time-step, unless CFL numbers of the order of unity are used. This is not sufficient to overcome the stability condition with time steps ∆t O(M). Thus, the implicit scheme was formulated with the matrix coefficients of the original Roe scheme, i.e. setting z = 1 in the matrix coefficients of the LU-SGS scheme. So, the fix is not accounted for in the resulting implicit scheme.

Following the above discussion, the most restrictive stability condition [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF] was considered for the standard CFL condition for the Roe scheme

∆t cell = CFL h cell boundary (|q n | + c) (30)
and the characteristic time-step matrix was formulated for the original eigenvalues of the flux Jacobian matrix

∆t + = CFL h |λ * + | , ∆t -= CFL h |λ * -| , ∆t 0 = CFL h |λ * 0 | , (31)
where the star symbol indicates as previously that the eigenvalues must be bounded away from zero. However, it was experienced that the characteristic time-step matrix does not improve the damping properties of the numerical procedure in this case where the modified dissipation vector could not be taken into account implicitly. Additionally, CFL numbers of the order CFL O( 1 M ) could not always be used without enforcing the LU-SGS method for the solution of the implicit scheme, with some additional under-relaxation. The underlying reason also lies in the fact that the stability condition leading to the usual time-step formulation [START_REF] Xie | An accurate and robust HLLC-type Riemann solver for the compressible Euler system at various Mach numbers[END_REF] or [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF] is not properly formulated.

Asymptotic Analysis in the Incompressible Limit

The asymptotic analysis presented next allows to get a better understanding of the different behaviours characterizing both schemes in the very incompressible limit.

Preliminary discrete analysis

Expressions ( 7)-(Rossow) and ( 15)-(Rieper) for the dissipation vector can be easily rearranged according to the Weiss-Smith decomposition (4) and recast using the interface Mach number, with:

Rossow: δq n = (1 -|M 0 |) ∆p ρc + M 0 ∆q n , δp = M 0 ∆p + ρc (1 -|M 0 |)∆q n ;
Rieper:

δq n = (1 -|M 0 |) ∆p ρc + M 0 z∆q n , δp = M 0 ∆p + ρc(z -|M 0 |)∆q n ;
(Roe):

δq n = (1 -|M 0 |) ∆p ρc + M 0 ∆q n , δp = M 0 ∆p + ρc(1 -|M 0 |)∆q n . (32)
In the asymptotic limit, with c = O(1) and z = O(M) we get for the modified normal velocity:

Rossow: δq n 1 ρc ∆p = O(∆p), Rieper: δq n 1 ρc ∆p = O(M∆p). (33)
This coefficient for the pressure-difference term actually enforce the pressure-velocity coupling in low Mach-number flow regions, following a similar mechanism of the AUSM+(P) scheme developed for all-speed flows in [START_REF] Edwards | Low-Diffusion Flux-Splitting Methods for Flows at All Speeds[END_REF][START_REF] Mary | Deville An Algorithm for Unsteady Viscous Flows at all Speeds[END_REF][START_REF] Mary | Large Eddy Simulation of Flow Around an Airfoil Near Stall[END_REF], with the numerical flux

h = 1 2 U i (Q L + Q R ) + P ress - 1 2 |U i |(Q R -Q L ), (34)
where P ress denotes the centred pressure term and U i is an interface fluid velocity defined with a pressure stabilization term according to Rhie and Chow [START_REF] Rhie | Numerical Study of the Turbulent Flow past an Airfoil with Trailing Edge Separation[END_REF]

U i = 1 2 [(q n ) L + (q n ) R ] - c 2 ρ re f V re f (p R -p L ) (35)
with c 2 being a small constant and ρ re f , V re f are some reference density and velocity, respectively, for homogeneity purposes. The AUSM+(P) scheme is especially adapted to large eddy simulations, by monitoring the scalar dissipation, which can be locally switched-off in the absence of flow oscillations [START_REF] Mary | Large Eddy Simulation of Flow Around an Airfoil Near Stall[END_REF]. The numerical flux (34)-( 35) can be reformulated as the simple centred scheme with dissipation

h = 1 2 [F(W R ) + F(W L )] - 1 2 d
and the dissipation vector across the cell interface given by

d = |U d |(W R -W L ) + δq n Q + ∆p           0 0 |U d |           , where |U d | = |U i | + 1 2 ∆q n , δq n = c 2 ρ re f V re f ∆p and Q = 1 2 (Q R + Q L ) .
This is very similar to decomposition (4) for the Roe scheme. Note that for the AUSM+(P) scheme, δq n O(∆p) independently of the local Mach number. It has been especially shown in [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF]-(section 4) that the coefficient δq n is strongly related to the ability of the numerical scheme to damp out checkerboard pressure modes. As it can be seen from ( 33), the Rossow's artificial speed of sound has a similar behaviour to the AUSM+(P) scheme in the incompressible limit, since δq n 1 ρV re f ∆p with the asymptotic behaviour (9) of the artificial speed of sound. A similar behaviour for the Roe-Turkel scheme was also found in the analysis carried out in [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF]. On the other hand, δq n becomes very small for the Rieper's fix as indicated by [START_REF] Darmofal | The Importance of Eigenvectors for Local Preconditioners of the Euler Equations[END_REF], with the decreasing local Mach number. Actually the low Mach number fix does not introduced any modification to the δq n coefficient in the low-speed limit, and the asymptotic behaviour O(M∆p) is identical to the original Roe scheme. This explains why the Rieper's low Mach number fix may have transient undamped spurious pressure modes in the incompressible limit.

This mechanism is illustrated in Fig. 1 for the effect of a vanishing c 2 constant for the AUSM+(P) scheme, in the case of the non-lifting NACA0012 airfoil at M ∞ = 10 -3 , enforcing δq n → 0. The pressure checkerboard issue with decreasing dissipation when lowering values for c 2 is clearly visible for the Cp coefficient and the spurious entropy error, near the trailing edge, where the mesh is characterized by a strong refinement.

Regarding the pressure coefficient δp in [START_REF] Klein | Semi-Implicit Extension of Godunov-Type Scheme Based on Low Mach Number Asymptotics I: One-Dimensional Flow[END_REF], as observed by Li and Gu, the modifications in the dissipation matrix have to be found only in the contribution ρc∆q n when M 0 → 0, with effects on the non-physical behaviour limited to the momentum equation. For the Rieper's fix, this term is modified by the low Mach number fix and is vanishing with ρcz∆q n → 0 when M → 0. This is in accordance with the reduced dissipation characterizing the Rieper's fix in the incompressible limit, as mentioned earlier. On the other hand, this term is still present and modified with ρc ∆q n by the introduction of an artificial speed of sound. This is a common feature of most of the modified Roe-type schemes investigated in [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF].

These preliminary trends given by the discrete analysis are confirmed by the following more general framework analysis, considering the continuous problem and by the truncation error analysis performed in section 6.

Normalized Euler equations with matrix dissipation

In the following, we shall go into some more details of the analysis to get an insight into the effects of the dissipation matrix in the incompressible limit. We will consider the normalized Euler equations in the more appropriate primitive variables U for the analysis, with explicit artificial dissipation terms. We suppose the flow regular in the low-speed limit, with all flow variables sufficiently differentiable in space. Then the two-dimensional Euler equations can be formulated in quasi-linear form assuming locally constant matrix coefficients

∂U ∂t + A ∂U ∂x + B ∂U ∂y = h x 2 D x ∂ 2 U ∂x 2 + h y 2 D y ∂ 2 U ∂y 2 (36)
where A, B and D x , D y are the flux Jacobian matrix and the matrix-valued dissipation, respectively, in the (x, y) space directions, formulated for the primitives variables. The space increments h x , h y have been explicitly introduced as scale factors for homogeneity purposes. The normalized quantities (time, space, flow variables) used for the incompressible normalization are defined with

t = t v re f l re f , x = x l re f , ỹ = y l re f , ρ = ρ ρ re f , ũ = u v re f , ṽ = v v re f , c = c c re f , p = p p re f ,
where l re f is a reference length scale, ρ re f a reference density, v re f the reference speed, c re f a reference speed of sound and p re f = ρ re f c 2 re f a reference pressure. The reference velocity in independent of the reference speed of sound p re f /ρ re f . This ensures that the normalized speed velocity remains of order O(1) in the limit of a vanishing reference Mach number, defined in all the following with = v re f c re f .

So, in the normalization process, all above quantities are of the same order of magnitude, around unity. We will see that most of the asymptotic properties of the fully discrete schemes can be recovered formally considering explicitly the matrix valued dissipation in the continuous case. Therefore, the following analysis does not assume that a specific type of mesh is used (structured Cartesian/curvilinear or unstructured) or that the scheme depends upon a specific variable arrangement. System ( 36) is then expanded as follows

∂ρ ∂t + V • ∇ρ + ρ ∇.V = d ρ ∂V ∂t + ∇VV + ∇p ρ = d V ∂p ∂t + ρc 2 ∇.V + V • ∇p = d p (37)
where

d ρ , d V = [d u , d v ]
T and d p represent the dissipation vector for the mass, velocity and pressure equations, respectively. The resulting non-homogeneous normalized equations and following asymptotic analysis are discussed next for both schemes.

Rossow's artificial speed of sound

In the case of the artificial speed-of-sound approach, we get from [START_REF] Turkel | Assessment of Preconditioning Methods For Multidimensional[END_REF] the following matrix coefficients in primitive variables for the space directions (x, y)

D x =                                           |u| ρu c 0 r x c (1 - |u| c ) 0 c r x (1 - |u| c ) + |u| 0 u ρc 0 0 |u| 0 0 ρuc 0 cr x (1 - |u| c ) + |u|                                           , D y =                                           |v| 0 ρv c r y c (1 - |v| c ) 0 |v| 0 0 0 0 c r y (1 - |v| c ) + |v| v ρc 0 0 ρvc cr y (1 - |v| c ) + |v|                                          
with the directional ratio (12) in the (x, y) space directions

r x = c c x = 1 α 2 |u| c 2 + β 2 , r y = c c y = 1 α 2 |v| c 2 + β 2 .
The dissipation vector is then explicitly expanded with

d = [d ρ , d u , d v , d p ] T = h x 2 D x ∂ 2 U ∂x 2 + h y 2 D y ∂ 2 U ∂y 2 . ( 38
)
We shall recast system (37)- [START_REF] Wong | The Solution of the Euler Equations at Low Mach Number Using a Stabilized Finite Element Algorithm[END_REF] for the normalized quantities, with r x = 1 rx and r y = 1 ry , where rx = O(1) and ry = O(1) are the normalized ratios c c in the (x, y) space directions. Omitting the tilde symbol for sake of clarity, we then have successively for the normalized equations with dissipation terms:

Density: ∂ρ ∂t + V • ∇ρ + ρ ∇.V = h x 2 |u| ∂ 2 ρ ∂x 2 + ρu c ∂ 2 u ∂x 2 + 1 2 r x c (1 - |u| c ) ∂ 2 p ∂x 2 + h y 2 |v| ∂ 2 ρ ∂y 2 + ρv c ∂ 2 v ∂y 2 + 1 2 r y c (1 - |v| c ) ∂ 2 p ∂y 2 (39)
Velocity components:

ρ ∂u ∂t + ρu ∂u ∂x + ρv ∂u ∂y + 1 2 ∂p ∂x = h x 2 ρc 1 r x (1 - |u| c ) + |u| c ∂ 2 u ∂x 2 + 1 u c ∂ 2 p ∂x 2 + h y 2 ρ|v| ∂ 2 u ∂y 2 (40)
ρ ∂v ∂t + ρu ∂v ∂x + ρv ∂v ∂y + 1 2 ∂p ∂y = h x 2 ρ|u| ∂ 2 v ∂x 2 + h y 2 ρc 1 r y (1 - |v| c ) + |v| c ∂ 2 v ∂y 2 + 1 v c ∂ 2 p ∂y 2 (41)
Pressure:

∂p ∂t + ρc 2 ∇.V + V • ∇p = h x 2 ρuc ∂ 2 u ∂x 2 + c 1 2 r x (1 - |u| c ) + |u| c ∂ 2 p ∂x 2 + h y 2 ρvc ∂ 2 v ∂y 2 + c 1 2 r y (1 - |v| c ) + |v| c ∂ 2 p ∂y 2 .
(42)

Therefore, the solution for the normalized primitive variables are some functions of the reference Mach number and we look for asymptotic expansions in the incompressible limit

U = U 0 + U 1 + 2 U 2 + O( 3 ) as -→ 0, (43)
with the U j = U j (x, t), where x are the space coordinates and for the normalized variables, U j = O(1). In the following, we assume constant coefficients for the leading order "0-state" in the normalized system (39)-(42), according to the assumption of the quasi-linear form [START_REF] Rhie | Numerical Study of the Turbulent Flow past an Airfoil with Trailing Edge Separation[END_REF]. Introducing expansions (43) in equations ( 39)-(42) and collecting terms with equal power for the parameter , we find that for the leading order equations Order 1 2 : ∇p 0 = 0 (44)

Order 1 : ∇p 1 = 0 (45)
Order 1:

∂ρ 0 ∂t + V 0 • ∇ρ 0 + ρ 0 ∇.V 0 = h x 2 |u 0 | ∂ 2 ρ 0 ∂x 2 + r x c 0 ∂ 2 p 2 ∂x 2 + h y 2 |v 0 | ∂ 2 ρ 0 ∂y 2 + r y c 0 ∂ 2 p 2 ∂y 2 (46) ρ 0 ∂u 0 ∂t + (ρu) 0 ∂u 0 ∂x + (ρv) 0 ∂u 0 ∂y + ∂p 2 ∂x = h x 2 ρc 1 r x + |u| c 0 ∂ 2 u 0 ∂x 2 + h y 2 (ρ|v|) 0 ∂ 2 u 0 ∂y 2 (47) ρ 0 ∂v 0 ∂t + (ρu) 0 ∂v 0 ∂x + (ρv) 0 ∂v 0 ∂y + ∂p 2 ∂y = h x 2 (ρ|u|) 0 ∂ 2 v 0 ∂x 2 + h y 2 ρc 1 r y + |v| c 0 ∂ 2 v 0 ∂y 2 (48) d p 0 dt + (ρc 2 ) 0 ∇.V 0 = h x 2 (cr x ) 0 ∂ 2 p 2 ∂x 2 + h y 2 (cr y ) 0 ∂ 2 p 2 ∂y 2 . ( 49
)
So, with (44) and (45), this continuous asymptotic analysis shows that for the artificial speed-of-sound approach, the pressure is constant in space up to a fluctuation of order 2 and neglecting the higher order terms p(x, y, t) = p 0 (t) + p 1 (t) + 2 p 2 (x, y, t) = P 0 (t)

+ 2 p 2 (x, y, t), (50)
where the pressure disturbance field p 2 can be interpreted as the "incompressible" pressure with system (46)-( 49).

Taking into account that the leading order "0-state" is the surrounding thermodynamic state satisfying the ideal gaz law, then (ρc 2 ) 0 = γp O in the pressure equation (49). Therefore, the divergence of the leading-order velocity can be eliminated in the continuity equation ( 46), inserting its expression retrieved from (49), to finally give

D Dt ln(ρ 0 ) = 1 γ d dt ln(p 0 ) + h x 2 |u| ρ 0 ∂ 2 ρ 0 ∂x 2 + h y 2 |v| ρ 0 ∂ 2 ρ 0 ∂y 2 (51)
where D Dt = ∂ ∂t + V 0 .∇ is the material derivative along particle paths convected with velocity V 0 . The continuity equation (51) has a physical interpretation and describes a quasi-adiabatic compression of mass elements along particle paths, with artificial dissipation terms only depending on the surrounding density gradients. This result also has the noticeable consequence pointed out by Klein [32] that the divergence constraint derives from the pressure equation and not from the continuity equation. However, in the case of the artificial speed-of-sound approach, no divergence constraint can be found from the pressure equation (49). Even assuming no global compression for incompressible flows -1 γp 0 d p 0 dt = 0 equation (49) yields the following elliptic equation for the incompressible pressure

h x 2 1 ρc x 0 ∂ 2 p 2 ∂x 2 + h y 2 1 ρc y 0 ∂ 2 p 2 ∂y 2 = ∇.V 0 . (52)
This equation may be though as a local balance of incompressible pressure forces according to surrounding local flow compression or expansion. We see from the Order 1 system (46)-(49) that the rescaling of the Roe scheme following Rossow introduces stabilization terms and especially in the leading-order density and pressure equations, tightly coupled with the gradient of the pressure field p 2 . The situation is very different with the Rieper's fix.

Rieper's low Mach number fix

Similarly, in the case of the low Mach number fix, we get from ( 16) the following dissipation coefficients in primitive variables for the space directions (x, y)

D x =                                      |u| ρzu c 0 1 c (1 - |u| c ) 0 cz 0 u ρc 0 0 |u| 0 0 ρzuc 0 c                                      , D y =                                      |v| 0 ρzv c 1 c (1 - |v| c ) 0 |v| 0 0 0 0 cz v ρc 0 0 ρzvc c                                      .
As done previously, the non-homogeneous equations are recast for the normalized quantities, with for the low Mach number fix z = z, where z = O(1) is the normalized fix. We get (tilde symbol dropped):

Density: ∂ρ ∂t + V • ∇ρ + ρ ∇.V = h x 2 |u| ∂ 2 ρ ∂x 2 + 2 ρzu c ∂ 2 u ∂x 2 + 1 1 c (1 - |u| c ) ∂ 2 p ∂x 2 + h y 2 |v| ∂ 2 ρ ∂y 2 + 2 ρzv c ∂ 2 v ∂y 2 + 1 1 c (1 - |v| c ) ∂ 2 p ∂y 2 (53)
Velocity components:

ρ ∂u ∂t + ρu ∂u ∂x + ρv ∂u ∂y + 1 2 ∂p ∂x = h x 2 ρcz ∂ 2 u ∂x 2 + 1 u c ∂ 2 p ∂x 2 + h y 2 ρ|v| ∂ 2 u ∂y 2 (54) ρ ∂v ∂t + ρu ∂v ∂x + ρv ∂v ∂y + 1 2 ∂p ∂y = h x 2 ρ|u| ∂ 2 v ∂x 2 + h y 2 ρcz ∂ 2 v ∂y 2 + 1 v c ∂ 2 p ∂y 2 (55)
Pressure:

∂p ∂t + ρc 2 ∇.V + V • ∇p = h x 2 2 ρzuc ∂ 2 u ∂x 2 + 1 c ∂ 2 p ∂x 2 + h y 2 2 ρzvc ∂ 2 v ∂y 2 + 1 c ∂ 2 p ∂y 2 . ( 56
)
Looking for expansions (43) for the solution of system (53)-( 56), we find the following leading order equations Order 1 2 : ∇p 0 = 0 (57)

Order 1 : ∇p 1 = 0 (58)
Order 1:

∂ρ 0 ∂t + V 0 • ∇ρ 0 + ρ 0 ∇.V 0 = h x 2 |u 0 | ∂ 2 ρ 0 ∂x 2 + h y 2 |v 0 | ∂ 2 ρ 0 ∂y 2 (59)
ρ 0 ∂u 0 ∂t + (ρu) 0 ∂u 0 ∂x + (ρv) 0 ∂u 0 ∂y + ∂p 2 ∂x = h x 2 (ρcz) 0 ∂ 2 u 0 ∂x 2 + h y 2 (ρ|v|) 0 ∂ 2 u 0 ∂y 2 (60)
ρ 0 ∂v 0 ∂t + (ρu) 0 ∂v 0 ∂x + (ρv) 0 ∂v 0 ∂y + ∂p 2 ∂y = h x 2 (ρ|u|) 0 ∂ 2 v 0 ∂x 2 + h y 2 (ρcz) 0 ∂ 2 v 0 ∂y 2 (61)
d p 0 dt + (ρc 2 ) 0 ∇.V 0 = 0. ( 62
)
Then, for the Rieper's fix, the pressure is also constant in space up to a fluctuation of order 2 and is given by expression (50). As mentioned previously in the case of the artificial speed of sound, we see that the divergence of the leadingorder velocity can be readily eliminated from the pressure equation (62) and we find that the continuity equation ( 59) can be interpreted exactly with the same equation (51). However, there are major differences with previous results obtained for the artificial speed-of-sound approach. First, we see that the pressure equation (62) has vanishing artificial dissipation terms. This is a consequence of result (58) for the first-order pressure p 1 . Thus, assuming no global compression, equation (62) enforces a divergence-free leading-order velocity ∇.V 0 = 0. This result was obtained by Rieper within a pure discrete analysis [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF], assuming a Cartesian mesh. Second, comparing the dissipation terms of Order 1 systems (46)-( 49) and ( 59)-(62), we see that in the case of the low-Mach number fix, no additional stabilization terms involving the gradient of the incompressible pressure field p 2 are present in the Order 1 system. Then, this system can be formulated in a condensed form for the set of independent primitive variables [ρ 0 , u 0 , v 0 , p 2 ] T as follows

             ∂V 0 ∂t + V 0 • ∇V 0 = - ∇p 2 ρ 0 + ν 0 ∇2 V 0 ∇.V 0 = 0 (63) using the shorthand notation for convenience ∇2 V 0 = [ ∂ 2 V 0 ∂x 2 , ∂ 2 V 0 ∂y 2 ]
T and where ν 0 is a dissipation operator acting as a numerical viscosity, defined from the RHS of equations ( 60 13) for the normalized fix, we find rearranging terms

∇ 2 p 2 = h x 2 (ρ|v|) 0 ∂ 3 u 0 ∂x 3 + h y 2 (ρ|u|) 0 ∂ 3 v 0 ∂y 3 . (64)
Thus, system (63)-( 64) is a mixed elliptic (in space)-parabolic (in time) system, similar to the unsteady incompressible Navier-Stokes equations. It is well known from the incompressible theory that using a collocated arrangement for V 0 and p 2 variables on unstaggered grids, without specific stabilization, the resulting pressure field may be prone to checkerboard pressure modes. In other words, spurious pressure modes in the discrete solution cannot be ruled out. This is very much dependent on both the mesh definition and the reference Mach number. In [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF], undampted spurious pressure modes for the Rieper's fix are triggered for the cylinder problem, using a two-dimensional cellcentred triangular mesh with a collocated variable arrangement. It is demonstrated that this issue is the consequence of a pure discrete problem, when the number of edges of the computational cells is identical to the number of degrees of freedom. In our computations, spurious pressure modes exhibiting a checkerboard issue are also visible in the very low-Mach number limit, using highly stretched curvilinear structured meshes. On the other hand, using Cartesian meshes, no checkerboard issue was observed for the Rieper's fix, has discussed in [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF].

Some considerations for the acoustic time-scale

The acoustic normalization is defined with a different time scale, based on a reference speed of sound t = t c re f l re f . With this acoustic time scale, the time derivatives in the normalized systems ( 39)-( 42) and ( 53)-( 56) are modified with the derivative of the fast time variable t so that ∂ ∂t → 1 ∂ ∂t for all normalized primitive variables. As done previously, looking for expansions (43) for the solution of the modified normalized systems, we find that in both cases, the Order 1 2 system is unchanged with ∇p 0 = 0. However, the Order 1 and Order 1 systems are different from systems obtained for the incompressible scale (45)-( 49) or ( 58)-(62). The resulting systems are characterized by non-constant first-order pressure for both schemes. For the acoustic time-scale

p(x, y, t) = p 0 (t) + p 1 (x, y, t) + O( 2 ) (65)
and therefore, the acoustic pressure (65) is one order larger than the incompressible pressure (50), which is consistent with the Euler equations. It is well known that for the normalized Euler equations, p 1 is solution of a second-order hyperbolic equation describing the propagation of linear waves with speed c 0 , corresponding to the leading-order speed of sound [START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF][START_REF] Klein | Semi-Implicit Extension of Godunov-Type Scheme Based on Low Mach Number Asymptotics I: One-Dimensional Flow[END_REF].

Artificial Speed of Sound

In the case of the artificial speed of sound, a closed 5-equations system for (ρ 0 , u 0 , v 0 , p 0 , p 1 ) cannot be exactly found. We get from the Order 1 system the following relationship

∂p 0 ∂t = γ p 0 ρ 0 ∂ρ 0 ∂t
obtained by eliminating the dissipation terms function of the non uniform first-order pressure p 1 in the equations for ∂ρ 0 ∂t and ∂p 0 ∂t . This relationship can also be interpreted as ∂ ∂t ln(S 0 ) = 0 for the surrounding entropy S 0 , constant in time for this time scale. On the other hand, the Order 1 and Order 1 systems can be further combined to give the following non-homogeneous wave equation for p 1

∂ 2 p 1 ∂t 2 -c 2 0 ∇ 2 p 1 = ψ 1 + ψ 2 , (66)
with the forcing terms

ψ 1 = - h x 2 ∂ ∂t (r x |u|) 0 ∂ 2 p 1 ∂x 2 - h y 2 ∂ ∂t (r y |v|) 0 ∂ 2 p 1 ∂y 2 and ψ 2 = h x 2 ∂ ∂t (cr x ) 0 ∂ 2 p 2 ∂x 2 + h y 2 ∂ ∂t (cr y ) 0 ∂ 2 p 2 ∂y 2 . ( 67
)
This wave equation can be formulated for the pressure fluctuation p = pp 0 , and it can be readily seen that ψ 2 depending on p 2 is actually a one-order higher O( ) term in equation (66).

Low Mach Number Fix

Regarding the low Mach number fix, a resulting closed 5-equations system for (ρ 0 , u 0 , v 0 , p 0 , p 1 ) can be derived for the acoustic scale. We find that the Order 1 system has vanishing dissipation terms for all equations, with ∇p 0 = 0. This system especially yields

∂p 0 ∂t = 0 => p(x, y, t) = p 0 + p 1 (x, y, t) + O( 2 )
, and ∂ρ 0 ∂t = 0 for the acoustic scale, as there is no global pressure change on this time scale, which is in line with the Euler equations.

In the case of Rieper's fix, the non-homogeneous wave equation for p 1 does not introduce higher order terms involving p 2 in the forcing term, and combining the Order 1 and Order 1 systems, we find that

∂ 2 p 1 ∂t 2 -c 2 0 ∇ 2 p 1 = ψ 1 with ψ 1 = h x 2 c 0 ∂ ∂t ∂ 2 p 1 ∂x 2 + h y 2 c 0 ∂ ∂t ∂ 2 p 1 ∂y 2 . ( 68
)
Then, for the acoustic scale, the asymptotic analysis indicates that both schemes will also have a quite different acoustic component in the discrete solution. For the continuous Euler equations, the wave equation for p 1 is homogeneous and the forcing terms present in equations ( 66)-( 67) and (68) correspond to the contribution of the dissipation vector in the asymptotic systems. With this continuous analysis, it is difficult to compare the effects of different Order 1 systems and different forcing terms in the wave equation for the solution accuracy and the damping properties of the corresponding schemes. The solution of the Riemann problem commented in the following for the discrete problem gives more insight into the differences between the artificial speed of sound and the low Mach number fix.

Comments

Two incompressible limits exist for the pressure field, very different from a mathematical point of view, with p 2 solution of an elliptic equation for the slow incompressible time-scale, while p 1 is solution of an hyperbolic equation for the fast acoustic scale. As shown by Guillard et al within a functional analysis assuming isentropic [START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF] and barotropic flow [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF], in the theoretical framework used by Dellacherie [START_REF] Dellacherie | Analysis of Godunov Type Schemes Applied to The Compressible Euler System at Low Mach Number[END_REF] and Klein [START_REF] Klein | Semi-Implicit Extension of Godunov-Type Scheme Based on Low Mach Number Asymptotics I: One-Dimensional Flow[END_REF] in the development of a low Mach number extension of an explicit higher order shock-capturing scheme, both incompressible and acoustic limits are relevant and co-exist for the Euler equations in the incompressible limit. It has been especially demonstrated in [START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF] that solving the exact Riemann problem in the low speed limit, the pressure p (int) across artificial jumps at a cell interface between the right and left states has the following form for the isentropic flow

p (int) ρ re f c 2 re f = p 0 - 1 2 √ γρ 0 p 0 ∆u 1 + O( 2 ) = p 0 + p (int) 1 + O( 2 ) with ∆u 1 = (u 1 ) R -(u 1 ) L ,
even when the initial right and left pressure in the Riemann problem have the correct scaling (50). Therefore, artificial jumps of the velocity ∆u 1 in the discrete solution are responsible of pressure fluctuations of order one p (int) 1 , which has a pure acoustic origin, and these acoustic components remain present in the solution during the iterative process. In canceling out the jump of the normal velocity in the dissipation vector, the low Mach number fix aims at eliminating these pressure fluctuations of order one in the low Mach number limit, following the modified Godunov scheme formulated in [START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF], with

∆u 1 -→ 0 and p (int) ρ re f c 2 re f -→ p 0 + O( 2 ).
This result was a guideline reported by Rieper in the development the low Mach number fix. Thus, part of the solution corresponding the acoustic time-scale is damped out using the low Mach number fix. The artificial speed of sound does not modify the jump of the normal velocity and then, the discrete solution should contain visible acoustic components in the numerical procedure. On the other hand, the incompressible component in the discrete solution may lead to well known drawbacks explained previously, when enforcing the discrete divergence constraint. Practical aerodynamic computations in the very low-speed limit show that the flowfield undergoes a very slow inertial motion corresponding to a slow incompressible component, while the fastest wave speed travels at almost infinite speed. The solution is also characterized by fast acoustic components, as mentioned above. This has consequences on the condition number of the Jacobian matrix and the Von Neumann stability condition, but also for the solution accuracy. The contribution of the incompressible pressure p 2 in the discrete solution also indicates a trend of the hyperbolic system of conservation laws in the incompressible limit, characterized by an almost instantaneous wave propagation from the boundaries into the computational domain, which is close to an elliptic behaviour: in a pure elliptic problem, the solution is a function depending explicitly on the Neumann/Dirichlet boundary conditions everywhere inside the bounded domain of integration.

Truncation Error Analysis

The following truncation error analysis is only presented for the one-dimensional system (1) in divergence form, discretized in a uniform Cartesian grid with nodes x i , and grid spacing δx. Same conclusions are obtained for the multi-dimensional system. Assuming subsonic flow conditions, the eigenvalues of the Roe matrix keep the same sign from x i-1 to x i and from x i to x i+1 , and the truncation error in space for the Roe scheme at point x i is given by

i = δx 2 | A(W)| i ∂ 2 W ∂x 2 i + O(δx 2 ).
In other words, the first-order Roe scheme solves at the second order in space the "modified" equations

∂W ∂t + ∂F(W) ∂x = δx 2 | A(W)| ∂ 2 W ∂x 2 .
Considering now the case of the modified dissipation matrix according to Rossow and Rieper, the truncation error depends on the entries of modified matrix-valued dissipation D(W) formulated in conservative variables, with 1. For the first-order rescaled Roe schemes explored in this paper, the truncation error only increases in the energy equation for a fixed mesh as M → 0, with the specific total enthalpy becoming very large H O( 1 M 2 ) in the incompressible limit. So, the proper scaling of the pressure field requires grid refinement or grid adaption with δx → 0, and possibly the use of some MUSCL extrapolation, thus reducing the truncation error, as shown in the last section 7; 2. The Rieper's low Mach number fix has a smaller truncation error when M → 0 for a fixed mesh and fixed Mach number than the Roe-Turkel scheme and the Rossow's artificial speed of sound approach. As the fix does not modify the pressure jumps, the truncation error in the energy equation is of the same order of magnitude as the Roe scheme.

Results for the Steady-State Problem

Two baseline test-cases using structured meshes were considered to illustrate the above analysis, the cylinder problem and the NACA0012 airfoil for non-lifting low Mach numbers flows. For both test-cases, a grid convergence study was carried out and effects of grid stretching on the discrete solutions in the very low speed limit were also investigated. Results are discussed in detailed and are compared to the potential flow theory.

In the incompressible limit, we shall especially consider the proper behaviour of the pressure field

p(x, y) = p 0 + p 2 (x, y)M 2 ∞ as M ∞ -→ 0, ( 69
)
where p 0 is a reference surrounding pressure and M ∞ is the prescribed inflow Mach number. Therefore, reformulating the normalized pressure field with the following quantity p(x, y) = p(x, y)p min p max -

p min = p 2 (x, y) -p 2 min p 2 max -p 2 min = p2 (x, y) ∈ [0, 1] as M ∞ -→ 0, ( 70
)
this expression returns the normalized disturbance pressure field p2 (x, y), which should be independent on the inflow Mach number. The proper behaviour of the pressure field in the incompressible limit can also be displayed with the amplitude of the pressure disturbance versus the inflow Mach number

δ p = p max -p min p max p 2 max -p 2 min p 0 M 2 ∞ as M ∞ -→ 0, ( 71
)
which should be a quadratic function δ p(M ∞ ) of the inflow Mach number. The above quantity δ p only indicates the consistency of the pressure field, i.e. with the correct amplitude of the disturbance pressure, not the overall accuracy of the computations for the pressure field. As indicated in the following, looking at in more details the pressure coefficient, it will be shown that the solution accuracy depends on both the mesh definition and the inflow Mach number in the very incompressible limit, as indicated by the truncation error analysis. A simple indicator for accuracy is the pressure coefficient at the wall, which returns, according to the asymptotic analysis

C p (x, y) = p 2 (x, y) 1 2 ρ 0 c 2 0 as M ∞ -→ 0. ( 72
)
Therefore, the C p coefficient is the non-dimensionalized disturbance pressure at the wall, which should be also independent on the inflow Mach number in the asymptotic limit. Some authors have reported the global cut-off issue when the preconditioned scheme requires some reference Mach number to avoid a singular dissipation matrix with a vanishing Mach number. It should be reminded that in the case of the Roe-Turkel scheme formulated in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF] and the preconditioning used in [START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF], the use of a reference Mach number is not mentioned in the global preconditioning parameter and in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF], the reference Mach number of the Roe-Turkel scheme was set to the zero-level machine. This also the case in all low Mach-number flow computations presented in this paper, where in the definition (8) of the artificial speed of sound, β 2 = min(max(M 2 , 2 ), 1), where 2 was set to the zero-level machine. The Rieper's fix does not use any cut-off parameter, although in our implementation, the fix ( 13) was formulated with z = min(1, max(M, )) to avoid a vanishing acoustic eigenvalue at the stagnation point, with also set to the zero-level machine.

Non-lifting cylinder

For this elementary test-case, different types of mesh were considered in order to illustrate the asymptotic behaviour and the solution accuracy for a decreasing inflow Mach number M ∞ in the range [10 -2 , 10 -3 , 10 -4 , 10 -5 ]. Converged solutions with such low Mach number flow conditions could only be obtained with the stability conditions discussed in section 4, using the implicit LU-SGS scheme and high CFL numbers. The meshes used are described in table 1, for basically 3 types of structured mesh: Cartesian, an irregular curvilinear mesh and an adapted curvilinear mesh, represented in Fig. 2. The irregular mesh is characterized by a grid refinement in the circumferential direction, near the stagnation point downstream of the cylinder only. This yields a highly stretched mesh designed to trigger checkerboard pressure modes. An adapted mesh for solution accuracy was also generated with grid refinements in the radial and circumferential directions, at the upstream and downstream stagnation points. This adapted mesh has also a larger extension for the outer boundary. For all meshes, a slip condition with zero-order pressure extrapolation is applied at the solid wall. At the outer boundary, a non-reflecting boundary condition was used. Note that the solution of the potential theory was not imposed as external state at the outer boundary.

Disturbance pressure consistency

A grid convergence study was carried out for the Cartesian mesh, with 4 different mesh densities indicated in table 1. Such Cartesian meshes with almost quadrilateral cells and an aspect ratio of one within the mesh are not adapted at the wall to inviscid flow computations, with spurious entropy error generated near the stagnation points. However, the mesh refinement effect provides an indicator of the consistency of the pressure field with the inflow Mach number, for the amplitude of the disturbance pressure (71). The behaviour of δ p with the inflow Mach number and the mesh density is represented in the Fig. 3 in logarithmic scales, comparing the results for the modified Roe schemes according to Rieper, Rossow and the Roe-Turkel scheme for the finest mesh. From these plots, it can be observed that:

1. For all schemes, the coarse and medium mesh 1 are sufficient to reproduce the proper quadratic behaviour of the pressure field (69) up to Mach number M ∞ = 10 -3 . For lower Mach numbers, these meshes are not sufficiently refined. This is especially the case of the Rieper's fix, for which even the finest mesh is not sufficient to capture the proper scaling of the pressure disturbances. On the other hand, for the modified Roe scheme according to Rossow, the medium mesh 2 is fine enough to recover the proper scaling up to Mach number M ∞ = 10 -4 and only the finest mesh allows to get a consistent behaviour for δ p at M ∞ = 10 -5 (green symbols). Identical result is obtained in the finest grid for the Roe-Turkel scheme (bottom figure). Note that for the analytical potential flow solution, we get δ p = 2.799M 2 ∞ . This constant slope in the logarithmic scale is almost exactly reproduced in the finest Cartesian grid by the artificial speed of sound and the Roe-Turkel scheme. 2. The artificial speed of sound has a faster mesh convergence for the consistency of the pressure field, but these results indicate nothing else on the solution accuracy. This will be illustrated next. In particular, this result only indicates that the Rieper's fix would required finer meshes of some adapted mesh to recover the proper scaling of the pressure field in the very incompressible limit of Mach numbers smaller than 10 -3 . This is a consistency issue in computing the disturbance pressure field p 2 when the Mach number is lowered to the incompressible limit. Using the adapted mesh, the pressure field for the Rieper's fix almost recovers the proper scaling up to M ∞ = 10 -5 , as indicated with the "+" symbols for this adapted mesh, and the proper scaling is exactly reproduced using a third-order MUSCL extrapolation in this mesh ("circle" symbols). The artificial speed of sound also returns the proper behaviour of pressure field in the very low speed limit in this adapted mesh, with or without MUSCL extrapolation.

So, this first result for the mesh convergence study shows that the mesh must be in some way adapted or refined in the very incompressible limit, to get the consistent quadratic behaviour of the pressure field.

Checkerboard pressure modes

Many preliminary tests conducted during this work have clearly shown the occurrence of checkerboard pressure modes in highly stretched grids for the Rieper's low Mach number fix. This is shown in Fig. 4, with results illustrating the behaviour of the normalized pressure field (70) with the inflow Mach number, in the irregular curvilinear mesh. As explained by the asymptotic analysis in section 5, the Rieper's fix suffers from pressure-velocity decoupling, by enforcing a divergence-free leading-order velocity ∇.V 0 = 0, with a collocated arrangement of pressure and velocity components in the finite-volume scheme. This result was demonstrated first by Guillard et al. [31]. This issue is especially visible for the lowest Mach numbers M ∞ < 10 -3 , in areas where the mesh undergoes a higher grid stretching, and in the case M ∞ = 10 -5 no converged solution could be obtained using this mesh. This also results in inconsistent contours for the normalized pressure, which should be independent on the inflow Mach number. In contrast, the artificial speed of sound approach is not sensitive to the mesh definition, with no checkerboard issues observed, and consistent contours for the normalized pressure identical for all Mach numbers. This is also the case of the Roe-Turkel scheme, being free of any checkerboard pressure modes, as also shown in [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF].

Again, these results do not show that the solution obtained with the artificial speed of sound approach is accurate in this mesh. It only tells us that the pressure field is consistently reproduced. This can be seen with the comparison to the analytical solution of the potential theory (using the Bernoulli theorem for the pressure), plotted in the same mesh with the grey solid line in the background of figures 4 for the normalized pressure field (Theory in the legend). It must be observed that the solution for the pressure in Fig. 4-Right using the artificial speed of sound has large discrepancies with the potential theory, with a significant loss of symmetry upstream and downstream in this unadapted mesh, while the same solutions in Fig. 4-Left using the Rieper's fix are globally more symmetrical and match somewhat better the potential theory in the field, although strongly oscillating. So the consistency of the pressure field according to the asymptotic analysis is only a necessary condition, and not in any case a sufficient condition, for accuracy. This important point is discussed in more details in the following. with the C p distributions indicated in previous Fig. 6. Both original first-order schemes still have inconsistencies, with different contours of p2 according to the Mach number. The solutions for the pressure field are not perfectly symmetrical and discrepancies can still be seen with the analytical potential flow solution. Using the third-order MUSCL extrapolation leads to a dramatic improvement on the solution accuracy for the pressure field, for both schemes. With the Rieper fix, except the extreme case M ∞ = 10 -5 not shown in the figures, the computed p2 is nearly identical to the potential flow solution (Fig. 7b) and inconsistencies are no longer visible, with identical p2 contours for M ∞ = 10 -3 and M ∞ = 10 -4 . With the MUSCL extrapolation applied to the artificial speed of sound, contours of p2 are also symmetrical and very close to the potential theory for this test-case, for all Mach number up to M ∞ = 10 -5 (Fig. 7d). Numerical solutions obtained with uniform incoming flow conditions when M ∞ → 0 correspond to isentropic flows for this test-case. However, spurious entropy in the discrete solutions is still visible downstream of the cylinder at M ∞ = 10 -3 . This is illustrated in the next figures 8 for the effect of the MUSCL extrapolation in the adapted mesh. Comparing the left figures 8a and 8c, corresponding to the solution of the first-order schemes, it can be seen that the low Mach number fix has much less error than the artificial speed of sound. This was shown by the truncation error analysis. The expected effect of the MUSCL extrapolation has also a dramatic effect on the reduction of spurious entropy shown in Fig. 8b and Fig. 8d, as it could be observed with the previous plots for p2 . In the case of the low Mach number fix, the spurious entropy is suppressed within the entire flowfield, while the artificial speed of sound still has visible residual spurious entropy only generated at the downstream stagnation point.

Non-lifting NACA0012 airfoil

Two types of mesh were considered for this test-case, with a decreasing inflow Mach number M ∞ in the range [10 -1 , 10 -2 , 10 -3 , 10 -4 ]. The meshes used are described in table 2 and represented in Fig. 9 (identical to the meshes used in a previous work [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF] for the Roe-Turkel scheme). The same boundary conditions as previously are applied to this test-case. The family of Vasseberg-Jameson Cartesian-like meshes, defined with O-meshes generated from a conformal transformation developed in [START_REF] Vassberg | Jameson In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF], were first considered for a grid convergence study. The finest mesh used is a high-density grid with more that 1 million grid points. The C-mesh is a standard curvilinear mesh with grid refinement at the leading edge. This standard curvilinear mesh triggers spurious checkerboard pressure modes for the Rieper's fix arising from inflow conditions corresponding to very low Mach numbers. This test-case also demonstrates that transient errors persist in the solutions at the stagnation point for both schemes, when the inflow Mach number is lower than M ∞ = 10 -3 . In all the following figures, solutions obtained with the O-meshes will be referred to "VJ Mesh" and "Stretched Mesh" when the C-mesh was used. Solutions for the C p coefficient are compared to the potential flow solution, which was computed using a panel method. As done previously for the cylinder problem, the consistency of the solutions in the incompressible limit with expressions (70), ( 71) and ( 72

Disturbance pressure consistency

The behaviour of δ p with the inflow Mach number is represented in the next Fig. 10 for both schemes and both types of mesh. Note that the coarse O-mesh is sufficient to recover the proper quadratic behaviour of the pressure field according to Eq. (71), up to M ∞ = 10 -4 . Identical results are obtained for the low Mach number fix and the artificial speed of sound. Again, this only indicates that the consistency of the disturbance field is achieved in these meshes for both schemes, and for all inflow Mach numbers tested. This is further confirmed with the next Fig. 11, illustrating the behaviour of the normalized pressure field with the inflow Mach number M ∞ → 0 , in the curvilinear stretched mesh. We see that both schemes are in perfect accordance with expression (70) for the disturbance pressure p2 , with identical contours in the incompressible limit up to M ∞ = 10 -4 . Additionally, in opposite to the previous test-case using the highly stretched irregular mesh, no pressure-velocity decoupling is visible in the disturbance pressure field for the low Mach number fix. With the standard C-mesh used for this test-case, the grid stretching is not that high than the irregular mesh used in the case of the cylinder. However, the pressure-velocity decoupling is still present in the solution and especially visible at the solid wall, as it is shown next.

Checkerboard pressure modes and solution accuracy

We first look at the behaviour of the pressure coefficient for the fixed C-mesh with M ∞ → 0 for both schemes. In the case of the artificial speed of sound, plots in Fig. 12-Left shows that the C p distributions are identical when M ∞ = 10 -1 to M ∞ = 10 -3 , with close results compared to the panel method. This is expected from result (72) in the incompressible limit. When M ∞ = 10 -4 , the C p distribution is shifted, with a large error at the stagnation point, as explained for the cylinder test-case. In any case, the solution is free of any checkerboard pressure modes, which is again in line with the asymptotic analysis. For the low Mach number fix in Fig. 12-Right, similar results for the C p distribution are found, except that first, a larger error at the stagnation point can be observed when M ∞ = 10 -4 , corresponding to a larger shift of the pressure coefficient and second, a pressure-velocity decoupling issue is visible at M ∞ = 10 -3 , and is further amplified at M ∞ = 10 -4 . Here, pressure checkerboard modes are triggered at the trailing edge, where the grid is characterized by a strong refinement. This is very similar to the behaviour of the AUSM+(P) scheme in the same mesh, where the mechanism of checkerboard pressure modes was illustrated in Fig. 1. On the other hand, looking at in more details these results, it can be observed that, except the extreme case with M ∞ = 10 -4 , the C p distribution obtained with the low Mach number fix is more accurate and better fits the potential theory than the artificial speed of sound. This is especially noticeable at the leading edge, where the stagnation pressure returns C * p = 1 with the low Mach number fix, which is not exactly achieved with the artificial speed of sound. This is again in accordance with the truncation error analysis. However, considering the extreme case with M ∞ = 10 -4 , a much larger error at the stagnation point is found using the low Mach number fix compared to the artificial speed of sound (red lines in the figures). This is due to lack a dissipation characterizing the low Mach number fix in the very incompressible limit, as illustrated in the following for a grid convergence study.

The stagnation pressure error arising for the lowest Mach number M ∞ = 10 -4 tested in these computations, could be only suppressed upon grid refinement. This is especially the case with the artificial speed sound. The effect of grid refinement using the Vasseberg-Jameson meshes was considered for the fixed inflow Mach number M ∞ = 10 -4 , as shown in Fig. 13. In the case of the artificial speed of sound in Fig. 13-Left, the shift in the C p distribution is reduced as the mesh density is increasing, with a diminishing spurious entropy error. Using the fine 513 × 513 and high density 1025 × 1025 meshes, the pressure coefficient matches much better the potential theory, reducing the error at the stagnation point. With the low Mach number fix, the grid convergence in Fig. 13-Right is slower, as observed with the cylinder test-case. However, the stagnation point error in only reduced up to the fine mesh. Using the high-density mesh, it could be observed that the stagnation pressure error increases with iterations, although the spurious entropy is further reduced at the leading edge in this mesh and the computation has no convergence issue. This shows that the discrete solution has growing errors at the stagnation point at large times using extremely refined meshes in the incompressible limit. Errors are also characterized in the transient phase with spurious lift produced by the scheme, even using coarser meshes. According to numerous numerical experiences carried out by the author, using different types of meshes, this transient spurious lift can usually only be damped out after a very large number of iterations, indicating that the Low Mach number fix has insufficient dissipation to efficiently overcome the mechanism of transient error growth pointed out by Darmofal. It was shown with the cylinder test-case that using a third-order MUSCL extrapolation could improve the accuracy of the original first-order schemes by drastically reducing the truncation error. This was also observed for this testcase using the stretched C-mesh. Results for the C p distribution are shown in the next Fig. 14 and have to compared with same results in Fig. 12 obtained for the first-order schemes. For both schemes, the pressure coefficient perfectly agrees with the potential solution up to M ∞ = 10 -3 and therefore, the residual stagnation point errors are significantly reduced. The corresponding C p distributions even better match the potential theory than using the high density Omesh. However, when M ∞ = 10 -4 , errors at the stagnation are not sufficiently reduced, and the pressure coefficient is only shifted downwards with the artificial speed of sound and is just reduced with the low Mach number fix. This further confirms that the critical issue of the accuracy at the stagnation point can only be overcome upon using extremely high-density grids in the very incompressible limit. Also note that checkerboard pressure modes are reduced for the low Mach number fix using a MUSCL extrapolation, but not suppressed.

Acoustic-incompressible interactions for the non-lifting NACA0012 airfoil

The following considerations regarding acoustic-incompressible interactions in the discrete solutions are obviously not configuration dependent and are just discussed here for the NACA0012 airfoil. The presence or not of the unsteady acoustic components p 1 , can be detected, even for the steady-state problem, looking at first the behaviour of the convergence history in the incompressible limit. Wave equations ( 66)-( 67) and (68) previously derived for the continuous problem clearly indicate that the acoustic pressure is present in the discrete solutions, for both schemes.

The convergence history for the non-lifting NACA0012 airfoil is compared for both schemes in the next figures, for M ∞ = 10 -1 , M ∞ = 10 -2 and M ∞ = 10 -4 , using the usual stretched C-mesh (Fig. 15(a)-(b)-(c)) and the fine O-mesh 513 × 513 (Fig. 15(d)). With the implicit scheme, large CFL numbers were used as indicated in the figures, following the stability analysis described in section 4. Von Neumann conditions [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF]

(using CFL = 1 M 2 ∞
) and ( 29)

(with CFL = 1
M ∞ ) were applied to the artificial speed of sound and the low Mach number fix, respectively. Residuals are normalized with the maximum residual found during iterations. As illustrated in the figures, a significant shift with the zero-level machine can be observed, with a noticeable increasing shift as M ∞ → 0. All solutions actually reach a limit cycle, which is usually not observed for compressible or nearly incompressible flows, where residuals for the steady-state problem can be dropped to the zero-level machine. Comparing Fig. 15(c) and Fig. 15(d), it can be observed that the same limit cycle is reached M ∞ = 10 -4 in the C-mesh and the fine O-mesh, which has a much higher grid density. Limit cycles were also found using lower values for the CFL number. This is not a matter of numerical stability issue, mesh definition or insufficient damping of the numerical dissipation. The reason for this limit cycle was not very well understood in a previous work [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF]. It was believed in the case of the Roe-Turkel scheme that the numerical procedure could not damp pressure disturbances of the order of O(M 2 ), when such disturbances become round-off errors for the machine precision. This could also be explained for the artificial speed of sound, with the presence of the disturbance pressure p 2 in the Order 1 system (46)-(49) at the incompressible time scale.

To illustrate this, we give in table 3 typical outputs for the min and max non-dimensional pressure p p 0 in the entire mesh (identical outputs could be extracted for the cylinder test-case). Both schemes clearly capture a pressure disturbance corresponding to the incompressible pressure at the second order O(M 2 ∞ ) when M ∞ → 0. The min and max non-dimensional pressure are represented with 15 digits, corresponding to the machine precision for double precision computations. With the artificial speed of sound, a perfect correlation can be found for all decreasing Mach numbers between the remaining significant figures for the representation of the second-order term p 2 p 0 M 2 ∞ , and the order of magnitude of the residual drop where the limit cycle is reached. However, it turns out that this is no longer the case with the low Mach number fix, for which the residual drop in much more significant. No straightforward correlation between the remaining significant figures and the limit cycle can be found as indicated in the table and also illustrated in Fig. 15. Furthermore, the low Mach number fix has less dissipation than the artificial speed of sound. This suggests that the occurrence of a limit cycle for the steady-state problem originates from a higher frequency content in the discrete solution computed with the artificial speed of sound, as it is also the case with the Roe-Turkel scheme. Such high frequencies in the discrete solutions, not damped by the implicit scheme, may result from spurious acoustic pressure waves, which can be only be displayed from unsteady computations, considering the dilatation rate, as shown next.

Unsteady computations were performed using a simple implicit Backward Euler scheme, first-order accurate in time, in the fine O-mesh 513 × 513 for the inflow Mach number M ∞ = 10 -4 . Note that the following unsteady results do not have any physical meaning for the displayed acoustic field present in the discrete solutions. For the steady-state problem, acoustic sources of pure numerical origin may only by justified by an oscillating pressure at the stagnation point or accross a contact discontinuity at the trailing edge. Physically relevant acoustic waves can only result from the interaction of turbulent boundary layers occurring at the trailing edge for incompressible flows. Therefore, an arbitrary time-step ∆t = 0.01s was used, since no physical frequency can be retrieved from these inviscid flow simulations assuming a steady-state solution. This arbitrary time-step is by far much larger than the minimum of the local time steps allowed by the explicit scheme. However, the purpose of these unsteady computations is to only illustrate the different acoustic components captured by the artificial speed of sound and the low Mach number fix, using a limited number of unsteady iterations. From the unsteady computations, it can first be noted that the incompressible disturbance pressure p2 is well

M ∞ p min /p 0 p max /p 0 Residual drop Residual drop (Rossow & Rieper) (Rossow & Rieper ) Rossow
Rieper M ∞ = 10 -1 0.990430049136692 0.999924485661354 13 orders of magnitude 14 orders of magnitude M ∞ = 10 -2 0.999904126570691 0.999999293524798 11 orders of magnitude 13 orders of magnitude M ∞ = 10 -3 0.999999043260276 0.999999994827128 9 orders of magnitude 12 orders of magnitude M ∞ = 10 -4 0.999999990618559 1.000000000134090 7 orders of magnitude 11 orders of magnitude reproduced, for both schemes. Plots represented in the right Fig. 16-Left show that the p2 contours are almost identical for the steady-state and the unsteady problems, with the proper scaling as indicated in Fig. 16-Right. The p2 contours become independent of the unsteady iterations in the large times, indicating that the incompressible component in the discrete solutions is also consistently reproduced at the second order O(M 2 ∞ ) for the unsteady problem. Although already visible for the original first-order schemes, the presence of an unsteady component in the discrete solutions is better displayed increasing the accuracy in space of both schemes formulated with the thirdorder MUSCL extrapolation. This is illustrated in the next plots in Fig. 17 for the dilation rate in the farfield, at different unsteady time steps, using high saturation in the definition of the ∇.V contour levels. The acoustic component corresponding to the artificial speed of sound in the left column is characterized by an acoustic pulse for the short times, propagating in the farfield (t=30s), and then reflected back from the farfield boundary into the computational domain (t=60s). This acoustic pulse may be related to oscillating spurious pressure forces found in the short times in the discrete solution for both schemes, due to numerical instabilities. In addition, the farfield boundary being not an absorbing boundary condition for acoustic waves, an unsteady cycle can be observed until the acoustic pulse is dissipated in the large times. As expected, with the cancellation of the jump of the normal velocity, such an acoustic pulse is not reproduced with the low Mach number fix for the short times (t < 60s), as shown by the same plots in the right column.

However, in the large times, the discrete solution computed with the low Mach number fix exhibits a spurious acoustic component downstream, very slowly damped out by this low dissipation scheme, and not observed with the artificial speed of sound. It is expected that this acoustic component should be dissipated in the time limit where the non-lifting steady-state solution should be recovered. Nevertheless, looking at in details the patterns for the dilatation rate in the farfield at large times, it can be seen that both schemes have permanent acoustic disturbances in their discrete solutions, although with a very small intensity. Actually, in the large times, the history of these acoustic disturbances is clearly visible in the farfield, as it can be seen in the bottom figures. This may explained the limit cycle characterizing the convergence history, previously illustrated in Fig. 15 for both schemes and the consistency issues encountered with the low Mach number fix in very fine meshes, which require a very large number of iterations to damp spurious pressure modes.

These results illustrate a complex mechanism with acoustic-incompressible interactions, although of pure numerical origin for this inviscid flow computation, when a compressible scheme is used in the very incompressible limit.

Conclusions

The framework analysis developed in this paper considering the inviscid flow assumption was thoroughly supported by grid convergence studies and grid stretching effects in the low speed limit, for very low inflow Mach numbers M ∞ ≤ 10 -3 . Up to M ∞ = 10 -3 , the artificial speed-of-sound according to Rossow and the Rieper's low Mach number fix provide a consistent rescaling of the Roe scheme. Solutions were found in good agreement with the potential flow theory, for the two test-cases considered, provided that an adapted mesh is used, especially near the stagnation points. It has been shown numerically that using a third-order MUSCL extrapolation, an almost perfect agreement with the potential flow theory could be obtained. However, important differences supported by the asymptotic analysis could be illustrated numerically in the very incompressible limit, at M ∞ = 10 -4 for the NACA0012 airfoil and up to M ∞ = 10 -5 for the cylinder problem.

The low Mach number fix was designed to reduce the numerical dissipation by canceling the jump of the normal velocity at cell interfaces, corresponding to a reduced scheme truncation error. This may be an issue for inviscid flow computations, with transient errors insufficiently damped by the numerical procedure. This was especially observed using very fine meshes and Mach numbers M ∞ < 10 -3 , with a numerical solution characterized by large errors at stagnation points. Additionally, with a vanishing acoustic speed at a stagnation point, the Von Neumann stability condition is not well posed, although close to the standard Von Neumann condition for compressible flows. On the other hand, lowering the dissipation is an advantage for the simulation of more complex turbulent flows, as illustrated with the L2 Roe scheme described in [START_REF] Oßwald | a Low Dissipation Version of Roe's Approximate Riemann Solver for Low Mach Numbers[END_REF], where the dissipation was even further reduced, with the fix also applied to the jump of the tangential velocity. However, an asymptotic analysis conducted in this paper for the continuous problem, explicitly accounting for the matrix dissipation, clearly shows that the Rieper's fix is intrinsically prone to pressure-velocity decoupling, regardless of the mesh used, structured or non-structured. This is a consequence of enforcing the divergence-free constraint. Schemes enforcing the divergence-free constraint of the leading-order velocity, using a collocated arrangement of the velocity components and the pressure in a cell-centred finite-volume discretization, are not consistent with the incompressible theory without specific requirements. This main drawback of the low Mach number fix was illustrated for the cylinder problem and the NACA0012 airfoil considered in this work, at very low inflow Mach numbers, using highly stretched structured meshes.

The concept of artificial speed of sound was thoroughly validated by Rossow in [START_REF] Rossow | Efficient Computation of Compressible and Incompressible Flows[END_REF], addressing the analytical stiffness for low Mach number flow, the discrete stiffness for Reynolds number effects and mesh convergence. Our contribution with this paper was to emphasize the asymptotic properties, not given by Rossow, for this rescaling of the Roe scheme aiming at increasing the amplitude of the pressure jump in the matrix-valued dissipation. It has been pointed out that the artificial speed of sound is very similar to the Roe-Turkel scheme in terms of asymptotic properties, with however a reduced numerical dissipation in the incompressible limit. It has been particularly demonstrated that the artificial speed of sound does not enforce the divergence-free constraint of the leading-order velocity. Numerous numerical experiments using different structured meshes and different low Mach numbers have shown that the resulting scheme is insensitive to the mesh definition and is free of any pressure checkerboard issue. This also the case of the Roe-Turkel scheme, as demonstrated in [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF]. The Von Neumann analysis for stability using the Birken and Meister theorem [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF] shows that a much more restrictive stability condition must be applied in the low speed limit, than the suggested stability condition given by Rossow. A stringent stability condition also applies to the artificial speed of sound with ∆t ≤ O(M 2 ) as M → 0, which is asymptotically identical to the Roe-Turkel scheme, as demonstrated in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF]. However, it is worth mentioning that the rescaling of the Roe scheme using the artificial speed of sound is much more simple to implement than the complex matrix dissipation resulting from the Roe-Turkel scheme when formulated in a general framework [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF]. On the other hand, the artificial speed of sound approach also suffers from large errors at a stagnation point in the very incompressible limit, which could only be circumvented in the case of the NACA0012 airfoil using an extremely high density mesh at M ∞ = 10 -4 . This is a typical weakness in formulating schemes defined with a rescaling of the matrix-valued dissipation. This stems from the fact that no specific low speed preconditioning of the flux Jacobian matrix is introduced, as indicated by the mathematical framework used by Darmofal et al [START_REF] Darmofal | The Importance of Eigenvectors for Local Preconditioners of the Euler Equations[END_REF].

As shown by Klein [START_REF] Klein | Semi-Implicit Extension of Godunov-Type Scheme Based on Low Mach Number Asymptotics I: One-Dimensional Flow[END_REF], Guillard et al [START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF][START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review Rémi[END_REF], Dellacherie [START_REF] Dellacherie | Analysis of Godunov Type Schemes Applied to The Compressible Euler System at Low Mach Number[END_REF], the Euler equations in the incompressible limit are characterized by incompressible-acoustic interactions. This has major consequences on the stability and the accuracy of numerical schemes designed to handle incompressible flow using a density-based approach. One of main guideline in designing the low Mach number fix was also to cancel out spurious acoustics pressure waves arising from the jump of the normal velocity, as shown with the low Mach number flow extension of the Godunov scheme proposed in [START_REF] Guillard | Murrone On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II[END_REF]. Unsteady computations carried out in this paper for the NACA0012 airfoil clearly illustrate with the dilatation rate that the artificial speed of sound contains an unsteady acoustic component in the discrete solution of the steady-state problem, characterized by spurious acoustic pulses at short times, quickly dissipated at large times. This was correlated with typical limit cycles found in the convergence history for the steady-state problem, also observed in the case of the Roe-Turkel scheme [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF], and also visible in some convergence histories given by Rossow (see figures 2, 3 in [START_REF] Rossow | Efficient Computation of Compressible and Incompressible Flows[END_REF]). Such acoustic pulses are not visible as expected with the low Mach number fix at short times, although some spurious acoustic modes are still present in the discrete solution at large times, which should be dissipated in the time limit where the steady-state solution is recovered. However, both schemes have permanent residual acoustic disturbances in the discrete solution, as indicated by the asymptotic analysis for the acoustic time-scale. Acoustic components found in both schemes originate from pure numerical instabilities and do not have any physical meaning for the steady-state problem. This behaviour of the discrete solutions illustrates the complex acoustic-incompressible interaction occurring when a compressible flow solver is used in the very incompressible limit. A recent work especially points out that usual schemes developed for incompressible flows, including the Roe-Turkel scheme and the Rieper's fix, are not accurate enough for pure acoustic computations [START_REF] Bruel | A low Mach correction able to deal with low Mach acoustics[END_REF].

It might be tempting to combine both schemes, what can be easily formulated, since different coefficients of the matrix dissipation are affected by the low Mach number fix and the artificial speed of sound. For instance, the formulation of the artificial speed of sound could be corrected using the fix applied to the jump of the normal velocity. However, this strategy has to be considered with care and the formulation of such a blended scheme may not be as straightforward. From preliminary investigations carried out by the author, the resulting scheme may also combine the main drawbacks of both approaches.      γ -1 2

|V| 2 ρc (1 -γ) u ρc (1 -γ) v ρc (γ -1) ρc -u ρ 1 ρ 0 0 -v ρ 0 1 ρ 0 1 ρc γ -1 2 |V| 2 -c 2 (1 -γ) u ρc (1 -γ) v ρc (γ -1) ρc      .
In practice, the dissipation vector is formulated for the conservative variables within the change of variables

d = ∂W ∂ W0 D0 ∂ W0 ∂W (W R -W L ) ,
and is computed in our implementation using the algebraic algorithm described in [START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows AIAA paper[END_REF]. Note that the modified dissipation matrix [START_REF] Turkel | Assessment of Preconditioning Methods For Multidimensional[END_REF] for the Rossow's artificial speed of sound or [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF] for the Rieper's fix must be diagonalized for the transonic regime, as some entropy fix must be applied to eigenvalues approaching zero, with

D0 = R0 Λ L0 ,
where Λ = diag(µ + , µ -, µ 0 , µ 0 ) is the diagonal matrix of the modified eigenvalues and the left and right eigenvector matrix given explicitly in the following. The formulation of the entropy fix must be considered with care as indicated in [START_REF] Boniface | Rescaling of the Roe Scheme in Low Mach-Number Flow Regions[END_REF] for the case of the Roe-Turkel scheme. Identical treatments were considered for the modified eigenvalues ( 17)

Q + = 1 r (c -|q n |) + |q n | -µ + , Q -= 1 r (c -|q n |) + |q n | -µ -, with r = c c .
Note that with the identity Q + Q -= -q 2 n , it can be readily shown that the right and left eigenvector matrix are orthogonal. Thus, the rescaling of the Roe scheme according to Rossow preserves the properties of the original flux Jacobian matrix in symmetrizing variables in the subsonic range.

Low Mach number fix

The following right and left eigenvector matrix are derived for the dissipation matrix [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF] in d W0 variables:

R0 =                               Q + Q - 0 0
-n x q n -n x q n -n y 0

-n y q n -n y q n n x 0 and the sign of µ + -|q n | is also the sign of We know from the previous result that ∆ -M 2 0 (z + 1) 2 ≥ 0 and therefore |λ + | -µ -≥ 0 if -1 ≤ M 0 ≤ 0.

Figure 1 :

 1 Figure 1: Effect of vanishing δq n coefficient for the AUSM+(P) scheme on the pressure checkerboard for the non-lifting NACA0012 airfoil at M ∞ = 10 -3 .

  )-(61). For incompressible flows, a Poisson equation for the pressure is derived taking the divergence of the momentum equations and applying the divergence constraint. Using equations (60)-(61) for the leading-order velocity components, and with z = |u| c + |v| c from (

Figure 2 :

 2 Figure 2: Types of structured meshes used for the cylinder problem. Left: Cartesian, Center: Curvilinear irregular, Right: Curvilinear adapted.

Figure 3 :

 3 Figure 3: Amplitude of the normalized disturbance pressure δ p with the inflow Mach number M ∞ → 0 computed in the Cartesian and adapted meshes.

Figure 8 :

 8 Figure 8: Effect of the third-order MUSCL extrapolation for the spurious entropy in the adapted curvilinear mesh for M ∞ = 10 -3 .

Table 2 :

 2 ) are discussed for the non-lifting NACA0012 airfoil.Type of meshMesh definition (nodes) Mesh density (cells) Mesh extension (chord) O-mesh coarse 129 × Definition of structured meshes used for the NACA0012 airfoil.

Figure 9 :

 9 Figure 9: Types of structured meshes used for the NACA0012 airfoil. Left: Cartesian-like Vasseberg-Jameson O-meshes, Right: Curvilinear stretched C-mesh.

Figure 10 :

 10 Figure 10: Amplitude of the normalized disturbance pressure δ p with the inflow Mach number M ∞ → 0 computed in the coarse O-mesh and in the curvilinear stretched C-mesh.

Figure 11 :

 11 Figure 11: Contours of the disturbance pressure p2 with the inflow Mach number M ∞ → 0 in the curvilinear stretched C-mesh. Left: Rossow, Right: Rieper.

Figure 12 :

 12 Figure 12: C p distribution in the curvilinear stretched C-mesh with the inflow Mach number M ∞ → 0. Left: Rossow, Right: Rieper.

Figure 13 :Figure 14 :

 1314 Figure 13: Mesh convergence for the C p and entropy distributions using the Vasseberg-Jameson O-meshes with inflow Mach number M ∞ = 10 -4 . Left: Rossow, Right: Rieper.

2 Iteration 4 Figure 15 :

 2415 Figure 15: Comparison of the convergence history with M ∞ → 0 for the artificial speed of sound and the low Mach number fix.

Figure 16 :

 16 Figure 16: Comparison of the steady and unsteady results for the disturbance pressure p2 (Left) and the amplitude of the disturbance pressure δ p (Right), computed by the artificial speed of sound and the low Mach number fix in the fine O-mesh at the inflow Mach number M ∞ = 10 -4 .

Figure 17 :

 17 Figure 17: Farfield view of spurious acoustic pressure waves at M ∞ = 10 -4 computed by the artificial speed of sound (Left column) and the low Mach number fix (Right column) in the fine O-mesh 513 × 513 (∇.V contour levels in [-10 -7 , 10 -7 ]).

= c 2 ( 1 +

 21 (µ + -|q n |) (µ -+ |q n |) ≥0 [ √ ∆ + (z + 1)M n ] ≥0 = c 2 z [ √ ∆ + (z + 1)M n ] ≥0 +c 2 M n [∆ -(z + 1) 2 M 2 n ]. With ∆ -(z + 1) 2 M 2 n = (1 -M 2 n )(1z) 2 ≥ 0 if M n ≤ 1, then µ + -|q n | ≥ 0. 2. Similarly [µ + -ρ( A)][µ -+ ρ( A)] = c 2 (1 + M n )M n [ √ ∆ -(z + 1)]c 2 (1 + M n )(1z)M n ,with the last term -c 2 (1+ M n )(1z)M n ≤ 0.Then the sign of µ + -ρ( A) is also given by[µ + -ρ( A)] [µ -+ ρ(A)] M n )M n ≥0 [∆ -(z + 1) 2 ]c 2 (1 + M n )(1z)M n [ √ ∆ + (z + 1)] ≥0 , with ∆ -(z + 1) 2 = -4z(1 -M 2 n ) ≤ 0. Then we conclude that µ + -ρ( A) ≤ 0 ∀M n ≤ 1.3. With the previous result, we then haveµ -≤ µ + ≤ ρ(A) = max(|λ -|, |λ + |), with |λ -| = cq n = c(1 -M 0 ) and |λ + | = c + q n = c(1 + M 0 ), and therefore if -1 ≤ M 0 ≤ 0, then ρ( A) = |λ -| and if 0 ≤ M 0 ≤ 1, then ρ( A) = |λ + |.Then we just need to show thatµ -≤ |λ -| if 0 ≤ M 0 ≤ 1 and µ -≤ |λ + | if -1 ≤ M 0 ≤ 0.Let first assume 0 ≤ M 0 ≤ 1. We proceed as previously to derive the following expression(|λ -| -µ -)(|λ + | + µ + ) = c 2 [(1 -M 2 0 )(1z) + √ ∆ -M 0 (z + 1)].So the sign of |λ -| -µ -is also the sign of(|λ -| -µ -) (|λ + | + µ + ) ≥0 [ √ ∆ + M 0 (z + 1)] ≥0 = c 2 (1 -M 2 0 )(1z) ≥0 [ √ ∆ + M 0 (z + 1)] ≥0 +c 2 [∆ -M 2 0 (z + 1) 2 ]. With ∆ -M 2 0 (z + 1) 2 = (1 -M 2 0 )(z -1) 2 ≥ 0, then |λ -| -µ -≥ 0 if 0 ≤ M 0 ≤ 1. Now assume -1 ≤ M 0 ≤ 0 or else M 0 = -|M 0 |. Then (|λ + | -µ -)(|λ -| + µ + ) = c 2 [(1 -M 2 0 )(1z) + √ ∆ -|M 0 |(z + 1)].So with a similar argument as in the previous case(|λ + | -µ -) (|λ -| + µ + )

Table 1 :

 1 Definition of structured meshes used for the cylinder problem (mesh definition = circumferential × radial).

	Type of mesh	Mesh definition (nodes) Mesh density (cells) Mesh extension (diameter)
	Cartesian mesh coarse	75 × 25	1 776	5
	Cartesian mesh medium 1	150 × 50	7 301	5
	Cartesian mesh medium 2	300 × 100	29 601	5
	Cartesian mesh fine	600 × 200	119 201	5
	Curvilinear irregular mesh	140 × 40	5 421	5
	Curvilinear adapted mesh	449 × 60	26 432	20

Table 3 :

 3 Typical double precision outputs for the min & max non-dimensional pressure p p 0 = 1 + p 2 p 0 M 2 ∞ with M ∞ → 0 in conjunction with the limit cycle of both schemes.

Figure 4: Contours of the disturbance pressure p2 with the inflow Mach number M ∞ → 0 in the irregular curvilinear mesh. Theory in grey solid lines correspond to the analytical solution for the potential flow, using the Bernoulli theorem.

Overall solution accuracy

It is interesting to look at the behaviour of the C p distribution at the solid wall while decreasing the inflow Mach number to very low values. Let consider first the case of previous results obtained in the irregular mesh (Fig. 4) and the corresponding pressure coefficient at the wall plotted in Fig. 5. The C p coefficient is compared with the potential theory ("square" grey symbols). The low Mach number fix has a large pressure-velocity decoupling, which is amplified when M ∞ is becoming very low (≤ 10 -3 ). The solution could not be plotted at M ∞ = 10 -5 with the scale used in the figures. This is a consequence of errors at the stagnation point, as discussed next. However, for M ∞ = 10 -2 , the C p distribution has almost no pressure-velocity decoupling and compares very well with the potential theory (green curve in the Fig. 5-Left), what does not achieve the artificial speed of sound in this irregular mesh.

With the artificial speed of sound, no pressure-velocity decoupling can be observed at the wall and a consistent pressure coefficient is obtained for all Mach numbers M ∞ ≥ 10 -4 with identical plots, according to (72). The loss of symmetry w.r.t. the C p distribution computed from the potential theory is due to the non-symmetrical mesh used and a larger truncation error. At M ∞ = 10 -5 , it can be noted that the C p distribution is shifted as indicated with the dashed line in the right figure. This shift is visible at the upstream and downstream stagnation points for both schemes in the figures and is a consequence of larger errors occurring in the very incompressible limit in computing the stagnation pressure.

Errors at the stagnation points are amplified in the adapted mesh, characterized by a strong grid refinement at the upstream stagnation point, compared to the irregular grid. This is illustrated next in Fig. 6-Left comparing the C p distribution for both schemes at inflow Mach numbers M ∞ = 10 -3 and M ∞ = 10 -4 . On the other hand, we see that the low Mach number fix does not exhibit any pressure-velocity decoupling issue in this adapted mesh. The pressure coefficient matches with the potential solution at M ∞ = 10 -3 (solid lines). At M ∞ = 10 -4 (dashed lines), a large shift previously observed in the C p distribution is again found for both schemes. This can be readily explained with the stagnation pressure p * in the incompressible regime, that should be sum of the surrounding constant pressure p 0 and the dynamic pressure:

Any error in the computation of p * propagates all along the solid wall, being the streamline containing the upstream stagnation point. The wall pressure p w for the incompressible flow is given by p w = p * -ρ w

, where ρ w and V w denote the density and velocity vector at the solid wall, respectively. Therefore, a shift in the computation of p * yields a constant shift of the pressure all along the solid wall. This is exactly what can be observed for both schemes when M ∞ = 10 -4 . For the lower Mach number M ∞ = 10 -5 tested in these computations, the shift becomes so large that the C p could not be plotted with the scale used in these figures. Similar errors, with even much larger errors for the low Mach number fix, were also found the finest Cartesian mesh used in the previous mesh convergence study, for inflow Mach numbers < 10 -2 . Large errors for the stagnation pressure were also observed with the adapted mesh when moving the outer boundary up to

Adapted mesh

Figure 6: Cp distribution in the adapted curvilinear mesh for M ∞ = 10 -3 (solid lines) and M ∞ = 10 -4 (dashed lines). Left: Original first-order schemes. Right: Schemes with a third-order MUSCL extrapolation.

150 diameters. Identical computations in the adapted mesh were performed using a third-order MUSCL extrapolation for the reconstructed primitive variables, without limiter. The corresponding schemes are indicated with "Rieper-O3" and "Rossow-O3" in the next figures. Note that the MUSCL extrapolation does not modify the asymptotic behaviour of the truncation error described in section 6. However, a slight reduction of the truncation error was expected. The effect of the MUSCL extrapolation is indicated in Fig. 6-Right, compared to the same results obtained with the original firstorder schemes in the left figure. The pressure coefficient even better matches with the potential solution at M ∞ = 10 -3 (solid lines). With the artificial speed of sound, we get exactly C p = 1 at the stagnation points. When M ∞ = 10 -4 (dashed lines), the large shift in the C p distribution found for the first-order schemes is reduced, but not canceled. This effect is especially noticeable for the low Mach number fix.

The underlying mechanism for errors at stagnation points has to be found in the eigenvector structure of the flux Jacobian matrix, with a critical loss of orthogonality of the left and right eigenvectors as M ∞ → 0, when the Euler equations are formulated for the primitive or conservative variables. The mathematical framework developed by Darmofal et. al in [START_REF] Darmofal | The Importance of Eigenvectors for Local Preconditioners of the Euler Equations[END_REF] has demonstrated a transient error growth due to the non-normality of the Jacobian matrix, resulting in a lack of robustness near stagnation points. Significant improvements could be obtained in designing specific preconditioners for the Roe scheme or incorporating the symmetrizing entropy variables in the SUPG scheme developed in [START_REF] Wong | The Solution of the Euler Equations at Low Mach Number Using a Stabilized Finite Element Algorithm[END_REF] for low Mach number flows. The Rossow's artificial speed of sound and the Rieper's low Mach number fix are not defined from a low-speed preconditioning of the flux Jacobian matrix and only consistently rescale the matrix-valued dissipation. Numerical evidence exists that these schemes based on a rescaling of the matrix dissipation have accuracy issues at stagnation points, for very low Mach numbers < 10 -3 . This is also the case of the Roe-Turkel scheme. This will be also shown with the next test-case for the NACA0012 airfoil. On the other hand, modifying the flux Jacobian matrix requires to reformulate all boundary conditions based on characteristic variables, even for viscous flows [START_REF] Turkel | Assessment of Preconditioning Methods For Multidimensional[END_REF][START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows AIAA paper[END_REF]. This is one the main drawback of low-speed preconditioning techniques, especially for large aerodynamics codes, and for time-dependent problems, the accuracy may be lost as the preconditioning changes the time-dependent behaviour of the conservation laws [START_REF] Turkel | Vatsa Choice of Variables and Preconditioning for Time Dependent Problems AIAA paper[END_REF].

The effect of the MUSCL extrapolation in the adapted mesh is illustrated in the next Fig. 7 for the pressure disturbance field, with same plots as in Fig. In this appendix, it is demonstrated that for the eigenvalues (17) of the modified dissipation matrix according to the Rossow's artificial speed of sound [START_REF] Turkel | Assessment of Preconditioning Methods For Multidimensional[END_REF], the following inequalities hold in the subsonic range M n ≤ 1 (with c ≤ c or else r ≥ 1):

1. µ + ≥ |q n | (Then µ + is the spectral radius of the dissipation matrix); 2. µ + ≥ |q n | + c = ρ( A) the spectral radius of the flux Jacobian matrix; 1. Taking advantage of having both

, we can derive the following expression:

So the sign of µ + -|q n | is also the sign of

2. Following the same idea, we also derive

And the sign of µ + -ρ( A) is also the sign of

With ∆ -4r 2 M 2 n = (r 2 -1) 2 (1 -M n ) 2 ≥ 0, then all terms in the RHS of the previous expression are non-negative for M n ≤ 1 Therefore µ + -ρ( A ≥ 0 ∀M n ≤ 1, ∀r ≥ 1.

Appendix C.

In this appendix, it is demonstrated that for the eigenvalues (25) of the modified dissipation matrix with the Rieper's fix [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF], the following inequalities hold in the subsonic range M n ≤ 1 (with 0 ≤ z ≤ 1): 1. We proceed as in the previous appendix, using identity µ + µ -= zc 2 (1 -M 2 n ) with µ ± = c 2 z + 1 ± √ ∆ ≥ 0. We can derive the following expression: