Jean-Christophe Boniface

Rescaling of the Roe Scheme in Low Mach-Number Flow Regions

Keywords: Compressible flow, Euler equations, Low Mach number, Preconditioning, Roe scheme, Roe-Turkel scheme, Characteristic time-step

A rescaled matrix-valued dissipation is reformulated for the Roe scheme in low Mach-number flow regions from a well known family of local low-speed preconditioners popularized by Turkel.

The rescaling is obtained by suppressing the iterative preconditioning and by deriving explicitly the full set of eigenspaces of the Roe-Turkel matrix dissipation. This formulation preserves the time consistency and does not require to reformulate the boundary conditions based on the characteristic theory. The dissipation matrix achieves by construction the proper scaling in low-speed flow regions and returns the original Roe scheme at the sonic line. We find that all eigenvalues are nonnegative in the subsonic regime. By removing the iterative preconditioning, it becomes necessary to formulate a stringent stability condition to the explicit scheme in the low-speed flow regions based on the spectral radius of the rescaled matrix dissipation. This formulation also raises a two-timescale problem for the acoustic waves, which is circumvented for a steady-state iterative procedure by the development of a robust implicit characteristic matrix time-stepping scheme. The behaviour of the modified eigenvalues in the incompressible limit and at the sonic line also suggest applying the entropy correction carefully for complex non-linear flows.

Introduction

Iterative preconditioning based on Chorin artificial compressibility has become widely used for computing low-speed flow configurations with numerical schemes developed for compressible flows. This iterative preconditioning technique was designed to achieve an optimal conditioning of the iterative procedure and to guaranty the proper scaling of the artificial dissipation Email address: jean-christophe.boniface@onera.fr () Preprint submitted to Journal of Computational Physics April 1, 2016 when the Mach number approaches zero. The low-speed preconditioning approach has proved to be very efficient to overcome the accuracy issue of compressible flow solvers in the incompressible limit. Actually, the low-speed preconditioning should always be used since many industrial applications are characterized by mixed compressible and incompressible flows, over a wide range of Reynolds numbers.

However, the standard iterative preconditioning suffers from the complexity of its practical implementation. Since the local preconditioner modifies the characteristic relations, all boundary conditions based on characteristic variables or Riemann invariants must be reformulated accordingly. For large aerodynamics codes in which a large number of boundary conditions may be implemented, it is then necessary to reformulate most of the boundary conditions. Furthermore, the extension of iterative preconditioning to unsteady flows is not trivial and without a special treatment, time-accuracy may be lost.

The problem of the accuracy in the asymptotic limit of the incompressible flow for the discretization of the normalized Euler equations was first addressed by Guillard and Viozat [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF][START_REF]Viozat Implicit Upwind Schemes for Low Mach Number Compressible Flows INRIA Report[END_REF].

For their non-iterative preconditioning formulated with the Roe scheme as defined in [START_REF]Viozat Implicit Upwind Schemes for Low Mach Number Compressible Flows INRIA Report[END_REF], it was shown that the checkerboard modes for the leading and second-order pressure fields are cancelled out by the rescaled matrix dissipation and that the pressure field should be constant in space up to a fluctuation in space of order two. Furthermore, the authors have clearly pointed out a lack of dissipation of the standard conservative schemes in the incompressible limit.

A number of authors have further considered the discrete analysis based on the normalized equations for the asymptotic behaviour of the pressure and velocity fields. An "all-speed Roe scheme" has been developed by Li and Gu [START_REF] Li | An All-Speed Roe-Type Scheme and its Asymptotic Analysis of Low Mach Number Behaviour[END_REF] in order to recover at the discrete level the divergence constraint of the leading order velocity and the Poisson equation for the second-order pressure, which are not satisfied by the rescaling of the Roe scheme formulated in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF][START_REF]Viozat Implicit Upwind Schemes for Low Mach Number Compressible Flows INRIA Report[END_REF]. However, checkerboards modes are not automatically suppressed by their numerical flux and a low Mach number fix proposed by Rieper for the Roe's approximate Riemann solver [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF] seems attractive, as combining the advantages in the incompressible limit of both approaches investigated in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF] and [START_REF] Li | An All-Speed Roe-Type Scheme and its Asymptotic Analysis of Low Mach Number Behaviour[END_REF].

Over the last years, many modified Roe-type [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF], AUSM-type [START_REF] Mary | Deville An Algorithm for Unsteady Viscous Flows at all Speeds[END_REF] or Godunov-type schemes [START_REF] Guillard | On the Behaviour of Upwind Schemes in the Low Mach Number Limit: II. Godunov Type Schemes INRIA report[END_REF] have been formulated to apply conservative finite-difference schemes to low-speed flows. It has been found necessary to propose a unified theoretical framework to analyse their respective discrete properties and to understand why they fail to be accurate in the incompressible limit without specific corrections. It is worth mentioning the work of Li and Gu for the analysis of Roe-type schemes based on expansions in the freestream mach number of solutions of the normalized Euler equations [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF] and the contribution of Dellacherie for Godunov-type schemes [START_REF] Dellacherie | Analysis of Godunov Type Schemes Applied to The Compressible Euler System at Low Mach Number[END_REF] using the Hodge decomposition for solutions derived from the one-dimensional barotropic Euler equations.

On the other hand, few contributions have addressed the numerical stability of shock-capturing schemes adapted for low-speed flows. A Fourier Analysis is carried out by Dellacherie for the one-dimensional wave equations using the low Mach Godunov scheme and an explicit CFL condition is formulated for both the explicit and the implicit scheme [START_REF] Dellacherie | Analysis of Godunov Type Schemes Applied to The Compressible Euler System at Low Mach Number[END_REF]. For a non-iterative lowspeed preconditioning formulation of the compressible Euler equations, the issue of the stability for a matrix-valued dissipation formulated from the Lax-Friedrich scheme is pointed out for the first time by Birken and Meister [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF] on the basis of the asymptotic behaviour of the largest eigenvalue in the incompressible limit. Formulating a stability criteria is also an essential feature when designing numerical schemes for complex flows, especially when an "all-speed scheme" is being developed. Results obtained by Birken and Meister clearly show for the Euler equations that the standard CFL condition used for the computation of compressible flows is no longer valid in the incompressible limit and that a stringent stability condition for the time step with ∆t ≃ O(M 2 )

when M → 0 must be accounted for when a non-iterative preconditioning is considered. Nevertheless, the eigenspaces of the matrix-valued dissipation are not derived and a practical CFL condition for the local time step is not formulated explicitly for the steady-state problem.

Thus, a practical point of view has essentially motivated this contribution. For the compressible Euler equations, it was found interesting to investigate the effect of removing the premultiplication of the preconditioning matrix with the time-derivative of the independent flow variables. This formulation called "non-iterative preconditioning" and corresponding to a rescaling of the matrix-valued dissipation has been investigated particularly in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF][START_REF]Viozat Implicit Upwind Schemes for Low Mach Number Compressible Flows INRIA Report[END_REF] and [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF]. With this simplification, the explicit scheme recovers a basic structure with the centred scheme and stabilization terms. Then it becomes no longer necessary to reformulate the characteristic curves for the preconditioned Jacobian matrix and the time accuracy is preserved. Some recent attempts of improving the accuracy of conservative schemes in the low speed limit as mentioned above may have been also motivated by a drastic simplification of the implementation of the low-speed preconditioning. This is especially the case of the Rieper low Mach-number fix proposed in [START_REF] Rieper | A Low Mach Number Fix for Roe's Approximate Riemann Solver[END_REF],

and the Thornber et al. "Low Mach" LMRoe scheme [START_REF] Thornber | An Improved Reconstruction Method for Compressible Flows with Low Mach Number Features[END_REF] modifying jumps of the discrete velocities, which was further developed by Oßwald et al. with the "Low dissipation Low Mach"

L 2 Roe scheme [START_REF] Oßwald | a Low Dissipation Version of Roe's Approximate Riemann Solver for Low Mach Numbers[END_REF]. This may also explain the success of the AUSM-familly schemes and their modification for low Mach-number flows.

The stabilization terms must be then interpreted as a rescaling of the original dissipation in the low-speed flow regions and therefore the necessary Von Neumann criteria for the linear stability must be reconsidered completely. This can be achieved only by deriving the eigenvalues and the full set of the right and left eigenvectors of the matrix-valued dissipation. The diagonalization is performed step by step, by considering first as working variables the symmetrising entropic variables for the stream-aligned formulation of the Euler equations. The formulation for the conservative variables is obtained by successive similarity transformations. We have considered the Roe scheme as baseline formulation for the matrix dissipation. The necessary rescaling of the Roe scheme in the low speed flow regions is formulated from a family of preconditioners popularized by Turkel [START_REF] Turkel | Review of Preconditioning Methods For Fluid Dynamics NASA CR[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 family of preconditioners include the Choi-Merkle preconditioner [START_REF] Choi | The Application of Preconditioning To Viscous Flows[END_REF], with the particular case considered in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF]. The scheme is sometimes termed as the Roe-Turkel scheme [START_REF]Viozat Implicit Upwind Schemes for Low Mach Number Compressible Flows INRIA Report[END_REF][START_REF] Dellacherie | Analysis of Godunov Type Schemes Applied to The Compressible Euler System at Low Mach Number[END_REF]. The diagonalization of the rescaled matrix-valued dissipation must be achieved for the computation of complex flows, also because somehow an entropy fix may be used to prevent eigenvalues from approaching zero and to select the relevant physical solution satisfying the entropy condition across shocks.

By removing the iterative preconditioning, we loose the optimal conditioning but we save the essential feature of the accuracy in the incompressible limit. Without iterative preconditioning, we see that by forcing the proper scaling of the dissipation matrix, we cannot avoid the two-timescale issue, since the isentropic solution breaks up into a fast component that should disappear when the Mach number goes to zero, and a slow component solution of the incompressible equations [START_REF] Guillard -A | On the Behaviour of Upwind Schemes in the Low Mach Number Limit : II. Godunov Type Schemes INRIA Report[END_REF]. The system then becomes very stiff and numerically difficult to solve with the explicit scheme and a time step defined by the stability condition for the fastest acoustic wave. For steady-flow problems, a robust implicit scheme must be developed to improve both the stability and the damping properties of the iterative procedure. However, we must stress that our final objective is not to develop an optimal approach specialized to handle almost incompressible flows, but only a more accurate compressible flow solver in low speed flow regions, without reformulating the boundary conditions and keeping the consistency in time for unsteady flows.

In the following, we shall first consider the iterative and non-iterative preconditioning following Turkel's analysis. This allows a smooth transition to the formulation of the rescaled Roe scheme and the corresponding stability condition. For steady-state problems, a robust implicit scheme has been developed in order to circumvent the severe stability bound occurring for low Mach number flows. The implicit scheme combines the introduction of a characteristic time-step matrix and under-relaxation. Preliminary results for low-speed and transonic flows are then presented. It is experienced that higher accuracy is achieved with the rescaled matrix dissipation not only in the incompressible limit but also at the stagnation point. It is especially shown in the transonic case that the known spurious entropy produced by the Roe scheme at the stagnation point is strongly reduced by the rescaling of the Roe matrix dissipation.

Iterative and non-iterative preconditioning

Iterative preconditioning

When considering the iterative preconditioning, the time derivatives are premultiplied by a matrix which slows down the speed of the acoustic waves down the fluid speed. Hence, at low Mach numbers all the waves will travel at the same speed and therefore, the equations will be again well conditioned. The preconditioned Euler equations for the iterative preconditioning are formulated as follows

P -1 ∂W ∂t = - ( ∂F ∂x + ∂G ∂y ) = - ( A ∂W ∂x + B ∂W ∂y ) ,
where P is a preconditioning matrix. The expression of the Jacobian matrices A, B depends on the choice of the variables W. However, the preconditioning can be viewed as an artificial compressibility effect of the incompressible Euler equations and was first formulated for the Euler equations formulated for the primitive variables. The formulation of the preconditioner described below is the generalization of previous preconditioners due to Turkel [13], considering the in-

dependent variables W 0 = [p, u, v, S ]
T where p is the pressure u, v are the velocity components and S is the entropy. This leads to a non-conservative formulation for ease of presentation. The Euler equations are then formulated in quasi-linear form

∂W 0 ∂t = - [ A 0 ∂W 0 ∂x + B 0 ∂W 0 ∂y ] , (1)
with the corresponding Jacobian matrices in the (x, y) directions

A 0 =                            u ρa 2 0 0 1 ρ u 0 0 0 0 u 0 0 0 0 u                            , B 0 =                            v 0 ρa 2 0 0 v 0 0 1 ρ 0 v 0 0 0 0 v                            . ( 2
)
The family of preconditioners we shall consider contain a free parameter β as formulated by

Choi-Merkle and is defined as

P 0 -1 ∂W 0 ∂t = - [ A 0 ∂W 0 ∂x + B 0 ∂W 0 ∂y ] , (3)
with

P 0 -1 =                                                       a 2 β 2 0 0 δ αu ρβ 2 1 0 0 αv ρβ 2 0 1 0 0 0 0 1                                                       , P 0 =                                                       β 2 a 2 0 0 - β 2 a 2 δ - αu ρa 2 1 0 αu ρa 2 δ - αv ρa 2 0 1 αv ρa 2 δ 0 0 0 1                                                       . (4)
For δ = 0 we get the family of preconditioners including the Weiss-Smith formulation as suggested by Turkel [START_REF] Turkel | Preconditioned Methods For Solving the Incompressible And Low Speed Compressible Equations[END_REF], while for δ = 1 and α = 0 we recover the preconditioner introduced by Choi and Merkle [START_REF] Choi | The Application of Preconditioning To Viscous Flows[END_REF]. Setting α = 0 corresponds to the original artificial compressibility method formulated for the primitive variables while α = -1 corresponds to the artificial compressibility method formulated for the conservatives variables [START_REF] Turkel | Review of Preconditioning Methods For Fluid Dynamics NASA CR[END_REF]. In all the following the free parameter is formulated as a square of the local Mach number with

β 2 a 2 = min ( max ( k 1 M 2 , ϵ 2 ) , 1 ) (5)
where usually

ϵ 2 ≃ k 2 M 2 ∞ .
Computationally, k 2 is found to depend on the number of grid nodes near the stagnation points. For viscous flows, k2 also depend on α and may depend on the local cell Reynolds number [START_REF] Choi | The Application of Preconditioning To Viscous Flows[END_REF]. Note that in the next sections, the preconditioning parameter homogeneous to a speed |β| > 0 is introduced, with the following definition

|β| = a √ min ( max ( k 1 M 2 , ϵ 2 ) , 1 ) . (6)
To transition smoothly to transonic flow, we demand that β 2 a 2 → 1 with α = 0 as M → 1. In that case the preconditioning matrix returns the identity matrix.

We consider now the steady-state problem for the preconditioned Euler equations. We reformulate system (3) as follows

∂W 0 ∂t = - [ (P 0 A 0 ) ∂W 0 ∂x + (P 0 B 0 ) ∂W 0 ∂y ] . (7)
The eigenvalues of P 0 A 0 are independent of the parameter δ, with

             λ 0 = u, λ ± = zu ± √ (zu) 2 + β 2 a 2 (a 2 -u 2 ), z = 0.5 ( 1 -α + β 2 a 2 ) (8)
In a straightforward implementation of this procedure, the preconditioning algorithm can be converted from any set of variables to the appropriate conditioning matrix for the conservation variables. Let W 1 represent some set of independent variables and W 0 the solution of the preconditioned equations [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF], expressed as a function of the variables W 1 with W 0 (W 1 ). Then we also

have in variables W 1 ∂W 1 ∂t = - [ (P 1 A 1 ) ∂W 1 ∂x + (P 1 B 1 ) ∂W 1 ∂y ] . (9)
with the preconditioner and the flux Jacobian matrices obtained by similarity transformation

P 1 = ∂W 1 ∂W 0 P 0 ∂W 0 ∂W 1 , A 1 = ∂W 1 ∂W 0 A 0 ∂W 0 ∂W 1 , B 1 = ∂W 1 ∂W 0 B 0 ∂W 0 ∂W 1 . (10)
If W 1 = W represents the conservative variables, with A = ∂F ∂W , B = ∂G ∂W , we recover the preconditioned Euler equations for the conservative variables

∂W ∂t = -P ( ∂F ∂x + ∂G ∂y ) . (11)
Confronting expressions [START_REF] Li | Mechanism of Roe-type Schemes for All-Speed Flows and its Application[END_REF] and [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF] for the matrix coefficients, the preconditioned Jacobian matrices (P 0 A 0 ), (P 0 B 0 ) and (P 1 A 1 ), (P 1 B 1 ) are similar by any change of independent variables and therefore have the same eigenvalues. The same transformations apply to the rescaled Roe scheme derived in the next sections.

Extension to time-dependent flows

Turkel has proposed a straightforward extension of low speed preconditioning for time dependent problems, based on the dual-time stepping method [START_REF] Turkel | Vatsa Choice of Variables and Preconditioning for Time Dependent Problems AIAA paper[END_REF]. The unsteady equations are then replaced by

∂W ∂τ + [ ∂W ∂t + ∂F ∂x + ∂G ∂y ] = 0, ( 12
)
where τ is the dual time. The time derivative is discretized by some first or second-order differencing formula. A typical implicit difference scheme for the preconditioned system in pseudotime with generic artificial dissipation terms can be expressed with the following conservative formulation

∆W = -∆t * P         1 ∆t         c 0 W n + ∑ k∈{-1,1} c k W n+k         + ( ∂F ∂x + ∂G ∂y ) - ∂ ∂x ( P -1 |P A| ∂W ∂x ) - ∂ ∂y ( P -1 |PB| ∂W ∂y )] , (13)
where ∆t * is a pseudo time-step. The numerical treatment for the artificial dissipation is then unchanged by the introduction of the backward differencing formula considered for the approximation of the time derivative.

Non-iterative preconditioning

A second approach to handle both steady and time-dependent flows is described in [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF]. A time-consistent method is considered with the preconditioning used as a way of controlling the effects of the artificial dissipation in case of low speed flow. The proposed method is formulated for the conservative variables by removing the pre-multiplication of the residual with the preconditioning matrix

∆W = -∆t [( ∂F ∂x + ∂G ∂y ) - ∂ ∂x ( P -1 |P A| ∂W ∂x ) - ∂ ∂y ( P -1 |PB| ∂W ∂y )] . (14)
We can see that system ( 14) is a time-consistent discretization of the non-preconditioned Euler equations based on a backward Euler approximation of the time derivative. In the framework of a steady state problem, the time step ∆t is replaced by a local time step. This non-iterative preconditioning algorithm is interpreted by Turkel as a preconditioning of the stabilization term only. Following the same idea, we may also consider the following dual time-stepping scheme for unsteady problems

∆W = -∆t *         1 ∆t         c 0 W n + ∑ k∈{-1,1} c k W n+k         + ( ∂F ∂x + ∂G ∂y ) - ∂ ∂x ( P -1 |P A| ∂W ∂x ) - ∂ ∂y ( P -1 |PB| ∂W ∂y )] . (15)
We clearly see that both schemes can be interpreted as a modification of the matrix-valued dissipation scheme. Thus, we will prefer to consider the formulations ( 14) or [START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows NASA Contractor Report 201605[END_REF] as a rescaling of the dissipation for low-speed flows, since the two approaches are quite different. The noniterative preconditioning achieves by construction the proper scaling of the matrix dissipation in the incompressible limit, but doesn't improve the condition number of the iterative procedure and the stability bound is drastically modified in the low-speed flow regions. Additionally formulations [START_REF] Turkel | Preconditioned Methods For Solving the Incompressible And Low Speed Compressible Equations[END_REF] or [START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows NASA Contractor Report 201605[END_REF] alleviate its implementation in a CFD code since once removed the iterative preconditioning, there is no need to reformulate the boundary conditions using the characteristic theory or the Riemann invariants. As a final comment, for the steady-state problem or the time-dependent problem, the time or pseudo-time derivative vanishes and since the preconditioning matrix is not singular, both iterative and non-iterative preconditioning methods fall into the same steady-state or unsteady equations. At the asymptotic convergence of the iterative process it is worth noticing that in both cases the governing matrix dissipation is given by coefficients (P -1 |P A|, P -1 |PB|) in the (x, y) space directions.

However, for the linear problem, the stability of the iterative preconditioning is characterized by |P A| and |PB|, while removing the pre-multiplication of the residual with P, the stability must be reconsidered from matrices P -1 |P A| and P -1 |PB|. It becomes then essential to reformulate the Von Neumann condition by deriving algebraically the corresponding eigenspaces and eigenvalues. This can be achieved step by step as indicated in the two next sections, starting from the stream-aligned formulation of the Euler equations for the symmetrizing entropic variables.

3. Formulation of the rescaled Roe scheme 3.1. One-dimensional framework

Iterative preconditioning

In order to investigate the limit to incompressible equations as the Mach number goes to zero, we follow Turkel [START_REF] Turkel | Preconditioned Techniques in Computational Fluid Dynamics[END_REF] and explicitly describe the scheme for the one-dimensional equations, including a second-order artificial viscosity term. For the equations termed in entropic variables, we shall consider first a streamwise two-dimensional coordinate system with the x-axis aligned to the flow. Then setting ∂v ∂y = 0, ∂p ∂y = 0 in the basic system (1) and taking as independent variables for the pressure dΦ = d p ρa and dS = d p-a 2 dρ ρa proportional to the entropy for the symmetrizing variables, system (1) can be recast as follows

                                                     ∂Φ ∂t + u ∂Φ ∂x + a ∂u ∂x = 0, ∂u ∂t + a ∂Φ ∂x + u ∂u ∂x = 0, ∂v ∂t + u ∂v ∂x = 0, ∂S ∂t + u ∂S ∂x = 0. ( 16
)
The equations for v and S decouple and so we shall consider the two acoustic equations. Let introduce some artificial dissipation terms with constant coefficients. We then get the following

scheme                      ∂Φ ∂t + u ∂Φ ∂x + a ∂u ∂x = Q 1,1 ∂ 2 Φ ∂x 2 + Q 1,2 ∂ 2 u ∂x 2 , ∂u ∂t + a ∂Φ ∂x + u ∂u ∂x = Q 2,1 ∂ 2 Φ 2 ∂x 2 + Q 2,2 ∂ 2 u ∂x 2 . (17)
Let d W0 = [dΦ, du] T . In matrix form, the previous system reads

∂ W0 ∂t + Ã0 ∂ W0 ∂x = Q ∂ 2 W0 ∂x 2 , (18)
with the symmetric Jacobian matrix

Ã0 =                  u a a u                  . (19)
The matrix Ã0 has the two eigenvalues λ ± = u ± a corresponding to the acoustic waves. The limit to incompressible equations for the scheme [START_REF] Turkel | Preconditioned Techniques in Computational Fluid Dynamics[END_REF] has been investigated in [START_REF] Turkel | Preconditioned Techniques in Computational Fluid Dynamics[END_REF][START_REF] Turkel | Leer Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes[END_REF]. From now on, all quantities ρ, u, v, p, S are nondimensionalized by some reference quantities ρ r , u r , p r = ρ r u 2 r , S r = r r , where r r is the reference perfect gas constant. Let introduce ϵ = u r a r = M. In the incompressible limit, the fluid is assumed isentropic. Then

p = ρ γ γϵ 2 ≃ O( 1 M 2 ), a = ρ γ-1 2 ϵ ≃ O( 1 M ) as M → 0. ( 20
)
Within a formal asymptotic analysis, the mesh size being kept constant, it was found that in the incompressible limit, the proper scaling of the matrix-valued dissipation requires

Q 1,1 = O( 1 M 2 ), Q 1,2 = O( 1 M ), Q 2,1 = O( 1 M ), Q 2,2 = O(1). (21)
For the system (17), we assume that the artificial dissipation is some function of the matrix | Ã0 |.

In the case of subsonic flow |u| ≤ a, we have the following diagonalization

| Ã0 | =                  a a a -a                                   a + u 0 0 a -u                                   1 2a 1 2a 1 2a -1 2a                  =                  a u u a                  . ( 22
)
Thus, as M → 0:

| Ã0 | ≃                    O( 1 M ) O(1) O(1) O( 1 M )                    . ( 23
)
This implies that the momentum equation has too much artificial viscosity, while the the divergence equation has too little viscosity which will cause oscillations. For the system in variables d W0 , the preconditioner becomes

P0 -1 =                      a 2 β 2 0 αua β 2 1                      , P0 =                    β 2 a 2 0 - αu a 1                    . ( 24
)
The preconditioned Jacobian matrix

P0 Ã0 =                       β 2 a 2 u β 2 a a ( 1 -α u 2 a 2 ) (1 -α)u                       , (25)
has eigenvalues λ ± given by ( 8). We have the following diagonalization

| P0 Ã0 | = R|Λ| R-1 , (26)
where

|Λ| = diag(|λ + |, |λ | ), R and
R-1 are respectively the matrices of the right and left eigenvectors

R =                     β 2 a β 2 a λ + - β 2 a 2 u λ -- β 2 a 2 u                     , R-1 = 1 (λ + -λ -)                     a β 2 [λ + -(1 -α)u] 1 - a β 2 [λ --(1 -α)u] -1                     . ( 27
)

Rescaled Roe scheme

By removing the iterative preconditioning, we introduce the rescaled Roe scheme

∂ W0 ∂t + Ã0 ∂ W0 ∂x = Q0 ∂ 2 W0 ∂x 2 , (28)
with the matrix-valued dissipation

Q0 = P0 -1 | P0 Ã0 | = P0 -1 R|Λ| R-1 .
In the following, we will assume the flow subsonic. Let χ = (zu

) 2 + β 2 a 2 (a 2 -u 2 ). We always have χ > 0 and χ ≥ (zu) 2 ⇒ ( zu - √ χ ) ( zu + √ χ ) ≤ 0.
Hence, eigenvalues λ + and λ -have opposite signs for a subsonic flow with

if u ≤ 0 → λ -≤ 0 ⇒ λ + ≥ 0, if u ≥ 0 → λ + ≥ 0 ⇒ λ -≤ 0.
Thus |λ + | = λ + and |λ -| = -λ -and we get an explicit expression of the dissipation matrix

P0 -1 | P0 Ã0 | = 1 √ χ                     [ zM 2 + ( 1 -M 2 )] a 2 z u a [ z + α ( 1 -M 2 )] u a [ zM 2 + β 2 a 2 ( 1 -M 2 ) ] a 2                     , (29)
where M = |u| a . Note that for the two acoustic equations

P0 -1 | P0 Ã0 | = 1 √ χ [ zu Ã0 + β 2 a 2 (a 2 -u 2 ) P0 -1
] as a consequence of the Cayley-Hamilton theorem. Therefore the dissipation matrix can be interpreted as a weighted average between the Jacobian matrix and the preconditioner. At a sonic point the second coefficient vanishes and it can be easily verified that for both possible values u = ±a the rescaled dissipation matrix returns | Ã0 | for α = 0. Let be Q i j the coefficients of matrix [START_REF] Turkel | Robust Low Speed Preconditioning for Viscous High Lift Flows AIAA paper[END_REF]. In the limit to the incompressible equations, assuming M ∞ << 1 and

β 2 ≃ u 2 = O(1),we
successively have

β 2 a 2 ≃ u 2 a 2 = O(M 2 ), χ = z 2 u 2 + β 2 ( 1 -M 2 ) ≃ u 2 = O(1),
and

Q 11 = 1 √ χ [ zM 2 + ( 1 -M 2 )] a 2 ≃ a 2 |u| = O( 1 M 2 ), Q 12 = 1 √ χ z u a ≃ a = O( 1 M ) Q 21 = 1 √ χ [ z + α ( 1 -M 2 )] u a ≃ a = O( 1 M ), , Q 22 = 1 √ χ [ zM 2 + β 2 a 2 ( 1 -M 2 ) ] a 2 ≃ |u| = O(1).
Then, as M → 0

P0 -1 | P0 Ã0 | ≃                    O( 1 M 2 ) O( 1 M ) O( 1 M ) O (1)
                   , (30)
and the scheme gives the correct asymptotic order for the artificial dissipation terms. It is worth noticing that the the rescaling (30) corresponds to a net increase of numerical dissipation in the pressure equation by a factor 1 M compared to the unpreconditioned scheme [START_REF] Mulder | Leer Experiments with Implicit Upwind Methods for the Euler Equations[END_REF] whereas in the equation for the stream-aligned velocity, the coefficients of the stabilization terms are consistently rescaled by the preconditioner.

We also see that matrix P0

-1 | P0 Ã0 | has two real eigenvalues given by

µ ± = a 2 √ χ [ 2z + α ( 1 -M 2 ) ± √ ∆ ] , (31)
with

∆ = ( 1 - β 2 a 2 ) 2 ( 1 -M 2 ) 2 + 4z [ z + α ( 1 -M 2 )] M 2 , χ = z 2 M 2 + β 2 a 2 (1 -M 2 ).
We find that both eigenvalues µ ± are positive as long as the flow is subsonic with

µ + µ -= Q 11 Q 22 -Q 12 Q 21 = det ( P0 -1 | P0 Ã0 | ) ,
and det

( P0 -1 | P0 Ã0 | ) = [ ( z -β 2 a 2 ) 2 -α β 2 a 2 ] u 2 + β 2 [ z 2 -β 2 a 2 ] u 2 + β 2 ( a 2 -u 2 ) = a 2 -u 2 = det ( | Ã0 | ) > 0.
Note that rescaled matrix dissipation is symmetric for α = 0 and therefore is positive definite with the two positive eigenvalues for in case of subsonic flow. This property is lost by any change of variables that does not preserve the symmetry of P0 -1 | P0 Ã0 |. In the incompressible limit, the eigenvalues behave as

µ ± ≃ a 2 2|β| [ 1 + β 2 a 2 ± |1 - β 2 a 2 | ] .
Since

β 2 a 2 ≤ 1 µ + ≃ a 2 |β| = O( 1 M 2 ), µ -≃ |β| = O(1).
Thus, in the asymptotic limit acoustics waves associated to µ + will travel at infinite speed while acoustic waves associated to µ -are slowed down to the flow velocity. This is a general behaviour of the eigenvalues µ ± in the low speed limit, because (30) holds necessarily with the coefficients Q i j for a proper scaling of the artificial dissipation, whatever the definition of the preconditioner P0 . This raises the well known singularity of the two time-scale problem. The optimal conditioning of the iterative preconditioning is lost and is even worse that the original system. The behaviour of eigenvalues µ ± in the low Mach number limit will become critical for the stability of the rescaled Roe scheme and a CFL condition must be formulated accordingly. For the other extreme case, at a sonic line, for α = 0, we see that we recover the values of the unpreconditioned system

µ + = 2a, µ -= 0.
When the preconditioning is not activated, setting α = 0 and β 2 a 2 = 1 in expressions (31) gives µ ± = a(1 ± M) = a ± |u|, the eigenvalues of | Ã0 | recast for the velocity module.

We can also easily derive the following diagonalization of the rescaled matrix dissipation

P0 -1 | P0 Ã0 | = R0 Λ R0 -1 , (32)
where Λ = diag(µ + , µ -), R0 is the matrix of the right eigenvectors for the acoustic eigenvalues

µ ± and R0 -1 is obtained by direct inversion of R0 R0 =                  [ µ + -Q 22 ] Q 12 Q 21 [ µ --Q 11 ]                  , R0 -1 = 1 det ( R0 )                  [ µ --Q 11 ] -Q 12 -Q 21 [ µ + -Q 22 ]                 
, with the determinant of matrix R0 given by det

( R0 ) = [ µ --Q 11 ] [ µ + -Q 22 ] -Q 12 Q 21 = µ + [ µ --Q 11 ] + µ - [ µ + -Q 22 ] = [ µ + -Q 22 ] [ µ --µ + ] = [ Q 11 -µ - ] [ µ --µ + ] . (33)
At a stagnation point, Q 12 = 0, Q 21 = 0, and

µ + = Q 11 = a 2 |β| , µ -= Q 22 = |β|, µ + -µ -= |β| ( a 2 β 2 -1
) .

Thus with µ + -µ -0, we always have µ --Q 11 0, µ + -Q 22 0, and the decomposition (32) can never become singular in the range |u| < a. Additionally, with α = 0, since by construction

µ + + µ -= Q 11 + Q 22 , we have (µ + -Q 22 )Q 12 + (µ --Q 11 )Q 21 = 0 with Q 12 = Q 21 , (34)
and we can see that the right and left eigenvectors are orthogonal. In this special case det

( R0 ) = - [ (µ + -Q 22 ) 2 + Q 2 12
] .

Then defining

M0 = 1 √ (µ + -Q 22 ) 2 + Q 2 12 R0 (35)
the symmetric matrix P0 -1 | P0 Ã0 | is diagonalized by the following unitary congruence

P0 -1 | P0 Ã0 | = M0 Λ M0 T .
(36)

Multi-dimensional extension

Iterative preconditioning

We derive the multidimensional extension of the previous analysis, assuming the flow subsonic. Let be d W0 = [dΦ, du, dv, dS ] T the set of entropy variables for the two-dimensional Euler equations. In matrix form, the two-dimensional scheme reads

∂ W0 ∂t + Ã0 ∂ W0 ∂x + B0 ∂ W0 ∂y = Q x ∂ 2 W0 ∂x 2 + Q y ∂ 2 W0 ∂y 2 , ( 37
)
with the expressions of the Jacobian matrix in the (x, y) directions andQ x and Q y are some matrix-valued artificial dissipation coefficients in the (x, y) directions.

Ã0 =                                                 u a 0 0 a u 0 0 0 0 u 0 0 0 0 u                                                 , B0 =                                                 v 0 a 0 0 v 0 0 a 0 v 0 0 0 0 v                                                 ,
For a general direction n = [n x , n y ] T , we will consider a matrix-valued dissipation in the general

form Q(n) = n x Q x + n y Q y .
As the Mach number goes to zero, a proper scaling of system (37) is achieved when the coefficients of the matrix-valued dissipation satisfy [START_REF] Turkel | Preconditioned Techniques in Computational Fluid Dynamics[END_REF][START_REF] Turkel | Leer Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes[END_REF]]

Q(n) ≃                                                 O( 1 M 2 ) O( 1 M ) O( 1 M ) 0 O( 1 M ) O(1) O(1) 0
O( 1 M ) O(1) O(1) 0 0 0 0 O(1)                                                 as M -→ 0. ( 38
)
A standard dissipation for the scheme (37) will be some function of the Jacobian matrix Ã0 (n) = n x Ã0 + n y B0 . Note that we can always formulate Ã0 (n) for any unitary vector n 1 = n |n| since Ã0 (n) = Ã0 (n 1 )|n| independently of the choice of variables. This is obviously also the case for the dissipation matrices | Ã0 (n)| and P0

-1 | P0 Ã0 (n)| that we will consider in the following. In the case of the Roe scheme, the matrix-valued dissipation is characterized by | Ã0 (n)|. The matrix Ã0 (n) has eigenvalues λ ± = q ± a|n| and double λ 0 = q, with q = un x + vn y . Assuming subsonic flow, we have |q| ≤ a|n| and |λ

+ | = a|n| + q, |λ -| = a|n| -q. Hence |Λ| = diag(a|n| + q, a|n| -q, |q|, |q|).
Thus, we get the explicit expression

| Ã0 (n)| =                                                     a|n| n x |n| q n y |n| q 0 n x |n| q a |n| ( n 2 x + n 2 y M n ) n x n y |n| a (1 -M n ) 0 n y |n| q n x n y |n| a (1 -M n ) a |n| ( n 2 y + n 2 x M n ) 0 0 0 0 |q|                                                     , where M n = |q| a|n| is the directional Mach number. Let M = √ u 2 + v 2 a be the local Mach number.
It is actually important in the following to differentiate the local Mach number defining the preconditioning parameter β 2 a 2 and the directional (normal) Mach number arising in the formulation of the Jacobian matrices. We have M n ≤ M. Then as M -→ 0, M n -→ 0 and

| Ã0 (n)| ≃                                                 O( 1 M ) O(1) O(1) 0
O(1) O( 1 M ) O( 1 M ) 0 O(1) O( 1 M ) O( 1 M ) 0 0 0 0 O(1)                                                 as M -→ 0.
This is to be compared to result (38). As for the one-dimensional case, we can see that without preconditioning in the incompressible limit, the pressure equation has to little numerical viscosity while the momentum equation has too much numerical viscosity. We can see that the eigenvalues of the preconditioned matrix P0 Ã0 are independent of the parameter δ in the preconditioner. This is also true for the rescaled matrix dissipation P0 -1 | P0 Ã0 |. Thus, in the following, we will set δ = 0 in the family of preconditioners, since we are basically interested in finding out a stability criteria for the rescaled Roe scheme.

For the set of variables d W0 , the preconditioning matrix in the two-dimensional case is then given by

P0 -1 =                                                       a 2 β 2 0 0 0 αua β 2 1 0 0 αva β 2 0 1 0 0 0 0 1                                                       , P0 =                                                    β 2 a 2 0 0 - αu a 1 0 - αv a 0 1 0 0 0                                                    .
In direction n, the Jacobian matrix of the preconditioned system reads

P0 Ã0 (n) =                                                     β 2 a 2 q β 2 a n x β 2 a n y 0 a ( n x - αuq a 2 ) q -αun x -αun y 0 a ( n y - αvq a 2 ) -αvn x q -αvn y 0 0 0 0 q                                                     .
We have the following diagonalization [START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows NASA Contractor Report 201605[END_REF]

| P0 Ã0 (n)| = R|Λ| R-1 , with eigenvalues              λ 0 = q, λ ± = zq ± √ (zq) 2 + β 2 a 2 ( a 2 |n| 2 -q 2 ) . (39)
The matrix of the right eigenvectors given by

R =                                                         β 2 a β 2 a 0 0 β 2 n x -αuλ + λ + -q β 2 n x -αuλ - λ --q -n y 0 β 2 n y -αvλ + λ + -q β 2 n y -αvλ - λ --q n x 0 0 0 0 1                                                        
, and the matrix of the left eigenvectors can be expressed as follows

R-1 =                                                        a β 2 (λ + -(1 -α)q) λ + -λ - n x λ + -λ - n y λ + -λ - 0 - a β 2 (λ --(1 -α)q) λ + -λ - - n x λ + -λ - - n y λ + -λ - 0 αa (un y -vn x ) (λ + -q)(λ --q) (β 2 n y -αvq) (λ + -q)(λ --q) - (β 2 n x -αuq) (λ + -q)(λ --q) 0 0 0 0 1                                                        .

Rescaled Roe scheme

We now consider the non-iterative preconditioning scheme

∂ W0 ∂t + Ã0 ∂ W0 ∂x + B0 ∂ W0 ∂y = P0 -1 | P0 Ã0 | ∂ 2 W0 ∂x 2 + P0 -1 | P0 B0 | ∂ 2 W0 ∂y 2 .
In case of subsonic flow, we know the sign of the two acoustic eigenvalues, with λ + ≥ 0 and λ -≤ 0. Therefore

|Λ| = diag(λ + , -λ -, |q|, |q|).
We proceed as for the one-dimensional case to find an explicit expression of the matrix valued dissipation

P0 -1 | P0 Ã0 (n)| = P0 -1 | R|Λ| R-1 .
Let first introduce the following coefficients

Q 11 = 1 √ χ |n| 2 [ zM 2 n + ( 1 -M 2 n )] a 2 , Q 12 = 1 √ χ z q a, Q 21 = 1 √ χ |n| 2 [ z + α ( 1 -M 2 n )] q a Q 22 = 1 √ χ |n| 2 [ zM 2 n + β 2 a 2 ( 1 -M 2 n ) ] a 2 ,
with z = 0.5

( 1 -α + β 2 a 2
) and

χ = (zq) 2 + β 2 a 2 ( a 2 |n| 2 -q 2 ) = a 2 |n| 2 [ z 2 M 2 n + β 2 a 2 (1 -M 2 n )
] .

We also define

ξ x = αuq -β 2 n x (λ + -q)(λ --q) , ξ y = αvq -β 2 n y (λ + -q)(λ --q) ,
and

ζ x = n y (un y -vn x ) (λ + -q)(λ --q) , ζ y = - n x (un y -vn x ) (λ + -q)(λ --q)
, with λ ± given by (39). Using the identity (λ +q)(λ -q) = αq 2 -β 2 |n| 2 , we have

ξ x n x + ξ y n y = 1. (40)
We also see that

ζ x n x + ζ y n y = 0. ( 41
)
The rescaled dissipation matrix can be expressed as follows

P0 -1 | P0 Ã0 (n)| =                                                 Q 11 n x Q 12 n y Q 12 0 ξ x Q 21 -αζ x a|q| n x ξ x (Q 22 -|q|) + |q| n y ξ x (Q 22 -|q|) 0 ξ y Q 21 -αζ y a|q| n x ξ y (Q 22 -|q|) n y ξ y (Q 22 -|q|) + |q| 0 0 0 0 |q|                                                 . (42)
In the multidimensional case, we find the following decomposition

P0 -1 | P0 Ã0 (n)| = 1 √ χ [ zq Ã0 (n) + β 2 a 2 (a 2 |n| 2 -q 2 ) P0 -1 ] + (|q| -Q 22 ) L0 (n) (43)
with

L0 (n) =                                                 0 0 0 0 -αζ x a n y ξ y 0 0 -αζ y a 0 n x ξ x 0 0 0 0 1                                                 .
The matrix (|q| -Q 22 ) L0 (n) has essentially |q| -Q 22 as non-zero eigenvalue with |q| ≃ Q 22 in the whole subsonic range except in the low speed limit. It doesn't affect the pressure equation. On the other hand, the resulting matrix of the two first contributions has the two acoustic eigenvalues µ ± > 0 given in the following and double eigenvalue Q 22 > 0. So, this resulting matrix is close to the rescaled matrix dissipation (42). In particular, this "acoustic" contribution has the same behaviour as the rescaled matrix dissipation in the low speed limit. Additionally, when α = 0 we see that P0 -1 and L0 (n) have only diagonal entries and therefore the rescaled dissipation matrix is symmetric. At a sonic point, i.e. |q| = a|n|, if α = 0 we also get Q 22 = a|n|. Therefore the two last contributions vanish and again for both possible values q = ±a|n| we check that

P0 -1 | P0 Ã0 (n)| = | Ã0 (n)|.
In the incompressible limit, as M -→ 0, we have successively

χ ≃ |q| 2 ,
and for the coefficients

Q i j Q 11 ≃ a 2 |q| |n| 2 = O( 1 M 2 ), Q 12 ≃ a = O( 1 M ), Q 21 ≃ a|n| 2 = O( 1 M ), Q 22 ≃ |q| = O(1).
We also have

ξ x ≃ O(1), ξ y ≃ O(1), ζ x ≃ O(1), ζ y ≃ O(1). Finally P0 -1 | P0 Ã0 (n)| ≃                                                 O( 1 M 2 ) O( 1 M ) O( 1 M ) 0 O( 1 M ) O(1) O(1) 0
O( 1 M ) O(1) O(1) 0 0 0 0 O(1)                                                 as M -→ 0,
which gives by construction the proper scaling in the incompressible limit. Again, the rescaling can be interpreted as a net increase of numerical dissipation in the pressure equation by a factor 1 M and a consistent rescaling for the flow velocity.

From expression (42), we can easily find the eigenvalues µ of the matrix Q0 . Using the identities (40) and (41), we find for the linear waves

µ 0 = |q| (double),
and for the acoustic waves

µ ± = 1 2 [ Q 11 + Q 22 ± √ (Q 11 -Q 22 ) 2 + 4 Q 12 Q 21 ] .
Upon substitution with the expressions of the Q i j

µ ± = a|n| 2 √ χ [ 2z + α(1 -M 2 n ) ± √ ∆ ] with ∆ = ( 1 - β 2 a 2 ) 2 ( 1 -M 2 n ) 2 + 4z [ z + α(1 -M 2 n ) ] M 2 n , χ = z 2 M 2 n + β 2 a 2 (1 -M 2 n ).
We find as expected from the one-dimensional analysis that as long as the flow remains subsonic, the eigenvalues µ ± are positive with the following identity

µ + µ -= a 2 |n| 2 -q 2 . ( 44
)
Thus, with α = 0 and µ 0 ≥ 0, the rescaled dissipation matrix formulated for the set of independent variables d W0 is positive semidefinite in the subsonic range. For the two extreme cases

as M n -→ 0 ⇒ µ + ≃ a 2 |β| |n| = O( 1 M 2 ), µ -≃ |β||n| = O(1), µ 0 ≃ O(1)
as

M n -→ 1 ⇒ µ + ≃ 2a|n|, µ -≃ 0, µ 0 ≃ a|n|.
As anticipated in the one-dimensional case, we can see that severe stability constraints will occur in the incompressible limit, as it is discussed later.

In the multi-dimensional case, we derive the following diagonalization of the rescaled matrix dissipation

P0 -1 | P0 Ã0 (n)| = R0 Λ R0 -1 ,
with Λ = diag(µ + , µ -, |q|, |q|), R0 is the matrix of the right eigenvectors

R0 =                                                 (µ + -Q 22 ) Q 12 0 0 [ ξ x Q 21 -αζ x a|q|R 1 ] [ ξ x (µ --Q 11 ) -αζ x a|q|R 2 ] -n y 0 [ ξ y Q 21 -αζ y a|q|R 1 ] [ ξ y (µ --Q 11 ) -αζ y a|q|R 2 ] n x 0 0 0 0 1                                                 , (45)
where

R 1 = µ + -Q 22 µ + -|q| , R 2 = Q 12 µ --|q| ,
and R0 -1 is the matrix of the left eigenvectors

R0 -1 =                                                   (µ --Q 11 ) D - Q 12 D n x - Q 12 D n y 0 - Q 21 D (µ + -Q 22 ) D n x (µ + -Q 22 ) D n y 0 -αζa|q|L 1 -ξ y -αζ y a|q|L 2 ξ x + αζ x a|q|L 2 0 0 0 0 1                                                   , (46)
with

ζ = ζ x n y = - ζ y n x = (un y -vn x ) (λ + -q) (λ --q) ,
and

L 1 = Q 22 -|q| (µ + -|q|) (µ --|q|) , L 2 = Q 12 (µ + -|q|) (µ --|q|)
,

where D = det ( R0 )
given by (33).

The above scaling of the eigenvectors corresponds to the most straightforward extension of the one-dimensional formulation. It is worth noticing that the apparent singularity µ --|q| = 0 which may occur if α 0 for a specific directional Mach number M n < 1, root of the non-linear equation

M 2 n + µ + a|n| M n -1 = 0,
doesn't create any numerical difficulty. We actually find by a deep inspection of the eigenspaces corresponding the second right eigenvector and the third left eigenvector that if µ -= |q|, then necessarily the velocity vector is aligned with the actual direction n and therefore

ζ x = ζ y = ζ = 0.
Thus, for this specific state, the second right eigenvector and the third left eigenvector recover non-singular expressions. It will be also shown in the stability analysis that µ + > |q| for all possible values of the preconditioning parameter α. So the above diagonalization can never become singular in practice in the subsonic range.

We also remind that by disactivating the preconditioning matrix with In the special case α = 0, we also see that the "shear" and "entropic" eigenvectors are both orthogonal to the acoustic eigenvectors. Taking |n| = 1, with (34) the acoustic eigenvectors are also orthogonal to each others and the diagonalization with unitary congruence (35)-(36) can be applied. So the symmetry properties of the original unpreconditioned system are preserved for the rescaled matrix dissipation. However, this is no longer true when the rescaled matrix dissipation is formulated from the primitive or conservative variables since the eigenvectors are modified by any change of variables. We know from [START_REF] Darmofal | The Importance of Eigenvectors for Local Preconditioners of the Euler Equations[END_REF] that the non normality of the preconditioned system leads to significant robustness problems at stagnation points.

The formulation of the rescaled matrix dissipation expressed for the conservative variables

Stability analysis for the explicit scheme

We now consider the explicit scheme with rescaling of the dissipation matrix derived in the previous section, expressed for the conservative variables

W ∂W ∂t = - [ ∂F ∂x + ∂G ∂y - ∂ ∂x ( Q x ∂W ∂x ) - ∂ ∂y ( Q y ∂W ∂y )] .
In the two-dimensional case, the explicit time-discretization at mesh point x i j = (iδx, jδy), where δx and δy are the constant mesh spacing in the (x, y) directions, can be formulated in conservation

form ∆W i j ∆t + h i+ 1 2 j -h i-1 2 j δx + g i j+ 1 2 -g i j-1 2 δy = 0, (47)
with the time increment ∆W i j = W n+1 i j -W n i j . For the first-order Roe scheme, the numerical flux writes in the x-direction

h i+ 1 2 j = 1 2 [ F(W n i j ) + F(W n i+1 j ) ] - 1 2 | A roe | i+ 1 2 j ( W n i+1 j -W n i j ) . (48)
For the rescaled Roe scheme, we substitute (48) with

h i+ 1 2 j = 1 2 [ F(W n i j ) + F(W n i+1 j ) ] - 1 2 P -1 |P A| i+ 1 2 j ( W n i+1 j -W n i j ) , (49)
with matrix P -1 |P A| formulated for the Roe average. There is no need to explicitly formulate the matrix dissipation for the conservative variables. In practice, the dissipation vector can always be computed using similarity transformations following the algebraic algorithm described in [START_REF] Turkel | Radespiel Preconditioning Methods For Low Speed Flows NASA Contractor Report 201605[END_REF] for the iterative preconditioning

P -1 |P A| i+ 1 2 j ( W n i+1 j -W n i j ) = [ ∂W ∂ W0 R0 Λ R0 -1 ∂ W0 ∂W ] i+ 1 2 j ( W i+1 j -W i j
) .

The stability constraint is strongly concerned with the spectral radius ρ(P -1 |P A|) of the matrixvalued dissipation. We have ρ(P -1 |P A|) = max(µ + , |q|) since µ + > µ -. We found that µ + > |q| when M n → 0 and M n → 1. This is actually true ∀M n ≤ 1 and for α = 0, α = 1 and α = -1 as demonstrated in the following.

Using identity (44) and taking advantage of having both µ ± ≥ 0, we can derive the following expression

(µ + -|q|)(µ -+ |q|) = a 2 |n| 2         1 + M n √ ∆ χ -2M 2 n         . (50)
From (50), we have the sufficient condition that as long as M n ≤ √ 2 2 , then µ + -|q| ≥ 0 and ρ(P -1 |P A|) = µ + , without condition on the free parameter α. We also find that if ∆χ ≥ 0, then

√ ∆ χ ≥ 1 and µ + -|q| ≥ 0 with 1 + M n √ ∆ χ -2M 2 n ≥ 1 + M n -2M 2 n ≥ 0 for all 0 ≤ M n ≤ 1.
We can formulate the quantity ∆χ as follows

∆ -χ = P 0 (1 -M 2 n ) + P 1 M 2 n ,
with

P 0 = (1 - β 2 a 2 ) 2 (1 -M 2 n ) + 2α(1 -α + β 2 a 2 )M 2 n - β 2 a 2 and P 1 = 3 4 (1 -α + β 2 a 2 ) 2
Since P 1 > 0, a sufficient condition for ∆χ ≥ 0 is P 0 ≥ 0. Taking α = 1, we see that as long as M n lies in the range √ 2

2 ≤ M n ≤ 1, P 0 ≥ 0. Thus setting α = 1 ensures µ + -|q| > 0 and then ρ(P -1 |P A|) = µ + for all M n ≤ 1. On the other hand, considering now the cases α = 0 and α = -1, we cannot conclude on the sign of µ + -|q| by using the expression ∆χ. However, using identity (50), we see that

(µ + -|q|)(µ -+ |q|)         1 + M n √ ∆ χ + 2M 2 n         = a 2 |n| 2         1 + 2M n √ ∆ χ + M 2 n ( ∆ -4M 2 n χ χ )         ,
with χ > 0. Following the same idea, we also have

(µ + -ρ(A))(µ -+ ρ(A))         √ ∆ χ + 2M n         = a 2 |n| 2 (1 + M n ) ( ∆ -4M 2 n χ χ ) ,
where ρ(A) = |q| + a|n| is the spectral radius for the unpreconditioned scheme.

With α = 0

∆ -4M 2 n χ = (1 - β 2 a 2 ) 2 (1 -M 4 n ) ≥ 0 for all 0 ≤ M n ≤ 1. With α = -1 ∆ -4M 2 n χ = ( (1 - β 2 a 2 ) 2 -M 2 n ) (1 -M 2 n ) => 1 + M 2 n ( ∆ -4M 2 n χ χ ) > 0.
We then conclude that

ρ(P -1 |P A|) = µ + ∀M n ≤ 1,
for all possible values of the parameter α. Furthermore, if α = 0, we also have µ + ≥ ρ( A) for all M n ≤ 1 which is not always the case with α = 1 or α = -1.

In the framework of a linear Von-Neumann analysis, W n i j is replaced by its Fourier expansion. The Jacobian matrices A and P -1 |P A| are assumed constant matrices and the amplification matrix corresponding to the explicit scheme (47) with the numerical flux ( 49) is given for one space-dimension by

G(ϕ, σ) = Id + (cos(ϕ) -1)σP -1 |P A| -i sin(ϕ)σA,
where ϕ is the wave number (normalized by 1/δx), σ = ∆t δx and i 2 = -1. This amplification matrix is the same as for the preconditioned Lax-Friedrich type schemes studied in [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF]. Let ϵ > 0. Then, there exits some subordinate matrix-norm ||.|| G,ϵ depending on G and ϵ such that

ρ(G) ≥ ||Id + (cos(ϕ) -1)σP -1 |P A| -i sin(ϕ)σ A|| G,ϵ -ϵ.
We also have for any matrix-norm

||A|| ≥ ρ( A) ≃ O( 1 M ) and ||P -1 |P A||| ≥ ρ(P -1 |P A|) ≃ O( 1 M 2 ) >> ρ( A) as M → 0.
Therefore, the Von-Neumann stability condition gives a stringent constraint on the time-step in the incompressible limit as formulated in [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF], with for a fixed δx

1 ≥ ρ(G) ≥ ||G|| G,ϵ -ϵ ≥ ∆t δx ρ(P -1 |P A|) -ϵ ≃ O( ∆t M 2 ) =⇒ ∆t ≃ O(M 2 ) as M → 0. ( 51
)
This proof can be extended to the multi-dimensional formulation in a straightforward way. From (51), we conclude that the stability condition for the explicit rescaled Roe scheme is formulated for the fasted wave speed as follows

∆t ≤ h ρ(P -1 |P A|) , (52)
where h represents some characteristic cell distance. In the two extreme cases, the stability condition yields

∆t ≤ h µ + ≃ O(M 2 ) as M → 0,
and

∆t ≤ h |q| + a|n| as M → 1. ( 53
)
So for steady-state problems, the local time step will be very small ≃ O(hM 2 ) in the lowspeed flow regions only. In supersonic flow regions, the "standard" stability condition (53) must be applied. Therefore, with the strong limitation on the time step, an efficient implicit stage must be added to the baseline explicit rescaled Roe scheme as described in the next section for the steady-state problem.

It is interesting to compare the acoustic eigenvalues µ ± of the rescaled dissipation matrix and the original eigenvalues |λ ± | in the whole subsonic range M n ≤ 1. In the next figures, all eigenvalues are normalized by the local speed of sound and therefore can be expressed only as functions of the dimensionless normal velocity q a|n| ∈ [-1, 1] (= u a in the one-dimensional case). The eigenvalues µ ± strongly depend on the local Mach number M within the parameter β 2 a 2 with a noticeable effect of the second parameter α for low Mach numbers. In Fig. 123we have compared the effect of α = -1, α = 0 and α = 1, considering different values of the local Mach number. As demonstrated above, we see that µ + is the largest eigenvalue in the entire subsonic range with µ + > |q|, ∀α, µ + ≥ ρ(A) for α = 0 and that µ -takes small values especially when M ≤ 0.1. The value of α has little impact on µ -but a large effect on µ + in the incompressible limit, as illustrated in Fig. 1 for M = 0.1. As it can be clearly seen, the stiffness of the rescaled dissipation matrix increases dramatically as the local Mach number is becoming smaller and smaller. The worst situation arises for α = 1, which could not be tested in our computations of low Mach-number flows. Globally α = -1 gives the better conditioning except when q a|n| → 0

where the asymptotic value of µ + is reached. This asymptotic value given by a 2 |β| |n| is independent of α but strongly depends on the local Mach number. When M → 1, the effect of α is almost negligible and the stiffness of the rescaled Roe matrix becomes identical to the original Roe matrix (see Fig. 3 for the case M = 1 with µ + = ρ(A) if α = 0). 5. Implicit matrix time-stepping scheme for the steady-state problem

Formulation

The Jacobi preconditioning or preconditioning-squared is known to smooth out high frequency errors and is especially efficient with the multigrid technique to accelerate the convergence. It basically doesn't help to improve the stiffness of the rescaled Roe scheme in the low speed limit. On the other hand the low-speed preconditioning is not designed to treat high frequency errors. Turkel suggested to combine both approaches. However, here we apply a block Jacobi preconditioner to the rescaled Roe scheme in order to formulate an implicit matrix timestepping scheme for the steady-state problem in the sens that each characteristic variables is updated with its own ∆t (k) [START_REF] Van Leer | Characteristic Time-Stepping or Local Preconditioning of the Euler Equations[END_REF]. Actually the matrix time-stepping approach provides a mechanism enforcing the damping properties of the scheme as it is described next within a linear stability analysis. In the case of the rescaled Roe scheme, the formulation of a time-step matrix is not obvious, since the dissipation matrix is not an explicit function of the original Jacobian matrix.

For the fully discrete scheme in a two-dimensional Cartesian mesh (47) with the numerical flux corresponding to the rescaled Roe scheme (49), the Jacobi preconditioning may be formulated as follows

( A p + B p ρ(P -1 |P A|) + ρ(P -1 |PB|) ∆W ) i j = -∆t        h i+ 1 2 j -h i-1 2 j δx + g i j+ 1 2 -g i j-1 2 δy        , ( 54
)
where the matrix coefficients A p and B p must fulfill the following conditions:

1. The eigenvectors must be the same as the original Jacobian matrix in order to preserve the formulation of the characteristics variables;

2. The eigenvalues must be non negative in order to let unchanged the signs of the characteristic speeds and bounded away from zero (thus positive) for the formulation of a characteristic time-step matrix;

3. The stability condition (52) must be satisfied for the characteristic variable corresponding to the fastest wave speed;

4. The matrix coefficients A p and B p should achieve a smooth transition with |A| and |B| at the sonic line.

Matrix coefficients A p = P -1 |P A| and B p = P -1 |PB| perfectly match conditions 2., 3. and 4.

but must be rejected as not satisfying the first condition. Setting

A p = | A| = R A |Λ A |R A -1
and

B p = |B| = R B |Λ B |R B -1
satisfy conditions 1., 2. and 4., provided that the vanishing eigenvalues |λ (k) | are corrected by some entropy fix. This latter choice reduces the spread of the characteristic speeds from a factor 1/M with respect to the formulation with the rescaled dissipation matrix.

However, condition 3. is lost with ρ(A) ≤ µ + and in practice, this formulation is unstable and can be used only with small or moderate values of the CFL number. A intermediate choice is to build matrix coefficients from the eigenvectors of the original Jacobian matrix and the eigenvalues µ (k) of the rescaled matrix dissipation supposed to be all bounded away from zero. In that case, the eigenvalues µ (k) must be corrected to ensure the transition at the sonic line with |λ (k) | the eigenvalues of the original Roe matrix dissipation. This intermediate formulation provides optimal stability properties and yields a robust scheme but let the spread of the characteristic speeds unchanged with respect to the rescaled Roe scheme. We introduce a characteristic time-step matrix

Hence, let consider

A p = R A ΛA R A -1 and B p = R B ΛB R B -1 , with Λ = diag(µ (k) ),
∆t c = CFL h ( A p + B p ) -1 = [ (R A ∆t A -1 R A -1 ) + (R B ∆t B -1 R B -1 ) ] -1 ,
where ∆t is the diagonal matrix ∆t = diag(∆t + , ∆t -, ∆t 0 , ∆t 0 ) with

∆t + = CFL h µ + * , ∆t -= CFL h µ - , ∆t 0 = CFL h µ 0 * .
The star symbol indicates that the eigenvalues may be bounded away from zero and especially that µ ± must be modified to ensure that the above condition 4 is satisfied. This point is discussed later, in the supersonic flow extension of the rescaled Roe scheme.

Thus, the Jacobi preconditioner amounts to reformulate the original conservative scheme (47) with a time-step matrix

∆W i j = -(∆t c ) i j       h i+ 1 2 j -h i-1 2 j δx + g i j+ 1 2 -g i j-1 2 δy       .
This scheme has optimal properties in damping high frequency errors, and is well adapted to transonic flows. By setting P to the identity matrix Id in (54), we can also formulate a matrix time-stepping scheme for the original Roe scheme with

∆t + = CFL h |λ + | , ∆t -= CFL h |λ -| , ∆t 0 = CFL h |λ 0 | .
and λ ± = q ± a|n|, λ 0 = q. Note that if ∆t -and ∆t 0 are set to the maximum value allowed by the CFL condition, we recover the usual local time step with ∆t c = ∆tId. In the matrix timestep approach, we would like to let the time steps associated to the slowest waves unbounded.

Intuitively, since in the incompressible limit we have for a fixed mesh

∆t + ≃ h CFL O(M 2 ), ∆t -≃ h CFL O(1), ∆t 0 ≃ h CFL O(1),
the original idea behind the formulation of a time-step matrix was to modify the characteristic time-steps so that the solution advance in time along the characteristic curves with time steps of the same order of magnitude. This is actually not suitable, since it is essential to keep the numerical stability of the rescaled Roe scheme. Equalizing all characteristic time steps would give the additional constraint that ∆t -and ∆t 0 must be scaled to the smallest time step ∆t + , and we recover the behaviour of the numerical scheme formulated for the scalar time step. In order to overcome the issue of using very small time steps in the incompressible limit, we see that the matrix time-stepping scheme must be combined with a robust implicit scheme allowing the use of large CFL numbers, as recommended in [START_REF] Birken | Meister Stability of Preconditioned Finite-Volume Schemes at Low Mach Number[END_REF].

A usual way to define an implicit scheme for the Roe scheme, is to consider the following backward Euler scheme as formulated in [START_REF] Mulder | Leer Experiments with Implicit Upwind Methods for the Euler Equations[END_REF] for an approximated differentiation of the Roe flux

-σ x A + i-1 2 j ∆W i-1 j + σ x A - i+ 1 2 j ∆W i+1 j + ( Id -σ x A - i+ 1 2 j + σ x A + i-1 2 j -σ y B - i j+ 1 2 + σ y B + i j-1 2 ) ∆W i j -σ y B + i j-1 2 ∆W i j-1 + σ y B - i j+ 1 2 ∆W i j+1 = ∆W exp i j (55) G(ϕ, θ) and respectively λ (k) A , λ (k) B and µ (k) A , µ (k)
B the eigenvalues of the matrix coefficients A, B and A p , B p , P -1 |P A| and P -1 |PB|. Under this assumption, the time-step matrix has eigenvalues

∆t (k) = CFL h |λ (k) A | + |λ (k) B | (Roe scheme), ∆t (k) = CFL h µ (k) A + µ (k) B
(rescaled Roe scheme).

Important for differentiating the different formulations is to have a close upper estimate of the spectral radius of the amplification matrix. Without loss of generality, we have for the Roe scheme in one space-dimension

ρ(G) = 1 min k (|1 + λ (k) K |) , with |1 + λ K | = √ (1 + 2(1 -cos(ϕ))σ|λ|(1 + σ|λ|)) ≥ 1, ∀σ.
Denoting ν = σ|λ|, we see that ν -→ ρ(G)(ϕ, ν) is a strictly diminishing function and therefore

ρ(G) = 1 √ (1 + 2(1 -cos(ϕ))ν min (1 + ν min )) ,
where

ν min = σmin k (|λ (k) |) if no matrix time-stepping scheme is formulated, ν min = min k (σ (k) |λ (k) |), σ (k) = ∆t (k)
δx with the matrix time-stepping scheme. The optimal rate of convergence to a steady state is then obtained for large CFL numbers with lim CFL→∞ ρ(G) = 0, ∀ϕ ∈]0, 2π[, ∀λ 0. The same results hold for the two-dimensional problem with the upper-bound estimate

ρ(G) ≤ 1 √ (1 + 2(1 -cos(ϕ))(ν min ) x + 2(1 -cos(θ))(ν min ) y ≤ 1, ∀σ x , σ y .
So the higher are (ν min ) x or (ν min ) y , the lower is ρ(G) and the better is the convergence rate to a steady state.

In the case of the rescaled Roe scheme, using similar arguments, we get in the one-dimensional case

ρ(G) = 1 √ (1 + (1 -cos(ϕ))[2 + (1 -cos(ϕ))(ν µ ) min ](ν µ ) min + sin 2 (ϕ)(ν λ ) min , where (ν µ ) min = σmin k (µ (k) ), (ν λ ) min = σmin k (|λ (k) |) if no matrix time-stepping scheme is formu- lated, (ν µ ) min = min k (σ (k) µ (k) ), (ν λ ) min = min k (σ (k) |λ (k) |)
with the matrix time-stepping scheme.

Now comparing the effect of the matrix time-stepping scheme on the non-dimensional quantities (ν min ) x and (ν min ) y for both the Roe scheme and the rescaled Roe scheme, we get -For the Roe scheme:

σ x |λ (k) A | = ∆t δx |λ (k) A | = CFL h δx |λ (k) A | ρ(A) + ρ(B) ≤ CFL h δx |λ (k) A | |λ (k) A | + |λ (k) B | = ∆t (k) δx |λ (k) A | = σ (k) x |λ (k) A | ∀k,
and similarly

σ y |λ (k) B | ≤ σ (k) y |λ (k) B |, ∀k. Assuming that (ν min ) x is reached for k = k 0 , we see that (ν min ) x = σ x |λ (k 0 ) A | ≤ σ x |λ (k) A | ≤ σ (k) x |λ (k) A |, ∀k ⇒ (ν min ) x ≤ min k (σ (k) x |λ (k) A |).
-For the rescaled Roe scheme, the previous result writes

[(ν µ ) min ] x = σ x µ (k 0 ) A ≤ min k (σ (k) x µ (k) A ) and [(ν λ ) min ] x = σ x |λ (k 0 ) A | ≤ min k (σ (k) x |λ (k) A )|.
Therefore the matrix time-stepping scheme always yields higher (ν min ) x and (ν min ) y and thus may improve the damping properties of the iterative procedure. It is also interesting to look at the situation in the incompressible limit. Since for the Roe scheme in a fixed mesh

|λ + | ≃ |λ -| ≃ O( 1 M ) and |λ 0 | ≃ O(1) as M → 0, with find that (ν min ) x ≃ O(M) < min k (σ (k) x |λ (k) |) ≃ O(1) as M → 0.
For the rescaled Roe scheme in the same conditions

µ + ≃ O( 1 M 2 ) and µ -≃ µ 0 ≃ O(1) as M → 0, thus [(ν µ ) min ] x ≃ O(M 2 ) << min k (σ (k) x µ (k) ) ≃ O(1) as M → 0.
So the optimal damping properties of the matrix time-stepping scheme should be enforced in the incompressible limit.

It can also be demonstrated that the rescaled Roe scheme is also second-order dissipative in the sens of Kreiss [START_REF] Lerat | Difference Schemes for Nonlinear Hyperbolic Systems -A General Framework Lecture Notes in Math[END_REF] and can be used for the computation of non-linear flows.

Under relaxation of the implicit scheme

In our implementation of the rescaled Roe scheme, we have considered the finite-volume method in structured grids. The implicit scheme (55) is solved iteratively using a LU-SGS method, with successive forward and backward diagonal sweeps within the mesh. Numerical experiments for low Mach-number flow simulations have shown that a large number of diagonal sweeps should be considered to get an efficient inversion of the linear system. Furthermore, considering that all boundary conditions were not linearized implicitly (only the implicit treatment of the wall and symmetry boundary conditions was implemented) and the error introduced the iterative procedure, this practically limits the use of large CFL numbers required by the stability condition. As shown in the previous section, the efficiency of rescaled Roe scheme to handle low-speed flow computations depends on first the use of large CFL numbers and second the use of a characteristic time-step matrix. This will be illustrated next for the case of the NACA0012 airfoil. Then, it has been found necessary to improve the efficiency the LU-SGS algorithm and a simple under-relaxation of the implicit scheme turned out to very efficient to overcome these limitations.

Let first reformulate the implicit scheme (55) in condensed form

(Id -∆t M) ∆W = ∆t Res exp ,
where ∆t = ∆tId if a scalar time step is used or ∆t = ∆t c for the time-step matrix. We modify this with

(Id -ϕ∆t M) ∆W = ∆t Res exp ,
and we take ϕ > 1, which implies under-relaxation. So large values of ϕ will slow down the convergence rate for the steady problem, but will permit the use large CFL numbers. We experienced that, depending upon the low-Mach number inflow conditions, the relaxation coefficient ϕ should be increased from ϕ = 1.5 at M ∞ = 10 -1 to ϕ = 150 at M ∞ = 10 -3 , in order to obtain a stable algorithm supporting very large CFL numbers, with typically CFL ≃ O(1/M 2 ∞ ).

Practical implementation of the rescaled Roe scheme

In our implementation, the Choi-Merkle preconditioning parameter will assume the following formulation

β 2 a 2 = min(max(M 2 , ϵ 2 ), 1),
where ϵ 2 is used as a small cut-off value. It is worth mentioning that for high subcritical or transonic flows, it is essential to set the small parameter ϵ 2 to a much smaller value than M 2 ∞ to preserve the higher accuracy of the rescaled Roe scheme in the low-speed flow regions. For lowspeed flow computations, the small parameter ϵ 2 was set to a cut-off value for all inflow Mach numbers considered. In [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF], the use of the small parameter ϵ 2 was not mentioned. This doesn't seem to create any numerical difficulties. On the other hand, no spurious dissipation effects due to corrected eigenvalues are introduced. For low Mach-number flow computations, the Van Leer formulation was also considered

β 2 a 2 = min( 2M 2 1 -2M 2 , 1), so β 2 a 2 = 1 as M ≥ 0.5.

Entropy fix

The so-called entropy fix is considered as a mean to select the relevant physical solution satisfying the entropy condition to the hyperbolic problem, which may be violated by the Roe's approximate Riemann solver. The formulations proposed by Harten have been widely accepted by the CFD community, as being very simple to implement in existing CFD codes based on the Roe scheme. These formulations act on all the eigenvalues approaching zero and therefore also correct the contact discontinuity with consequently additional numerical dissipation and loss of accuracy. The behaviour of the modified eigenvalues µ ± and µ 0 = |q| suggests to apply selectively the entropy fix. With µ + > 0, a correction has to be found only for µ 0 at the stagnation point and for µ -at the sonic line, since the formulation of the rescaled matrix prevents µ -from being too small near a stagnation point. With µ -≃ |q| near the stagnation point, this suggests to equalize µ 0 and µ -when M → 0. We suppose the Mach number to be small enough so that the preconditioning parameter recovers its cut-off value β 2 a 2 = ϵ 2 . Then, a possible correction will assume the following form

µ 0 = max(|q|, ϵa|n|) with ϵa|n| = lim M→0 µ -≃ O(1). ( 58
)
This correction might be interpreted as formally similar to the first entropy correction proposed by Harten and Hyman [START_REF] Harten | Self Adjusting Grid Methods for one-dimensional hyperbolic conservation laws[END_REF]. However, notice that the quantity ϵa|n| is actually a bit smaller than the threshold employed for the entropy fix, which is usually based on some fraction of the spectral radius |q| + a|n| (even for the linear waves) and not on the local speed of sound only. For viscous or turbulent flows, as it is discussed below, the formulation of the cut-off value ϵ 2 must be reconsidered, and the correction (58) may assume a different formulation.

Regarding the acoustic eigenvalue µ -at the sonic line when M n = 1, it seems reasonable to be consistent with the formulation of the original matrix dissipation, by considering the same formulation of the entropy fix as used for the Roe scheme. This choice is motivated by the extension to a supersonic flow described in the following. Considering again the symmetrizing entropic variables, in the case of subsonic flow |u| ≤ a we have the following diagonalization

| Ã0 | =                  a a a -a                                   |u + a| 0 0 |u -a|                                   1 2a 1 2a 1 2a -1 2a                  .
Although in case of transonic flow the entropic variables are not appropriate, suppose for the analysis that u → a. Then the acoustic eigenvalue |λ -| = |u -a| may be corrected by some entropy fix and therefore |λ -| → 0 is modified with

| Ã0 | formulated for diag(|u + a|, Ψ (|u -a|)),
where Ψ is some given cut-off function preventing |λ -| from reaching zero. On the other hand, assuming β 2 a 2 → 1 with α = 0 in the expression of the rescaled dissipation matrix as |u| → a, we get

P0 -1 | P0 Ã0 | =                  |u| u u -|u|                                   a + |u| 0 0 a -|u|                                   1 2|u| 1 2u 1 2u -1 2|u|                  .
Thus at the sonic line, µ -= a -|u| → 0 and P0 -1 | P0 Ã0 | may be also reformulated with diag(a + |u|, Ψ p (a-|u|)). The cut-off function Ψ p (µ -) must be defined in order to assure P0

-1 | P0 Ã0 | → | Ã0 |
when |u| → a. In the limit |u| → a, when applying the entropy fix, we can see that we get the identity P0

-1 | P0 Ã0 | = | Ã0 | provided that          Ψ p (µ -) = Ψ (|λ + |) and µ + = |λ -| if u → -a Ψ p (µ -) = Ψ (|λ -|) and µ + = |λ + | if u → a (59)
If no entropy fix is applied

| Ã0 | = P0 -1 | P0 Ã0 | =                  a ±a ±a a                  when u = ±a.
The formulation (59) of the entropy fix for the rescaled Roe scheme is unchanged by similarity transformations and can be also extended to the Euler equations in conservation form. The definition of Ψ p is also formally unchanged in the multidimensional case, as indicated below.

For complex non-linear flows, we may consider the classical Harten's entropy fix formulated in

[27]              Ψ (|λ|) = λ 2 + δ 2 h 2δ h if |λ| < δ h Ψ (|λ|) = |λ| if |λ| ≥ δ h (60)
with δ h being the threshold for the correction of the eigenvalues approaching zero which is usually a fraction of ρ(A).

Supersonic flow

If the flow becomes locally supersonic, the formulation of the rescaled Roe scheme is no longer valid and the original Roe scheme must be selected. It becomes then necessary to introduce a "switch" between the subsonic and supersonic flow conditions. The matrix valued dissipation Q must be formulated as follows

if a|n| -|q| ≥ 0 then Q = P -1 |P A| else Q = | A roe |, with P -1 |P A| = | A roe | if |q| =
a|n| and α = 0 provided that the rescaled dissipation matrix is formulated for the Roe's average. In the supersonic flow regions, the CFL condition must be selected accordingly with (53). As aforementioned, the entropy fix must be applied selectively in order to assure a smooth transition with the Roe scheme near the sonic line. For the rescaled dissipation matrix, the entropy fix must be applied according to (59) with

Ψ p (µ -) = Ψ (|λ ± |) and µ + = |λ ∓ | if |q ± a|n|| ≤ δ h ,
whenever the Harten's entropy fix (60) is considered.

Since µ + → ρ( A) as M → 1, no special treatment has to be considered for the formulation of the scalar time-step. However, this is not true for the eigenvalues as illustrated in Fig. 123.

So in the formulation of the time-step matrix, a smooth transition

A p = R A ΛA R A -1 → | A roe |
can only be achieved by modifying the eigenvalues µ (k) . A possible reformulation of the acoustic speeds µ * ± has been found by considering the "acoustic portion" of decomposition (43). Upon substitution of matrices A and P -1 by their respective eigenvalues, we get

µ * + = Q 11 + a|n| zq √ χ , µ * -= |Q 22 -a|n| zq √ χ |,
where the absolute value on µ * -ensures positivity. The behaviour of the corrected eigenvalues as function of q a|n| is plotted in Fig. 4 for M = 0.1 and α = 0. It can be especially seen that In practice, |q -a|n|| and |q| may locally vanish or be very small and must be bounded away from zero with the entropy fix.

High-order extension

The rescaled Roe scheme can be easily extended to higher orders of space accuracy assuring the TVD property by using the MUSCL extrapolation, since it doesn't affect the formulation of the matrix valued dissipation.

Extension to viscous flows

The local preconditioning approach used to construct the rescaled dissipation matrix (whatever it is iterative or not) may be further developed to handle viscous and turbulent flows. So far, the extension to the Navier-Stokes/RANS equations has been considered by extending the definition of the preconditioning parameter β 2 a 2 in order to enforce the accuracy and especially the robustness of the preconditioned scheme in stagnation point regions and in the boundary layer 40 with high-aspect ratio grids. It has been especially demonstrated in [START_REF] Colin | A Robust Low-Speed Preconditioning Formulation[END_REF] that a better control of the free parameter β 2 a 2 has to be achieved in the low Reynolds number flow regions. The authors suggest the cut-off value ϵ 2 to be related to the isentropic Mach number inside the boundary layer. In [START_REF] Turkel | Robust Low Speed Preconditioning for Viscous High Lift Flows AIAA paper[END_REF], the "inviscid" preconditioning parameter is modified for high lift configurations and a cell Reynolds number is considered as a viscous correction. The necessary corrections of the original formulation are especially relevant for the computation of some industrial configurations characterized by co-existing transonic flows and large three-dimensional separated flows. These modifications obviously do not modify the derivation of the eigenspaces and the formulation of stability condition for the rescaled dissipation matrix, the low-speed preconditioner remaining formally the same.

Computational results

We were basically interested in testing numerically the stability and the accuracy of the rescaled matrix dissipation for the Roe scheme in the two extreme cases M → 0 and M → 1.

Since it was necessary to assess the effect of removing the iterative preconditioning on both the numerical stability and the accuracy of the solution, flow configurations were selected assuming steady conditions. In order to frame the assumption of inviscid flow considered previously and to demonstrate numerically the properties of the rescaled Roe scheme, we have considered a NACA0012 airfoil configuration widely used for the assessment of CFD solvers. The meshes used are the Vasseberg-Jameson meshes, a family of O-meshes generated from a conformal transformation, yielding quadrilateral cells with an aspect ratio of one, with gridlines essentially orthogonal at each vertex within the mesh [START_REF] Vassberg | Jameson In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF]. This mesh topology allows to avoid introducing poor grid quality or highly grid stretching effects in the numerical solutions. Note that the original formulation in the CFD solver of the boundary conditions using characteristic variables has been left unchanged when the rescaled Roe scheme was selected.

Low Mach number flows about the NACA0012 airfoil

The experiment of reference [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF] was reproduced for the reformulation of the rescaled Roe matrix dissipation. In [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF] the rescaled dissipation matrix was computed from the diagonalization of preconditioned Jacobian matrix |P A|, while in our implementation, the diagonalization was formulated from the full matrix P -1 |P A| and the implicit scheme (55) with matrix coefficients (56) was used. The primary interest of this test case is to demonstrate that the proper scaling of the pressure field is recovered at the discrete level in the incompressible limit by the rescaled Roe scheme. The test-case corresponds to a nonlifting solution at zero angle-of-attack. The grid considered for this preliminary assessment has N c = 128 cells in both directions.

For the normalized equations, we know that the discrete solutions of the first-order Roe scheme support pressure fluctuations in space of order M ∞ , with p(x, t) = p 0 (t) + p 1 (x, t)M ∞ , while in the continuous case, the pressure fluctuations scale as M 2

∞ . It was also demonstrated that the rescaled Roe scheme (with α = 0) support pressure fluctuations in space of secondorder p(x, t) = p 0 (t) + p 2 (x, t)M 2 ∞ . Thus, without preconditioning or rescaling of the dissipation operator, a compressible flow solver usually fails to compute low Mach number flows.

A rigorous way to prove this result numerically consists in considering a sequence of computations with decreasing inflow Mach number using the same mesh and investigating whether the pressure fluctuations undergo a linear or quadratic behaviour with the Mach number. Looking for steady solutions, the pressure field will then assume the form p(x) = p 0 + δp(x), with δp(x) = p f (x)M n ∞ , where p 0 is the leading order pressure corresponding to the thermodynamic static pressure, p f is the magnitude of the pressure fluctuations and n = 1 or 2 is the order of the fluctuations. Following [START_REF] Guillard | Viozat On The Behaviour of Upwind Schemes in The Low Mach Number[END_REF], the accuracy of the computations in the incompressible limit can be characterized by the evolution of a normalized pressure field p ∈ [0, 1] and an indicator characterizing the order of the pressure fluctuations as the Mach number goes to zero. The normalized pressure used for the comparison of the pressure fields at different inflow Mach numbers is defined as p

(x) = p(x) -p min p max -p min = p f (x) -(p f ) max (p f ) max -(p f ) min = p f (x),
and then also corresponds to a normalized field for pressure fluctuations. Note that the normalized pressure p is independent from the inflow Mach number and therefore as M ∞ → 0, numerical solutions should also converge to a consistent approximation of the incompressible solution. A normalized pressure fluctuation δ p is also introduced with

δ p = p max -p min p max ≃ [ (p f ) max -(p f ) min p 0 ] M n ∞ for M ∞ << 1.
Since in the normalization process all flow variables are of the same order of magnitude (around unity), the above definition of δ p gives the variation of the pressure fluctuations with the Mach number.

The behaviour p = O(1/M 2 ) for the pressure field was reproduced for a fixed mesh within a sequence of computations by increasing the inflow pressure while keeping constant the inflow velocity and the density, for decreasing Mach numbers M ∞ in the range 10 -1 to 10 -3 . The entropy correction was deactivated for all inflow Mach numbers considered. Numerical experiments are represented in the next figures for the Van Leer formulation of the preconditioning parameter and α = -1. As expected, solutions obtained with the rescaled matrix dissipation converge to a unique isentropic solution with constant density, corresponding to a decoupling of the velocity and pressure fields, which cannot be achieved with the Roe scheme, as illustrated in Fig. 56. The behaviour of the pressure fluctuations with the inflow Mach number is represented in figure 7. With the rescaled Roe scheme, the pressure fluctuations exactly scale with M 2 ∞ . As also indicated, the same experience was conducted for the Choi-Merkle formulation of the preconditioning parameter and α = 0, with exactly the same quadratic behaviour reproduced for the pressure fluctuations. As shown by many authors, we see that the Roe scheme supports fluctuations of order M ∞ in the incompressible limit. Additionally, for inflow Mach numbers lower than 0.1, a grid-convergence study carried out in [START_REF]Viozat Implicit Upwind Schemes for Low Mach Number Compressible Flows INRIA Report[END_REF] rigorously demonstrates that the effort in grid density required to the Roe scheme to consistently converge to the correct pressure field is inversely proportional to the Mach number. This is a consequence of the behaviour of the leading order term of the truncation error in the incompressible limit.

The effects of the CFL number and of the use of the matrix time-stepping scheme are illustrated in Fig 8 for the rescaled Roe scheme. As demonstrated in the stability analysis, the convergence rate is drastically increased when large CFL numbers are used, and the time-step matrix proved to be more efficient than the standard scalar time for any tested CFL number.

The use of large CFL numbers was made possible by under relaxation of the implicit scheme.

Largest CFL numbers correspond to ∆t + ≃ O(h), so typically CFL ≃ O(1/M 2 ∞ ), with CFL = 100 (ϕ = 1) or larger for M ∞ = 10 -1 , CFL = 10 000 (ϕ = 15) for M ∞ = 10 -2 and CFL = 1 000 000 (ϕ = 150) for M ∞ = 10 -3 . In addition, as indicated in the figures with the dashed line, the zero level machine is reached for the pressure fluctuations of order O(M -2 ). Therefore, the zero level machine corresponding to the limit cycle is shifted according to the inflow Mach number, from 13 orders of convergence when M ∞ = 10 -1 for pressure fluctuations of order O(10 -2 ), to 9 orders of convergence when M ∞ = 10 -3 , with pressure fluctuations of order O(10 -6 ).

The effect of a usual grid stretching was also investigated for the rescaled Roe scheme. The However, the effect of the time-step matrix on the convergence rate turned out to be much more important. The convergence rate is more than 30 times faster with the characteristic time-step matrix for M ∞ = 10 -1 , and practically, the zero machine level could not be reached with a reasonable number of time steps using the scalar time step, when M ∞ ≤ 10 -2 .

Transonic flow about the NACA0012 airfoil

A transonic flow condition for M ∞ = 0.8, α = 0 0 was also considered in order to assess the smooth transition of the rescaled Roe scheme with the original Roe scheme at the sonic line, and the effect of the rescaled matrix dissipation at the stagnation point. This was achieved ∞ and α = 0. For both schemes, the entropy fix was activated, with the special treatment described previously for the rescaled Roe scheme, and a second-order MUSCL extrapolation was considered. For a transonic flow, both schemes have a similar conditioning and therefore similar convergence rates were obtained, with zero-level machine convergence obtained for each mesh. Note that with the matrix time-stepping scheme, the Roe scheme also achieves a faster convergence to the steady state compared to the standard scalar time step.

The comparison for the pressure coefficient obtained in the finest mesh is plotted in Fig. 10left. The respective solutions have a very similar discrete shock structure, with the same computed shock location. The effect of the proper scaling of the matrix dissipation has to be found near the stagnation point, where a slightly higher peak pressure is predicted with the rescaled Roe scheme. Then, the suction effect on the leading edge is slightly reduced and a larger pressure drag coefficient is found with the rescaled Roe scheme, for all meshes used in the grid-convergence study. A consequence of a better control of the numerical dissipation can be readily pointed out with the entropy distribution on the airfoil. The effect of the grid refinement on the entropy

Conclusions

A rescaled matrix-valued dissipation was formulated by deriving explicitly the eigenspaces of the Roe-Turkel matrix dissipation. This formulation is interpreted as a rescaling of the Roe scheme for low Mach number flows. When the entropic variables are considered as independent variables, it is demonstrated that the rescaled matrix dissipation is positive semidefinite and can be diagonalized by unitary congruence. So, robustness problems could be significantly alleviated especially at stagnation points. However, our final objective is to developed a more accurate compressible flow solver in low speed flow regions, and the optimal normality properties of the rescaled matrix dissipation are lost by any change of variables. By construction, the rescaled Roe scheme gives the correct asymptotic order for the artificial dissipation terms, and in particular returns the proper quadratic behaviour of the pressure fluctuations in the incompressible limit.

However, as the Mach number goes to zero, the system becomes very stiff and very small time steps must be used, independently of the definition of the local preconditioner used for the formulation of the Roe-Turkel scheme. With the rescaling of the Roe scheme, it is then necessary to completely reformulate the stability condition. A CFL condition was formulated explicitly, based on the spectral radius of the rescaled matrix dissipation. For steady-state problems, it is shown that a robust implicit scheme can be formulated in order to circumvent the severe stability bound occurring for low Mach number flows. The implicit scheme is characterized by the introduction of a characteristic time-step matrix, enforcing the damping properties of the scheme, and under-relaxation, allowing the use of very large CFL numbers. However, as indicated in the introduction, the rescaled Roe scheme is intended to improve locally the prediction of complex high Reynolds number flows characterized by the coexistence of high gradients of the solution in mixed incompressible and compressible regions, arising typically in industrial configurations.

For such complex flows, it is anticipated that the overall stability of the numerical procedure will not be dominated by the limited number of computational cells where the flowfield can be retained as incompressible. In transonic flow conditions, the rescaled Roe scheme ensures a smooth transition with the original Roe scheme, provided that the entropy fix is carefully formulated. A grid-convergence study demonstrates that the discrete shock structure is preserved with the TVD property enforced by the MUSCL extrapolation. As a consequence of the proper scaling of the matrix dissipation, a significant effect can be especially observed for the inviscid flow, with a drastic reduction of the spurious entropy generation at the stagnation point.

β 2 a 2 =

 2 1 with α = 0, expressions (45)-(46) do not return the right and left eigenvector matrix of | Ã0 (n)|, associated to eigenvalues |λ + | = a|n| + q and |λ -| = a|n|q. They actually correspond to the diagonalization of matrix | Ã0 (n)| with the acoustic eigenvalues recast for the module of the normal velocity a|n| + |q|and a|n| -|q|. In the case q -→ -a|n|, we especially see that |λ + | -→ 0 and |λ -| -→ 2a|n| while µ + -→ 2a|n| and µ --→ 0 respectively. This justifies to formulate carefully a possible entropy fix at the sonic line as described later.

W

  is achieved by additional similarity transformations, using the intermediate primitive variables W 0 = [p, u, v, S ] T and applying the change of variables with the similarity matrix ∂W

Figure 1 : 1 Figure 2 :

 112 Figure 1: Dimensionless eigenvalues for the original matrix-valued dissipation (Roe) and for the rescaled Roe scheme when M = 0.1

Figure 3 :

 3 Figure 3: Dimensionless eigenvalues for the original matrix-valued dissipation (Roe) and for the rescaled Roe scheme when M = 1

  R A and R B being the left eigenvector matrices of the original Jacobian matrices A and B. In a Cartesian mesh, the CFL condition (52) can be expressed with a more restrictive formulation of the local time-step ∆t = CFL h ρ(P -1 |P A|) + ρ(P -1 |PB|).

  ± | the eigenvalues of |A roe | as q → ±a|n| and that the asymptotic behaviour of µ * + ≃ µ + as |q| → 0 is reproduced by the correction. Note that µ * -≃ |λ -| except when |q| → 0 where µ * -becomes singular. So the final expression of µ * -is further simplified in order to avoid any singularity, and we simply take µ * -= |λ -| in the whole subsonic range |q| ≤ a|n|. So for the rescaled Roe scheme, the time-step matrix will assume the following expression for all inflow Mach numbers ∆t + = CFL h Q 11 + a|n| zq √ χ , ∆t -= CFL h |q -a|n|| , ∆t 0 = CFL h |q| .

Figure 4 :

 4 Figure 4: Dimensionless corrected eigenvalues µ * ± in the definition of time-step matrix when M = 0.1 -Grey delta symbols: µ * + . Black diamond symbols: µ * -

Figure 5 :Figure 6 :

 56 Figure 5: Velocity contours with the inflow Mach number in the Vasseberg-Jameson mesh. Left: rescaled Roe scheme. Right: Roe scheme (V ∞ = 270m/s)

Figure 7 :

 7 Figure 7: Pressure fluctuations with the inflow Mach number for the Roe scheme and the rescaled Roe scheme

Figure 10 :

 10 Figure 10: Comparison of the flow solution for the Roe scheme and the rescaled Roe scheme at M ∞ = 0.8, α = 0 0 . Left: pressure coefficient. Right: entropy distribution

Acknowledgements

The author is grateful to his colleagues Lionel Castillon, Christelle Wervaecke and Antoine Dumont for providing the 2D CFD code together with many useful test-cases. The author would also like to thank Daniel Destarac for providing the Vasseberg-Jameson conformal grids.

where ∆W exp i j is the time increment for the explicit scheme, σ x = ∆t δx , σ y = ∆t δy and A + , A -, B + , B -are respectively the positive and negative parts of the flux Jacobian matrices A and B in the (x, y) space directions. For the rescaled Roe scheme, a baseline implicit scheme is then formulated by replacing in (55) respectively

The implicit scheme (55) can be also formulated for the characteristic matrix time-stepping scheme, upon the following substitutions

Stability of the implicit scheme

In the following, the stability of the implicit rescaled Roe scheme with scalar and matrix time-steps is considered. Both versions correspond actually to small perturbations of the original implicit Roe scheme (55) in the transonic regime. Both are especially unconditionally stable for the linear problem as shown below. But as the free-stream Mach number is decreasing, their dissipative properties substantially differ from the original Roe scheme. This can be pointed out within a Fourier analysis. Let (ϕ, θ) be the normalized wave numbers in the two space directions.

Then, for the linear problem the amplification matrix takes the general form

with for the implicit scheme H(ϕ, θ) = Id + K(ϕ, θ). We get successively -For the original Roe scheme

-For the rescaled Roe scheme

-For the the matrix time-stepping schemes, we replace σ x and σ y according to (57).

Stability and damping properties of the four schemes can be characterized by the following estimates, assuming that all matrix coefficients simultaneously diagonalizable according to the Lerat's general framework [START_REF] Lerat | Multidimensional Centered Schemes of the Lax-Wendroff Type Computational Fluid Dynamics Review[END_REF]. Let introduce λ (k) G the eigenvalues of the amplification matrix Right: Convergence history for the scalar (dashed lines) and the matrix (solid lines) time steps