On the Behaviour of Upwind Schemes in the Low Mach Number Limit

Chapter 8

On the Behaviour of Upwind Schemes in the Low Mach Number Limit: A Review H. Guillard and B. Nkonga Universit
e C^ ote d’Azur, Inria, CNRS, LJAD, France

Chapter Outline 1 Introduction 2 The Multiple Low Mach Number Limits of the Compressible Euler Equations 2.1 Incompressible Limit 2.2 Acoustic Limit 2.3 Acoustic–Incompressible Interactions 2.4 Finite Volume Schemes 2.5 The Diagnosis 2.6 The Remedies

204

205 205 207 208 213 215 217

3 Numerical Illustrations 3.1 Order 1: Quadrangular Cartesian Grids 3.2 Order 1: Vertex-Centred Triangular Meshes 3.3 Order 1: Cell-Centred Triangular Meshes 3.4 Order 2: Vertex-Centred Triangular Meshes 4 Conclusion References

221 221 223 224 226 228 228

ABSTRACT This work is devoted to a review of different modifications proposed to enable compressible flow solvers to compute accurately flows near the incompressible limit. First the reasons of the failure of upwind solvers to obtain accurate solutions in the low Mach number regime are explained. Then different correction methods proposed in the literature are reviewed and discussed. This work concludes by some numerical experiments to illustrate the behaviour of the different methods. Keywords: Low Mach number flows, Upwind schemes, Finite volume methods MSC2010 Classification Codes: 76B03, 76G25, 76M12, 65M08

Handbook of Numerical Analysis, Vol. 18. http://dx.doi.org/10.1016/bs.hna.2016.09.002 © 2017 Elsevier B.V. All rights reserved.

203

204 Handbook of Numerical Analysis

1 INTRODUCTION Computational fluid dynamics (CFD) is now-a-days a mature discipline with a wide range of applications ranging from traditional hydro- or aerodynamics to biological flows, environmental studies or even pedestrian flows (Helbing, 1992) or social network modelling. CFD thus has to deal with a huge set of different types of flows that can be characterized by nondimensionless quantities like the Froude, Ekmann, Strouhal, Reynolds, etc., numbers expressing the relative importance of the forces acting on a fluid element. In this collection of dimensionless numbers, one of these, the Mach number plays a particular role from the point of view of numerical methods. For low Mach number flows close to the incompressible limit, staggered grids or mixed finite element methods are a common choice and the pressure gradient is usually computed in a centred and implicit way. On the other side of the Mach number range, the important points are the use of the conservative form of the equations, a preference for explicit schemes and a rare recourse to staggered arrangement of the variables. Of course, the design of a numerical method that would be equally valid and efficient on the whole range of Mach number is a kind of Graal for CFD researchers and numerous attempts in this direction have been reported. Some of these works extend to the compressible regime the numerical methods used for the computation of incompressible flows. Examples of these type of methods to cite but a few are Bijl and Wesseling (1998) or Zienkiewicz et al. (1990) and Munz et al. (2003). A large number of works, also approach the computation of low Mach number flows from the point of view of time discretization and proposed pressure correction methods (Degond and Tang, 2011; Gallouet et al., 2008; van der Heul et al., 2003), the origin of these can be traced back to the ICE method of Harlow and Amsden in the 1970s (Harlow and Amsden, 1971). From the compressible side, since incompressible flows are a particular case of compressible ones, one can argue that in principle, a compressible flow solver should be able to compute these flows. Unfortunately, there are experimental evidences showing that on a fixed mesh, the solutions of the compressible flow discretized equations are not an accurate approximation of the solutions of the incompressible model. One of the first documented article concerned with this problem is a paper by Volpe in 1993 who in a very lucid way made the observation that compressible flow solvers fail to compute near incompressible flows and ask: “why ?” This question has generated an important literature and we will not try to give here an exhaustive report of it. It appears that the answer is deeply linked to the complex behaviour of the compressible models in the low Mach number regime. The crucial point to understand is that in general the strong limit solution of the compressible Euler (or Navier–Stokes) model is not described by the corresponding incompressible equations. Once, this point has been cleared, the answer to the question by Volpe is relatively easy and consequently the design of compressible flow solvers able to compute low Mach number flows can be undertaken. In this paper, we have selected three basic strategies that have been proposed along the years to enhance the capabilities

On the Behaviour of Upwind Schemes Chapter

8 205

of shock capturing techniques to compute low Mach flows as these three types of techniques are representative of the main works in this area. The plan of this work is as follows. In Section 2, we recall some results on the behaviour of the solutions of the Euler equations in the low Mach number limit. This section will show that this behaviour cannot be identified to a simple strong convergence towards the solutions of the incompressible flow system. In Section 3, we will use these results to understand why compressible flow solvers do not deliver accurate approximations of the solutions of the incompressible models, and we will describe some remedies proposed in the literature to correct the deficiencies of the usual upwind schemes. The last section will contain some numerical experiments illustrating the behaviour of these different methods.

2 THE MULTIPLE LOW MACH NUMBER LIMITS OF THE COMPRESSIBLE EULER EQUATIONS 2.1

Incompressible Limit

We begin in this section to recall how to derive in a formal way the incompressible system as a singular limit of the compressible Euler equations. As we will see in the following sections, this simple analysis is incomplete and does not consider the multiple scale character of the Euler system. The compressible Euler system writes @ r + div ru ¼ 0 @t

(1a)

@ ru + div ru  u + rp ¼ 0 @t

(1b)

@ p + u  rp + ra2 div u ¼ 0 @t

(1c)

where r is the fluid density, u the velocity and p the pressure while pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ @p=@rjs is the sound velocity. Let us introduce the following normalization variables:  r¼

r , r



u ¼

u , u



p ¼

p r ða Þ

, 2



x ¼

x , d



t ¼

tu d

(2)

where a* is a reference speed of sound and d* is an arbitrary length scale. We remark that the normalization of the pressure follows directly from the definition of the speed of sound. We then obtain the system of nondimensionalized equations: @t r + div ru ¼ 0 @t ru + div ru  u +

1 rp ¼ 0 e2

@t p + u  rp + ra2 div u + ¼ 0

(3a) (3b) (3c)

206 Handbook of Numerical Analysis

where e ¼ u*/a* is the reference Mach number. We now look for solutions of the system (3) in the form of asymptotic expansion in power of the Mach number:  



















 ðr , u, p Þ ¼ ðr 0 , u 0 , p 0 Þ + eðr 1 , u 1, p 1 Þ + e2 ðr 2 , u 2 , p 2 Þ + ⋯

(4)

Introducing these expressions in Eq. (3) and collecting the terms with equal power of e we obtain (we have dropped the  for convenience): 1. Order 1/e2 rp0 ¼ 0

(5)

rp1 ¼ 0

(6)

@ r + div ðr0 u0 Þ ¼ 0 @t 0

(7a)

@ r u0 + div ðr0 u0  u0 Þ + rp2 ¼ 0 @t 0

(7b)

@ p0 + u0  rp0 + r0 a20 div u0 ¼ 0 @t

(7c)

2. Order 1/e

3. Order 1

The order 1/e2 and 1/e equations imply that the pressure is constant in space up to fluctuations of order e2. Thus we may write: pðx, tÞ ¼ p0 ðtÞ + e2 p2 ðx, tÞ In the presence of open boundaries, the thermodynamic pressure p0 will be imposed and be equal to the exterior pressure. For the sake of simplicity, we assume that the exterior pressure do not change with time and thus the thermodynamic pressure p0 will be a constant in space and time: dPext dp0 ¼ ¼0 dt dt and we will have pðx,tÞ ¼ p0 + e2 p2 ðx, tÞ

p0 ¼ cte

(8)

thus Eq. (7c) implies that the order 0 velocity is divergence free. Note that in the context of nonbarotropic flow, it is the pressure (or energy equation) that implies the incompressible constraint. Some authors (e.g. Klein, 1995) have emphasized this point to stress that incompressible flows can support nonconstant densities. Introducing the divergence constraint into the continuity Eq. (7a), we see that this implies that the material derivative of the density

On the Behaviour of Upwind Schemes Chapter

8 207

is zero. Here for simplicity, we assume that all particle paths come from regions with the same density and conclude that: r0 ¼ Cte The order 1 system thus reduces to incompressible constant density Euler system:  r0

r0 ¼ Cte

(9a) 

@ u0 + div ðu0  u0 Þ + rp2 ¼ 0 @t div ðu0 Þ ¼ 0

2.2

(9b) (9c)

Acoustic Limit

As already mentioned the previous analysis is incomplete. Actually, a fundamental assumption is hidden in the choice of the reference scalings (2). By choosing t ¼ d =u we have implicitely assumed that the reference time of interest is associated with phenomena moving at the material velocity u . However, another choice is possible and one can instead consider the alternate time scale: t ¼ d =a where emphasis is put on phenomena occurring at the speed of sound. The ratio of these two time scales is precisely the Mach number. Using the new dimensionless time a*t/d, one obtains instead of (3) the new nondimensional system: 8 1 > > @t r + div ru ¼ 0 (10a) > > e > > > < 1 1 (10b) @t ru + div ru  u + 2 rp ¼ 0 > e e > > > > > 1 > : @t p + u  rp + ra2 div u ¼ 0 (10c) e The analysis of this system can proceed as the analysis of (3). Introducing the asymptotic expansion in power of the Mach number (4) we will get as previously at order 1/e2 rp0 ¼ 0 but the order 1/e equations now give 8 @t r0 ¼ 0 > < @t r0 u0 + rp1 ¼ 0 > : @t p0 ¼ 0

(11)

(12a) (12b) (12c)

208 Handbook of Numerical Analysis

therefore for this time scale, the order 0 pressure is a constant (in time and space) and the density is constant in time. For simplicity, we will assume that the density at t ¼ 0 is constant in space (note that acoustic can be defined for a space varying background density: this assumption is therefore not mandatory but largely simplifies the analysis). In order to close the system, we have to consider the order 1 pressure equation that gives: @t p1 + a20 div r0 u0 ¼ 0

(13)

where a0pisffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the constant in space and time background speed of sound (i.e a0 ¼ @p=@rjs ðr0 , p0 Þ). Eqs. (12b) and (13) constitute a closed system of equations that in the sequel, we will called the acoustic system. Actually, taking the divergence of (12b) and using (13), we obtain the wave equation @tt p1  a20 5 p1 ¼ 0

(14)

a second-order hyperbolic equation that describes the propagation of linear waves with speed a0

2.3 Acoustic–Incompressible Interactions The two limits that we have just established represent two very different systems from a mathematical point of view: the incompressible system is elliptic in nature (and this why this limit is called a singular limit: the mathematical nature of the system changes from hyperbolic to elliptic) while the acoustic one is hyperbolic and describes the propagation of waves of velocity a0. The problem we face now is to determine which of the two low Mach limits, the incompressible, or the acoustic one is relevant. The answer is that both limits are relevant: Our common experience shows us that we do not have incompressible phenomena or acoustic ones but that we have at the same time incompressible and acoustic phenomena. In the physical world, the two limits coexist and superpose. For instance, the computation of the aerodynamical parameters (drag, lift, etc.) of a slow moving body, say the reader’s car for instance,a can be done with great accuracy with an incompressible system of equations. But at the same time, this moving object is accompanied by sound emission and transmission.b The question is therefore to describe the superposition of the two phenomena and their interactions (if any). This is not a simple question and the answer depends on many factors: the physical model can be the Euler system of some dissipative system (Navier– Stokes equations), the state law may be isentropic or not, the boundary conditions and the domain that can be bounded, periodic or unbounded play also a crucial importance in the results as well as the nature of the initial data that

a

Specially if the reader respects the legal speed limits. Specially if the reader does not respect the legal decibel limits.

b

On the Behaviour of Upwind Schemes Chapter

8 209

can be closed to an incompressible solution or not. From a mathematical point of view, this question has be examined in many works including but not restricted to Klainerman and Majda (1982), Ukai (1986), Schochet (1988), Lions (1993), Schochet (1994), Grenier (1997), Gallagher (1998) and Lions and Masmoudi (1998). Chapter 2 of Majda (2012) is a good introduction to these questions. The article of Schochet (2005) and the whole volume of Abgrall and Guillard (2005) are devoted to the analysis of several mathematical and numerical problems related to low Mach number flows while a review of the mathematical questions related to the convergence towards the incompressible limit can be found in Alazard (2008). In this article to simplify the problem while retaining its most salient feature, we consider the periodic case where the spatial domain O ¼ d , where d is the two- or three-dimensional torus and we assume that the flow is barotropic, that is, there exists a direct relation between pressure and density p ¼ p(r). The case of nonbarotropic flows is more complex and the establishment of general results much more technical. For an overview of the mathematical problems encountered in the nonbarotropic case, we refer to Metivier and Schochet (2001), M
etivier and Schochet (2001) and again to Alazard (2008). For the sake of simplicity, we will in addition assume that the flow is isothermal: p ¼ ra2 with a ¼ constant. This assumption is by no means essential, and all the results established in this section are equally valid for barotropic flows. However, this assumption largely simplifies the algebra and makes clearer the essence of the results. With these assumptions, the nondimensional system that we want to analyze writes 8 @r > > < + div m ¼ 0 @t (15) > > : @m + div ðu  mÞ + 1 rr ¼ 0 @t e2 where we have here chosen to scale the time variable by the slow time scale d*/u*. We remark that in the incompressible limit, pressure fluctuations are of order e2 (see Eq. (8)) p ¼ p0 + Oðe2 Þ while in the acoustic limit, pressure fluctuations are an order of magnitude larger: p ¼ p0 + OðeÞ but in the two case, the order 0 pressure is constant, this justifies the following change of variables: r ¼ 1 + ec

210 Handbook of Numerical Analysis

where p0 ¼ 1 with an appropriate scaling of the pressure. Thanks to this change of variable, system (15) becomes: 8 @ e 1 > e > > (16a) <@t c + e div m ¼ 0 > 1 >@me > (16b) : + div ðm  mÞ + rce ¼ OðeÞ @t e since u ¼ m + OðeÞ. Let us denote ve ¼ (ce, me)t the vector of the dependent variables, we also introduce the notations:   div m v ¼ rc while Qðv, vÞ is a quadratic operator defined by   0 Qðv, vÞ ¼ mm and the system (16) can be written: @ve 1 + Qðve , ve Þ + ve ¼ OðeÞ @t e to which is associated the initial data:  t ve ðt ¼ 0Þ ¼ ce0 , me0 on O

(17)

(18)

An important remark is that the operator  is a linear skew symmetric operator with constant coefficients, thus we have for any j: ¼ (j1, …, jd) Z   jd @xj ve  @xj ve ¼ 0 where @xj : ¼ @xj11,…, ,…, xd O

for all sufficiently smooth ve. This is the mathematical translation of the fact that the acoustic operator does not create or dissipate energy but only transport it across the domain. Thus despite the factor 1/e that appears in (17) in front of this operator, the Hs(O) energy estimates obtained by taking Z   @xj ve  @xj ½equation ð17Þ O

do not depend on e. This has been used in Klainerman and Majda (1981, 1982) (see also Kreiss et al., 1991 and Chapter 2 of Majda, 2012) to establish the existence of uniformly bounded in C([0, T], Hs(O)) solutions of (17) on a time T independent of e by using classical iteration techniques relying on Friedrich inequalities for hyperbolic symmetric systems.

On the Behaviour of Upwind Schemes Chapter

8 211

However, for general initial data, this is not sufficient to obtain a limit on ve when e ! 0. The basic reason is that we have no uniform bound on the time derivative @ tve since Eq. (17) precisely says that the time derivative is of order 1/e. An important exception concerns the case where the initial data are “wellprepared”: Assume that the initial conditions can be written ve ðt ¼ 0Þ ¼ v0 + eve1 where v0 and ve1 are bounded in Hs(O) and where v0 2 Ker, then the e term will balance with the 1/e factor in front of the operator ve , giving a bound independent of e on the time derivative of ve. This allows to establish the convergence of the solution of (17) to the solution of the incompressible system with initial data v0. For rigorous results concerning the case of well-prepared initial data, one can see Klainerman and Majda (1981), Klainerman and Majda (1982), Majda (2012), Schochet (1988) as well as the interesting works of Kreiss and his coworkers on the “bounded derivative” principle (Browning and Kreiss, 1982; Gustafsson, 1908; Kreiss, 1980; Kreiss et al., 1991). However, in the case of general initial data, where the time derivative at t ¼ 0 is genuinely of order 1/e, one has to use a different strategy. Such a strategy has been introduced in Schochet (1994) and then used with different variations in other works (Gallagher, 1998; Grenier, 1997; Lions and Masmoudi, 1998). The key idea behind this strategy is that in the presence of ill-prepared data, the acoustic operator  will be responsible of the generation of fast oscillatory waves on time scales of order e. It is therefore necessary to “filter” these components of the solution in order to bound the time derivative. This is done as follows: Let us introduce the group LðtÞ ¼ expðtÞ

(19)

of the solutions of the acoustic equation. In other words vðt, xÞ ¼ LðtÞv0 ðxÞ

if

@v + v ¼ 0 @t

with vðt ¼ 0, xÞ ¼ v0 ðxÞ

Now consider the “filtered variable” we defined by we ¼ Lðt=eÞve

(20)

From the definition of L, we deduce that @we  e @ve ¼ w + Lðt=eÞ @t @t e and then with (17) that we satisfies: @we + Lðt=eÞQðLðt=eÞwe , Lðt=eÞwe Þ ¼ OðeÞ @t

(21)

with initial data ve(t ¼ 0) as Lð0Þ is the identity. In contrast with Eq. (17), (21) does not contain any unbounded 1/e term. Indeed using again that  is a skew-symmetric linear operator with constant coefficient, L is an isometry:

212 Handbook of Numerical Analysis

k LðtÞvkHs ¼k vkHs and thus (21) shows that the time derivative of we is uniformly bounded independently of e. We can then establish that we converge to a function w solution of the system @w + Qðw, wÞ ¼ 0 @t

(22)

where Qðw, wÞ is an averaged operator defined as the limit of Lðt=eÞQðLðt=eÞw, Lðt=eÞwÞ The original variable v can then be obtained by “un-filtering” back the solution of the filtered equations and we have for all e sufficiently small ve ¼ Lðt=eÞw + oð1Þ

(23)

The system (22) describe completely the behaviour of the solution for sufficiently small Mach number. To obtain an explicit description of the solution, we can reformulate this system by splitting its solution w into a slow part ws and a fast one wf defined as ws ¼ Pw

wf ¼ ðI  PÞw

where P is the L2 projection into the kernel of , that is, for any w ¼ (f, u)t, we have  t Pw ¼ f, udiv where f is the spatial mean of f while udiv ¼ Pu ¼ u rD1ru is the L2 projection of the velocity on the space of divergence-free vector fields. Applying the projection P to (22) we get @ws + PQðws , ws Þ + PQðws , wf Þ + PQðwf , ws Þ + PQðwf , wf Þ ¼ 0 @t

(24)

the first line of this system simply gives that the pressure fluctuation cs is zero. To establish the equation on the velocity, we will use that PQðwf , ws Þ ¼ PQðws , wf Þ ¼ 0. The quadratic term involving the fast variable Qðwf , ws Þ can be shown to be written as a gradient (see Schochet, 1994 or Grenier, 1997) then it also disappears under the action of the projection P and the second line reduces to the Leray formulation of incompressible system: @udiv + Pr  ðudiv  udiv Þ ¼ 0 @t

(25)

Finally, using that L is the identity on the kernel of , we have from (23)   0 + Lðt=eÞwf + oð1Þ (26) ve ¼ udiv where udiv solves the incompressible equation.

On the Behaviour of Upwind Schemes Chapter

8 213

Let us summarize the results of this section. In the low Mach number limit, the solution can be written as the sum (26) of a slow part ws(x, t) where the velocity is divergence-free plus a fast component wf(x, t, t/e) that depends on the fast time t/e. Moreover, the slow part is solution of the incompressible equations:  r0

divðws Þ ¼ 0

(27a) 

@ udiv + div ðudiv  udiv Þ + rp ¼ 0 @t

(27b)

In general, udiv(x, t) is only a weak limit of the solutions except for well prepared initial data close to the kernel of . This special class of initial data is characterized by:  cðx, 0Þ ¼ er2 ðx, 0Þ (28) uðx, 0Þ ¼ u01 ðx, 0Þ + eu2 ðx, 0Þ with divu01 ¼ 0 Going back to the original variables (p, u) in dimensional units, we see that this implies that the initial data are such that:  pðx, 0Þ ¼ r a2 ðConstant + e2 p2 ðx, 0ÞÞ (29) uðx, 0Þ ¼ a ð0 + eu1 ðx,0Þ + e2 u2 ðx, 0ÞÞ with divu1 ¼ 0 where p2(x, 0), u1(x, 0), u2(x, 0) are smooth, regular functions.

2.4

Finite Volume Schemes

We consider in this section the application of standard compressible flow solvers to the computation of steady or unsteady solutions of the Euler equations of gas dynamics in the low Mach limit. More precisely, given a fixed mesh, we want to understand the behaviour of finite volume schemes when the reference Mach number (for instance the inflow Mach number for an aerodynamical computation) goes to zero, and in particular why this particular type of approximation fails to compute the incompressible limit. In order to simplify the presentation and to point out the difference with the continuous case, we again restrict ourselves to the isothermal case. Now, let O be a 2D or 3D flow domain that for simplicity, we consider polygonal and let Th be a tessellation of this domain in nonoverlapping control volumes Ci, i 2{1, …, nh}. The present analysis does not assume any particular properties of the mesh Th. As is quite common with compressible solvers, the approximation, which we are dealing with, localizes the degrees of freedom (density, pressure and velocity) in a colocated manner at the centre of the cell Ci. Note that in contrast, incompressible flow solvers use very often staggered grids and that if colocalized grids are used for the computations of incompressible flows, some specific procedures have to be used to overcome the inf-sup criterion. We will come back to this point later. We use the notations of Guillard and Viozat (1999): The set of neighbours of cell

214 Handbook of Numerical Analysis

i is denoted by VðiÞ. The area (volume in 3-D) of cell i is jCij. The symbol il denotes an edge (face) between cell i and cell l whose length (surface) is dil. The vector n il denotes the unit normal vector to edge il oriented from i to l. For any variable f, we define the jump operator Dilf as Dil f ¼ fi  fl

(30)

and we adopt the specific notation DilU to denote the normal component of the velocity jump at the interface between two cells, denoted by Dil U ¼ ðu i  u l Þ:n il

(31)

Finally, the notation fil will denote some average (here for simplicity the Roe average) of the variable f between fi and fl. With these notations, the application of a first-order upwind finite volume method yields the following set of semidiscrete equations:   X d ri + jCi j dil Fðvi , vl , n il Þ ¼ 0 dt mi l2VðiÞ

where Fðvi ,vl ,n il Þ is the numerical flux. It is customary to express the numerical flux as a sum of a centred part and a numerical dissipation: 1 jAil j Dil v Fðvi , vl , n il Þ ¼ ðFðvi Þ  n il + Fðvl Þ  n il Þ  2 2 where Ail is some approximation of the Jacobian matrix @FðvÞ  n il =@v. From the hyperbolicity of the system, this matrix can be diagonalized and the numerical flux can be written as: 1 1 Fðvi , vl , n il Þ ¼ ðFðvi Þ  n il + Fðvl Þ  n il Þ  Ril jLil jR1 il Dil v 2 2

(32)

where Ril is the matrix composed of the right eigenvectors and Lil, the diagonal matrix formed by the corresponding eigenvalues. In the present case for instance, the application of Roe scheme results in: Fðv8 i , vl , n il Þ Uil > > ðrl u l :n il + ri u i :n il Þ + a Dil r + ril Dil U > > a >     1< Uil2 ¼ Dil U n il ðr u :n u + r u :n u + r n Þ + 2a Uil Dil r + ril a + 2> > > l l il l i i il i l il a > > : +jUil jðril Dil V + VDil rÞn?il where in this expression, a* is the constant speed of sound and the terms proportional to a jump Dil represent a numerical dissipation.

On the Behaviour of Upwind Schemes Chapter

2.5

8 215

The Diagnosis

The analysis of this scheme will be done using exactly the same strategy than in the continuous case. First we normalize the variables using the reference quantities given in (2). As in the continuous case, the reference Mach number e ¼ u*/a* appears explicitly in the equations and we obtained the following set of scaled semidiscrete equations Mass conservation d 1X dil rl ðu l  n il Þ jCi j ri + dt 2 l2VðiÞ 1 X eX + dil Dil r + dil ril Uil Dil U ¼ 0 2e l2VðiÞ 2 l2VðiÞ

(33a)

Momentum conservation d 1X 1 X jCi j ri u i + dil rl ðu l  n il Þu l + 2 dil rl n il dt 2 l2VðiÞ 2e l2VðiÞ  1X 1 X eX + ðdil Uil Dil rÞn il + ðril Dil UÞn il + r U2 Dil U n il e l2VðiÞ 2e l2VðiÞ 2 l2VðiÞ il il 1X + dil jUil jðril Dil V + Vil Dil rÞn?il ¼ 0 2 l2VðiÞ (33b) then we again proceed as in the continuous case and perform the change of variablecr ¼ 1 + ec to obtain: 1Xd  X dil  d il jCi j ci + u l :n il + Dil c ¼ OðeÞ cl ðu l  n il Þ + Uil Dil U + 2 dt e l2VðiÞ 2 l2VðiÞ (34a)

X dil  d jCi j u i + ðu l  n il Þu l + cil ðDil UÞn il + 2Uil Dil c + jUil jðDil VÞn?il 2 dt l2VðiÞ 1 X dil  + cl n il + Dil Un il ¼ OðeÞ e l2VðiÞ 2

(34b)

c

This requires, however, to prove that the 0-order pressure is constant, we assume this result here and refer to Guillard and Viozat (1999) for the proof.

216 Handbook of Numerical Analysis

This change of variable “symmetrize” the system and as in the continuous case, we can write the previous equations in a compact form similar to (17): @veh 1 + Qh ðveh , veh Þ + h veh ¼ OðeÞ @t e

(35)

where as in the continuous case, the discrete acoustic operator h defined by:   1 X Dil c + u l :n il dil ½h vh i ¼ (36) cl n il + Dil Un il 2jCi j l2VðiÞ is a linear constant coefficient operator. The terms u l :n il and cl n il obviously represent approximations of div u and rc. However, in contrast with the continuous case where the acoustic operator is defined by   div u v ¼ rc the discrete acoustic operator contains also artificial viscosity terms Dilc and Dil Un il that are approximation of some second-order dissipative terms with diffusion coefficient proportional to the mesh size. In view of the results of Section 2, the origin of these terms is clear. The low Mach behaviour of the continuous Euler equations contains two potential limits: the incompressible one and the acoustic one. But the pressure fluctuations associated to the acoustic limit are an order of magnitude larger than the ones associated with the incompressible limit (see Section 2.2). Consequently, these terms appear as the dominant ones in the artificial viscosity of upwind schemes. Indeed, if one interpret Fðvi , vl , n il Þ the numerical flux at the interface between two cells as the flux coming from the solution of a Riemann problem between the states vi, vl, it has been shown in Guillard and Murrone (2004) that the interface state is actually the solution of a linear acoustic problem defined by the background state. From a practical point of view, the immediate consequence is that on a fixed mesh, the discrete solutions of compressible upwind solvers cannot converge to a reasonable approximation of the incompressible solution. Instead the discrete solutions will converge to the nonzero solution of some Stokes like problem defined by:

To see that, subtract u0i 

X l2V

½h vh  ¼ 0 dil n il ¼ 0 to the first line of (36) and ci 

X l2V

(37) dil n il ¼ 0

to the second line, it is seen that the system (37) can be written as X

ðDil c  Dil UÞdil ¼ 0

(38a)

l2VðiÞ

X

l2VðiÞ

ðDil c  Dil UÞdil n il ¼ 0

(38b)

On the Behaviour of Upwind Schemes Chapter

8 217

This can be considered as a linear system of 3  nh equations (nh being the number of cells), for the ne/2 unknowns Dilc DilU where ne is the number of edges.d In general, this system have nonzero solutions implying that the discrete solution will contain nonconstant c corresponding to pressure fluctuations of order e, while the pressure fluctuations associated with incompressible phenomena are of order e2. The analysis performed here for the barotropic equations extends easily to the complete Euler system with the same results. From now on, we consider again the full complete Euler system since this system is the one used in practical applications.

2.6

The Remedies

With this analysis at hand, one can propose different cures to recover a proper convergence towards the incompressible limit. The common point of these approaches will be to modify the artificial viscosity to suppress the unwanted acoustic fluctuations.

2.6.1 Preconditioned Dissipation: The Roe–Turkel Scheme This strategy is the oldest and was inspired from the preconditioned approach introduced in Turkel (1993, 1999), Turkel et al. (1994, 1997) and Choi and Merkle (1993) to overcome the convergence problems encountered by compressible solvers in the computation of steady states in the low Mach limit. In these approaches the standard scheme X d dil Fðvi , vl , n il Þ ¼ 0 jCi j vi + dt l2VðiÞ is modified by the introduction of a preconditioning matrix P multiplying the time derivative X d dil Fðvi , vl , n il Þ ¼ 0 P1 jCi j vi + dt l2VðiÞ the matrix P being chosen such that the eigenvalues of P@FðvÞ  n=@v are well balanced and bounded when e ! 0 in order to accelerate the convergence properties towards steady state of the time advancing scheme. The construction of such matrices has been studied in several works. Turkel (1993) presents an overview of the different possibilities, while the preconditioner proposed in Weiss and Smith (1995) is extremely popular in applications. The basic idea of the preconditioned dissipation schemes is to use the good properties of the preconditioning matrices to modify the artificial viscosity matrices in order to recover the convergence towards the incompressible d

The precise number of equations and unknowns have actually to take into account the boundary conditions.

218 Handbook of Numerical Analysis

solution that is destroyed by upwind schemes. In the preconditioned dissipation schemes, the numerical flux is thus changed into: !

Fðqi , ql , n!il Þ ¼

!

F ðqi Þ + F ðql Þ 1  n il + Pðqil Þ1 j Pðqil ÞAðqil Þ j Dil q 2 2

(39)

where P(qil) is a preconditioning matrix constructed on each edges. With respect to the original Roe scheme, only the dissipative terms are altered and therefore the numerical scheme remains a consistent approximation of the time dependent compressible Euler equations. The term Roe–Turkel scheme was coined to designate this particular type of low Mach scheme. In Guillard and Viozat (1999), an asymptotic analysis of why this approach allows to obtain a correct scaling of the pressure fluctuation with the square of the Mach number was presented. Actually, this scheme modifies in a very important way the dissipation matrix. In a Mach number expansion, the dominant term of the dissipation matrix becomes instead of (36) l

Continuity equation X

Dil c dil pffiffiffiffiffi ¼ 0 Yil l2VðiÞ l

Momentum equation X l2VðiÞ

 dil

 ðUn + 2uÞil pffiffiffiffiffi Dil c ¼ 0 cl n il + Yil

(40)

(41)

where Y is an interface quantity positive and of order 1. Comparing the case of the preconditioned dissipation with the one of the original upwind scheme, we note that there is an important difference in the structure of the dissipation terms. In the original Roe scheme, the 1/e equations link the pressure jumps with the normal velocity jumps. A consequence is that the order e pressure c cannot be constant. On the other hand, for the preconditioned dissipation, the equations involve only the pressure jumps and force c to be a constant. It must also be noted that the preconditioned dissipation schemes realize the convergence towards incompressible solution by augmenting the numerical dissipation on the pressure. On the opposite, the next low Mach scheme to be described try to attain this goal by reducing the numerical viscosity.

2.6.2 The All-Speed Scheme This approach has been introduced in a series of works (Li and Gu, 2008, 2013; Li et al., 2009) and stems from the observation that AUSM+ schemes used for large eddy simulation have good low dissipation properties for low Mach flows (Mary and Sagaut, 2002). While in the Roe–Turkel scheme, both the eigenvalues and the eigenvectors of the numerical dissipation are

On the Behaviour of Upwind Schemes Chapter

8 219

modified, this technique proposes to change only the eigenvalues. Thus the expression of the numerical flux becomes: 1 1 Fðvi , vl , n il Þ ¼ ðFðvi Þ  n il + Fðvl Þ  n il Þ  Ril jLAil jR1 il Dil v 2 2

(42)

where the u  n il a eigenvalues of Lil are modified into u  n il f ðeÞa with f(e) a smoothly varying function that goes to zero and 1 with e. This scheme reveals itself extremely oscillatory and then it was proposed to modify the centred part of the flux using a momentum interpolation method as in Rhie and Chow (1983) and Mary and Sagaut (2002). In Li and Gu (2008), the fact that the Roe–Turkel scheme of Section 2.6.1 does not satisfy a divergence-free constraint at Mach ¼ 0 was also presented as an argument to prefer the approach (42) since in this latter approach, the discrete velocity exactly verifies at e ¼ 0 X dil u l :n il ¼ 0 (43) l2VðiÞ

This argument has to be used with some caution. Actually and as recalled earlier, compressible flow solvers of the type studied here use a colocated arrangement of the variables and therefore this will be also the case for the limit scheme obtained at Mach ¼ 0. But from finite element theory applied to incompressible flows, we know that the choice of the pressure and velocity approximation spaces cannot be arbitrary and have to satisfy the inf-sup criterion to yield an unique solution. We also know that colocalized arrangement of variable does not satisfy the inf-sup criterion and therefore that these approximations suffer from the presence of spurious pressure modes. To overcome this difficulty, approximations of incompressible models use therefore mixed finite element spaces or staggered grids. Another way to overcome this difficulty is to use a stabilized finite element formulation as introduced in Hughes et al. (1986). In this approach, to stabilize the Galerkin formulation of the Stokes (or Darcy) problem that yield a discrete problem of the form:    ch 0 rh  ¼F Dh vh rh stabilization matrices are added to the diagonal giving a discrete problem with the structure:    rh  ch Mh ¼F Dh vh rh and it is shown in Hughes et al. (1986) that a proper choice of the stabilization matrices preserves the convergence rate of the approximation. We observe that the Roe–Turkel scheme precisely adds in the discrete mass conservation equation a similar stabilization term proportional to the pressure jumps while

220 Handbook of Numerical Analysis

this term is absent from (42). The momentum interpolation technique used in Li and Gu (2008), Li et al. (2009) and Li and Gu (2013) is another way to reintroduce this stabilization term to avoid the appearance of spurious pressure modes. We do not elaborate on the subject anymore but emphasize that there are sound theoretical reasons to actually prefer a scheme that does not enforce the discrete divergence-free constraint (43) over those where this constraint is respected. The numerical experiments section will give a clear illustration of this point.

2.6.3 Rieper Fix This technique introduced in Rieper (2011) represents a further step in attempts to simplify the preconditioned dissipation technique. In this last method, both the eigenvectors and the eigenvalues of the dissipation matrix are modified. In the All-Mach scheme (Li and Gu, 2008, 2013; Li et al., 2009), the modification concerns only the eigenvalues.e Here, the coefficient of the dissipation matrix is directly modified. The analysis conducting to the Rieper fix comes from the observation in Thornber et al. (2008a) that the unphysical entropy production of upwind schemes at low Mach numbers is related to the jumps of the normal velocity component at the cell interfaces. Then in Thornber et al. (2008b) an improved reconstruction method of MUSCL type was presented to minimize these jumps to achieve high accuracy at low Mach number. The same idea in a quite extreme manner was also used in Dellacherie et al. (2010) where the normal velocity jumps DilU were set exactly to zero in the artificial viscosity matrix, leading to a centring of the pressure gradient in the momentum equation. In Rieper (2011), the same idea is exploited in a smoothest way. Going back to the expression of the dominant term (36) of upwind schemes:   1 X Dil c + u l :n il dil ½h vh i ¼ (44) cl n il + Dil Un il 2jCi j l2VðiÞ

the numerical dissipation in the momentum equation associated to the DilU term in the second line is interpreted as preventing the pressure variable c to be constant. Rieper (2011) thus proposes to change this term into: Dil U ! min ðe, 1ÞDil U. At the dominant term in an expansion in Mach number the dissipation term will thus be:   1 X Dil c + u l :n il dil ½h vh i ¼ (45) cl n il 2jCi j l2VðiÞ

and the presence of nonconstant pressure c will be suppressed. Later on, this idea was pursue in Oswald et al. (2016) by reducing also the jump in the transverse velocity in an attempt to use as little numerical viscosity as possible e

In principle, however, in practice, use of momentum interpolation techniques is required.

On the Behaviour of Upwind Schemes Chapter

8 221

while preserving the stability. However, we remark from (45) that if the firstorder pressure c is a constant then the scheme will verify the divergence-free constraint (43) and according to the discussion in the previous section, the presence of spurious pressure fluctuations cannot be ruled out a priori. In Rieper (2011), arguments are presented to show that the scheme avoids the presence of checkerboard pressure modes. These argument, however, are only valid for regular cartesian meshes, and it is not sure that they can be applied for other type of meshes. Actually, in the following numerical experiment section, we will exhibit a situation where the Rieper fix produces pressure fluctuation of the wrong order while the results obtained with the classical standard Roe scheme are perfectly consistent with the incompressible limit.

3

NUMERICAL ILLUSTRATIONS

We use as numerical test the simple case of the computation at a low Mach number of the flow around a cylinder of radius r0. This test case is useful since an analytical reference solution of the incompressible solution is known. The flow O domain is an annulus [r0, r1]  [0, 2p]. Here, we have used r0 ¼ 0.5 and r1 ¼ 5. The initial data are uniform and set equal to r0 ¼ 1, u0 ¼ ðu0 , 0Þt , p0 ¼ 1

pffiffiffi with u0 ¼ gM where g ¼ 1.4 is the adiabatic index and M* is the mach number at infinity set here equal to 0.1. The exact solution at infinity is uniform pffiffiffi r∞ ¼ 1, u0 ¼ ð gM , 0Þt , p∞ ¼ 1 In all the computation presented here, the boundary conditions are implemented by considering a layer of ghost cells whose value is set to the exact solution at infinity. The computations are time-advanced until a steady state characterized by a residual of 107 is reached. The time advancing scheme is an implicit linearized method using Jacobian matrices corresponding to each of the four studied numerical flux (Roe, Roe–Turkel, All-Speed, Rieper). As noted in Birken and Meister (2005) (see also Dellacherie, 2010) the Roe–Turkel scheme suffers from a severe stability restriction Dt ¼ OðMa2 Þ that forbids the use of time explicit integration for this method. The other three schemes are stable under the standard CFL stability criterion Dt ¼ OðMaÞ and in principle explicit schemes can be used with these schemes although explicit time advancing becomes more and more costly as the Mach number is decreased.

3.1

Order 1: Quadrangular Cartesian Grids

The first set of numerical experiments is done using a polar quadrangular grid with a resolution of 50 (radius) by 150 (angle). A first-order accurate scheme is used.

222 Handbook of Numerical Analysis Contour Var: Pressure fluctuations

4.0

Contour Var: Pressure fluctuations

4.0

Max: 0.01350 Min: —0.01730

2.0

Y-axis

Y-axis

2.0

0.0

0.0

–2.0

–2.0

–4.0

–4.0 –4

–2

0 X-axis

2

4

–4

Contour Var: Pressure fluctuations

4.0

–2

0 X-axis

2

4

0 X-axis

2

4

Contour Var: Pressure fluctuations

4.0

Max: 0.007786 Min: —0.02019

2.0

Max: 0.008425 Min: —0.01984

2.0

Y-axis

Y-axis

Max: 0.007633 Min: —0.01797

0.0

0.0

–2.0

–2.0

–4.0

–4.0 –4

–2

0 X-axis

2

4

–4

–2

FIG. 1 Pressure fluctuations. Top left: Roe scheme. Top right: Roe–Turkel scheme. Bottom left: Rieper fix. Bottom right: all-speed scheme.

Fig. 1 displays the pressure fluctuations defined as dp ¼ (p  p0)/p0. 20 isocontours between  0.007 and 0.007 (0:007 ¼ gM2 =2) are displayed. As expected, the Roe scheme gives results quite far from the incompressible solution while the three modified low Mach schemes present a good agreement with the incompressible solution. The Roe–Turkel (Guillard and Viozat, 1999) and Rieper fix (Rieper, 2011) are extremely close to each others. The All-Speed scheme (Li and Gu, 2008) presents some small oscillations near the inflow boundary but is also very close to the other two low Mach modified schemes. For this last scheme, we have used the momentum interpolation method (the equation (8) and (9) of Li and Gu, 2008) with a value c2 ¼ 0.05, slightly larger than the recommended value of 0.04 in an attempt to reduce the possible apparition of spurious pressure modes (see the discussion in Section 2.6.2).

On the Behaviour of Upwind Schemes Chapter

3.2

8 223

Order 1: Vertex-Centred Triangular Meshes

In our next set of experiments, we construct a triangular mesh by cutting each quadrangle of the previous mesh in two by the diagonal in an alternate way (the so-called English flag mesh). Then around each node of the triangulation, a dual control volume is constructed by means of the barycentres of the triangles and the median of the edge. The number of degrees of freedom of this mesh is thus equal to the number of nodes of the triangulation and thus almost equal to the number of degrees of freedom of the previous quadrangular mesh (the difference is due to the boundary condition). However, in contrast to the previous quadrangle mesh, the number of neighbours of one node is not constant and is alternatively equal to 4 or 8. This kind is mesh is thus a good test to judge of the capability of the scheme to support or suppress spurious pressure modes (Fig. 2). Contour Var: Pressure fluctuations

4.0

Contour Var: Pressure fluctuations

4.0

Max: 0.01302 Min: —0.01415

2.0

Y-axis

Y-axis

2.0

0.0

0.0

–2.0

–2.0

–4.0

–4.0 –4

–2

Contour Var: Pressure fluctuations

4.0

0 X-axis

2

4

–4

–2

Contour Var: Pressure fluctuations

4.0

Max: 0.007547 Min: —0.01732

2.0

0 X-axis

2

4

0 X-axis

2

4

Max: 0.008198 Min: —0.01651

2.0

Y-axis

Y-axis

Max: 0.007523 Min: —0.01718

0.0

0.0

–2.0

–2.0

–4.0

–4.0 –4

–2

0 X-axis

2

4

–4

–2

FIG. 2 Pressure fluctuations. Top left: Roe scheme. Top right: Roe–Turkel scheme. Bottom left: Rieper fix. Bottom right: all-speed scheme.

224 Handbook of Numerical Analysis

On this kind of meshes, large differences between the schemes appear. The solution obtained with the classical Roe scheme is clearly deteriorated with respect to the one obtained on the quadrangular mesh. The Roe–Turkel does not seem to be much sensible to the mesh and give a result almost identical to the one obtained on a quadrangular much. The solutions obtained with the All-Speed scheme and Rieper fix are closed together and are clearly worse than on the quadrangular mesh. Moreover, oscillations appear in the results. The method of Rieper (2011) seems to have a better control of these oscillations but does not seem to escape totally to the presence of spurious pressure modes. As explained in Section 2.6.2, the origin of these oscillations can be traced back to the fact that these methods verify the discrete divergence-free constraint (43) while the Roe–Turkel scheme does not. The respect of this constraint is clearly not an advantage on the present mesh.

3.3 Order 1: Cell-Centred Triangular Meshes These experiments use the same triangle mesh as in the previous section but the discretization employs the cell-centred approach where the finite volume cells are the triangle themselves. Fig. 3 shows that the results are surprisingly good even for the classical Roe. At first glance, this seems to be in contradiction with the analysis of Section 2.5. This point has been noted for the first time in Rieper and Bader (2009) where is was shown that the classical Roe scheme do converge to the incompressible limit on triangular cell-centred meshes. Rieper and Bader (2009) present an analysis of this behaviour on special 2D triangular grids. Later on, Guillard (2009) extends this results on general triangulation in 2D or 3D. To understand this unexpected behaviour, let us go back to the expression (38) that gives the behaviour at the dominant order of the discrete solution X ðDil c  Dil UÞdil ¼ 0 (46a) l2VðiÞ

X

l2VðiÞ

ðDil c  Dil UÞdil n il ¼ 0

(46b)

For each cell, this expression is a system of d + 1 equations where d is the space dimension (d scalar equations for the momentum + 1 for the pressure). In Section 2.5, we conclude that in general this equation has nonzero solutions and that these nonzero solutions were responsible of the apparition of pressure fluctuations of order OðeÞ. However, in the special case of cell-centred simplicial meshes, considering xil ¼ (Dilc DilU) as the unknowns of (46), we remark that the number of unknowns is also equal

On the Behaviour of Upwind Schemes Chapter

4.0

4.0 –

–

2.0

Y-axis

2.0

Y-axis

8 225

0.0

0.0

–2.0

–2.0

–4.0

–4.0 –4

–2

0 X-axis

2

–4

4

–2

0 X-axis

2

4

–2

0 X-axis

2

4

4.0 –

Y-axis

2.0

0.0

–2.0

–4.0 –4

FIG. 3 Pressure fluctuations. Top left: Roe scheme. Top right: Roe–Turkel scheme. Bottom right: all-speed scheme.

to d + 1 since by definition a simplex has exactly d + 1 neighbours. Thus (46) can be viewed as a square linear system for the xil. It can be proved (Guillard, 2009) that the determinant of the corresponding matrix is nonzero. From this results, it can be deduced that the only possible solution of (46) is constant pressure fields and velocities satisfying a divergence-free constraint. Thus, in contrast with the other type of meshes, the solutions computed by Roe schemes on cell vertex triangular grid converge to the incompressible limit and upwind schemes in this context are particularly efficient to enforce a discrete divergence-free constraint with no spurious pressure modes. It can be seen in Fig. 3 that this property seems to be preserved by the Roe– Turkel and All-Speed schemes. This is, however, not the case for Rieper fix.

226 Handbook of Numerical Analysis

The computations done with this scheme lead to numerical blow-up (negative densities appear). It is easy to understand why. Rieper fix replaces the expression (36) of the large operator by:   1 X Dil c + u l :n il dil ½h vh i ¼ (47) cl n il 2jCi j l2VðiÞ

looking at the kernel of this operator, we can write it in the form (46) as: X X

ðDil c  Dil UÞdil ¼ 0

(48a)

ðDil c  Dil UÞdil n il ¼ Dil Udil n il

(48b)

l2VðiÞ

l2VðiÞ

this shows, for the same reasons proving that (46) has only zero solutions, that (48) has nonzero solutions. On cell-centred triangular meshes, the Rieper fix thus does not reduces pressure fluctuations but on the opposite, create them.

3.4 Order 2: Vertex-Centred Triangular Meshes We conclude this experimental section by some results on second-order schemes that at the present time are the standard of compressible solvers. The largest differences between the schemes appear on the vertex-centred triangular meshes already used in Section 3.2. The conditions are exactly the same than for the order 1 computations. For all the schemes, the same MUSCL procedure using a minmod limiter and a multislope reconstruction is applied. Fig. 4 displays 15 isocontours of pressure fluctuations comprised between 0.007. In addition, this figure is also superposed pseudocolour plots of the fluctuations of entropy comprised between 105 (i.e. white colour if < 105) and 103 (magenta if > 103). The Roe scheme as well as the Roe– Turkel schemes converge to a steady state. The production of entropy by the Roe scheme is important and affect all the downwind part of the cylinder and the wake. This unphysical production of entropy is greatly reduced by the Roe–Turkel scheme and as shown in Fig. 5 the solution presents a very good agreement with the incompressible solution. As regards the other two low Mach schemes, it has not been possible to obtain a steady state. The All-Speed and Rieper schemes exhibit an unsteady behaviour characterized by pseudo von Karman alleys. In addition, some oscillations in the pressure profile are noticeable that again seem to be better controlled by the Rieper scheme. The numerical viscosities that are greatly reduced by these two schemes do not seem to be sufficient to obtain a steady solution in this case.

4.0

4.0 –

–

2.0

Y-axis

Y-axis

2.0

0.0

0.0

–2.0

–2.0

–4.0

–4.0 –4

–2

0 X-axis

2

4

–4

0 X-axis

2

4

–2

0 X-axis

2

4

4.0

4.0 –

–

2.0

Y-axis

2.0 Y-axis

–2

0.0

0.0

–2.0

–2.0

–4.0

–4.0 –4

–2

0 X-axis

2

–4

4

FIG. 4 Pressure (isocontour) and entropy (colour). Top left: Roe scheme. Top Right: Roe–Turkel scheme. Bottom left: Rieper fix. Bottom right: all-speed scheme.

1.5 1.0

Y-axis

0.5 0.0 –0.5 –1.0 –1.5

–1.5

–1.0

–0.5

0.0

0.5

1.0

1.5

X-axis

FIG. 5 Pressure fluctuation isolines: Comparison of the incompressible potential solution (red, dotted line) with the numerical solution (blue) obtained by the order 2 Roe–Turkel scheme.

228 Handbook of Numerical Analysis

4 CONCLUSION Low Mach number flows appear in a wide variety of situations and the design of appropriate discretizations motivates an important activity. A first type of approaches used for this purpose starts from numerical methods originally designed for the approximation of incompressible flows. In this work, we have considered the opposite point of view and have examined how to adapt upwind compressible flow solvers to compute solutions close to the incompressible limit. A large number of experimental evidences show that upwind schemes are unable to compute accurately these flows. We have seen that this failure of compressible flows solvers is deeply connected to the fact that the incompressible limit of the compressible system is not the only possible one. The acoustic limit coexits with the incompressible part of the solution. Moreover, the pressure fluctuations associated to the acoustic limit, dominate by one order of magnitude the pressure variation due to incompressible effects. The consequence is then that the numerical viscosity of upwind scheme is mainly due to acoustic effects and not to incompressible ones. To overcome this difficulty, several different strategies have been proposed in the past. We have selected three of them as representative examples of the different proposed approaches. The rationale behind these approaches has been explained, and we have presented numerical experiments on a simple test case to compare them. The results show that the performance of the different method is highly dependent on the type of mesh used. Overall the preconditioned dissipation (Roe–Turkel) (Guillard and Viozat, 1999) is the less sensible to the type of meshes. The other two methods (Li and Gu, 2008; Rieper, 2011) have good performance on cartesian grids but presents oscillatory results on triangular vertex-centred grids. On the “pathological” case of cell-centred triangular meshes, the low Mach methods preserve the accuracy of the classical Roe scheme except the method of Rieper (2011) that creates huge oscillation leading to numerical divergence of the computations. Of course, it would be comfortable to identify a “best” method among the studied ones. This is unfortunately not the case and all the methods have their weaknesses. The preconditioned dissipation (Guillard and Viozat, 1999) method must be used with time implicit schemes. The All-Speed (Li and Gu, 2008) and Rieper (2011) methods enjoy a standard CFL stability constraint but give oscillatory results that can prevent the convergence to steady state. It must be emphasized that the test-case studied here is a simple inviscid one. Care must be taken before extrapolating the conclusions of this chapter for DNS or LES simulations where the reduction of numerical viscosity is an important issue.

REFERENCES Abgrall, R., Guillard, H., 2005. Special issue on low Mach number flows conference. ESAIM: M2AN 39 (3), 437–621. http://dx.doi.org/10.1051/m2an:2005018. Alazard, T., 2008. A minicourse on the low Mach number limit. Discrete Contin. Dyn. Syst. Ser. S 1 (3), 365–404.

On the Behaviour of Upwind Schemes Chapter

8 229

Bijl, H., Wesseling, P., 1998. A unified method for computing incompressible and compressible flows in boundary fitted coordinate. J. Comput. Phys. 141, 153–173. Birken, P., Meister, A., 2005. Stability of preconditioned finite volume schemes at low Mach numbers. BIT Numer. Math. 45 (3), 463–480. Browning, G., Kreiss, H.O., 1982. Problems with different time scales for non-linear partial differential equations. SIAM J. Appl. Math. 42, 704–718. Choi, Y.H., Merkle, C.L., 1993. The application of preconditioning in viscous flows. J. Comput. Phys. 105, 207–223. Degond, P., Tang, M., 2011. All speed method for the Euler equation in the low Mach number limit. Commun. Comput. Phys. 10, 1–31. Dellacherie, S., 2010. Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number. J. Comput. Phys. 229 (4), 978–1016. ISSN 0021-9991. http://dx.doi.org/10.1016/j. jcp.2009.09.044. http://www.sciencedirect.com/science/article/pii/S0021999109005361. Dellacherie, S., Omnes, P., Rieper, F., 2010. The influence of cell geometry on the Godunov scheme applied to the linear wave equation. J. Comput. Phys. 229 (14), 5315–5338. Gallagher, I., 1998. Asymptotic of the solutions of hyperbolic equations with a skew-symmetric perturbation. J. Diff. Equat. 150 (2), 363–384. ISSN 0022-0396. http://dx.doi.org/10.1006/ jdeq.1998.3487. http://www.sciencedirect.com/science/article/pii/S0022039698934878. Gallouet, T., Gastaldo, L., Herbin, R., Latch
e, J.C., 2008, 3. An unconditionally stable pressure correction scheme for the compressible barotropic Navier-Stokes equations. ESAIM: Math. Model. Numer. Anal. 42 (2), 303–331. http://eudml.org/doc/250414. Grenier, E., 1997. Oscillatory perturbations of the Navier-Stokes equations. J. Math. Pures. Appl. 76, 477–498. Guillard, H., 2009. On the behavior of upwind schemes in the low Mach number limit. IV: P0 approximation on triangular and tetrahedral cells. Comput. Fluids 38 (10), 1969–1972. ISSN 0045-7930. http://dx.doi.org/10.1016/j.compfluid.2009.06.003. http://www.sciencedirect. com/science/article/pii/S0045793009000838. Guillard, H., Murrone, A., 2004. On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes. Comput. Fluids 33 (4), 655–675. Guillard, H., Viozat, C., 1999. On the behavior of upwind schemes in the low Mach number limit. Comput. Fluids 28, 63–96. Gustafsson, B., 1908. Asymptotic expansions for hyperbolic problems with different time scales. SIAM J. Numer. Anal. 17, 623–634. Harlow, F.H., Amsden, A.A., 1971. A numerical fluid dynamics calculation method for all flow speeds. J. Comput. Phys. 8 (2), 197–213. ISSN 0021-9991. http://dx.doi.org/10.1016/00219991(71)90002-7. http://www.sciencedirect.com/science/article/pii/0021999171900027. Helbing, D., 1992. A fluid-dynamic model for the movement of pedestrians. Complex Syst. 6, 391–415. Hughes, T.J.R., Franca, L.P., Balestra, M., 1986. A new finite element formulation for computational fluid dynamics: V. Circumventing the Babusˇka-Brezzi condition: a stable Petrov-Galerkin formulation of the stokes problem accommodating equal-order interpolations. Comput. Methods Appl. Mech. Eng. 59 (1), 85–99. ISSN 0045-7825. http://dx.doi.org/10.1016/0045-7825(86)90025-3. http://www.sciencedirect.com/science/article/pii/0045782586900253. Klainerman, S., Majda, A., 1981. Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids. Commun. Pure Appl. Math. 34 (4), 481–524. ISSN 1097-0312. http://dx.doi.org/10.1002/cpa. 3160340405. Klainerman, S., Majda, A., 1982. Compressible and incompressible fluids. Commun. Pure Appl. Math. 35, 629–653.

230 Handbook of Numerical Analysis Klein, R., 1995. Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: one-dimensional flow. J. Comput. Phys. 121, 213–237. Kreiss, H.O., 1980. Problems with different time scales for partial differential equations. Commun. Pure Appl. Math. 33, 399–440. Kreiss, H.O., Lorenz, J., Naughton, M.J., 1991. Convergence of solutions of the compressible to the solutions of the incompressible Navier-Stokes equations. Adv. Appl. Math. 12, 187–214. Li, X.S., Gu, C.W., 2008. An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour. J. Comput. Phys. 227 (10), 5144–5159. ISSN 0021-9991. http://dx.doi.org/ 10.1016/j.jcp.2008.01.037. http://www.sciencedirect.com/science/article/pii/S0021999108000697. Li, X.S., Gu, C.W., 2013. Mechanism of Roe-type schemes for all-speed flows and its application. Comput. Fluids 86, 56–70. ISSN 0045-7930. http://dx.doi.org/10.1016/j.compfluid.2013.07.004. http://www.sciencedirect.com/science/article/pii/S0045793013002752. Li, X.S., Gu, C.W., Xu, J.Z., 2009. Development of Roe-type scheme for all-speed flows based on preconditioning method. Comput. Fluids 38 (4), 810–817. ISSN 0045-7930. http://dx.doi.org/10.1016/j.compfluid.2008.08.002. http://www.sciencedirect.com/science/ article/pii/S0045793008001928. Lions, J.L., 1993. Limite incompressible et acoustique pour les fluides visqueux, compressibles et isentropiques. CRAS 317, 1197–1220. Serie I. Lions, P.L., Masmoudi, N., 1998. Incompressible limit for a viscous compressible fluid. J. Math. Pures Appl. 77 (6), 585–627. ISSN 0021-7824. http://dx.doi.org/10.1016/S0021-7824(98) 80139-6. http://www.sciencedirect.com/science/article/pii/S0021782498801396. Majda, A., 2012. Compressible Fluid Flow and Systems of Conservation Laws in Several Space Variables. Applied Mathematical Sciences, Springer, New York, NY. https://books.google. fr/books?id¼U4zfBwAAQBAJ. ISBN 9781461211167. Mary, I., Sagaut, P., 2002. Large Eddy simulation of flow around an airfoil near stall. AIAA J. 40 (6), 1139–1145. M
etivier, G., Schochet, S., 2001. The incompressible limit of the non-isentropic Euler equations. Arch. Ration. Mech. Anal. 158 (1), 61–90. ISSN 1432-0673. http://dx.doi.org/10.1007/ PL00004241. Metivier, G., Schochet, S., 2001. Limite incompressible des
equations d’Euler non-isentropiques. http://www.maths.univ-rennes1.fr/metivier/preprints.html. S
eminaire
equations aux d
eriv
ees partielles de l’
ecole Polytechnique, to be found at. Munz, C.D., Roller, S., Klein, R., Geratz, K.J., 2003. The extension of incompressible flow solvers to the weakly compressible regime. Comput. Fluids 32 (2), 173–196. ISSN 00457930. http://dx.doi.org/10.1016/S0045-7930(02)00010-5. http://www.sciencedirect.com/science/ article/pii/S0045793002000105. Oswald, K., Siegmund, A., Birken, P., Hannemann, V., Meister, A., 2016. L2Roe: a low dissipation version of Roe’s approximate Riemann solver for low Mach numbers. Int. J. Numer. Methods Fluids 81 (2), 71–86. ISSN 1097-0363. http://dx.doi.org/10.1002/fld.4175. Fld.4175. Rhie, C.M., Chow, W.L., 1983. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA J. 21, 1525–1532. Rieper, F., 2011. A low-Mach number fix for Roe’s approximate Riemann solver. J. Comput. Phys. 230 (13), 5263–5287. ISSN 0021-9991. http://dx.doi.org/10.1016/j.jcp.2011.03.025. http://www.sciencedirect.com/science/article/pii/S0021999111001689. Rieper, F., Bader, G., 2009. The influence of cell geometry on the accuracy of upwind schemes in the low Mach number regime. J. Comput. Phys. 228 (8), 2918–2933. ISSN 0021-9991. http://dx.doi.org/10.1016/j.jcp.2009.01.002. http://www.sciencedirect.com/ science/article/pii/S0021999109000096.

On the Behaviour of Upwind Schemes Chapter

8 231

Schochet, S., 1988. Asymptotics for symmetric hyperbolic systems with a large parameter. J. Diff. Equat. 75, 1–27. Schochet, S., 1994. Fast singular limits of hyperbolic PDEs. J. Diff. Equat. 114, 476–512. Schochet, S., 2005. The mathematical theory of low Mach number flows. ESAIM: M2AN 39 (3), 441–458. http://dx.doi.org/10.1051/m2an:2005017. Thornber, B., Drikakis, D., Williams, R.J.R., Youngs, D., 2008a. On entropy generation and dissipation of kinetic energy in high-resolution shock-capturing schemes. J. Comput. Phys. 227, 4853–4872. Thornber, B., Mosedale, A., Drikakis, D., Youngs, D., Williams, R.J.R., 2008b. An improved reconstruction method for compressible flows with low Mach number features. J. Comput. Phys. 227 (10), 4873–4894. ISSN 0021-9991. http://dx.doi.org/10.1016/j.jcp.2008.01.036. http://www.sciencedirect.com/science/article/pii/S0021999108000429. Turkel, E., 1993. Review of preconditioning methods for fluid dynamics. Appl. Numer. Math. 12, 257–284. Turkel, E., 1999. Preconditional techniques in computational fluid dynamics. Annu. Rev. Fluid Mech. 31, 385. Turkel, E., Fiterman, A., van Leer, B., 1994. Preconditioning and the limit of the compressible to the incompressible flow equations for finite difference schemes. In: Caughey, D.A., Hafez, M.M. (Eds.), Frontiers of Computational Fluid Dynamics. John Wiley and Sons, Chichester, pp. 215–234. Turkel, E., Radespiel, R., Kroll, N., 1997. Assessment of preconditioning methods for multidimensional aerodynamics. Comput. Fluids 26, 613–634. Ukai, S., 1986. The incompressible limit and the initial layer of the compressible Euler equations. J. Math. Kyoto Univ. 26, 323–331. van der Heul, D.R., Vuik, C., Wesseling, P., 2003. A conservative pressure-correction method for flow at all speeds. Comput. Fluids 32 (8), 1113–1132. ISSN 0045-7930. http://dx.doi.org/10.1016/ S0045-7930(02)00086-5. http://www.sciencedirect.com/science/article/pii/S0045793002000865. Volpe, G., 1993. Performance of compressible flow codes at low Mach number. AIAA J. 31, 49–56. Weiss, J.M., Smith, W.A., 1995. Preconditioning applied to variable and constant density flows. AIAA J. 33 (11), 2050–2057. Zienkiewicz, O.C., Szmelter, J., Perraire, J., 1990. Compressible and incompressible flows: an algorithm for all seasons. Comput. Methods Appl. Mech. Eng. 78 (1), 105–121.