Yann Moguen
email: yann.moguen@univ-pau.fr

Erik Dick

Diffusion and dissipation in acoustic propagation simulation by convection-pressure split algorithms in all Mach number form

Keywords: Convection-pressure split algorithms, All Mach number schemes, Co-located finite volume method, Convective transport, Acoustic propagation

The topic of the paper is accuracy analysis of acoustic propagation simulation in low Mach number flows, by finite volume co-located discretisation methods of the time-dependent compressible fluid Euler equations that use the concept of convection-pressure splitting (CPS). These are algorithms that split the flux vectors into a part associated to the convection by the fluid particles, and a part associated to the propagation of the pressure waves. For the convection part, the appropriate space discretisation is the upwind one. For the pressure part, there are alternatives. We discern five types of algorithms that all are adapted for use in low Mach number flows, and thus are considered as all Mach number algorithms. We study the behaviour of the different types for the propagation of small pressure perturbations, of discontinuous or smooth shape, in low Mach number flows. We demonstrate that four of the proposed algorithms of convection-pressure split type are dissipative for such applications, although they are designed for low Mach number flows. The objective of the paper is to analyse why some algorithms are appropriate for acoustic propagation simulation and why some are not appropriate.

Introduction

Within the co-located finite volume method, there are two main categories of convectionpressure split algorithms. One is by construction of face variables by which the pressure flux vectors at the volume faces are composed. The other one is by direct construction of pressure flux vectors at the faces. Both methodologies belong to the family of approximate Riemann solvers. Based on specific choices, existing algorithms can be grouped into five types, which we detail in the paper. Because of their origin of approximate Riemann solvers, these algorithms are not a priori suited for low Mach number problems. An aspect with much attention in the literature is the adaptation for compressible fluid flow in the low Mach number limit towards incompressible fluid flow. The pressure wave propagation mechanism then disappears, which thus necessitates adaptation of the discretisation of the pressure flux part, in particular targeted to low diffusion. An aspect with much less attention in the literature is that the resulting adapted algorithms for use in all Mach number regimes also can be used for acoustic propagation in low Mach number flows. But there is no guarantee that these algorithms are accurate for such applications.

The paper starts with a short discussion of the transport mechanisms of the time-dependent Euler equations for a compressible fluid. This discussion is included because of its importance for understanding the algorithms, although it can be found in many text books on computational fluid mechanics, e.g. [START_REF] Wesseling | Principles of Computational Fluid Dynamics[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. Then, a short discussion follows on the behaviour of the equations in the low Mach number hydrodynamic and acoustic limits. Next is the concept of convection-pressure splitting, followed by the used implicit pressurecorrection time stepping method. The subsequent sections are on the choice of the convection velocity and the discretisation of the pressure part of the equations. The adaptation to low Mach number flows and the stability of the numerical solution of high Mach number flows with strong shocks are discussed. Five particular algorithms are chosen, which are representative for the possible choices. Then follow sections with numerical results. It is demonstrated by one-dimensional examples that all algorithms produce accurate results for shocks, contact discontinuities and expansion fans in steady and unsteady flows at medium-high Mach number. It is further demonstrated that all algorithms are accurate for low Mach number Riemann problems, but that only one is accurate for the propagation of a smooth small pressure perturbation in a low Mach number flow. An explanation is formulated for the dissipative solution by four of the algorithm types. The crucial ingredient is the pressure diffusion term in the expression of the velocity at the cell faces. The role of this diffusion is illustrated by a numerical one-dimensional experiment. Finally, the good functioning of one of the algorithms is illustrated for propagation of a pressure perturbation in a two-dimensional flow with very low Mach number.

The transport mechanisms in the Euler equations

We illustrate the propagation properties for one-dimensional flow. Most numerical results, later in this paper, are also for one-dimensional flow. The 1-D system of Euler equations is

∂ t        v E        + ∂ x        v vv + p Hv        = 0, (1)
where , v, p, E, H represent density, velocity, pressure, total energy and total enthalpy. E = e + v 2 /2, H = h + v 2 /2 and h = e + p, where e and h are the internal energy and the enthalpy. We further assume an ideal and perfect gas, satisfying as equations of state p = RT , h = c p T , and e = c v T , with γ = c p /c v , where R = c p -c v is the gas constant, T is the temperature, c p and c v are the constant specific heats and γ the ratio of the specific heats.

The conservation equations may be expanded into quasi-linear equations,

∂ t + v∂ x + (∂ x v) = 0, ∂ t v + v∂ x v + ∂ x p = 0, (2)
∂ t p + v∂ x p + γp(∂ x v) = 0.
By (∂ x v) is indicated that this derivative represents, in multiple dimensions, the divergence of the velocity. The other x-derivative is the gradient.

For an ideal and perfect gas, the velocity of sound is c = γp/ , so that γp = c 2 . With respect to the vector of variables W = [ , v, p] T , these equations form the system ∂ t W + A∂ x W = 0, with the system matrix

A =        v 0 0 v 1 0 γp v        =        v 0 0 v 1 0 c 2 v        . (3)
The eigenvalues are v -c, v, and v + c. The associated right and left eigenvector matrices are

R = L -1 =        1 -c 0 c c 2 0 c 2        , L = R -1 =        0 -1 2c 1 2 c 2 1 0 -1 c 2 0 1 2c 1 2 c 2        . (4)
The coefficients of the transported combinations of increments of the variables ( , v, p) are given by the left eigenvectors. The transported quantities are:

• with velocity v: δ -1 c 2 δp or δp p -γ δ , which means entropy ;

• with velocity v -c: -δv + δp c ;

• with velocity v + c: δv

+ δp c . (6)
The transport with velocity v is by the fluid particles. Such transport is called advection or convection. The quantities [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit: A review[END_REF] and ( 6) are combinations of velocity perturbations and pressure perturbations that are transported at the velocity of sound, superposed to the motion of the fluid particles. The corresponding transport is called acoustic propagation and it is said that the transported quantities move by waves.

In this paper, we study the behaviour of discrete forms of the Euler equations (1) at low Mach number. For a general discussion on low Mach number aspects, we refer to Guillard and Viozat [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF], Guillard and Murrone [START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes[END_REF], and Guillard and Nkonga [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit: A review[END_REF]. Hereafter, we illustrate the basic limit behaviour of the Euler equations, which is sufficient for the purposes in this paper.

For analysis of the limit behaviour at low Mach number, reference values of density r , velocity v r (thought of as convective), pressure p r are introduced, from which a reference Mach number is defined by M r = v r / p r / r . Reference values of length, l r , and duration, t r , are also used, allowing for the definition of a reference Strouhal number by St r = (l r /v r )/t r . Thus, the non-dimensional quasi-linear continuity, momentum and energy equations are

St r ∂ t + v∂ x + (∂ x v) = 0, St r ∂ t v + v∂ x v + 1 M 2 r ∂ x p = 0, (7)
St r ∂ t p + v∂ x p + γp(∂ x v) = 0.

To allow for a single length-scale, multiple time-scale analysis, a time variable able to capture the fast acoustic waves is introduced by

τ = t M r .
It is then assumed that density can be expanded as

(x, t, M r ) = K k=0 M k r (k) (x, t, τ ) + o(M K r ), K = 0, 1, 2, M r → 0 ,
and that similar expansions hold for velocity v and pressure p. Substitution of these ex-pansions in Eqs. [START_REF] Zha | Numerical solution of Euler equations by using a new flux vector splitting scheme[END_REF] gives, at the convective time scale,

St r ∂ t (0) + v (0) ∂ x (0) = 0, St r (0) ∂ t v (0) + (0) v (0) ∂ x v (0) + ∂ x p (2) = 0, (8) (∂ x v (0) ) = 0,
where it is assumed, for simplicity of exposition, that the zeroth-order pressure p (0) (t), which is the thermodynamic pressure, is constant in space. Eqs. [START_REF] Liou | A new flux splitting scheme[END_REF] are the equations of the hydrodynamic low Mach number limit of the Euler equations (2) in quasi-linear form. By this limit is meant that the flow is considered on a time scale adapted to the convection velocity at a low Mach number obtained by a large velocity of sound. In particular, an incompressible flow obeys these equations. There is only transport by convection of density and velocity and there is no propagation by waves; but pressure has to satisfy the compatibility imposed by the divergence condition on the velocity. So, it is clear that algorithms designed for the full set of Euler equations (2) cannot function without adaptations for the reduced set [START_REF] Liou | A new flux splitting scheme[END_REF] in the hydrodynamic low Mach number limit.

At the acoustic time scale, the zeroth-order density (0) (x) is constant in time, and the linear acoustic equations are obtained,

St r (0) ∂ τ v (0) + ∂ x p (1) = 0, (9a) St r ∂ τ p (1) + γp (0) (∂ x v (0) ) = 0. (9b)
By the acoustic limit is meant that the flow is considered on a time scale adapted to the acoustic propagation at a low Mach number obtained by a small convection velocity. In the hydrodynamic low Mach number limit, the pressure is the "hydrodynamic" pressure M 2 r p (2) , which thus changes, by a change of the velocity, proportional to the square of the velocity. In the acoustic low Mach number limit, the pressure is the "acoustic" pressure M r p (1) , which changes linearly proportional to a change of the velocity.

The transport mechanisms in the conservation equations

The transport mechanisms in the conservative form of the set of equations ( 1) may be identified by recomposing these equations from the expanded equations [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. But, they may also be identified directly on the set of conservation equations by an analysis in the style of the flux vector splitting analysis of Steger and Warming [START_REF] Steger | Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods[END_REF].

As intermediate variables we consider now

W = [ , m, p] T , where m = v. Then E = e + 1 2 v 2 = 1 γ-1 p + 1 2 m 2 and H = h + 1 2 v 2 = γ γ-1 p + 1 2 m 2 .
The vector of the conservative variables and the flux vector are

U =        m 1 2 m 2 + 1 γ-1 p        , F =        m m 2 + p 1 2 m 3 2 + γ γ-1 mp       
.

The Jacobians are T = ∂ W U and A = ∂ W F . From the homogeneity of degree one of both U and F with respect to the intermediate variables, it follows that U = T W and

F = AW . Thus F = AT -1 U = T T -1 AT -1 U = T BT -1 U with B = T -1 A. The matrices T and B are T =        1 0 0 0 1 0 -1 2 v 2 v 1 γ-1        , B =        0 1 0 -v 2 2v 1 -vc 2 c 2 v       
.

The eigenvalues of the matrix B are v -c, v, and v + c. The corresponding right and left eigenvector matrices are

R = L -1 =        1 1 1 v -c v v + c c 2 0 c 2        , L = R -1 =        v 2c -1 2c 1 2c 2 1 0 -1 c 2 -v 2c 1 2c 1 2c 2       
.

A part of the matrix B can now be associated to convection by the part (v, v, v) of the eigenvalues. Similarly, a part can be associated to acoustic propagation by the part (-c, 0, +c) of the eigenvalues. This second part is usually called the pressure part. The parts of B corresponding to convection and to pressure are

B C = R        v 0 0 0 v 0 0 0 v        L =        v 0 0 0 v 0 0 0 v        , B P = R        -c 0 0 0 0 0 0 0 c        L =        -v 1 0 -v 2 v 1 -vc 2 c 2 0       
.

The corresponding parts of the flux vector are

F C = T B C W =        v vv v 3 2 + vp γ-1        =        v E        v , F P = T B P W =        0 p pv        . ( 10
)
This splitting of the flux vector in convection and pressure parts is the one used by Zha and Bilgen [START_REF] Zha | Numerical solution of Euler equations by using a new flux vector splitting scheme[END_REF]. These authors claimed a new way of splitting, but, actually, it was already derived by Steger and Warming [START_REF] Steger | Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods[END_REF]. But Steger and Warming made only a side remark, because their emphasis was on splitting in positive and negative parts and not in convection and pressure parts. The practical use of convection-pressure splitting became only clear by the work of Liou and Steffen [START_REF] Liou | A new flux splitting scheme[END_REF]. The splitting by [START_REF] Moguen | Pressure-velocity coupling allowing acoustic calculation in low Mach number flow[END_REF] is commonly called the Zha-Bilgen splitting. We follow this practise and denote it by ZB-splitting. The system of Euler equations ( 1) with the ZB-splitting is

∂ t        v E        + ∂ x               v E        v        + ∂ x        0 p pv        = 0. (11)
There are two other ways of convection-pressure splitting (CPS) in the literature. One is the splitting by Liou and Steffen [START_REF] Liou | A new flux splitting scheme[END_REF], called AUSM, which stands for Advection Upstream Splitting Method, and one by Toro and Vázquez-Cendón [START_REF] Toro | Flux splitting schemes for the Euler equations[END_REF]. These are:

∂ t        v H        -∂ t        0 0 p        + ∂ x               v H        v        + ∂ x        0 p 0        = 0, (12)
∂ t        v 1 2 v 2        + ∂ t        0 0 1 γ-1 p        + ∂ x               v 1 2 v 2        v        + ∂ x        0 p γ γ-1 pv        = 0. ( 13
)
We denote these ways of splitting by LS-splitting and TV-splitting. They are somewhat more complex than ZB-splitting, because they involve also splitting of the time-derivative part. The properties of the convective parts may be analysed as for the ZB-splitting. But, by similarity with the convective part of the ZB-splitting, one sees that the three propagation speeds for the LS-splitting are (v, v, v). In the convective part of the TV-splitting, the third equation is dependent on the two other ones and, again by analogy, one sees that the propagation speeds of the system of the first two equations are (v, v). So, the convective parts of the three ways of splitting have a similar character. The experience of the present authors is that the three ways of splitting do not react very differently to specific discretisation choices in time and space. In this work, we use the ZB-splitting, because of its simplicity, and we do not study the effects of the differences of the ways of splitting. By the analogy of the convection parts in the three ways of splitting, we apply some ways of discretisation illustrated in the literature for the LS-splitting or the TV-splitting to the ZB-splitting.

4 The implicit pressure-correction time-stepping method

We use the implicit pressure-correction method from our previous work [START_REF] Moguen | Pressure-velocity coupling allowing acoustic calculation in low Mach number flow[END_REF][START_REF] Moguen | A combined momentum-interpolation and advection upstream splitting pressure-correction algorithm for simulation of convective and acoustic transport at all levels of Mach number[END_REF] on the combination of momentum interpolation and LS-splitting (AUSM), but adapted to the ZB-splitting. A cell-centre finite volume space discretisation is used. With the convectionpressure splitting [START_REF] Moguen | A combined momentum-interpolation and advection upstream splitting pressure-correction algorithm for simulation of convective and acoustic transport at all levels of Mach number[END_REF], the convection part is discretised upwind for the transported quantities, which requires the definition of a face velocity. Discretisation of the pressure part either requires the definition of face velocity and face pressure, i.e. determination of a face state, or an expression of the pressure flux vector at a face. We illustrate the pressurecorrection method applied to an algorithm with determination of a state, but adaptation to an algorithm with determination of a flux vector can easily be done. So, we assume that the velocity and the pressure at faces with subscript i + 1/2 on iteration level k are available. We will detail possible choices in subsequent sections. The algorithm is presented with time integration by first-order backward differencing (BDF1). This time integration is also used in the numerical tests shown later, except for the acoustic wave tests, where second-order backward differencing (BDF2) is used. Again, adaptation is easy. Secondorder space discretisation is obtained by second-order accurate extrapolated values on the left-hand side and right-hand side of a cell face, denoted by (•) L and (•) R . Specifically, for the face i + 1/2:

φ L = φ i + 1 2 ψ i (φ), φ R = φ i+1 - 1 2 ψ i+1 (φ),
with the symmetrically limited difference:

ψ i (φ) = MinMod(φ i -φ i-1 , φ i+1 -φ i ).
We assume positive convection velocity so that for a convectively transported quantity, the (•) L value is taken. We set τ = ∆t/∆x. A predicted quantity is denoted by the superscript and update values by the superscripts and . Time levels are denoted by the superscripts n and n + 1. The iteration level is by the superscript k. At the first iteration, values at level k = 0 are these of the time level n.

Prediction step

• p i = p k i (no update of the pressure).

• i from the continuity equation by

i -n i + τ i + 1 2 ψ i ( k ) v k i+1/2 -τ i-1 + 1 2 ψ i-1 ( k ) v k i-1/2 = 0. (14)
• ( v) i from the momentum equation by

( v) i -( v) n i + τ ( v) i + 1 2 ψ i [( v) k ] v k i+1/2 -τ ( v) i-1 + 1 2 ψ i-1 [( v) k ] v k i-1/2 + τ (p k i+1/2 -p k i-1/2 ) = 0. ( 15
)
• v i = ( v) i i , (
E) i = p k i γ -1 + 1 2 [( v) i ] 2 i .
The updates of density and momentum may be written as

i = i -k i and ( v) i = ( v) i -( v) k i .
The equations ( 14) and ( 15) may be written as

(1 + τ v k i+1/2 ) i -τ v k i-1/2 i-1 + ∆t Res (Cont) k i = 0, (16)
(1 + τ v k i+1/2 )( v) i -τ v k i-1/2 ( v) i-1 + ∆t Res (Mom) k i = 0, (17)
where Res (Cont) k i and Res (Mom) k i are the residuals of the discretised continuity and momentum equations at the iteration level k. Thus

∆t Res (Cont) k i = k i -n i + τ k i + 1 2 ψ i ( k ) v k i+1/2 -τ k i-1 + 1 2 ψ i-1 ( k ) v k i-1/2 , ∆t Res (Mom) k i = ( v) k i -( v) n i + τ ( v) k i + 1 2 ψ i [( v) k ] v k i+1/2 -τ ( v) k i-1 + 1 2 ψ i-1 [( v) k ] v k i-1/2 + τ (p k i+1/2 -p k i-1/2 ).

Momentum correction

With ( v) = ( v) k+1 -( v) and p = p k+1 -p k , we derive from the momentum equation ( 17) at cell i (combined update of momentum and pressure):

1 + τ v k i+1/2 ( v) i + τ (p i+1/2 -p i-1/2 ) = 0. Thus ( v) i = -α k i (p i+1/2 -p i-1/2 ), with α k i = τ 1 + τ v k i+1/2 . ( 18
)
At a face i + 1/2 we postulate

( v) i+1/2 = -α k i+1/2 (p i+1 -p i ), with α k i+1/2 = τ 1 + τ v k i+1 . ( 19
)

Correction step

In the energy equation, no correction is done for the kinetic energy. So,

( E) k+1 -( E) = ( e) k+1 -( e) k = 1 γ -1 p .
The convection flux term is written as

( vE) k+1 = 1 2 v 2 v + ( ev) k+1 .
The correction is derived by (double correction):

( ev) k+1 ( ev) + e ( v) k+1 -( v) + v ( e) k+1 -( e) k .
The correction of the pressure flux term is derived by (single correction):

(pv) k+1 (pv) + p k ( v) k+1 -( v) .
Thus, the equation for pressure correction, derived from the energy equation, is

1 ∆t 1 γ -1 (p k+1 -p k ) + ∂ x e ( v) k+1 -( v) + ∂ x v ( e) k+1 -( e) k + ∂ x p k ( v) k+1 -( v) + 1 ∆t [( E) -( E) n ] + ∂ x 1 2 v 2 v + ∂ x ( ev) + ∂ x (pv) = 0.
By setting the face values of velocity back to the level k, the pressure-correction equations are written as

C i-1 p i-1 + C i p i + C i+1 p i+1 + ∆t Res (Ener) i = 0, (20)
with

C i-1 = - τ γ -1 v k i-1/2 -τ α k i-1/2 e i-1 -τ α k i-1/2 p k i-1/2 i-1/2 , C i = 1 γ -1 + τ γ -1 v k i+1/2 + τ (α k i+1/2 e i + α k i-1/2 e i-1 ) + τ   α k i+1/2 p k i+1/2 i+1/2 + α k i-1/2 p k i-1/2 i-1/2   , C i+1 = -τ α k i+1/2 e i -τ α k i+1/2 p k i+1/2 i+1/2 , ∆t Res (Ener) i = [( E) i -( E) n i ] + τ 1 2 [( v) i ] 2 i + 1 2 ψ i 1 2 [( v) k ] 2 k v k i+1/2 -τ 1 2 [( v) i-1 ] 2 i-1 + 1 2 ψ i-1 1 2 [( v) k ] 2 k v k i-1/2 + τ ( e) k i + 1 2 ψ i ( e) k v k i+1/2 -τ ( e) k i-1 + 1 2 ψ i-1 ( e) k v k i-1/2 + τ (v k i+1/2 p k i+1/2 -v k i-1/2 p k i-1/2 )
,

where e i = ( e) k i i = 1 γ -1 p k i i .
The coefficients of the updates in the equations ( 16), [START_REF] Liou | A sequel to AUSM: AUSM +[END_REF] and [START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[END_REF] require the face velocities at the iteration level k. The face velocities are defined in any convection-pressure splitting algorithm. We use the values of the particular method. The C-coefficients of the pressure updates in the energy equation [START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[END_REF] require also face values of density and pressure. For these quantities, we use arithmetic averages, thus i+1/2 = 1 2 ( i + i+1 ) and p i+1/2 = 1 2 (p i + p i+1 ). The residuals of the equations are expressed according to the particular method.

Updates

p k+1 i = p k i + p i , k+1 i = i , ( v) k+1 i = ( v) i + ( v) i = ( v) i -α k i (p i+1/2 -p i-1/2 ), with p i+1/2 = 1 2 (p i + p i+1 ), v k+1 i = ( v) k+1 i k+1 i , ( E) k+1 i = ( E) i + p i γ -1 .

Determination of the face state with characteristic equations

With an approximate Riemann solver by determination of a face state, a possible approach is to determine the face velocity in the convection part of Eq. ( 11) and the face velocity and face pressure in the pressure part of Eq. ( 11) by the quantities in the star-region in between the characteristic paths by the extreme propagation velocities. The characteristic paths are sketched in Fig. 1. The sketch is made with positive convection speed and convection speed lower than the speed of sound. The primitive variables are used for state definition. The star-region is the region between the leftmost and rightmost characteristic paths. These are paths of pressure waves, either expansions or compressions, across which density, velocity and pressure change. The star-region consists of two parts with states W * L and W * R divided by a path across which a density jump is possible: The states W L and W * L can be connected by a line parallel to a characteristic path with speed v + c. Then, by taking frozen values of and c on the left-hand state,

* L = * R , but
p * + L c L v * = p L + L c L v L . (21)
Similarly, the states W R and W * R can be connected by a line parallel to a characteristic path with speed v -c. Then,

p * -R c R v * = p R -R c R v R . (22)
The solution of the set ( 21) and ( 22) is

v * = L c L v L + R c R v R L c L + R c R - p R -p L L c L + R c R , (23)
p * = R c R p L + L c L p R L c L + R c R - L c L R c R L c L + R c R (v R -v L ). (24)
By taking the sum of Eqs. ( 21) and ( 22), the expression of the star-pressure may also be written as

p * = 1 2 (p L + p R ) - 1 2 R c R (v R -v * ) + 1 2 L c L (v L -v * ). ( 25
)
6 Determination of the face state with Rankine-Hugoniot equations By integration of the Euler equations on the control volume at the left-hand side of the t-axis, the balance is

(U * L -U L )(0 -S L ) + F * L -F L = 0 or F * L -F L = S L (U * L -U L ). ( 26
)
The balance on the rightmost control volume is

(U * R -U R )(S R -S * ) + 1 - S * S R (F R -F * R ) = 0 or F * R -F R = S R (U * R -U R ). ( 27
)
The balance on the middle control volume is

(U * L -U * R ) 1 - S * S R (S * -0) + 1 - S * S R (F * R -F * L ) = 0 or F * R -F * L = S R (U * R -U * L ). ( 28
)
The equations ( 26)-( 28) are called the Rankine-Hugoniot conditions across the paths. They are derived here for particular choices of the signs of v -c, v and v + c, but, clearly, they are generally valid.

Eq. ( 28) has an exact solution for velocity and pressure:

v * L = v * R = S * = v * and p * L = p * R = p * .
The set of equations ( 26) and ( 27) do not generally allow exact solutions for the states U * L and U * R with a priori chosen approximations of S L and S R . For exact solutions, S L and S R have to be eigenvalues of the Jacobian matrices relating the flux vector differences and the state vector differences. With K either L or R, the Jacobian matrices

A * K are obtained from F * K -F K = A * K (U * K -U K ).
The eigenvalues may be expressed with Roeaveraged state variables, but their expression needs knowledge of the state U * K , which is not available at the start of a discretisation. Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF][START_REF] Toro | The HLLC Riemann solver[END_REF] determines the star states from the continuity equation and the momentum equation, thus not satisfying the energy equation. These equations are

* K v * -K v K = S K ( * K -K ) or * K (S K -v * ) = K (S K -v K ), ( 29
) * K v * v * -K v K v K + p * -p K = S K ( * K v * -K v K ). (30)
The momentum balance may be written as

* K v * (S K -v * ) -K v K (S K -v K ) = p * -p K .
After substitution of the mass balance, this is

K (S K -v K )(v * -v K ) = p * -p K . ( 31
)
Because S L is an approximation of v L -c L and S R is an approximation of v R + c R , the equations ( 31) are similar to the characteristic equations ( 21) and [START_REF] Sun | A robust flux splitting method with low dissipation for all-speed flows[END_REF]. The solution of the set of the two equations ( 31) is similar to ( 23) and ( 24): 23) and ( 24) are recovered. Again, the pressure may be expressed by adding the two equations [START_REF] X.-D. Ren | Role of the momentum interpolation mechanism of the Roe scheme in shock instability[END_REF] into

v * = L (v L -S L )v L + R (S R -v R )v R L (v L -S L ) + R (S R -v R ) - p R -p L L (v L -S L ) + R (S R -v R ) , (32)
p * = R (S R -v R )p L + L (v L -S L )p R L (v L -S L ) + R (S R -v R ) - L (v L -S L ) R (S R -v R ) L (v L -S L ) + R (S R -v R ) (v R -v L ). (33) With S R = v R + c R and S L = v L -c L , the equations (
p * = 1 2 (p L + p R ) - 1 2 R (S R -v R )(v R -v * ) + 1 2 L (v L -S L )(v L -v * ). ( 34
)
The expressions of the star values of velocity and pressure are now in a more general form than with the method in the previous section, but the obtained equations are still characteristic equations. A possible choice of the propagation speeds of the pressure waves is

S R = max {v R + c R , v M + c M } , S L = min {v L -c L , v M -c M } , (35)
where v M and c M are quantities of an intermediate state between the states L and R. These allow taking into account the variation of the wave speeds across the wave paths. This way, it is avoided that an expansion fan is seen as a discontinuity, thus as an expansion shock. Typically, Roe-averages are used for the intermediate state, but AUSM-averages may be used as well.

In order to simplify the writing of the expressions, we define, as done by Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF],

α R = R (S R -v R ) > 0 and α L = L (v L -S L ) > 0. ( 36
)
The expressions of the star-values become

v * = α L v L + α R v R α L + α R - p R -p L α L + α R , (37)
p * = α R p L + α L p R α L + α R - α L α R α L + α R (v R -v L ), (38)
p * = 1 2 (p L + p R ) - 1 2 α R (v R -v * ) + 1 2 α L (v L -v * ). ( 39
)
7 Convection-pressure splitting with characteristic equations

An example of a method with convection-pressure splitting and face values of velocity and pressure by characteristic equations is the one by Toro and Vázquez-Cendón [START_REF] Toro | Flux splitting schemes for the Euler equations[END_REF]. They take the face velocity and face pressure equal to the star-quantities [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF] and [START_REF] Xie | On numerical instabilities of Godunov-type schemes for strong shocks[END_REF], with α L = L c L and α R = R c R . The convection part is expressed with the face velocity and and upwinded vector of conserved variables by choosing the state L or R depending on the sign of the transporting velocity. Hereafter, we use the term simple upwinding for this type of upwinding. With the face velocity and pressure, the pressure part of the equations is expressed. We take this method as a basic one, but for testing we use a more elaborated version, discussed in the next section.

Toro and Vázquez-Cendón maintain the expressions ( 37) and ( 38) for positive v -c or negative v + c. But, it seems more accurate to use face values from the leftmost state (v L and p L ) for v -c > 0 and from the rightmost state (v R and p R ) for v + c < 0. They do not consider this possibility. Their argument is that the propagation speeds in the pressure part of the equations are -c, 0, c, thus with fixed negative sign for the leftmost wave and fixed positive sign for the rightmost one. They interpret the guaranteed subsonic character of the pressure part of the equations as a justification for taking the star values as face values in a convection-pressure split algorithm. They even use a reduced value of the speed of sound, c = (γ -1)p/ , because -c , 0, c , are the propagation speeds of the system of equations with the convection part set to zero. We follow the choice of the face values by the star values, because it is consistent in the sense that these values are in between the left-hand (L) and right-hand (R) values at the face. We also follow the choice for simple upwinding, because, again, it is consistent. But, based on the analysis of the propagation speeds in section 3, we take the physical expression of the speed of sound.

The current study focuses on the qualities of discrete schemes for propagation of acoustic perturbations in flows with very low Mach number. For subsonic flows, the expressions (37)- [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF] for face values are certainly justified. For high Mach number flows, there may be an advantage in taking supersonic flow conditions into account. But, we will illustrate by examples that the expressions in the form (37)-( 39), together with simple upwinding, function perfectly for high Mach number flows.

In the current study, we call a method with face velocity and face pressure by the expressions ( 37) and [START_REF] Xie | On numerical instabilities of Godunov-type schemes for strong shocks[END_REF], or [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF] and [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF], with the coefficients by the expressions [START_REF] Chen | Affordable shock-stable item for Godunovtype schemes against carbuncle phenomenon[END_REF] and [START_REF] Simon | A simple cure for numerical shock instability in the HLLC Riemann solver[END_REF], or simplifications thereof, as a method using characteristic equations, and we denote it by CHAR-CPS-ZB. As an example, we take the recent method by Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF], although they call their method as being of HLLC-type (see section 12 on the HLL-method). The method uses the expressions ( 37) and (39) as basis. But adaptations are necessary for low Mach number flows and for high Mach number flows on structured grids with strong shocks aligned with cell faces of one type and other cell faces aligned with the flow. The term shock stability is often used for referring to the problem with high Mach number flows. We employ this terminology hereafter. We use the method of Kitamura and Shima as an example for the discussion on the necessary adaptations in the next sections. Afterwards, we apply the same principles to the other methods.

8 Adaptation of the characteristic expressions for the hydrodynamic low Mach number limit

As shown in section 2, the pressure should scale with the velocity squared in the hydrodynamic low Mach number limit. But in the expressions [START_REF] Xie | On numerical instabilities of Godunov-type schemes for strong shocks[END_REF] and ( 39), the pressure is linearly dependent on the velocity difference. In order to obtain the correct scaling, the coefficient of the velocity difference should be made proportional to the local velocity. Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF] realise this scaling by preconditioning of the speed of sound.

With K either L or R, they define

cK = f c ( MK )c K , MK = min 1, M 2 K + M 2 co and f c ( MK ) = MK (2-MK ), ( 40
)
where M co is a cut-off Mach number, and M K = V K /c K is a local Mach number defined with the modulus V K of the multi-dimensional velocity (in multiple dimensions). With

SR = max {v R + cR , v M + cM } , SL = min {v L -cL , v M -cM } , αR = R ( SR -v R ) > 0, αL = L (v L -SL ) > 0, (41)
the expressions used by Kitamura and Shima are

v * = αL v L + αR v R αL + αR - p R -p L α L + α R , (42)
p * = 1 2 (p L + p R ) - 1 2 αR (v R -v * ) + 1 2 αL (v L -v * ). (43)
These expressions are weighted averages with an added difference term which acts as a diffusor and dissipator.

The rescaling of the α L and α R coefficients in the pressure expression (43) is essential.

The velocity difference term added to the face pressure causes an artificial viscosity term in the momentum equation (or momentum equations in multiple dimensions), which has a large coefficient in low Mach number flows without the rescaling. With the rescaling, it is strongly reduced. The multi-dimensional Mach number ensures that the artificial viscosity stays high in the transversal direction of a high Mach number flow. This is meant to improve the shock stability. But, for ensuring shock stability, it looks more elegant to us to add specific terms with this purpose. Moreover, Kitamura and Shima do not obtain guaranteed shock-stability. We opt for the recently proposed shock-stability terms by Chen et al. for use with an AUSM [START_REF] Chen | An improved AUSM-family scheme with robustness and accuracy for all Mach number flows[END_REF]. We will also use this method as the example of an AUSM. The shock-stability terms by Chen et al. can be used with any algorithm that is prone to shock instability and are switched off in smooth low Mach number flows. So, for the analysis of propagation of smooth acoustic perturbations in flows with very low Mach number, we do not need the shock-stability terms, which is an important advantage.

With the role on shock stability removed from the face pressure expression (43), the multi-dimensionality of the Mach number definition is not needed anymore. Therefore, we change the Mach number M K in (40) into the one-dimensional Mach number on the face, M K = v K /c K . Further, we want to avoid a cut-off Mach number, but the expression of the face velocity (42) may then become singular. However, by using preconditioned values of the sound speed, the first part of the right-hand side of Eq. ( 42) is only marginally changed. Therefore, we take the unaltered expression [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF], which is, to our understanding, well suited as well, and we set the cut-off Mach number M co = 0. So, we consider the changes that we make to the method of Kitamura and Shima as minor ones.

The pressure difference term in the velocity expression [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF] or ( 42) is necessary for avoiding a checkerboard type oscillation in the pressure field by creating a diffusion term in the continuity equation and the energy equation, which are the pressure-determining equations. This necessity, which is known as the necessity for pressure-velocity coupling, comes from the central type discretisation of the pressure flux and becomes quite critical in a flow in the hydrodynamic low Mach number limit.

Li and Gu [START_REF] Li | Mechanism of Roe-type schemes for all-speed flows and its application[END_REF] studied the literature on coupling-formulations with upwind methods.

There is very broad experience with the AUSM-type flux vector splitting and with the Roe-type flux-difference splitting. Based on the available experience, Li and Gu formulated two rules. The first concerns the pressure-velocity coupling. The face velocity should be connected to the pressure difference across the face assuring a pressure-velocity coupling such that a checkerboard oscillation of pressure is avoided at the limit of small velocity. The order of magnitude of the coefficient of ∆p is ideally between (1/c) and (1). A coefficient of order (1/c) is marginal and may lead to a weak form of checkerboard oscillation. A coefficient of order (1), or larger, avoids the checkerboard oscillation but may cause a rather large artificial viscosity in the numerical result. The pressure-velocity coupling does not have to be present in all equations. Li and Gu judge that it is the easiest to realise it in the continuity equation.

The second rule of Li and Gu concerns the face pressure expression. Usually, a diffusion term appears by a term proportional to the velocity difference across the face. Very often, this term does not have the correct scaling for the hydrodynamic low Mach number limit. The coefficient of the velocity difference is often of order (c) and it should be reduced to order (v) or lower. This deficiency can mostly be repaired by multiplication of the difference term by a suitable function of the local Mach number. We conclude that the expressions ( 37) and ( 43) satisfy the rules of Li and Gu.

In the tests of the method denoted by CHAR-CPS-ZB, we use the star values as face values and use simple upwinding, as in the method of Toro and Vázquez-Cendón [START_REF] Toro | Flux splitting schemes for the Euler equations[END_REF]. Moreover, also following Toro and Vázquez-Cendón, we take α L = L c L and α R = R c R , and thus αL = L cL and αR = R cR , in the expressions (37) and (43) for face velocity and face pressure. There is ambiguity with the choice of the upwinding for face velocity equal to zero, but it has no effect because the convection flux is then zero. The pressure part of the flux vector is expressed with v 1/2 and p 1/2 . Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF] take more accurate values of the face velocity, the face pressure and the upwinded vector of conserved variables. For this goal, they use expressions that are also used in a HLLCmethod (see section 12 on the HLL-method). But their method is not a pure HLL-method in the sense of Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF][START_REF] Toro | The HLLC Riemann solver[END_REF].

With the star values of v * and p * , the corresponding vectors of conserved variables are

U * K = * K , * K v * , 1 2 * K v 2 * + 1 γ -1 p * T , K = L, R. (44)
The densities * K are determined by the equations [START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF]. The face values of velocity and pressure may then be defined by

v 1/2 = v L and p 1/2 = p L if 0 < S L , v 1/2 = v * and p 1/2 = p * if S L ≤ 0 ≤ S R , v 1/2 = v R and p 1/2 = p R if S R < 0. ( 45
)
The vector of conserved variables at the face may be defined by

U 1/2 = U L if 0 < S L , U 1/2 = U * L if S L ≤ 0 ≤ v * , U 1/2 = U * R if v * ≤ 0 ≤ S R , U 1/2 = U R if S R < 0. (46)
Again, the convection flux is zero for v * = 0. We do not consider the more detailed definitions (44)-( 46) used by Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF] and we take the simple definitions by Toro and Vázquez-Cendón [START_REF] Toro | Flux splitting schemes for the Euler equations[END_REF]. We will demonstrate by test examples that the simple definitions satisfy the basic requirements for approximate Riemann solvers. The more detailed definitions may improve the accuracy for high Mach number flows, but in the essential tests in the current research on acoustic propagation, they are anyhow not activated.

9 Adaptation of the characteristic expressions for the acoustic low Mach number limit

For small pressure and velocity differences, the limit equations of ( 37) and (43) are

v * = 1 2 (v L + v R ) - 1 2 p R -p L c , (47)
p * = 1 2 (p L + p R ) - 1 2 f c ( M ) c(v R -v L ), (48)
where M , and c are values between these at the positions L and R. These equations apply to a flow that is nearly incompressible because all differences over a small distance become then small.

In the acoustic low Mach number limit, the view is with the speed of sound as the unity value and the flow velocity as very small. By (47) one sees that the velocity expressions [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF] and (42) then introduce a pressure diffusion term with a rather large coefficient in the continuity equation and the energy equation. So, one may expect possibility of rather strong damping in the discrete solution of propagating pressure waves in a low Mach number flow. We will illustrate later that this is the case. Thus, for use in the acoustic low Mach number limit, the coefficient of the pressure difference in a face velocity expression should be lowered. However, this is in conflict with the necessity for a sufficiently large coefficient in the hydrodynamic low Mach number limit. It means that it is not possible to satisfy both requirements by a manipulation which is only based on a Mach number. So, with most methods, adaptation to the acoustic low Mach number limit is not targeted. The face value expressions that are adapted for the hydrodynamic low Mach number limit can be used for the simulation of acoustic propagation, but Eq. ( 47) makes clear that low-diffusion and low-dissipation solutions can only be obtained with smooth pressure perturbations, a large number of grid points along the perturbation, and high accuracy of the space discretisation, as demonstrated by Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF] and also by Shima and Kitamura [START_REF] Shima | New approaches for computation of low Mach number flows[END_REF] with a method called SLAU (see section 11 on SLAU).

v AUSM + 1/2 = c 1/2 M 1/2 , with c 1/2 = min {c L , cR } , ( 49
)
where cK = (c † K ) 2 / max{c † K ,|v K |}, (c † K ) 2 = 2 γ-1 γ+1 H K , M K = v K /c 1/2
, with K for L or R and H the total enthalpy. The interpolation polynomials are

M 1/2 = f + M (M L ) + f - M (M R ), p 1/2 = f + p (M L )p L + f - p (M R )p R . (50)
The polynomials of AUSM + are listed by Liou [START_REF] Liou | A sequel to AUSM: AUSM +[END_REF]. Their main behaviour is

f + M = 0 for M ≤ -1, f + M (M 0) = 3 8 + 1 2 M + O(M 4 ), f + M = M for M ≥ 1, f - M = M for M ≤ -1, f - M (M 0) = - 3 8 + 1 2 M + O(M 4 ), f - M = 0 for M ≥ 1, f + p = 0 for M ≤ -1, f + p (M 0) = 1 2 + 15 16 M + O(M 3 ), f + p = 1 for M ≥ 1, f - p = 1 for M ≤ -1, f - p (M 0) = 1 2 - 15 16 M + O(M 3 ), f - p = 0 for M ≥ 1.
For low Mach number, the interpolated Mach number is approximately the average Mach number. The face velocity is thus also approximately the average of the left-hand and right-hand values at the face. There is no contribution by a pressure difference. Thus, for use at low Mach number, a pressure difference has to be added. Examples of such terms are those of the AUSM + -up by Liou [START_REF] Liou | A sequel to AUSM, part II: AUSM+ -up for all speeds[END_REF] or those of the SLAU of Shima and Kitamura [START_REF] Shima | Parameter-free simple low-dissipation AUSM-family scheme for all speeds[END_REF] and SLAU2 by Kitamura and Shima [START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[END_REF].

For low Mach number, the interpolated pressure is approximately the average pressure with a term proportional to the velocity difference v L -v R added, thus of form [START_REF] Xie | On numerical instabilities of Godunov-type schemes for strong shocks[END_REF]. Insight in the structure of the pressure interpolation formula is gained by writing the expression (50), with the short notations f p + L and f p -R , as

p 1/2 = 1 2 (p L + p R ) + 1 2 (f p + L -f p - R )(p L -p R ) + 1 2 (f p + L + f p - R -1)(p L + p R ). (51)
For low Mach number, this expression is approximately

p 1/2 = 1 2 (p L + p R ) - 15 16 v c 1/2 (p R -p L ) - 15 16
p c 1/2 (v R -v L ), (52)
where v and p are values in between those at states L and R. The second term in the right hand-side of ( 51) is a pressure-diffusion term with a small coefficient. The third term in the right hand-side of ( 51) is a velocity-diffusion term with a coefficient of order c 1/2 , which is not small. Therefore, this term has to be damped for use in a low Mach number flow. Examples of damping methods are by manipulation of the polynomial coefficients with the AUSM + -up by Liou [START_REF] Liou | A sequel to AUSM, part II: AUSM+ -up for all speeds[END_REF] or by multiplication with a function which becomes small for small Mach number in the SLAU of Shima and Kitamura [START_REF] Shima | Parameter-free simple low-dissipation AUSM-family scheme for all speeds[END_REF] and SLAU2 by Kitamura and Shima [START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[END_REF].

As an example, we take here the recent method by Chen et al. [START_REF] Chen | An improved AUSM-family scheme with robustness and accuracy for all Mach number flows[END_REF]. The Mach number interpolation is

M 1/2 = f + M (M L ) + f - M (M R ) - 1 2 (1 -f M )(1 -g) p R -p L 1/2 c 2 1/2 . ( 53
)
The corresponding face velocity is

v 1/2 = v AUSM + 1/2 - 1 2 (1 -f M )(1 -g) p R -p L 1/2 c 1/2 , ( 54
)
where

1/2 = 1 2 ( L + R ), f M = 1 2 1 -cos(π M ) , M = min{1, max{|M L | ,|M R |}}. ( 55
)
The factor (1 -f M ) approaches unity in a low Mach number flow and it becomes zero for a high value of the Mach number associated to the velocity component normal to the face. The function g is a shock-detector function, which approaches unity at the position of a strong shock and approaches zero in smooth flows. We detail it in section 14 on shock stability. Thus, by the factor (1 -f M )(1 -g) the pressure-diffusion term is active in smooth low Mach number flows and the coefficient of the pressure-diffusion term has the correct magnitude according to the first rule by Li and Gu. The pressure diffusion term is near to zero in a high Mach number flow, also on a face aligned with the flow due to the multi-dimensional function g. This strong reduction of the pressure difference term improves the shock-stability (see section 14).

The pressure interpolation in the method of Chen et al. is the AUSM-formula (51) with the third term damped and an added velocity diffusion term for shock stability (p ds ):

p 1/2 = 1 2 (p L +p R )+ 1 2 (f p + L -f p - R )(p L -p R )+ 1 2 (f p )(f p + L +f p - R -1)(p L +p R )+p ds , (56)
where f p = f M . There is actually a cut-off Mach number in the expression of f p by Chen et al., but we do not use it in our later tests. We detail the shock-stability term in section 14. But in our later tests, it is not needed and we set it to zero.

We denote the method with the expressions (54)-( 56) by AUSM-CPS-ZB.

11 Convection-pressure splitting with SLAU

In the SLAU-method (Simple Low-dissipation Advection Upstream method) by Shima and Kitamura [START_REF] Shima | Parameter-free simple low-dissipation AUSM-family scheme for all speeds[END_REF] and the SLAU2 by Kitamura and Shima [START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[END_REF], the expression of the mass flux is derived from the Roe flux-difference splitting method:

ṁ1/2 = 1 2 ( L v L + R v R ) - 1 2 L |v L | + R |v R | L + R ( R -L ) -χ 1 2 p R -p L c , (57)
with c = 1 2 (c L + c R ).
There is actually a factor multiplying the first two terms in the right-hand side, but it is not written here because it is only different from unity for a very strong supersonic expansion fan and we do not have such phenomenon in the tests that we present later.

In one dimensional flow, the χ-factor is

χ = (1 -M1/2 ) 2 , with M1/2 = min    1, 1 c 1 2 (v 2 L + v 2 R )    . ( 58
)
For use with ZB-splitting, a face velocity is needed, which we may derive from (57) as

v 1/2 = L v L + R v R L + R - L |v L | + R |v R | L + R R -L R + L -χ p R -p L ( L + R )c 1/2 , ( 59
)
where χ = (1 -M1/2 ) 2 , with M1/2 = min 1, M 1/2 , (60)
and where we take the AUSM + expressions (49) and (50) of c 1/2 and M 1/2 .

In a low Mach number flow, the χ-factor (58) is near to unity. The pressure difference term in the face velocity expression (59) is thus similar to that of the AUSM-method (54).

The face pressure expression is (56):

p 1/2 = 1 2 (p L + p R ) + 1 2 (f p + L -f p - R )(p L -p R ) + 1 2 (f p )(f p + L + f p - R -1)(p L + p R ), (61) with f p (M 1/2 ) = 1-(1-M1/2 ) 2 = M1/2 (2-M1/2 ), and
M1/2 = min 1, M 1/2 . ( 62
)
The shock-stability term in Eq. ( 56) may be added for high Mach number flows. The face pressure expression of SLAU2 is slightly different:

p 1/2 = 1 2 (p L + p R ) + 1 2 (f p + L -f p - R )(p L -p R ) + 1 2 (v 2 R + v 2 L )(f p + L + f p - R -1) 1/2 c 1/2 .
We denote the method defined by the expressions (59)-( 62) by SLAU-CPS-ZB. The essential difference with the AUSM-method of the previous section is the expression of the face velocity (59).

12 Convection-pressure splitting with the HLL-method

In the basic HLL-method, a contact discontinuity is ignored. With a unique state in the star-region, U * = U * L = U * R , the Rankine-Hugoniot equations ( 26) and ( 27) become

F * -F L = S L (U * -U L ) and F * -F R = S R (U * -U R ). (63)
These equations determine uniquely the vectors U * and F * . The solution is

U * = S R U R -S L U L S R -S L - F R -F L S R -S L and F * = S R F L -S L F R S R -S L + S R S L (U R -U L ) S R -S L . (64)
In a HLL-method, the face flux vector is

F 1/2 = F L if 0 < S L ; F 1/2 = F * if S L ≤ 0 ≤ S R ; F 1/2 = F R if S R < 0.
The star-region state vector U * is not used for the discretisation. Moreover,

F * = F (U * ).
By ignoring the contact discontinuity, a HLL-algorithm is dissipative for a contact discontinuity. In multiple dimensions, there is also possible discontinuity of the velocity components in the tangential direction in the flux vector on a face, called shear discontinuities. These discontinuities are transported convectively. The basic HLL-method is also dissipative for the shear discontinuities.

The contact discontinuity is taken into account in the HLLC-method, and in multiple dimensions, also discontinuities of the tangential velocity components. With the state vectors (44), the flux vectors in the star region are constructed by the Rankine-Hugoniot equations ( 26) and ( 27):

F * L = F L + S L (U * L -U L ) and F * R = F R + S R (U * R -U R ).
This way, a consistent approximation of the Rankine-Hugoniot equations across the wave paths is assured. The choice of the face flux vector is by

F 1/2 = F L if 0 < S L ; F 1/2 = F * L if S L ≤ 0 < v * ; F 1/2 = F * R if v * < 0 ≤ S R ; F 1/2 = F R if S R < 0.
There is ambiguity for v * = 0. The average of F * L and F * R may then be taken.

With a convection-pressure split method, the HLL-procedure is only applied to the pressure flux part. We use the formulation proposed by Mandal and Panwar [START_REF] Mandal | Robust HLL-type Riemann solver capable of resolving contact discontinuity[END_REF], with a small change of the transporting velocity in the convection part. The convection part of the flux is defined by flux vector upwinding:

F C,1/2 =        L ( v) L 1 2 ( v 2 ) L + 1 γ-1 p L        v 1/2 for v 1/2 ≥ 0. ( 65
)
Instead of the expression of Mandal and Panwar, we take the AUSM + face velocity [START_REF] Liou | A sequel to AUSM: AUSM +[END_REF] as transporting velocity. The replacement of the face velocity is a detail, but we feel safer by a face velocity that is a smooth function of the Mach number at the face. The pressure part of the flux is defined by HLL, but with replaced with p/c 2 in the dissipation term proportional to U R -U L :

F p,1/2 = S R F p,L -S L F p,R S R -S L + S R S L S R -S L         1 c 2 (p R -p L ) 1 c 2 [(pv) R -(pv) L ] 1 2c 2 [(pv 2 ) R -(pv 2 ) L ] + 1 γ -1 (p R -p L )         , ( 66
) with c = (c L + c R )/2, S R = max{0, v R + c R , ṽ + c} and S L = min{0, v L -c L , ṽ -c},
where ṽ and c are Roe-averages.

We denote the method with the expressions ( 65) and (66) by HLLP-CPS-ZB. The letter P is added to refer to the replacement of density by pressure in the diffusive part of the pressure flux vector expression (66). This replacement is done for obtaining exact representation of a contact discontinuity. This way, the diffusion of the contact discontinuity of the original HLL-method is removed. In multiple dimensions, the diffusion of the shear discontinuities is maintained. This diffusion has a role in the shock-stability of the method, as we discuss in section 14.

Mandal and Panwar [START_REF] Mandal | Robust HLL-type Riemann solver capable of resolving contact discontinuity[END_REF] demonstrate good functioning of HLLP-CPS-ZB for 1-D Riemann problems and 2-D high Mach number problems. In the later tests, we use the method as described above, with S R = v R +c R and S L = v L -c L . We remark that the interpretation of Toro and Vázquez-Cendón [START_REF] Toro | Flux splitting schemes for the Euler equations[END_REF] is that the pressure wave propagation in a convectionpressure splitting method is superposed to the convection motion, and thus the choice of the propagation speeds may also be S R = c R , S L = -c L and then c 2 = c L c R . With the last choices, the expressions (66) become similar to these of a characteristic method ( 37)- [START_REF] Xie | On numerical instabilities of Godunov-type schemes for strong shocks[END_REF]. Thus, also these last choices function well.

By definition of four factors,

T R = S R S R -S L , T L = S L S R -S L , T M = S R S L S R -S L , T A = T M c 2 , F p,1/2 has as components        0 + T A (p R -p L ) T R p L -T L p R + T A [(pv) R -(pv) L ] T R (pv) L -T L (pv) R + T A 2 [(pv 2 ) R -(pv 2 ) L ] + T M γ-1 (p R -p L )        . ( 67
)
With the HLLP method, the stepping method is still as described in section 4, but with the residuals defined by

k i -n i + τ k i + 1 2 ψ i ( k ) v k i+1/2 -τ k i-1 + 1 2 ψ i-1 ( k ) v k i-1/2 + τ (T A ) k i+1/2 p k i+1 - 1 2 ψ i+1 (p k ) -p k i + 1 2 ψ i (p k ) -τ (T A ) k i-1/2 p k i - 1 2 ψ i (p k ) -p k i-1 + 1 2 ψ i-1 (p k ) , ( v) k i -( v) n i +τ ( v) k i + 1 2 ψ i [( v) k ] v k i+1/2 -τ ( v) k i-1 + 1 2 ψ i-1 [( v) k ] v k i-1/2 + τ (T R ) k i+1/2 p k i + 1 2 ψ i (p k ) -τ (T L ) k i+1/2 p k i+1 - 1 2 ψ i+1 (p k ) -τ (T R ) k i-1/2 p k i-1 + 1 2 ψ i-1 (p k ) + τ (T L ) k i-1/2 p k i - 1 2 ψ i (p k ) + τ (T A ) k i+1/2 (pv) k i+1 - 1 2 ψ i+1 [(pv) k ] -(pv) k i + 1 2 ψ i [(pv) k ] -τ (T A ) k i-1/2 (pv) k i - 1 2 ψ i [(pv) k ] -(pv) k i-1 + 1 2 ψ i-1 [(pv) k ] , ( E) i -( E) n i + τ 1 2 [( v) i ] 2 i + 1 2 ψ i 1 2 [( v) k ] 2 k v k i+1/2 -τ 1 2 [( v) i-1 ] 2 i-1 + 1 2 ψ i-1 1 2 [( v) k ] 2 k v k i-1/2 + τ ( e) k i + 1 2 ψ i ( e) k v k i+1/2 -τ ( e) k i-1 + 1 2 ψ i-1 ( e) k v k i-1/2 + τ (T R ) k i+1/2 (pv) k i + 1 2 ψ i (pv) k -τ (T L ) k i+1/2 (pv) k i+1 - 1 2 ψ i+1 (pv) k -τ (T R ) k i-1/2 (pv) k i-1 + 1 2 ψ i-1 (pv) k + τ (T L ) k i-1/2 (pv) k i - 1 2 ψ i (pv) k + τ (T A ) k i+1/2 2 (pv 2 ) k i+1 - 1 2 ψ i+1 (pv 2 ) k -(pv 2 ) k i + 1 2 ψ i (pv 2 ) k -τ (T A ) k i-1/2 2 (pv 2 ) k i - 1 2 ψ i (pv 2 ) k -(pv 2 ) k i-1 + 1 2 ψ i-1 (pv 2 ) k + τ (T M ) k i+1/2 γ -1 p k i+1 - 1 2 ψ i+1 (p k ) -p k i + 1 2 ψ i (p k ) -τ (T M ) k i-1/2 γ -1 p k i - 1 2 ψ i (p k ) -p k i-1 + 1 2 ψ i-1 (p k ) .
In the limit for small velocity, the effective mass flux in the continuity equation becomes v -∆p 2c . The dissipation term realises a weak pressure-velocity coupling. In principle, it is sufficient. The continuity equation term in the pressure flux may thus stay

T A (p R -p L ). ( 68
)
In the momentum equation, the momentum flux becomes vv + p -1 2c (p∆v + v∆p). The coupling of the pressure difference to the velocity is as in the continuity equation. The effective face pressure is p -p 2c ∆v. With p ∼ c 2 , this may also be written as p -1 2 c∆v. The coefficient of ∆v does not have the correct order of magnitude in the hydrodynamic low Mach number limit. This may be cured by p -f p (M ) p 2c ∆v. So, an adapted momentum equation term in the pressure flux is

T R p L -T L p R + T A v(p R -p L ) + f p (M 1/2 )T A p(v R -v L ).
Li and Gu [START_REF] Li | Mechanism of Roe-type schemes for all-speed flows and its application[END_REF] actually suggest (with the substitution of by p/c 2 , which is not done by Li and Gu):

T R p L -T L p R + T A f p (M 1/2 ) [(pv) R -(pv) L ] , (69)
because they do not consider splitting of the (pv) term. Our numerical tests show that both modifications function equally well. This means that the coupling of the pressure difference to the face velocity is not essential in the momentum equation. Therefore, we take the expression (69) by Li and Gu.

In the energy equation, the flux becomes

1 2 v 2 v + 1 γ -1 pv + pv - 1 2c 1 2 v 2 ∆p + p v∆v - 1 γ -1 c 2 ∆p.
The combination of the first term and the fourth term shows a mass flux similar to that of the continuity equation. The combination of the second and third term with the fifth and sixth term shows a face pressure of the form p -p 2c ∆v and a face velocity of the form v -c 2p ∆p. The face pressure does not have the correct scaling for pressure-velocity coupling. Thus, the fifth term has to be multiplied by f p (M ). The face velocity has the correct form because it may also be written as v -1 2 ∆p c . Again, if we do not consider splitting of the term ∆(pv 2 ), the adapted expression of the energy term in the pressure flux is

T R (pv) L -T L (pv) R + T A 2 f p (M 1/2 ) (pv 2 ) R -(pv 2 ) L + T M γ -1 (p R -p L ). (70)
Li and Gu [START_REF] Li | Mechanism of Roe-type schemes for all-speed flows and its application[END_REF] do not propose any change to the energy difference term ∆( E) in the diffusive part of the energy equation, but they do not use convection-pressure splitting and they do not consider splitting of differences in the diffusive parts in several components. The low Mach number adaptation by Eq. ( 70) in the energy equation is done here with the HLLP-CPS-ZB in order to be consistent with the low Mach number adaptation in the momentum equation (69).

In the later numerical tests, the HLLP-CPS-ZB is defined by the convection part (65) with v 1/2 according to the AUSM + definition, and the pressure flux part by the components (68), ( 69) and (70) with f p (M 1/2 ) and M 1/2 according to Eq. (62).

Li and Gu [START_REF] Li | Mechanism of Roe-type schemes for all-speed flows and its application[END_REF] demonstrate good functioning of the HLL, with their adaptation for low Mach number, for a Riemann problem and for steady state 2-D low Mach number flows. Sun et al. [START_REF] Sun | A robust flux splitting method with low dissipation for all-speed flows[END_REF] demonstrate good functioning of the HLLP-CPS-TV with the low Mach number adaptation of Li and Gu for Riemann problems, steady state 2-D low Mach number flows and steady state 2-D high Mach number flows.

13 Convection-pressure splitting with momentum interpolation

The characteristic expressions (37)-( 38) are the solutions of the momentum equations [START_REF] X.-D. Ren | Role of the momentum interpolation mechanism of the Roe scheme in shock instability[END_REF]. This means that one of the equations ( 37)-( 38) may be replaced by a discretised momentum equation. Because there is a problem with obtaining the correct scaling of the expression of the face velocity [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF] in the same way for the hydrodynamic and the acoustic low Mach number limits, a practical approach is to maintain the expression of the face pressure in the form [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF], or a variant thereof, and to replace the expression of the face velocity by a discretised time-dependent momentum equation. By discretisation of the momentum equation at a face, a relation is then obtained between the face velocity and the pressure difference across the face. This may then serve as expression of the face velocity. If the discretisation of the momentum equation is done in the same way as for the quantities at the cell centres, the discrete equation has then automatically the correct scaling for all flow types. This requires thus a time-accurate discretisation of the momentum equation at a face.

In the acoustic low Mach number limit, the Euler equations ( 9) are

∂ t v + ∂ x p = 0 and ∂ t p + γp(∂ x v) = 0. (71)
Thus, for time-accuracy in the acoustic low Mach number limit, the time-dependent term and the pressure gradient term in the momentum equation have to discretised locally at the face. The result is then that the face velocity components become dynamic quantities that evolve in time together with the quantities in the cell centres. There are some methods of this type, including our method [START_REF] Moguen | Pressure-velocity coupling allowing acoustic calculation in low Mach number flow[END_REF][START_REF] Moguen | A combined momentum-interpolation and advection upstream splitting pressure-correction algorithm for simulation of convective and acoustic transport at all levels of Mach number[END_REF]. Other examples are by Ong and Chan [START_REF] Ong | A pressure-based Mach-uniform method for viscous fluid flows[END_REF], and by Xiao et al., Denner et al. and Bartholomew et al. [START_REF] Xiao | Fully-coupled pressure-based finitevolume framework for the simulation of fluid flows at all speeds in complex geometries[END_REF][START_REF] Denner | Pressure-based algorithm for compressible interfacial flows with acoustically-conservative interface discretisation[END_REF][START_REF] Bartholomew | Unified formulation of the momentum-weighted interpolation for collocated variable arrangements[END_REF]. A method by Li and Gu [START_REF] Li | The momentum interpolation method based on the time-marching algorithm for all-speed flows[END_REF] is similar, but the face velocity at the old time level is interpolated with as consequence that the face velocity is not completely time accurate.

We use here our own methodology. The face pressure is

p 1/2 = 1 2 (p L + p R ) - 1 2 f p (M 1/2 ) 1/2 c 1/2 (v R -v L ), (72)
where M 1/2 and c 1/2 are the Mach number and the velocity of sound of the AUSM + , by Eqs. ( 49) and (50), and the scaling factor f p (M 1/2 ) is by Eq. ( 62). The face pressure expression is simpler than the one of the MIAU scheme presented in our previous work [START_REF] Moguen | A combined momentum-interpolation and advection upstream splitting pressure-correction algorithm for simulation of convective and acoustic transport at all levels of Mach number[END_REF].

For the construction of the face velocity, the momentum equation is discretised in a preliminary way by using the face velocity of the AUSM + scheme into

1 ∆t [( v) i -( v) n i ] + A k i ( v) i -B k i + 1 ∆x (p k i+1/2 -p k i-1/2 ) = 0, with A k i = 1 ∆x v AUSM + k i+1/2 , B k i = ( v) k i-1 + 1 2 ψ i-1 ( v) k 1 ∆x v AUSM + k i-1/2 - 1 2 ψ i ( v) k 1 ∆x v AUSM + k i+1/2 . (73)
The face velocity v k i+1/2 is then obtained by interpolation:

B k i+1/2 = A k i+1/2 ( v) i+1/2 + 1 ∆x (p k i+1 -p k i ) + 1 ∆t [( v) i+1/2 -( v) n i+1/2 ], (74)
with

A k i+1/2 = 1 2 A k i + A k i+1 , B k i+1/2 = 1 2 B k i + B k i+1 , v k i+1/2 = ( v) i+1/2 1 2 ( k L + k R ) i+1/2 .
The essence of the time-accurate interpolation (74) is that the pressure gradient term and the inertia term are not interpolated, but discretised directly on the face. We remark that for use of the preliminary discretisation (73) in the interpolated formula (74), expressions of the face pressure in the preliminary discretisation are not necessary. The method by momentum interpolation is denoted by TAMI-CPS-ZB, where the acronym TAMI means time-accurate momentum interpolation.

Shock stability

Low-diffusion upwind schemes are susceptible to anomalies in the flow downstream of a strong shock aligned with cell faces of one type on structured grids and other cell faces aligned with the flow. Common deficiencies are odd-even transversal oscillations downstream of a moving shock and the carbuncle perturbation of the shock in a steady hypersonic blunt body flow. These phenomena, called shock anomalies or shock instabilities, are discussed an analysed by many authors.

The causes of the anomalies are not yet fully understood, but recent research on the topic results in certitude about two aspects. First, the velocity difference term in the face pressure expression and the pressure difference term in the face velocity expression drive towards instability. The second aspect is the necessity for sufficient damping of coupled pressure and velocity perturbations in the transversal direction.

Recently, Fleischmann et al. [START_REF] Fleischmann | A low dissipation method to cure the grid-aligned shock[END_REF] obtained shock stability of the Roe-method by decreased velocity difference in the effective face pressure by preconditioning the velocity of sound. This is thus similar to the reduction of the velocity difference terms in the expression of the face pressure (43). A pure Roe-method is shock-unstable, but by reduction of the velocity difference, the method may thus become shock-stable. The proposed reduction by Fleischmann et al. is similar to that by Li and Gu [START_REF] Li | Mechanism of Roe-type schemes for all-speed flows and its application[END_REF][START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF] and that by Rieper [START_REF] Rieper | A low-Mach number fix for Roe's approximate Riemann solver[END_REF], but their motivation was not shock stability but high accuracy in the hydrodynamic low Mach number limit. Obviously, the remedy is the same for both aspects. The explanation is that a small velocity difference in the normal direction of a cell face, approximately aligned with the flow, can cause a large change of the cell centre pressure by a large coefficient of the velocity difference in a face pressure expression of form [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF]. The perturbation of the pressure may then cause a large change of the velocity component perpendicular to the face by a face velocity expression of form [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF]. The two coupling mechanisms create the possibility for a coupled pressure and velocity oscillation in a direction transversal to the flow. By reduction of one of the coupling coefficients, the tendency towards oscillations is reduced.

The explanation thus means that shock instability may also be cured by lowering the coefficient of the pressure difference term in a face velocity expression of form [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF]. It was shown by Ren et al. [START_REF] X.-D. Ren | Role of the momentum interpolation mechanism of the Roe scheme in shock instability[END_REF] that this is correct. Also the SLAU-method by Shima and Kitamura [START_REF] Shima | Parameter-free simple low-dissipation AUSM-family scheme for all speeds[END_REF], the SLAU2 by Kitamura and Shima [START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[END_REF] become shock-stable by lowering the coefficient of the pressure difference term in the face velocity expression.

Reduction of the coupling between the face velocity and the pressure difference across the face or reduction of the coupling between the face pressure and the velocity difference across the face is, however, not a guarantee for obtaining shock stability. E.g., the method of Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF] based on characteristic equations, which we take as an example, adapted for low Mach number flows by reduction of the velocity difference term in the face pressure expression, is not shock-stable for flows at very high Mach number. That lowering the coupling between the face velocity and the pressure difference does not guarantee shock-stability was demonstrated by Li et al. [START_REF] Li | Shock-stable Roe scheme combining entropy fix and rotated Riemann solver[END_REF] for the Roe-method. These observations demonstrate that also sufficient damping of the coupled velocity and pressure perturbations is necessary. Some additional damping may thus be necessary, which has to be tailored for minimisation of dissipation [START_REF] Li | Shock-stable Roe scheme combining entropy fix and rotated Riemann solver[END_REF]. A further aspect is that lowering the coupling between face velocity and pressure difference improves the shock-stability, but increases the tendency towards expansion shocks, as demonstrated by Li et al. [START_REF] Li | Cures for expansion shock and shock instability of Roe scheme based on momentum interpolation mechanism[END_REF] for the Roe-method. Thus, the way of reducing the coupling has to take this possible effect into account [START_REF] Li | Cures for expansion shock and shock instability of Roe scheme based on momentum interpolation mechanism[END_REF].

Insight into the necessary damping has been obtained by comparison of the shock-stable HLL-method and the shock-unstable HLLC-method. The difference between the flux expressions of both methods shows terms that diffuse the contact discontinuity and the shear discontinuities. The shear diffusion was found to be crucial. A shock unstable HLLCmethod may thus be stabilised by adding shear diffusion. Examples of such methods are those of Shen et al. [START_REF] Shen | A robust HLLC-type Riemann solver for strong shock[END_REF], Chen et al. [START_REF] Chen | Affordable shock-stable item for Godunovtype schemes against carbuncle phenomenon[END_REF], and Simon and Mandal [START_REF] Simon | A simple cure for numerical shock instability in the HLLC Riemann solver[END_REF]. The HLLEM may be stabilised in a similar way. The HLLEM is an accurate method by adding anti-diffusion terms for the contact discontinuity and the shear discontinuities to the basic HLL. By adding the full terms, the HLLEM is shock-unstable, but by limiting the anti-diffusion terms for shear in the vicinity of strong shocks, a shock-stable method is obtained. Examples are the methods by Xie et al. [START_REF] Xie | A low diffusion flux splitting method for inviscid compressible flows[END_REF][START_REF] Xie | On numerical instabilities of Godunov-type schemes for strong shocks[END_REF]. The method by Mandal and Panwar [START_REF] Mandal | Robust HLL-type Riemann solver capable of resolving contact discontinuity[END_REF], described in section 12, is shock-stable, because the replacement of density by pressure in the diffusive part of the pressure flux vector removes the diffusion of a contact discontinuity, but keeps the diffusion of shear discontinuities. A similar method is by Sun et al. [START_REF] Sun | A robust flux splitting method with low dissipation for all-speed flows[END_REF].

We detail here the formulation by Chen et al. [START_REF] Chen | Affordable shock-stable item for Godunovtype schemes against carbuncle phenomenon[END_REF], which they use to stabilise the AUSM [START_REF] Chen | An improved AUSM-family scheme with robustness and accuracy for all Mach number flows[END_REF]. The terms for addition to the momentum fluxes through a face perpendicular to the x-direction in 2-D are

p ds x = -γ 1/2 c 1/2 (f p + L f p - R )(g)(v xR -v xL ), p ds y = -γ 1/2 c 1/2 (f p + L f p - R )(g)(v y R -v y L ),
where v x and v y denote the xand y-components of the velocity, and p ds x and p ds y are added to the momentum equation in x-direction and y-direction, respectively. There is thus diffusion added in normal direction and in tangential direction on all faces. The diffusion in tangential direction is essential. The one in normal direction is added for robustness. The factor (f p

+ L f p - R )
ensures that the additions are only done for a subsonic velocity component perpendicular to the face and the g-function ensures that the terms 29 are only added in the vicinity of a strong shock. The g-function is

g = 1 2 [1 + cos(πh)], with h = min k {h k }, h k = min p L k p R k , p R k p L k ,
where k indicates all the faces of the cells adjacent to the face where the pressure is evaluated. In vicinity of shocks, the g-function approaches unity, but for smooth flow, it is very small. Further, by the factor (1 -f M )(1 -g) in the face velocity expression (54), the pressure difference term is set to a small value in vicinity of a strong shock, which reduces the mechanism of shock instability.

15 Numerical results for steady flows

In order to demonstrate that the proposed methods all function well for steady flows at medium-high and high Mach numbers, we show results of one-dimensional flows with a stationary contact discontinuity or a stationary shock.

Stationary contact discontinuity

The initial data are listed in Table 1. We take N = 200 cells and CFL(c L ) = 500.

L (kg/m 3 ) v L (m/s) p L (Pa) R (kg/m 3 ) v R (m/s) p R (Pa) t f (s)
1.0 0 100 000 0.5 0 100 000 10

Table 1 Settings for the stationary contact discontinuity test of Sec. 15.1. The result of CHAR-CPS-ZB is shown in Fig. 3. The contact discontinuity is exactly represented. The results of the other algorithms are identical (not shown).

Stationary shock

The initial data are listed in Table 2. We take N = 1 000 cells and CFL(v L ) = 0.5. The result of HLLP-CPS-ZB is shown in Fig. 4. The results are similar by the other algorithms (not shown). 16 Numerical results for Riemann problems at medium-high Mach number

L (kg/m 3 ) v L (m/s) p L (Pa) R (kg/m 3 ) v R (m/s) p R (Pa) t f (s) 1
We study two cases, chosen from the list of test cases by Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. The first is the modified Sod test, with a right travelling shock wave, a right travelling contact discontinuity and a left travelling expansion wave. The second is a similar case, but with a very strong right travelling shock wave, and a strong right travelling contact discontinuity.

Modified Sod test

The initial data are listed in Table 3. We take N = 1 000 cells and CFL(v L ) = 0.5.

The results of AUSM-CPS-ZB and TAMI-CPS-ZB are shown in Fig. 5 and Fig. 6. In the AUSM-result, there is a slight oscillation in the flow downstream of the contact discontinuity. The reason is that the Mach number has different values upstream and downstream of the moving contact discontinuity. This leads to an interpolated pressure by Eq. (56) that is not exactly equal to the upstream and downstream values of the pressure. The results of the two other algorithms are similar to the result of TAMI-CPS-ZB (not shown).

L (kg/m 3 ) v L (m/s) p L (Pa) R (kg/m 3 ) v R (m/s) p R (Pa) t f (s) 1 0.75 1 0.125 0 0.1 0.2

Strong shock tube

The initial data are listed in Table 4. We take N = 1 000 cells and CFL(v L ) = 0.5.

The result of SLAU-CPS-ZB is shown in Fig. 7. The results of CHAR-CPS-ZB, AUSM-CPS-ZB and TAMI-CPS-ZB are similar (not shown). The HLLP-CPS-ZB does not converge with all dissipation terms set at second order accuracy in Eq. (66). A Solution somewhat more diffusive than the one shown in Fig. 7 was obtained by setting the pressure dis- sipation term in the continuity equation and in the energy equation to first order accuracy, thus: T A (p i+1 -p i ) and T M γ-1 (p i+1 -p i ).

L (kg/m 3 ) v L (m/s) p L (Pa) R (kg/m 3 ) v R (m/s) p R (Pa) t f (s) 1 0 1 000 1 0 0.01 1.2 × 10 -2

Numerical results for a Riemann problem at low Mach number

The test case is a Riemann problem at a very low Mach number and with very weak expansion waves. It is similar to a case suggested by Guillard and Murrone [START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes[END_REF], but modified such that the pressure changes and velocity changes stay small. We opt for small changes because these are the most challenging for an approximate Riemann solver.

The initial data are listed in Table 5. We take N = 1 000 cells and CFL(v L + c L ) = 0.5 and 10 -2 . with low CFL number, oscillations occur around the expansion waves. The pressure is obtained from the energy equation [START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[END_REF]. There are two mechanisms of pressure diffusion in this equation. A first is by pressure-velocity coupling. This diffusion has the form of the terms with coefficients (τ α) in the C-factors of the system matrix of the energy equation. With the CHAR-, AUSM-, SLAU-and HLLP-methods, the diffusion coefficients are proportional to τ , because the coupling coefficients of the pressure difference in the face velocity expressions (37), ( 54) and (66) are of the order of unity. Because the energy balance by Eq. ( 20) is multiplied with the time step, the coefficient of the pressure diffusion due to pressure-velocity coupling is of order unity, independent of the time step. With the TAMI-method, the coupling coefficient of the pressure difference term in the expression of the face velocity by Eq. ( 74) is of the order τ . The pressure diffusion due to pressure-velocity coupling is thus proportional to the time step. It becomes thus small for small time step. The second mechanism of pressure diffusion is by the upwinding of the convection term. This diffusion has the form of the terms with coefficients (τ v) in the C-factors of the system matrix of the energy equation. This diffusion becomes thus small for small velocity. The very small pressure diffusion by the TAMI-method for small time step and small velocity thus makes that it is not meaningful to use TAMI-CPS-ZB with a small CFL number in a flow with small convection velocity. The results by SLAU-CPS-ZB and HLLP-CPS-ZB are shown in Figs. 12 and13. The results by CHAR-CPS-ZB and AUSM-CPS-ZB are comparable (not shown). The dissipation of the pressure pulse by these methods is obvious. The result by TAMI-CPS-ZB is shown in Fig. 14. There is no visible dissipation; only a small dispersion. The explanation of the dissipation by the CHAR-, AUSM-, SLAU-and HLLP-methods is that the face velocity expressions (37), (54) and (66) are of the form

L (kg/m 3 ) v L (m/s) p L (Pa) R (kg/m 3 ) v R (m/s) p R (Pa) t f (s) 25
v i+1/2 = v - 1 2 p R -p L c , (75)
where the overbar values are forms of averages of left-hand side and right-hand side values at the face. The coupling of the face velocity with the pressure difference causes the dissipation. With TAMI-CPS-ZB, the face velocity, for small convection velocity, is obtained from the time-dependent momentum equation on the face, which, in space-discretised form, is

∂ t ( v) i+1/2 + p i+1 -p i ∆x = 0. ( 76
)
The only dissipation associated with this equation comes from the time discretisation and the space discretisation. The result shown in Fig. 14 thus proves the advantage by obtaining the face velocity in a time-accurate way by the discretisation of a momentum equation on the face.

A countereffect of the construction of the pressure-velocity coupling by the time-accurate discretisation of the momentum equation at a face (76) is, as discussed in the previous section, that the coupling disappears for small time step. It is thus not possible to use TAMI-CPS-ZB with a small CFL number based on the velocity of sound. That the pressure difference term in the algebraic expression of the face velocity by the CHAR-, AUSM-, SLAU-and HLLP-methods is the cause of the dissipation is demonstrated by the result shown in Fig. 15, obtained with the SLAU-method, by setting the χ-factor in the expression of the face velocity (59) to 10 -2 . The reduction of the dissipation is obvious. But, clearly, the quality of the solution is not yet that by the time-accurate momentum interpolation method. Moreover, it is not possible to obtain a better result than shown in Fig. 15, because oscillations appear in the solution with a lower value of the χfactor. Another aspect is that lowering the coefficient of the pressure difference term in the face velocity expression is not a practical method for problems of acoustic propagation in low Mach number flows, because this coefficient has to be large enough for application to in Fig. 19. The solution by the TAMI-method is very accurate but the solution of the CHAR-method shows some diffusion. The comparison of both results proves the benefit of using the time-accurate momentum interpolation method. The observation in Fig. 17 is that CHAR-CPS-ZB is able to produce an accurate result, but that this requires sufficient smoothness of the acoustic perturbation and a sufficiently fine space discretisation of the perturbation. This is in accordance with the statement by Kitamura and Shima [START_REF] Kitamura | AUSM-like expression of HLLC and its all-speed extension[END_REF] on the accuracy of CHAR-CPS-ZB for acoustic propagation. With the result shown in Fig. 19, it is clear that the requirements on smoothness of the acoustic perturbation and the space resolution are much less severe with TAMI-CPS-ZB. momentum equation at the face. The accuracy of the methods by algebraic face velocity definition (CHAR-, AUSM-, SLAU-, HLLP-methods) is lower, although the face velocity expression is derived by considering a time-dependent Riemann problem and the obtained expressions are time-accurate for Riemann problems. Only the method with time-accurate momentum interpolation (TAMI-method) for the definition of the face velocity proves to be genuinely accurate for acoustic propagation in low Mach number flows. An essential difference between the momentum interpolation method and the other ones is that the face velocity is a dynamic quantity which evolves in time in the same way as the cell-centre quantities. This feature is the key for obtaining accurate solutions of acoustic propagation in low Mach number flows.

  no difference in velocity and pressure: v * L = v * R = v * , and p * L = p * R = p * (contact discontinuity). The characteristic equations along the pressure wave paths are dp -c dv = 0 along d t x = v -c ; dp + c dv = 0 along d t x = v + c.

Fig. 1 .

 1 Fig. 1. Solution for star-values with characteristic equations.

Fig. 2

 2 Fig. 2 is a sketch of the characteristic paths in the x -t plane with the propagation speeds v -c, v and v + c, estimated by S L , S * and S R . The star-region values of the conserved variables at the left-hand and right-hand sides of the contact discontinuity are denoted by U * L and U * R . For derivation of relations across the paths, control volumes are chosen that connect the states on each side.

Fig. 2 .

 2 Fig. 2. Solution for star-values with Rankine-Hugoniot equations.

10

  Convection-pressure splitting with AUSMThe basic expressions of face velocity and face pressure of an AUSM (Advection Upstream Splitting Method) are by interpolation between left-hand and right-hand values on a face. The interpolation is done with polynomials formed by factors m + 1, m, m -1, where m is a Mach number. The use of these factors is inspired by the propagation speeds v + c, v and v -c of the characteristic quantities. For the definition of the face Mach number, first, a common speed of sound c 1/2 at the face 1/2 is defined. The expressions of the AUSM + by Liou[START_REF] Liou | A sequel to AUSM: AUSM +[END_REF] are

Fig. 3 .

 3 Fig. 3. Stationary contact discontinuity: CHAR-CPS-ZB.

Fig. 4 .

 4 Fig. 4. Stationary shock: HLLP-CPS-ZB.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Modified Sod test: AUSM-CPS-ZB.

Fig. 7 .

 7 Fig. 7. Strong shock tube: SLAU-CPS-ZB.

Fig. 11 .

 11 Fig. 11. Low Mach number Riemann problem: TAMI-CPS-ZB. CFL(v L + c L ) = 10 -2 .18 Numerical results for acoustic propagation in one dimension

Fig. 12 .

 12 Fig. 12. One-dimensional acoustic pulse: SLAU-CPS-ZB.

Fig. 13 .Fig. 14 .Fig. 15 .

 131415 Fig. 13. One-dimensional acoustic pulse: HLLP-CPS-ZB.

  t f = 1 ms.

Fig. 16 .Fig. 17 .

 1617 Fig.[START_REF] Shima | New approaches for computation of low Mach number flows[END_REF]. Two-dimensional acoustic pulse test. Evolution of the pressure distribution (Pa), by TAMI-CPS-ZB, at t f = 0, 0.33, 0.66 and 1 ms; σ = 0.05.

Table 2

 2 Settings for the stationary shock test of Sec. 15.2.

	.4	1 400	19 600	8.0	245	2 283 400 0.05

Table 3

 3 Settings for the modified Sod test of Sec. 16.1.

Table 4

 4 Settings for the strong shock tube test of Sec. 16.2.

Table 5

 5 Settings for the low Mach number Riemann problem of Sec. 17.The results of CHAR-CPS-ZB are shown in Figs.8 and 9for two values of the CFL number. The algorithm remains stable for the low value of the CFL number. The results are almost identical for AUSM-CPS-ZB, SLAU-CPS-ZB and HLLP-CPS-ZB (not shown).The results of TAMI-CPS-ZB with the same CFL numbers are shown in Figs.10 and 11. For CFL = 0.5, the result is almost identical to the result of the other methods, but Fig. 8. Low Mach number Riemann problem: CHAR-CPS-ZB, with CFL(v L + c L ) = 0.5.

	0.200	10 000.00	25	0.202	10 000.85 0.01

Fig. 9. Low Mach number Riemann problem: CHAR-CPS-ZB, with CFL(v L + c L ) = 10 -2 .

  Fig. 10. Low Mach number Riemann problem: TAMI-CPS-ZB. CFL(v L + c L ) = 0.5.

			Density (kg/m 3 )					Mach number
	25					0.00854
	24.9995					0.00852
	24.999					0.0085
	24.9985					0.00848
	24.998					0.00846
	24.9975					0.00844
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
			Density (kg/m 3 )					Mach number
	25					0.00854
	24.9995					0.00852
	24.999					0.0085
	24.9985					0.00848
	24.998					0.00846
	24.9975					0.00844
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1

low Mach number flows in the hydrodynamic low Mach number limit.

19 Two-dimensional acoustic pulse propagation in a low Mach number flow A two-dimensional Gaussian acoustic pulse is given at t = 0 by 0 = 0 + (δ ) 0 , u 0 = u 0 , v 0 = v 0 , p 0 = p 0 + (δp) 0 , where

.0886 × 10 -3 m/s, p 0 = 101 300 Pa. The exact solution of the initial value problem given by the Euler equations linearized around the uniform mean flow is (see Ref. [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF]):

sin(ξc 0 t) j 1 (ξη) ξ dξ, p(x , y , t) = c 2 0 (x , y , t), where

and j 0 and j 1 are the first kind cylindrical Bessel functions of orders zero and one, respectively. Here, M x = M y = v 0 /c 0 = 9×10 -6 . The computational domain is a square of side [0, 1] (m), divided into 500 × 500 cells forming a regular Cartesian grid. The numerical method is the 2-D direct extension of the one detailed in section 4, with an alternate direction procedure for solving the pentadiagonal system obtained from the energy equation.

The pressure distributions obtained with TAMI-CPS-ZB at four time instants are shown in Fig. 16 for σ = 0.05 and CFL(v 0 + c 0 ) = 10. The velocity distributions at t f = 1 ms on the line y = 0.5 obtained by TAMI-CPS-ZB and CHAR-CPS-ZB, and the analytical solution, are shown in Fig. 17. The accuracy of both methods is good. In Fig. 18 is shown similarly the pressure distribution by TAMI-CPS-ZB at t f = 1 ms for σ = 0.02, thus with a narrower initial pulse. The corresponding velocity distributions at t f = 1 ms on the line y = 0.5 by TAMI-CPS-ZB and CHAR-CPS-ZB, and the analytical solution, are shown