Jean-François Coulombel
email: jean-francois.coulombel@math.univ-toulouse.fr.

Tomas Lundquist
email: tomas.lundquist@math.univ-toulouse.fr.

On some stable boundary closures of finite difference schemes for the transport equation

Keywords: AMS classification: 65M12, 65M06, 65M20 transport equation, numerical schemes, extrapolation boundary condition, energy, stability

We explore in this article the possibilities and limitations of the so-called energy method for analyzing the stability of finite difference approximations to the transport equation with extrapolation numerical boundary conditions at the outflow boundary. We first show that for the most simple schemes, namely the explicit schemes with a three point stencil, the energy method can be applied for proving stability estimates when the scheme is implemented with either the first or second order extrapolation boundary condition. We then examine the case of five point stencils and give several examples of schemes and second order extrapolation numerical boundary conditions for which the energy method produces stability estimates. However, we also show that for the standard first or second order translatory extrapolation boundary conditions, the energy method cannot be applied for proving stability of the classical fourth order scheme originally proposed by Strang. This gives a clear limitation of the energy method with respect to the more general approach based on the normal mode decomposition.

Introduction

A general approach for studying the stability of numerical boundary conditions for discretized hyperbolic equations has been initiated in the fundamental contribution [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF]. However, for technical reasons, the stability estimates in [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF] are restricted to zero initial data and can be obtained only after verifying the fulfillment of some "algebraic" condition which is commonly referred to as the Uniform Kreiss-Lopatinskii Condition. The theory in [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF] has been successfully applied to some well-known extrapolation procedures at outflow boundaries, see a preliminary announcement in [START_REF] Kreiss | Difference approximations for hyperbolic differential equations[END_REF] and the complete proof in [START_REF] Goldberg | On a boundary extrapolation theorem by Kreiss[END_REF]. In the recent work [START_REF] Coulombel | The neumann numerical boundary condition for transport equations[END_REF], one of the authors has revisited the stability estimates for the outflow extrapolation procedures in [START_REF] Kreiss | Difference approximations for hyperbolic differential equations[END_REF][START_REF] Goldberg | On a boundary extrapolation theorem by Kreiss[END_REF] and shown that a suitably devised energy argument could bypass the (technical and lenghty) arguments of [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF]. In this article, we examine in a systematic way the possible applications but also the limitations of the energy method for analyzing the stability of extrapolation procedures at an outflow boundary for the most simple onedimensional transport equation. We begin by reviewing some central concepts and notation that will be extensively used throughout the remainder of this paper.

Continuous problem

We consider the scalar advection equation in 1D on a semi-infinite interval with an outflow boundary at x = 0, u t + u x = 0, -∞ ≤ x ≤ 0 , t ≥ 0 .

(1)

The differential operator ∂ t + ∂ x is semi-bounded on this domain, i.e. for all 0 ≤ t 1 ≤ t 2 and integrable smooth functions Φ(x, t) where Φ(-∞, t) = 0, we have,

Φ t + Φ x = 0 ⇒ Φ 2 (t 2 ) = Φ 2 (t 1 ) - t 2 t 1 Φ(0, τ ) 2 dτ ≤ Φ 2 (t 1 ) , (2)
where we have used the standard definition of the L 2 inner product and norm (with respect to the space variable x) for real valued functions, (Φ, Ψ) := 0 -∞ Φ(x) Ψ(x) dx , Φ 2 := (Φ, Φ) .

When later discretizing (1) in time and space, we shall seek to mimic the energy balance (2).

Discretizing on the whole real line

We seek to discretize the problem (1) in time and space using a finite difference method. In the interior of the spatial domain we use a repeated interior stencil, which we analyze for the Cauchy version of (1).

We introduce an equispaced spatial grid over the whole real line,

x j := j ∆x, -∞ ≤ j ≤ ∞ ,
where ∆x > 0 is the grid spacing. Unless otherwise stated, sequences (Φ j ) on the grid are assumed to be real valued, and the norm on ℓ 2 (Z) is defined by,

Φ 2 ℓ 2 (Z) := j∈Z Φ 2 j .
Next, we introduce a time step ∆t > 0 (and accordingly, t n = n∆t) and consider the ratio z := ∆t/∆x (the so-called Courant-Friedrichs-Lewy parameter [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF]) as a constant. In what follows, u n j stands for an approximation of the solution u to the Cauchy version of (1) in the neighborhood of (x j , t n ). We consider explicit finite difference discretizations to (1) of the form:

u n+1 j = A(z) u n j , -∞ ≤ j ≤ ∞ , n ∈ N , (3)
where A(z) is a polynomial with respect to z = ∆t/∆x of difference operators. E.g, decomposing A(z) in terms of simple shift operators, A(z) then acts on sequences (Φ j ) indexed by j according to,

A(z) Φ j := ℓ + ℓ=-ℓ - a ℓ (z) Φ j+ℓ ,
where the integers ℓ ± mark the extent (or bandwidth) of the difference stencil, and where the coefficients a ℓ (z) are polynomial expressions with respect to z. Most of the examples we will consider below fall into the classes of three point stencils (ℓ -= ℓ + = 1) and five point stencils (ℓ -= ℓ + = 2). It is clear that A(z) acts boundedly on ℓ 2 (Z), and we shall say that A(z) is a contraction if the norm of A(z) as an operator on ℓ 2 (Z) is not larger than 1.

In what follows, we only consider finite difference schemes that are consistent with the transport equation (1). In other words, see for instance [START_REF] Richtmyer | Difference methods for initial-value problems[END_REF][START_REF] Gustafsson | Time dependent problems and difference methods[END_REF], we only consider operators A(z) in (3) that satisfy at least the first order accuracy relations, ℓ + ℓ=-ℓ - a ℓ (z) = 1 , and

ℓ + ℓ=-ℓ - ℓ a ℓ (z) = -z .
In order to prove stability estimates, be they on the whole space Z or on the half-space Z -, we aim to use the energy method in order to prove that the scheme is a contraction with respect to a certain norm. For the Cauchy problem, we thus need for A(z) to be a contraction on ℓ 2 (Z). Since the coefficients of the operator A(z) do not depend on j, we can achieve this by splitting each local energy balance term into a telescopic and a dissipative part. In other words, the goal will be to show that, for an arbitrary sequence

(Φ j ) j∈Z , Ψ j = A(z) Φ j ⇒ Ψ 2 j = Φ 2 j + T j -T j-1 + S j , S j ≤ 0 , (4)
for all values of z within a certain specified interval. In (4), T j , T j-1 and S j are quadratic quantities over the Φ k entries for j -ℓ -≤ k ≤ j + ℓ + . When summing over all whole numbers, the telescopic part cancels and we are left with

Ψ 2 ℓ 2 (Z) = Φ 2 ℓ 2 (Z) + j∈Z S j ≤ Φ 2 ℓ 2 (Z) ,
showing that A(z) is a contraction, just as the continuous problem (2).

The main goal of analyzing schemes of the form (3) will be to derive integration by parts decompositions of the form (4) with S j ≤ 0 for any j. There is a lot of freedom involved in splitting the difference Ψ 2 j -Φ 2 j into telescopic and non-telescopic parts as in (4), see for instance [START_REF] Coulombel | The neumann numerical boundary condition for transport equations[END_REF]Lemma A.1]. In this paper we will take advantage of this in order to expand and refine the energy method arguments used in [START_REF] Coulombel | The neumann numerical boundary condition for transport equations[END_REF].

Finite difference operators

The basic unit for constructing all finite difference stencils will be the backward difference operator,

D Φ j := Φ j -Φ j-1 . (5)
In order to define a parametric family of schemes with a three point stencil, we will use the standard first and second order derivative approximations (the "centered finite difference" and the "discrete Laplacian", respectively). On normalized form with respect to the spacing ∆x, we thus consider,

D 0 Φ j = 1 2 (D Φ j + D Φ j+1 ) = 1 2 (Φ j+1 -Φ j-1 ) , (6)
∆ Φ j = D Φ j+1 -D Φ j = Φ j-1 -2 Φ j + Φ j+1 . (7)
Note that both of these operators can easily be shown to be of accuracy order 2. To define a family of five point stencils, we also need to approximate the third and fourth derivatives, which we define using,

D 0 ∆ Φ j = 1 2 (-Φ j-2 + 2 Φ j-1 -2 Φ j+1 + Φ j+2 ) , (8)
∆ 2 Φ j = Φ j-2 -4 Φ j-1 + 6 Φ j -4 Φ j+1 + Φ j+2 . (9)
Since we will not consider schemes with more than a five point stencil, in this paper we will limit our attention to only the finite difference operators introduced above.

Outflow boundary conditions

We close the discrete domain at the outflow boundary x = 0, and define a general vector on the grid as,

Φ = . . . , Φ -2 , Φ -1 , Φ 0 .
At the points nearest to the outflow boundary at x = 0, extrapolation to ghost points outside the computational domain are needed in order to close the finite difference stencils and thus the numerical scheme (3) itself. Namely, we need a method to define the Φ 1 , . . . , Φ ℓ + values in order to compute Ψ j := A(z) Φ j for all j ∈ Z -. Note that this extrapolation step may or may not reduce the local accuracy order from that of the interior stencil itself. Next, we define an inner product on Z -in the form of a quadrature rule, such that in general only the last r quadrature weights h j > 0 are non-unity,

Φ 2 H := (Φ, Φ) H , (Φ, Ψ) H := j≤-r Φ j Ψ j + 0 j=-r+1 h j Φ j Ψ j . (10)
Now let's assume that contractivity (4) holds for all sequences (Φ j ) ∈ ℓ 2 (Z). Then, using an arbitrary set of extrapolation conditions to define the Φ j values with positive indices j = 1, . . . , ℓ + , we have,

∀ j ∈ Z -, Ψ j = A(z) Φ j ⇒ Ψ 2 H -Φ 2 H - j≤-r S j = 0 j=-r (h j -h j+1 ) T j + 0 j=-r+1 h j S j =: E , (11)
where we have defined h -r := 1 and h 1 := 0. Since the local dissipation S j is nonpositive for any j, we get from (11):

∀ j ∈ Z -, Ψ j = A(z) Φ j ⇒ Ψ 2 H ≤ Φ 2 H + E .
The goal at this point is to determine the extrapolation procedure and the quadrature weights such that E is, at least, non-positive. In that case, the combination of the operator A(z) with the extrapolation procedure yields a contraction on ℓ 2 (Z -) provided that the latter space is equipped with the norm (10). We formalize this notion below.

Definition 1. Together with suitable extrapolation conditions at an outflow boundary, we say that the operator A(z) in (3) is semi-bounded with respect to the inner product defined in (10) if, for a given z ≥ 0, the following two conditions are satified.

1. The scheme is contractive, that is A(z) admits a decomposition (4) with S j ≤ 0.

2. The boundary contribution is non-positive, i.e.

E := 0 j=-r (h j -h j+1 ) T j + 0 j=-r+1 h j S j ≤ 0 , (12)
where we use the convention h -r := 1 and h 1 := 0.

As a matter of fact, it is useful for the purpose of convergence proofs (see, e.g., [Gus75, GKO95, CL20]) to infer from (11) a trace estimate on the sequence Φ. This will be the case provided that E is a negative definite quadratic form of its arguments, as evidenced below on several examples. This energy argument bypasses, as in [Wu95, CG11, CL20], the lengthy verification of the fulfillment of the so-called Uniform Lopatinskii Condition. We refer to [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF][START_REF] Coulombel | Stability of finite difference schemes for hyperbolic initial boundary value problems[END_REF] for a general presentation of this topic.

The rest of this article is organized as follows. In Section 2, we examine the application of the energy method described above in the case of three point schemes (ℓ -= ℓ + = 1). Since such schemes have at most second order accuracy [START_REF] Gustafsson | Time dependent problems and difference methods[END_REF], we restrict our analysis to the case of first and second order extrapolation at the outflow boundary. We then examine the case of five point schemes (ℓ -= ℓ + = 2) in Section 3. Since the involved algebra becomes heavier, we restrict sometimes the analysis to three particular cases which are analogues of some well-known three point schemes. We give examples of extrapolation procedures for which we can construct a suitable energy that yields semi-boundedness. In contrast to what is more commonly used, these energy stable extrapolation procedures are not translation invariant. Finally, in section 5 we demonstrate the limitation of more standard extrapolation procedures, proving that for a fourth order five point scheme combined with the most common first and second order extrapolation procedures, the energy method actually fails to work, though the corresponding closure of the operator A(z) is known to be power-bounded. This shows that the normal mode analysis developed in [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF] and subsequent works is sometimes necessary to capture stability in general.

Three point schemes

The family of three point stencils of at least first order accuracy can be defined using a single free parameter ν as follows (see [START_REF] Richtmyer | Difference methods for initial-value problems[END_REF][START_REF] Gustafsson | Time dependent problems and difference methods[END_REF]),

A(z) = I -z D 0 + ν 2 ∆ . (13)
Famous examples include the Euler forward (ν = 0), Lax-Friedrichs (ν = 1), the two-point upwind (ν = z) as well as the second order accurate Lax-Wendroff (ν = z 2 ) schemes.

Contractivity on the whole real line

For the family of schemes considered in (13), our first result is the following.

Lemma 1. Let z, ν ∈ R, and consider the scheme (3) with A(z) as in (13). Then (4) holds with

T j := Φ j D Φ j+1 T Q Φ j D Φ j+1 , Q :=    -z ν -z 2 ν -z 2 ν (1 -z) 2    , (14)
and

S j = D Φ j D Φ j+1 T M D Φ j D Φ j+1 , M := a 1 1 0 0 1 + a 2 1 -1 -1 1 , (15)
where

a 1 := - ν -z 2 2 , a 2 := ν 2 -z 2 4 .
In particular, A(z) is a contraction on ℓ 2 (Z), i.e. S j ≤ 0, if z and ν satisfy the relation z 2 ≤ ν ≤ 1.

Proof. With A(z) as in (13), we have in (4),

Ψ 2 j -Φ 2 j = -2 z Φ j (D 0 Φ j ) + z 2 (D 0 Φ j ) 2 + ν Φ j (∆ Φ j ) -z ν (D 0 Φ j ) (∆ Φ j ) + ν 2 4 (∆ Φ j ) 2 . ( 16
)
In order to split this expression into telescopic and non-telescopic parts, we first note that the (D 0 Φ j ) 2 term can be expressed as the sum of squares (this is nothing but the parallelogram identity),

(D 0 Φ j ) 2 = 1 4 2 (D Φ j ) 2 + 2 (D Φ j+1 ) 2 -(∆ Φ j ) 2 . ( 17
)
Moreover, we will use the purely telescopic formulas,

2 Φ j D 0 Φ j = Φ j Φ j+1 -Φ j-1 Φ j = Φ 2 j + Φ j D Φ j+1 -Φ 2 j-1 + Φ j-1 D Φ j , (18)
2 D 0 Φ j ∆ Φ j = (D Φ j+1 ) 2 -(D Φ j ) 2 . (19)
Finally, the Φ j ∆ Φ j term can be written as the combination of non-telescopic (sum of squares) and telescopic parts 1 ,

2 Φ j ∆ Φ j = -(D Φ j ) 2 -(D Φ j+1 ) 2 + Φ 2 j+1 -2 Φ 2 j + Φ 2 j-1 = -(D Φ j ) 2 -(D Φ j+1 ) 2 + 2 Φ j D Φ j+1 + (D Φ j+1 ) 2 -2 Φ j-1 D Φ j + (D Φ j ) 2 . ( 20
)
Inserting the formulas (17) through (20) into (16), we get the desired form Ψ 2 j -Φ 2 j = T j -T j-1 + S j with T j directly as in (14) and S j given by

S j = a 1 (D Φ j ) 2 + (D Φ j+1 ) 2 + a 2 (∆ Φ j ) 2 .
From the definition of the discrete Laplacian ∆ in (7), we can further rewrite this into (15). This proves the first part of the Lemma.

Next, an orthogonal set of eigenvectors to M is given by (1 1) T and (1 -1) T , and the associated eigenvalues are,

λ 1 = a 1 = - ν -z 2 2 , λ 2 = a 1 + 2 a 2 = - ν (1 -ν) 2 .
We have,

λ 1 ≤ 0 ⇔ z 2 ≤ ν , λ 2 ≤ 0 ⇔ ν ∈ [0, 1] .
and we have thus shown that the scheme defined by A(z) is a contraction if z 2 ≤ ν ≤ 1.

1 The decomposition mimics the equality 2 u u ′′ = -2

(u ′ ) 2 + (u 2 ) ′′ .
It is worth to mention here that the same exact bounds on z, ν for ℓ 2 stability given in Lemma 1 above can also be derived using Fourier analysis (they are both necessary and sufficient for A(z) to be a contraction). However, for a wider stencil we can not in general guarantee that the energy method will yield the same bounds as obtained with the Fourier method. Since we will use the energy method to derive extrapolation boundary conditions leading to semi-boundedness, we shall in this paper only consider stencils for which contractivity can be shown using the energy method in a similar fashion to the one above.

To conclude this section, we write down the eigenvalues associated with the quadratic form S j for the important special cases mentioned at the beginning of this section. The Euler forward method (ν = 0) yields,

λ 1 = z 2 2 , λ 2 = 0 ,
showing that this scheme is not contractive for any z. Moreover, the Lax-Friedrichs scheme (ν = 1) yields,

λ 1 = - 1 -z 2 2 , λ 2 = 0 ,
the upwind scheme (ν = z) yields,

λ 1 = λ 2 = - z (1 -z) 2 ,
and finally the Lax-Wendroff (ν = z 2 ) scheme,

λ 1 = 0 , λ 2 = - z 2 (1 -z 2 ) 4 .
This shows that the Lax-Friedrichs and Lax-Wendroff schemes are both contractive for all |z| ≤ 1 (meaning that we can freely change the sign in front of u x in (1) without modifying the stability property of the scheme), while the upwind scheme is contractive for 0 ≤ z ≤ 1.

Outflow boundary conditions

For simplicity, in (10) we let r = 1 from the start, i.e. the norm of Φ ∈ ℓ 2 (Z -) in ( 10) is defined by

Φ 2 H = j≤-1 Φ 2 j + h 0 Φ 2 0 ,
and for now we consider the last quadrature weight h 0 as a free (positive) parameter. In (12) we then have, using Lemma 1 from the previous paragraph,

E = (1 -h 0 ) T -1 + h 0 T 0 + h 0 S 0 = (1 -h 0 ) Φ -1 D Φ 0 T Q Φ -1 D Φ 0 + h 0 Φ 0 D Φ 1 T Q Φ 0 D Φ 1 + h 0 D Φ 0 D Φ 1 T M D Φ 0 D Φ 1 . (21)
Inserting ( 14) and (15) into the above expression (21) for E gives us right away a stability estimate for the most popular first order extrapolation boundary closure.

Corollary 1. Let z ∈ (0, 1] and let ν satisfy z 2 ≤ ν ≤ 1. Then the scheme (13) together with the first order extrapolation condition Φ 1 = Φ 0 is semi-bounded. In particular, the energy balance is given by,

j≤0 Ψ 2 j ≤ j≤0 Φ 2 j -z Φ 2 0 ≤ j≤0 Φ 2 j , (22)
which exactly mimics (2) with added dissipation.

Proof. We choose h 0 = 1 for the last quadrature weight. From the extrapolation condition we can use

D Φ 1 = 0 in (21) to immediately get, E = Φ 0 0 T Q Φ 0 0 + D Φ 0 0 T M D Φ 0 0 = -z Φ 2 0 + S 0 .
The energy balance (22) then follows from (11) and the fact that S j ≤ 0 for any j ≤ 0 if z 2 ≤ ν ≤ 1.

We now examine the case of a second order extrapolation condition at the outflow boundary, that is we impose ∆ Φ 0 = 0. Our result is the following.

Corollary 2. Let z ∈ (0, 1] and let ν satisfy z 2 ≤ ν ≤ 1, and consider the scheme (13) together with the second order extrapolation condition ∆ Φ 0 = 0. Then E in ( 12) is given by

E = Φ 0 D Φ 0 T B Φ 0 D Φ 0 , B = -z z( 1 2 -h 0 ) + ν 2 z( 1 2 -h 0 ) + ν 2 h 0 z 2 -ν(1+z) 2 . ( 23
)
Moreover, B is negative definite (and thus the scheme is semi-bounded) for the choice h 0 = (1-z +ν/z)/2.

Proof. From the second order extrapolation condition, we have,

∆ Φ 0 = 0 ⇒ D Φ 1 = D Φ 0 .
We also have, from the definition of D in (5),

Φ -1 = Φ 0 -D Φ 0 .
Going back to (21), we now write E as,

E = (1 -h 0 ) Φ 0 -D Φ 0 D Φ 0 T Q Φ 0 -D Φ 0 D Φ 0 + h 0 Φ 0 D Φ 0 T Q Φ 0 D Φ 0 + h 0 D Φ 0 D Φ 0 T M D Φ 0 D Φ 0 ,
which, using ( 14) and (15), can be further simplified into (23).

In order to see whether B in (23) is negative definite, we rotate to diagonal form using the similarity transformation,

R :=   1 ν + z 2 -z h 0 /z 0 1   ⇒ R T B R = -z 0 0 γ , where γ := (ν + z -2 z h 0 ) 2 + 4 h 0 z 3 -2 z (1 + z) ν /(4 z) .
The minimum of γ with respect to h 0 is obtained for

∂γ ∂h 0 = z 2 -z -ν + 2 h 0 z = 0 ⇔ h 0 = ν z -z + 1 2 ≥ 1 2 ,
which in all cases leads to the same negative value of γ,

γ = -z 3 /4 .
With this choice of h 0 , the quantity (10) (with r = 1) defines a norm on ℓ 2 (-∞, 0) and the second condition of Definition 1 is always satisfied for z ∈ (0, 1] and ν ∈ [z 2 , 1]. Semi-boundedness is thus obtained as long as the interior scheme is contractive.

In particular, the second order accurate Lax-Wendroff scheme (ν = z 2 ) is associated with the likewise second order accurate Trapezoidal rule with the constant value h 0 = 1/2 (independently of z ∈ (0, 1]). For all other three-point schemes, the above value for h 0 as obtained with this method depends on z and grows to infinity as z approaches 0 (ν being fixed).

Five point schemes

In this preliminary work we do not aim for a complete theory encompassing finite difference schemes of any order. Instead we go on to consider a slightly more involved example than in the previous section, namely the family of schemes with a five point stencil and at least second order accuracy. Such schemes can be parametrized with two free parameters σ, τ as follows,

A(z) = I -z D 0 + z 2 2 ∆ + σ D 0 ∆ + τ ∆ 2 . ( 24
)
We first derive a general integration by parts decomposition formula for A(z) and thereby obtain three interesting examples for which A(z) is contractive.

Contractivity on the whole real line

Generalizing Lemma 1 to the family of operators in (24), our result is the following.

Lemma 2. Let z, σ, τ ∈ R, and consider the scheme (3) with A(z) as in (24). Then (4) holds with

T j :=     Φ j D Φ j+1 ∆ Φ j D ∆ Φ j+1     T Q     Φ j D Φ j+1 ∆ Φ j D ∆ Φ j+1     , Q :=      -z -z (1-z) 2 σ σ 2 + τ -z (1-z) 2 z 2 (1-z) 2 -σ σ 2 -z σ -τ -z σ 2 + τ σ σ 2 -z σ -τ z 2 σ 2 + z τ z (1+z) 2 σ 2 + τ -a 1 3 σ 2 + τ -z σ 2 + τ z (1+z) 2 σ 2 + τ -a 1 3 z σ 4 + z 2 τ 2 + σ τ -a 1 3      , (25)
and

S j :=   ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   T M   ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   , M := a 1 3   1 0 0 0 1 0 0 0 1   + a 2   1 -1 0 -1 2 -1 0 -1 1   + a 3   1 -2 1 -2 4 -2 1 -2 1   , (26)
where

a 1 := z 2 (z 2 -1) 4 + 2 z σ + 2 τ , a 2 := - z σ 4 + σ 2 2 - z 2 τ 2 , a 3 := - σ 2 4 + τ 2 .
In particular, A(z) is a contraction on ℓ 2 (Z) in each of the three following cases:

• σ = 0, τ = -(1 -z 2 )/12 and z ∈ [-1, 1], • σ = z (1 -z 2 )/6, τ = -σ/2 and z ∈ [0, 1], • σ = z (1 -z 2 )/6, τ = -z 2 (1 -z 2 )/24 and z ∈ [-1, 1].
Proof. Using (24) in (4), we now have,

Ψ 2 j -Φ 2 j = -2 z Φ j D 0 Φ j + z 2 (D 0 Φ j ) 2 + z 2 Φ j ∆ Φ j -z 3 (D 0 Φ j ) (∆ Φ j ) + z 4 4 (∆ Φ j ) 2 + 2 σ Φ j D 0 ∆ Φ j -2 z σ (D 0 Φ j ) (D 0 ∆ Φ j ) + z 2 σ (∆ Φ j ) (D 0 ∆ Φ j ) + σ 2 (D 0 ∆ Φ j ) 2 + 2 τ Φ j (∆ 2 Φ j ) -2 z τ (D 0 Φ j ) (∆ 2 Φ j ) + z 2 τ (∆ Φ j ) (∆ 2 Φ j ) + 2 σ τ (D 0 ∆ Φ j ) (∆ 2 Φ j ) + τ 2 (∆ 2 Φ j ) 2 . ( 27
)
The terms in the first line on the right hand side above are covered by the previous formulas (17), ( 18), ( 19) and (20). Further, by substituting Φ j for ∆ Φ j in the same formulas, we get,

(D 0 ∆ Φ j ) 2 = 1 4 2 (D ∆ Φ j ) 2 + 2 (D ∆ Φ j+1 ) 2 -(∆ 2 Φ j ) 2 , ( 28
)
as well as

2 (∆ Φ j ) (D 0 ∆ Φ j ) = (∆ Φ j ) 2 + ∆ Φ j D ∆ Φ j+1 -(∆ Φ j-1 ) 2 + ∆ Φ j-1 D ∆ Φ j , (29) 2
(D 0 ∆ Φ j ) (∆ 2 Φ j ) = (D ∆ Φ j+1 ) 2 -(D ∆ Φ j ) 2 , ( 30
) and 2 ∆ Φ j ∆ 2 Φ j = -(D ∆ Φ j ) 2 -(D ∆ Φ j+1 ) 2 + 2 ∆ Φ j D ∆ Φ j+1 + (D ∆ Φ j+1 ) 2 -2 ∆ Φ j-1 D ∆ Φ j + (D ∆ Φ j ) 2 . ( 31
)
For the remaining terms on the right hand side of (27), we need to introduce the following four new high order formulas, the proof of which can be found in Appendix,

2 D 0 Φ j D 0 ∆ Φ j = -2 (∆ Φ j ) 2 + 1 4 (D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2 + (D Φ j+1 - 1 4 D ∆ Φ j+1 ) (D ∆ Φ j+1 + 2 ∆ Φ j ) -(D Φ j - 1 4 D ∆ Φ j ) (D ∆ Φ j + 2 ∆ Φ j-1 ) , (32)
as well as,

2 Φ j D 0 ∆ Φ j = [Φ j (D ∆ Φ j+1 + 2 ∆ Φ j ) + (∆ Φ j -D Φ j+1 ) D Φ j+1 ] (33) -[Φ j-1 (D ∆ Φ j + 2 ∆ Φ j-1 ) + (∆ Φ j-1 -D Φ j ) D Φ j ] , 2 D 0 Φ j ∆ 2 Φ j = [-(D ∆ Φ j+1 + ∆ Φ j ) ∆ Φ j + 2 D Φ j+1 D ∆ Φ j+1 ] (34) -[-(D ∆ Φ j + ∆ Φ j-1 ) ∆ Φ j-1 + 2 D Φ j D ∆ Φ j ] ,
and

2 Φ j ∆ 2 Φ j = 2 (∆ Φ j ) 2 + [2 Φ j D ∆ Φ j+1 -2 D Φ j+1 ∆ Φ j ] -[2 Φ j-1 D ∆ Φ j -2 D Φ j ∆ Φ j-1 ] . ( 35
)
Note that all the non-telescopic parts above are expressed in terms of the three variables (∆ Φ j ) 2 , (∆ 2 Φ j ) 2 and (D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2 . The last two of these can be written as,

(∆ 2 Φ j ) 2 =   ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   T   1 -2 1   1 -2 1   ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   =   ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   T   1 -2 1 -2 4 -2 1 -2 1     ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   ,
and

(D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2 =   ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   T   1 -1 0 -1 2 -1 0 -1 1     ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   .
Thus, in order to analyze the sign of the complete non-telescopic part, it is convenient to expand (∆ Φ j ) 2 into a similar matrix form with a telescopic correction,

(∆ Φ j ) 2 = 1 3   ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   T   1 0 0 0 1 0 0 0 1     ∆ Φ j-1 ∆ Φ j ∆ Φ j+1   - 1 3 (D∆ Φ j + ∆ Φ j ) 2 + (∆ Φ j ) 2 + 1 3 (D∆ Φ j-1 + ∆ Φ j-1 ) 2 + (∆ Φ j-1 ) 2 .
Inserting (28) through (35) into (27) and then rewriting (∆ Φ j ) 2 , (∆ 2 Φ j ) 2 and (D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2 according to the above, the first part of the Lemma follows. Note that the term a 1 (∆ Φ j ) 2 has been split between between the telescopic and dissipative parts.

The same orthogonal set of eigenvectors to each of the three terms in (26) is given by (1, 1, 1) T , (-1, 0, 1) T and (1, -2, 1) T , with the associated eigenvalues of M ,

λ 1 := a 1 3 , λ 2 := a 1 3 + a 2 , λ 3 := a 1 3 + 3 a 2 + 6 a 3 .
The operator A(z) is a contraction on ℓ 2 (Z) if all three eigenvalues λ 1,2,3 of M are non-positive. In contrast to the three point stencil case, we can not easily rewrite the conditions for contractivity into explicit relations for σ and τ . However, three special cases seem to be worth to mention. First, if we let σ = 0 and τ = -(1 -z 2 )/12, then the eigenvalues of M simplify into,

λ 1 = - (1 -z 2 ) (2 + 3 z 2 ) 36 , λ 2 = - (1 -z 2 ) (4 + 3 z 2 ) 72 , λ 3 = - 1 -z 2 72 .
Hence A(z) is a contraction for any z ∈ [-1, 1] in that case. Recall that for the second order Lax-Wendroff scheme, the dissipation approaches 0 as z approaches 0, which is not the case here. In terms of dissipation properties, this new scheme can instead be seen in some sense to be analogous to the three-point Lax-Friedrichs scheme (though one eigenvalue for the Lax-Friedrichs scheme is zero, which is not the case here).

A skewed 4-point stencil of accuracy order 3, which can be seen as analogous to the two-point upwind scheme, is given by σ := z (1 -z 2 )/6 and τ := -σ/2 = -z(1 -z 2 )/12, yielding,

λ 1 = - z (1 -z 2 ) (2 -z) 36 , λ 2 = - z (1 -z 2 ) (1 + z) (2 -z) 2 72 , λ 3 = - z (1 -z 2 ) (2 -z) (2 + 3 z -3 z 2 ) 72 .
Again, A(z) is a contraction on ℓ 2 (Z) for any z ∈ [0, 1]. Finally, the so-called Strang scheme of fourth order accuracy [START_REF] Strang | Trigonometric polynomials and difference methods of maximum accuracy[END_REF] is given by σ := z (1 -z 2 )/6 and τ := -z 2 (1 -z 2 )/24, with the corresponding eigenvalues,

λ 1 = 0 , λ 2 = - z 2 (1 -z 2 ) (4 -z 2 ) 144 , λ 3 = - z 2 (1 -z 2 ) (3 -z 2 ) (4 -z 2 ) 96 .
The operator A(z) is again a contraction on ℓ 2 (Z) for any z ∈ [-1, 1]. This scheme is a high (i.e., fourth) order analogue to the second order accurate Lax-Wendroff scheme.

We note that the corresponding stability analysis for the fourth order Strang scheme in [START_REF] Strang | Trigonometric polynomials and difference methods of maximum accuracy[END_REF] is performed by means of the Fourier transform and is therefore not applicable for the outflow problem, which will be considered next.

Outflow boundary conditions

We examine in this paragraph two sets of outflow boundary conditions based on second order accurate extrapolation, and thus leave higher order boundary conditions as a topic for future work. As opposed to the three point stencil case, we have not found a single set of second order boundary conditions which is stable for all of the three example schemes listed in Lemma 2, which is why we consider two different alternatives below.

In what follows, we shall use the inner product (10) with r = 1 and h 0 = 1/2. In (12), we thus have,

E = 1 2     Φ -1 D Φ 0 ∆ Φ -1 D ∆ Φ 0     T Q     Φ -1 D Φ 0 ∆ Φ -1 D ∆ Φ 0     + 1 2     Φ 0 D Φ 1 ∆ Φ 0 D ∆ Φ 1     T Q     Φ 0 D Φ 1 ∆ Φ 0 D ∆ Φ 1     + 1 2   ∆ Φ -1 ∆ Φ 0 ∆ Φ 1   T M   ∆ Φ -1 ∆ Φ 0 ∆ Φ 1   , (36)
with Q given in (25) and M given in (26).

Second order extrapolation of type 1

We can prove Proposition 1. Let r = 1, h 0 = 1/2, and consider the set of second order extrapolation conditions ∆ Φ 0 = D 0 ∆ Φ 0 = 0. In (36), we then have,

E =   Φ 0 D Φ 0 D 2 Φ 0   T B   Φ 0 D Φ 0 D 2 Φ 0   , B :=        -z z 2 2 σ 2 z 2 2 - z 3 2 -σ -zσ 2 σ 2 -zσ 2 2 τ 3 + σ τ + 5 z σ 12 - z 2 τ 2 + 2 τ 2 - z 2 (1 -z 2 ) 12       
.

(37) In particular, the scheme is semi-bounded (i.e. B is negative semi-definite) for all z ∈ (0, 1) at least in the following two special cases,

• σ = 0 and τ = -(1 -z 2 )/12, • σ = z (1 -z 2 )/6 and τ = -σ/2.
Proof. The two extrapolation conditions readily yield,

D Φ 1 = D Φ 0 ∆ Φ 1 = ∆ Φ -1 .
By combining the two we also get D ∆ Φ 1 = ∆ Φ -1 , and by definition we also have Φ -1 = Φ 0 -D Φ 0 as well as ∆ Φ -1 = D 2 Φ 0 . With these formulas, we can now simplify E in (36) into,

E = 1 2     Φ 0 -D Φ 0 D Φ 0 D 2 Φ 0 -D 2 Φ 0     T Q     Φ 0 -D Φ 0 D Φ 0 D 2 Φ 0 -D 2 Φ 0     + 1 2     Φ 0 D Φ 0 0 D 2 Φ 0     T Q     Φ 0 D Φ 0 0 D 2 Φ 0     + 1 2   D 2 Φ 0 0 D 2 Φ 0   T M   D 2 Φ 0 0 D 2 Φ 0   ,
which in turn, after some straightforward algebra, leads to (37).

For the semi-boundedness part of the Proposition, let us first consider the case σ = 0 and τ = -(1 -z 2 )/12. The general expression for the matrix B in (37) reduces to,

B :=      -z z 2 /2 0 z 2 /2 - z 3 2 0 0 0 - (1 -z 2 ) (3 + 4 z 2 ) 72      .
It is a simple exercise to verify that the upper left 2 × 2 block is negative definite for z ∈ (0, 1), hence B is negative definite for z ∈ (0, 1). We get in that case the conclusion of Proposition 1 as an immediate consequence. We now focus on the case σ = z (1 -z 2 )/6 and τ = -σ/2, for which the matrix B reads:

B = -z      1 -z/2 -(1 -z 2 )/12 -z/2 1 + 2 z 2 6 z (1 -z 2 )/12 -(1 -z 2 )/12 z (1 -z 2 )/12 (1 -z 2 ) (1 + z) (4 -3 z) 72      .
We then compute:

  U V W   T B   U V W   = -z U - z 2 V - 1 -z 2 12 W 2 + z 2 4 V 2 + (1 -z 2 ) 6 V + z 4 W 2 + 14 + 4 z -8 z 2 48 W 2 ,
thereby showing that B is negative semi-definite for z ∈ (0, 1). This completes the proof of Proposition 1.

Outflow boundary: second order extrapolation of type 2

We examine in this paragraph a second set of outflow boundary conditions with second order accuracy.

Proposition 2. Let r = 1, h 0 = 1/2, and consider the set of second order extrapolation conditions ∆ Φ 0 = ∆ 2 Φ 0 = 0. Then we have

E =   Φ 0 D Φ 0 D 2 Φ 0   T B   Φ 0 D Φ 0 D 2 Φ 0   , B :=    -z z 2 /2 -τ z 2 /2 -z 3 2 -σ z τ -τ z τ 2 τ 3 + σ τ + 5 z σ 12 -z 2 τ 2 + σ 2 2 -z 2 (1-z 2 ) 12    .
(38) In particular, the scheme is semi-bounded (i.e. B is negative semi-definite) for all z ∈ (0, 1) at least in the following two special cases,

• σ = z (1 -z 2 )/6 and τ = -σ/2, • σ = z (1 -z 2 )/6 and τ = -z 2 (1 -z 2 )/24.
Proof. We use (36) and insert the second order extrapolation conditions ∆ Φ 0 = ∆ 2 Φ 0 = 0, from which we deduce the relations D Φ 1 = D Φ 0 , and D ∆ Φ 1 = ∆ Φ 1 = -∆ Φ -1 . The boundary contribution E in the energy balance (11) reads

E = 1 2     Φ 0 -D Φ 0 D Φ 0 D 2 Φ 0 -D 2 Φ 0     T Q     Φ 0 -D Φ 0 D Φ 0 D 2 Φ 0 -D 2 Φ 0     + 1 2     Φ 0 D Φ 0 0 -D 2 Φ 0     T Q     Φ 0 D Φ 0 0 -D 2 Φ 0     + 1 2   D 2 Φ 0 0 -D 2 Φ 0   T M   D 2 Φ 0 0 -D 2 Φ 0   ,

which leads to (38).

We first focus on the case σ = z (1 -z 2 )/6 and τ = -σ/2, for which the matrix B reads:

B = -z      1 -z/2 -(1 -z 2 )/12 -z/2 1 + 2 z 2 6 z (1 -z 2 )/12 -(1 -z 2 )/12 z (1 -z 2 )/12 (1 -z 2 ) (1 + z) (4 -3 z) 72      ,
which we have already analyzed in the proof of Proposition 1. In particular, we have already shown that B is negative semi-definite for any z ∈ (0, 1) and the result of Proposition 2 follows in that case. Finally, and most importantly, the fourth order Strang scheme σ = z (1 -z 2 )/6, τ = -z 2 (1 -z 2 )/24 yields the expression,

B = -z     1 -z/2 -z (1 -z 2 )/24 -z/2 1 + 2 z 2 6 z 2 (1 -z 2 )/24 -z (1 -z 2 )/24 z 2 (1 -z 2 )/24 z (1-z 2 ) (4-z 2 +z (1-z 2 )) 144     .
We thus compute:

  U V W   T B   U V W   = -z U - z 2 V - z (1 -z 2 ) 24 W 2 + 1 6 V 2 + z 2 12 V + 1 -z 2 4 W 2 + z (1 -z 2 ) (4 -z 2 ) 144 W 2 ,
showing that B is negative semi-definite for any z ∈ (0, 1). Hence the result of Proposition 2 follows in the case of the fourth order Strang scheme.

In view of all above results, the energy method seems to be a rather efficient tool to prove stability estimates for numerical schemes that are contractive in the whole space Z combined with some carefully selected extrapolation conditions at an outflow boundary. It is the purpose of the next section to illustrate the limitations of this energy approach, namely that it can not be applied to the more often considered translation invariant extrapolation conditions. For such boundary conditions and more general schemes (based for instance on multistep quadrature methods in time), the general theory for proving stability of numerical boundary conditions initiated in [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF] remains the only available one.

On the limitation of the energy method for analyzing numerical boundary conditions

In this paragraph, we consider the fourth order Strang scheme [START_REF] Strang | Trigonometric polynomials and difference methods of maximum accuracy[END_REF]:

A(z) := I -z D 0 + z 2 2 ∆ + z (1 -z 2 ) 6 D 0 ∆ - z 2 (1 -z 2 ) 24 ∆ 2 . ( 39
)
With a prescribed five point stencil, the finite difference scheme (39) is the only one that achieves fourth order accuracy (with respect to both space and time), just like the Lax-Wendroff scheme is the only three point scheme that achieves second order accuracy. Our goal is to study the semi-boundedness of (39) when implemented on a half line Z -with extrapolation numerical boundary conditions, be they for instance of order 1 or 2. The extrapolation boundary conditions considered in [Kre66, Gol77, CL20] are, in the terminology of [START_REF] Goldberg | Scheme-independent stability criteria for difference approximations of hyperbolic initial-boundary value problems[END_REF][START_REF] Goldberg | Scheme-independent stability criteria for difference approximations of hyperbolic initial-boundary value problems[END_REF], translatory, meaning that they are of the exact same form in each ghost cell (as opposed to the extrapolation conditions considered in Propositions 1 and 2). The main results in [START_REF] Kreiss | Difference approximations for hyperbolic differential equations[END_REF][START_REF] Kreiss | Stability theory for difference approximations of mixed initial boundary value problems[END_REF][START_REF] Goldberg | On a boundary extrapolation theorem by Kreiss[END_REF][START_REF] Coulombel | The neumann numerical boundary condition for transport equations[END_REF] show that, whatever the extrapolation order at the outflow boundary, the corresponding operator on ℓ 2 (-∞, 0) is power-bounded. The proof of this result in [Kre66, [START_REF] Kreiss | Stability theory for difference approximations of mixed initial boundary value problems[END_REF][START_REF] Goldberg | On a boundary extrapolation theorem by Kreiss[END_REF] relies on the normal mode decomposition and power-boundedness follows from the general result in [START_REF] Wu | The semigroup stability of the difference approximations for initial-boundary value problems[END_REF], while the more direct proof of the same result in [START_REF] Coulombel | The neumann numerical boundary condition for transport equations[END_REF] relies first on the energy method for the Dirichlet boundary condition and on an induction argument with respect to the extrapolation order at the outflow boundary. Our goal below is to determine whether, as in Corollary 2, such stability estimates can be achieved by means of a straightforward energy method without using "auxiliary problems". The answer is negative, see Theorems 1 and 2 below, which seems to indicate that the induction argument in [START_REF] Coulombel | The neumann numerical boundary condition for transport equations[END_REF] is more or less the shortest way to derive stability estimates for translatory extrapolation numerical boundary conditions. Examples of uniformly stable, though non dissipative, boundary conditions have been known for quite some time in the context of hyperbolic partial differential equations, see for instance the examples provided in [START_REF] Benzoni-Gavage | Generic types and transitions in hyperbolic initial-boundary-value problems[END_REF][START_REF] Benoit | Geometric optics expansions for linear hyperbolic boundary value problems and optimality of energy estimates for surface waves[END_REF]. As far as we are aware of, the example given by Theorem 1 seems to be the first in the fully discrete setting.

Our first main result for the Strang scheme (39) is the following.

Theorem 1. There does not exist a parameter z 0 > 0, an integer r ∈ N, and a continuous map H from [0, z 0 ] with values in M r (R) such that for all z ∈ [0, z 0 ], H(z) is a symmetric matrix satisfying the following property: for any sequence Φ ∈ ℓ 2 (-∞, 2) satisfying the first order extrapolation boundary condition Φ 2 = Φ 1 = Φ 0 , the following energy inequality holds:

j≤-r Ψ 2 j +    Ψ -r+1 . . . Ψ 0    T H(z)    Ψ -r+1 . . . Ψ 0    - j≤-r Φ 2 j -    Φ -r+1 . . . Φ 0    T H(z)    Φ -r+1 . . . Φ 0    ≤ 0 , ( 40
)
with the sequence Ψ ∈ ℓ 2 (Z -) being defined by:

∀ j ≤ 0 , Ψ j := A(z) Φ j = Φ j -z D 0 Φ j + z 2 2 ∆ Φ j + z (1 -z 2 ) 6 D 0 ∆ Φ j - z 2 (1 -z 2 ) 24 ∆ 2 Φ j .
Let us observe that we do not even assume the matrix H(z) to be positive definite, which would be necessary to make the quantity:

j≤-r Φ 2 j +    Φ -r+1 . . . Φ 0    T H(z)    Φ -r+1 . . . Φ 0    ,
the square of a norm on the space of sequences Φ ∈ ℓ 2 (-∞, 2) with the prescribed boundary conditions. The obstacle for the existence of H(z) in Theorem 1 comes from the small values of z. It could very well be that for some z ∈ (0, 1), one can construct a real symmetric positive definite matrix H(z) of size r ∈ N such that one has the optimal energy balance:

j≤-r Ψ 2 j +    Ψ -r+1 . . . Ψ 0    T H(z)    Ψ -r+1 . . . Ψ 0    - j≤-r Φ 2 j -    Φ -r+1 . . . Φ 0    T H(z)    Φ -r+1 . . . Φ 0    + c (Φ 2 -1 + Φ 2 0 ) ≤ 0 ,
with c a positive constant. Theorem 1 shows, however, that such a construction will not be possible for all values of z ∈ (0, 1) with the additional requirement that H(z) extends continuously to z = 0 (as was the case, for instance, for the Lax-Wendroff scheme with second order extrapolation by Corollary 2).

Proof of Theorem 1. The proof of Theorem 1 is based on an induction argument with respect to the integer r. As will follow from the argument below, it appears that the initial step of the induction argument corresponds to r = 2. The first two cases r = 0 and r = 1 are dealt with separately, though the argument is similar to the one for r = 2. Let us start with a general argument from which the conclusion of Theorem 1 will follow.

• The derivative of the energy balance with respect to z at 0. In all what follows, we assume that the opposite statement to Theorem 1 holds, meaning that we assume that there exist a parameter z 0 > 0, an integer r ∈ N, and a continuous map H from [0, z 0 ] with values in M r (R) such that for any z ∈ [0, z 0 ], H(z) is a real symmetric matrix with the previously stated property. Let us consider a sequence Φ ∈ ℓ 2 (-∞, 2) that satisfies Φ 2 = Φ 1 = Φ 0 , and let us then denote with f (z) the quantity on the left hand side of (40). Since the sequence Ψ(z) depends in a C 1 fashion on z in ℓ 2 , and since H is continuous at 0, the function f is differentiable at 0 and it holds that,

f ′ (0) = j≤-r 2 Ψ j (0) Ψ ′ j (0) + 2    Ψ ′ -r+1 (0) . . . Ψ ′ 0 (0)    T H(0)    Φ -r+1 . . . Φ 0    . (41)
Substituting the value of Ψ ′ j (0), (41) can be expressed in terms of the sequence Φ only, i.e.

f ′ (0) = j≤-r -2 Φ j D 0 Φ j + 1 3 Φ j D 0 ∆ Φ j +    -(Φ -r+2 -Φ -r ) + (Φ -r+3 -2 Φ -r+2 + 2 Φ -r -Φ -r-1 )/6 . . . -(Φ 1 -Φ -1 ) + (Φ 2 -2 Φ 1 + 2 Φ -1 -Φ -2 )/6    T H(0)    Φ -r+1 . . . Φ 0    = - 4 3 Φ -r Φ -r+1 + 1 6 (Φ -r-1 Φ -r+1 + Φ -r Φ -r+2 ) +    -(Φ -r+2 -Φ -r ) + (Φ -r+3 -2 Φ -r+2 + 2 Φ -r -Φ -r-1 )/6 . . . -(Φ 1 -Φ -1 ) + (Φ 2 -2 Φ 1 + 2 Φ -1 -Φ -2 )/6    T H(0)    Φ -r+1 . . . Φ 0    .
We observe that f vanishes at z = 0 because Ψ j (0) = Φ j for all j, and f takes nonpositive values on [0, z 0 ] by (40). Hence 0 is a maximum of f on [0, z 0 ] and the derivative f ′ (0) is nonpositive. In other words, we have just obtained that there exists a real symmetric matrix H of size r such that for any sequence Φ ∈ ℓ 2 (-∞, 2) satisfying Φ 2 = Φ 1 = Φ 0 , the following inequality holds:

- 4 3 Φ -r Φ -r+1 + 1 6 (Φ -r-1 Φ -r+1 + Φ -r Φ -r+2 ) +    -(Φ -r+2 -Φ -r ) + (Φ -r+3 -2 Φ -r+2 + 2 Φ -r -Φ -r-1 )/6 . . . -(Φ 1 -Φ -1 ) + (Φ 2 -2 Φ 1 + 2 Φ -1 -Φ -2 )/6    T H    Φ -r+1 . . . Φ 0    ≤ 0 . ( 42
)
It remains to examine the consequences of (42), which is where the specific value of r comes into play because the values Φ 1 and Φ 2 are not arbitrary due to the extrapolation boundary conditions.

• The case r = 0. In that case, the inequality (42) reduces to:

- 4 3 Φ 0 Φ 1 + 1 6 (Φ -1 Φ 1 + Φ 0 Φ 2 ) ≤ 0 .
Using the extrapolation boundary conditions Φ 2 = Φ 1 = Φ 0 , we end up with:

- 7 6 Φ 2 0 + 1 6 Φ 0 Φ -1 ≤ 0 ,
which is obviously impossible since the values Φ 0 , Φ -1 are arbitrary. This means that we cannot use the standard ℓ 2 norm on Z -for proving the stability of (39) with first order extrapolation condition. Let us now deal with the next case (r = 1) in the induction argument.

• The case r = 1. In that case, the inequality (42) reduces to:

- 4 3 Φ -1 Φ 0 + 1 6 (Φ -2 Φ 0 + Φ -1 Φ 1 ) -H Φ 0 (Φ 1 -Φ -1 ) + H 6 Φ 0 (Φ 2 -2 Φ 1 + 2 Φ -1 -Φ -2 ) ≤ 0 ,
where H is a real number. After using the boundary conditions Φ 2 = Φ 1 = Φ 0 , we get:

- 7 6 Φ -1 Φ 0 + 1 6 Φ -2 Φ 0 -H Φ 0 (Φ 0 -Φ -1 ) - H 6 Φ 0 (Φ 0 -2 Φ -1 + Φ -2 ) ≤ 0 , (43)
where now the three values Φ 0 , Φ -1 , Φ -2 are arbitrary. It is useful to introduce the new variables:

y 3 := Φ 0 , y 2 := Φ 0 -Φ -1 , y 1 := Φ 0 -2 Φ -1 + Φ -2 ,
with which (43) is rewritten as follows:

∀ y ∈ R 3 , -y 2 3 + 5 6 -H y 3 y 2 + 1 -H 6 y 3 y 1 ≤ 0 , (44)
The latter inequality is obviously impossible since the only available parameter H should equal both 1 and 5/6 to cancel the off-diagonal terms y 3 y 2 and y 3 y 1 . This completes the proof of Theorem 1 in the case r = 1. Let us now deal with the case r = 2.

• The case r = 2. This is really the starting point of the induction argument, and we shall borrow the methodology introduced for the case r = 1. Namely, when r equals 2, the inequality (42) reduces to:

- 4 3 Φ -2 Φ -1 + 1 6 (Φ -3 Φ -1 +Φ -2 Φ 0 )+ -(Φ 0 -Φ -2 ) + (Φ 1 -2 Φ 0 + 2 Φ -2 -Φ -3 )/6 -(Φ 1 -Φ -1 ) + (Φ 2 -2 Φ 1 + 2 Φ -1 -Φ -2 )/6 T H Φ -1 Φ 0 ≤ 0 .
(45) After using the boundary conditions Φ 2 = Φ 1 = Φ 0 , (45) reduces to:

- 4 3 Φ -2 Φ -1 + 1 6 (Φ -3 Φ -1 + Φ -2 Φ 0 ) - (Φ 0 -Φ -2 ) + (Φ 0 -2 Φ -2 + Φ -3 )/6 (Φ 0 -Φ -1 ) + (Φ 0 -2 Φ -1 + Φ -2 )/6 T H Φ -1 Φ 0 ≤ 0 . ( 46
)
Let us extend the strategy used in the case r = 1, and introduce the new variables:

y 4 := Φ 0 , y 3 := D Φ 0 , y 2 := D 2 Φ 0 , y 1 := D 3 Φ 0 .
The inequality (46) can be equivalently rewritten as:

-y 2 4 + 3 y 4 y 3 -

2 3 y 4 y 2 - 1 6 y 4 y 1 - 13 6 y 2 3 + 5 6 y 3 y 2 + 1 6 y 3 y 1 - 1 6 13 y 3 -5 y 2 -y 1 6 y 3 + y 2 T H y 4 -y 3 y 4 ≤ 0 , (47)
where (47) holds for all y ∈ R 4 , since Φ 0 , Φ -1 , Φ -2 , Φ -3 in (46) are arbitrary. It is useful at this stage to introduce the coefficients of the symmetric matrix H , and we thus write:

H = h 11 h 12 h 12 h 22 .
The quadratic form in y ∈ R 4 on the left hand side of (47) is nonpositive, and furthermore it has no y 2 1 term. This implies that the coefficients of the cross products y 4 y 1 and y 3 y 1 must vanish. Computing those coefficients, we get:

h 11 = 1 , h 12 = 0 ,
which means that the matrix H reads:

H = 1 0 0 h 22 ,
and then (47) reduces to:

-y 2 4 + 5 6 -h 22 y 4 y 3 + 1 -h 22 6 y 4 y 2 ≤ 0 ,
which is nothing else but the inequality (44) we had obtained in the analysis of the case r = 1 except for the shift in the indeces (one should only substitute (y 4 , y 3 , y 2 ) in place of (y 3 , y 2 , y 1 ) in (44), and h 22 in place of H ). As already observed in the analysis of the case r = 1, we are led to a contradiction, which completes the proof of Theorem 1 in the case r = 2.

• The general case r ≥ 3. We go back to (42) and assume r ≥ 3. In particular, the first line on the left hand side of (42) does not involve the ghost cell values Φ 1 , Φ 2 . Substituting the first order extrapolation boundary conditions Φ 2 = Φ 1 = Φ 0 in (42), we get the inequality:

- 4 3 Φ -r Φ -r+1 + 1 6 (Φ -r-1 Φ -r+1 + Φ -r Φ -r+2 ) +        -(Φ -r+2 -Φ -r ) + (Φ -r+3 -2 Φ -r+2 + 2 Φ -r -Φ -r-1 )/6 . . . -(Φ -1 -Φ -3 ) + (Φ 0 -2 Φ -1 + 2 Φ -3 -Φ -4 )/6 -(Φ 0 -Φ -2 ) -(Φ 0 -2 Φ -2 + Φ -3 )/6 -(Φ 0 -Φ -1 ) -(Φ 0 -2 Φ -1 + Φ -2 )/6        T H        Φ -r+1 . . . Φ -2 Φ -1 Φ 0        ≤ 0 . ( 48
)
We introduce the variable y ∈ R r+2 defined by:

∀ ℓ = 1, . . . , r + 2 , y ℓ := D r+2-ℓ Φ 0 ,
which, conversely, corresponds to:

∀ ℓ = 0, . . . , r + 1 , Φ -ℓ = D ℓ y r+2 .
Using from now on the coordinates of y ∈ R r+2 as free parameters, (48) reads:

- 4 3 D r y r+2 D r-1 y r+2 + 1 6 (D r+1 y r+2 D r-1 y r+2 + D r y r+2 D r-2 y r+2 ) +       
-2 D r-2 y r+1 + D r-2 y r + (2 D r-3 y r-1 -D r-3 y r-2 )/6 . . .

-2 D y r+1 + D y r + (2 y r-1 -y r-2 )/6 -(13 y r+1 -5 y r -y r-1 )/6 -y r+1 -y r /6

       T H        D r-1 y r+2 . . . D 2 y r+2 D y r+2 y r+2        ≤ 0 . ( 49
)
We are not going to compute all the coefficients of the quadratic form (in y) arising on the left hand side of (49). It is useful however to observe that the first two coordinates y 1 and y 2 of y do not appear in the expressions of D r-1 y r+2 , . . . , D y r+2 , y r+2 . Therefore, if we rewrite the quadratic form (in y) arising on the left hand side of (49) as y T S y, with S a real symmetric matrix of size r + 2, then not only S is nonpositive because of (49), but S also reads:

S =   S Υ 2 Υ 1 Υ T 2 0 0 Υ T 1 0 0   ,
with Υ 1 , Υ 1 ∈ R r , and S a real symmetric matrix of size r. Since S is nonnegative, we must necessarily have Υ 1 = Υ 2 = 0. In other words, this means that no cross product of the form y 1 y 3 , . . . , y 1 y r+2 or y 2 y 3 , . . . , y 2 y r+2 arises on the left hand side of (49), or, equivalently, that the quantity on the left hand side of (49) does not depend on (y 1 , y 2 ). Computing the partial derivative with respect to y 1 , we get the relation:

(-1) r+1 6 (1 -H 11 ) D r-1 y r+2 + r ℓ=2
H 1ℓ D r-ℓ y r+2 = 0 , from which we deduce that the first line of H should read:

1 0 • • • 0 .
Since H is symmetric, (49) reduces to:

- 4 3 D r y r+2 D r-1 y r+2 + 1 6 (D r+1 y r+2 D r-1 y r+2 + D r y r+2 D r-2 y r+2 ) + -2 D r-2 y r+1 + D r-2 y r + 1 3 D r-3 y r-1 - 1 6 D r-3 y r-2 D r-1 y r+2 (50) +        -2 D r-3 y r+1 + D r-3 y r + (2 D r-4 y r-1 -D r-4 y r-2 )/6 . . . -2 D y r+1 + D y r + (2 y r-1 -y r-2 )/6 -(13 y r+1 -5 y r -y r-1 )/6 -y r+1 -y r /6        T H ♯        D r-2 y r+2 . . . D 2 y r+2 D y r+2 y r+2        ≤ 0 ,
where the real symmetric matrix H ♯ of size r -1 corresponds to the block decomposition of H :

H = 1 0 0 H ♯ .
We can simplify the first two lines of (50) by using the relation:

-2 D r-2 y r+1 + D r-2 y r + 1 3 D r-3 y r-1 - 1 6 D r-3 y r-2 = -D r-2 y r+2 -D r y r+2 + 1 6 D r-3 y r+2 -2 D r-2 y r+2 + 2 D r y r+2 -D r+1 y r+2 ,
and (50) can be rewritten as:

- 4 3 D r-1 y r+2 D r-2 y r+2 + 1 6 (D r y r+2 D r-2 y r+2 + D r-1 y r+2 D r-3 y r+2 ) +        -2 D r-3 y r+1 + D r-3 y r + (2 D r-4 y r-1 -D r-4 y r-2 )/6 . . . -2 D y r+1 + D y r + (2 y r-1 -y r-2 )/6 -(13 y r+1 -5 y r -y r-1 )/6 -y r+1 -y r /6        T H ♯        D r-2 y r+2 . . . D 2 y r+2 D y r+2 y r+2        ≤ 0 . ( 51
)
Shfiting the indeces in the variables, that is introducing the vector:

(ỹ r+1 , . . . , ỹ1 ) := (y r+2 , . . . , y 2 ) , and forgetting about the tilde, we see that ( 51) is exactly the same as (49) with the integer r -1 in place of r. By a finite induction process, we can therefore show that the validity of (49) for some real symmetric matrix H of size r implies the validity of (47) (which is exactly (49) in the particular case r = 2), and we have already seen that this leads to a contradiction. The proof of Theorem 1 is now complete.

The exact same argument of proof can be used to deal with the case of second order extrapolation at the boundary. We shall not reproduce the proof here and leave the (minor) modifications to the interested reader. We thus only state the final result, which is entirely similar to Theorem 1 above except for the extrapolation conditions.

Theorem 2. There does not exist a parameter z 0 > 0, an integer r ∈ N, and a continuous map H from [0, z 0 ] with values in M r (R) such that for any z ∈ [0, z 0 ], H(z) is a symmetric matrix satisfying the following property: for any sequence Φ ∈ ℓ 2 (-∞, 2) verifying the second order extrapolation boundary conditions ∆ Φ 1 = ∆ Φ 0 = 0, with the sequence Ψ being defined by:

∀ j ≤ 0 , Ψ j := A(z) Φ j = Φ j -z D 0 Φ j + z 2 2 ∆ Φ j + z (1 -z 2 ) 6 D 0 ∆ Φ j - z 2 (1 -z 2 ) 24 ∆ 2 Φ j ,
then the following energy inequality holds:

j≤-r Ψ 2 j +    Ψ -r+1 . . . Ψ 0    T H(z)    Ψ -r+1 . . . Ψ 0    - j≤-r Φ 2 j -    Φ -r+1 . . . Φ 0    T H(z)    Φ -r+1 . . . Φ 0    ≤ 0 .
Theorems 1 and 2 show that for the fourth order Strang scheme (39), stability for the translatory first or second order extrapolation conditions:

(first order) Φ 2 = Φ 1 = Φ 0 , (second order) ∆ Φ 1 = ∆ Φ 0 = 0 , cannot be obtained by a 'straightforward' energy argument (at least for all relevant values of the CFL parameter z), meaning by the construction of an energy that is a finite rank perturbation of the identity that is non-increasing for the associated evolution operator. It is known nevertheless that these numerical boundary conditions satisfy the strong stability condition of [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF], see [START_REF] Goldberg | On a boundary extrapolation theorem by Kreiss[END_REF]. These two above examples clearly indicate that the theory initiated in [START_REF] Gustafsson | Stability theory of difference approximations for mixed initial boundary value problems[END_REF] is the only one able to capture stability for numerical boundary conditions in general.

A Higher order integration by parts decompositions

This appendix is devoted to the proof of the relations (32)-( 35) which we have used in the proof of Lemma 2. Let us start with the proof of formula (32), which we rewrite here for the reader's convenience:

2 D 0 Φ j D 0 ∆ Φ j = -2 (∆ Φ j ) 2 + 1 4 (D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2 + (D Φ j+1 - 1 4 D ∆ Φ j+1 ) (D ∆ Φ j+1 + 2 ∆ Φ j ) -(D Φ j - 1 4 D ∆ Φ j ) (D ∆ Φ j + 2 ∆ Φ j-1 ) . ( 52
)
We first apply the formula (20) to get:

2 D 0 Φ j ∆ D 0 Φ j = -(D D 0 Φ j ) 2 -(D D 0 Φ j+1 ) 2 + 2 D 0 Φ j D D 0 Φ j+1 + (D D 0 Φ j+1 ) 2 -2 D 0 Φ j-1 D D 0 Φ j + (D D 0 Φ j ) 2 . ( 53
)
and then rewrite the first line on the right hand side of (53) (that is, the symmetric terms) as:

(D D 0 Φ j ) 2 + (D D 0 Φ j+1 ) 2 = 1 4 (∆ Φ j-1 + ∆ Φ j ) 2 + 1 4 (∆ Φ j + ∆ Φ j+1 ) 2 = 1 4 (-D ∆ Φ j + 2 ∆ Φ j ) 2 + 1 4 (2 ∆ Φ j + D ∆ Φ j+1 ) 2 = 2 (∆ Φ j ) 2 + ∆ Φ j ∆ 2 Φ j + 1 4 (D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2
We use again (20) for the term ∆ Φ j ∆ 2 Φ j and get:

(D D 0 Φ j ) 2 + (D D 0 Φ j+1 ) 2 = 2 (∆ Φ j ) 2 - 1 4 (D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2 + 1 2 (∆ Φ j+1 ) 2 - 1 2 (∆ Φ j ) 2 - 1 2 (∆ Φ j ) 2 - 1 2 (∆ Φ j-1 ) 2 .
Substituting this relation in the first line on the right hand side of (53), we have thus obtained the decomposition:

2 D 0 Φ j ∆ D 0 Φ j = -2 (∆ Φ j ) 2 + 1 4 (D ∆ Φ j ) 2 + (D ∆ Φ j+1 ) 2 + T j -T j-1 , with
T j := 2 D 0 Φ j D D 0 Φ j+1 + (D D 0 Φ j+1 ) 2 - 1 2 (∆ Φ j+1 ) 2 + 1 2 (∆ Φ j ) 2 = (D 0 Φ j+1 ) 2 -(D 0 Φ j ) 2 - 1 2 (∆ Φ j+1 ) 2 + 1 2 (∆ Φ j ) 2 .
The decomposition (52) then follows by just rewriting the quantity T j in terms of D Φ j+1 , ∆ Φ j and D ∆ Φ j+1 . The (small) details are left to the reader.

We now wish to justify the purely telescopic formula (33), which we also rewrite here for convenience:

2 Φ j D 0 ∆ Φ j = [Φ j (D ∆ Φ j+1 + 2 ∆ Φ j ) + (∆ Φ j -D Φ j+1 ) D Φ j+1 ] -[Φ j-1 (D ∆ Φ j + 2 ∆ Φ j-1 ) + (∆ Φ j-1 -D Φ j ) D Φ j ] . ( 54
)
Let us write the left hand side of (54) as:

2 Φ j D 0 ∆ Φ j =       Φ j+2 Φ j+1 Φ j Φ j-1 Φ j-2       T       0 0 1/2 0 0 0 0 -1 0 0 1/2 -1 0 1 -1/2 0 0 1 0 0 0 0 -1/2 0 0             Φ j+2 Φ j+1 Φ j Φ j-1 Φ j-2      
, and decompose the corresponding symmetric matrix in a telescopic way:

      0 0 1/2 0 0 0 0 -1 0 0 1/2 -1 0 1 -1/2 0 0 1 0 0 0 0 -1/2 0 0       =       0 0 1/2 0 0 0 0 -1 1/2 0 1/2 -1 0 0 0 0 1/2 0 0 0 0 0 0 0 0       -       0 0 0 0 0 0 0 0 1/2 0 0 0 0 -1 1/2 0 1/2 -1 0 0 0 0 1/2 0 0      
.

At this stage, we have obtained the telescopic decomposition:

2 Φ j D 0 ∆ Φ j = [Φ j Φ j+2 -2 Φ j Φ j+1 + Φ j-1 Φ j+1 ] -[Φ j-1 Φ j+1 -2 Φ j-1 Φ j + Φ j-2 Φ j ] ,
and the proof of (54) follows by rewriting the telescopic term as:

Φ j Φ j+2 -2 Φ j Φ j+1 + Φ j-1 Φ j+1 = Φ j ∆ Φ j+1 + Φ j ∆ Φ j -D Φ j D Φ j+1 = Φ j (D ∆ Φ j+1 + 2 ∆ Φ j ) + (∆ Φ j -D Φ j+1 ) D Φ j+1 .
This completes the proof of (54).

We now turn to the proof of the telescopic formula (34), which we also rewrite here for convenience:

2 D 0 Φ j ∆ 2 Φ j = [ -(D ∆ Φ j+1 + ∆ Φ j ) ∆ Φ j + 2 D Φ j+1 D ∆ Φ j+1 ] -[-(D ∆ Φ j + ∆ Φ j-1 ) ∆ Φ j-1 + 2 D Φ j D ∆ Φ j ] . ( 55
)
We expand the left hand side of (55) as follows:

2 D 0 Φ j ∆ 2 Φ j = (Φ j+1 -Φ j-1 ) (Φ j+2 -4 Φ j+1 + 6 Φ j -4 Φ j-1 + Φ j-2 ) = 6 Φ j (Φ j+1 -Φ j-1 ) -4 (Φ 2 j+1 -Φ 2 j-1 ) + (Φ j+1 -Φ j-1 ) (Φ j+2 + Φ j-2 ) = 6 Φ j Φ j+1 -4 Φ 2 j -4 Φ 2 j+1 -6 Φ j-1 Φ j -4 Φ 2 j-1 -4 Φ 2 j + (Φ j+1 -Φ j-1 ) (Φ j+2 + Φ j-2 ) = T j -T j-1 , with T j := 7 Φ j Φ j+1 -4 Φ 2 j -4 Φ 2 j+1 + Φ j+1 Φ j+2 -Φ j-1 Φ j+2 + Φ j-1 Φ j = -4 (D Φ j+1 ) 2 + (Φ j+2 -Φ j ) (Φ j+1 -Φ j-1 )

= -4 (D Φ j+1 ) 2 + (∆ Φ j+1 + 2 D Φ j+1 ) (2 D Φ j+1 -∆ Φ j ) = -∆ Φ j ∆ Φ j+1 + 2 D Φ j+1 D ∆ Φ j+1 = -(D ∆ Φ j+1 + ∆ Φ j ) ∆ Φ j + 2 D Φ j+1 D ∆ Φ j+1 .
This completes the proof of (55).

It only remains to prove the formula (35), that is:

Φ j ∆ 2 Φ j = (∆ Φ j ) 2 + [Φ j D ∆ Φ j+1 -D Φ j+1 ∆ Φ j ] -[Φ j-1 D ∆ Φ j -D Φ j ∆ Φ j-1 ] , (56)
which is the discrete counterpart of the relation u u ′′′′ = (u ′′ ) 2 + (u u ′′′ -u ′ u ′′ ) ′ . We compute: We have thus obtained the telescopic decomposition:

Φ j ∆ 2 Φ j -(∆ Φ j ) 2 =       Φ j+2 Φ j+1 Φ j Φ j-1 Φ j-2       T       0 0 1/2 0 0 0 -1 0 -1 0 1/2 0 2 0 1/2 0 -1 0 -1 0 0 0 1/2 0 0             Φ j+2 Φ j+1 Φ j Φ j-1 Φ j-2      
Φ j ∆ 2 Φ j -(∆ Φ j ) 2 = Φ j Φ j+2 -(Φ j+1 ) 2 + (Φ j ) 2 -Φ j-1 Φ j+1 -Φ j-1 Φ j+1 -(Φ j ) 2 + (Φ j-1 ) 2 -Φ j-2 Φ j ,
and (56) eventually follows from rewriting the telescopic quantity in terms of Φ j , D Φ j+1 , ∆ Φ j and D ∆ Φ j+1 :

Φ j Φ j+2 -(Φ j+1 ) 2 +(Φ j ) 2 -Φ j-1 Φ j+1 = Φ j (Φ j+2 -3 Φ j+1 +3 Φ j -Φ j-1 )-(Φ j+1 -Φ j ) (Φ j+1 -2 Φ j +Φ j-1 ) .

,

  and decompose the corresponding symetric matrix in a telescopic way:

Research of the author was supported by ANR project NABUCO, ANR-17-CE40-0025. Research of the author was supported by ANR project NABUCO, ANR-17-CE40-0025.


Evaluation of the Turbulence Model Influence on the Numerical Simulations of Unsteady Cavitation

Unsteady cavitation in a Venturi-type section was simulated by two-dimensional compu-tations of viscous, compressible, and turbulent cavitating flows. The numerical model used an implicit finite volume scheme (based on the SIMPLE algorithm) to solve Reynoldsaveraged Navier-Stokes equations, associated with a barotropic vapor/liquid state law that strongly links the density variations to the pressure evolution. To simulate turbulence effects on cavitating flows, four different models were implemented (standard k-RNG; modified k-RNG; k-with and without compressibility effects), and numeri-cal results obtained were compared to experimental ones. The standard models k-RNG and k-without compressibility effects lead to a poor description of the self-oscillation behavior of the cavitating flow. To improve numerical simulations by taking into account the influence of the compressibility of the two-phase medium on turbulence, two other models were implemented in the numerical code: a modified k-model and the k-model including compressibility effects. Results obtained concerning void ratio, velocity fields, and cavitation unsteady behavior were found in good agreement with experimental ones. The role of the compressibility effects on turbulent two-phase flow modeling was analyzed, and it seemed to be of primary importance in numerical simulations.

Introduction

Cavitating flows in turbomachinery lead to performance losses and modifications of the blades load. These consequences can be related to the quasi-steady cavitation effects, mainly depending on the time-averaged shapes of the vaporized structures. Moreover, turbomachinery under cavitating conditions is also submitted to transient phenomena, such as compressibility effects, flow rate fluctuations, noise, vibration, erosion, in which the unsteady behavior and the two-phase structure of the cavitating flow are influent. The study and modeling of the unsteady behavior of cavitation are thus essential to estimate the hydraulic unsteadiness in turbomachinery and the associated effects.

In this context, some numerical models taking into account the two-phase structure and the unsteady behavior of cavitation have been developed. Most of them are based on a single fluid approach: the relative motion between liquid and vapor phases is neglected and the liquid vapor mixture is treated as a homogeneous medium with variable density. The mixture density is related directly to the local void fraction and is managed either by a state law ͑Delannoy and Kueny ͓1͔, Song and He ͓2͔, and Merkle et al. ͓3͔͒, using a supplementary equation relating the void fraction to the dynamic evolution of bubble cluster ͑Kubota et al. ͓4͔, and Chen and Heister ͓5͔͒, or by a multiple species approach with a mass transfer law between liquid and vapor ͑Kunz et al. ͓6͔͒. The last model can be used in a full two-fluid approach, including relative motion between the phases ͑Alajbegovic et al. ͓7͔͒.

Main numerical difficulties are related to the strong coupling between the pressure field and the void ratio and to the coexistence of the strong compressibility of the two-phase medium with the quasi-incompressible behavior of the pure liquid flow. Moreover, the influence of the turbulence on unsteady two-phase compressible flows is not yet well known.

This paper presents unsteady cavitating flow simulations performed with a two-dimensional code. It is based on the model developed by Delannoy and Kueny ͓1͔ for inviscid fluids. Several physical modifications have been investigated by Reboud and Delannoy ͓8͔,Reboud et al. ͓9͔, to increase the range of applications and improve the physical modeling.

In the numerical code, the two-phase aspects of cavitation are treated by introducing a barotropic law that strongly links the fluid density to the pressure variations. From the numerical point of view, this physical approach is associated with a pressurecorrection scheme derived from the SIMPLE algorithm, slightly modified to take into account the cavitation process. The model has been validated on numerous cases, such as Venturi, ͓9,10͔, foils, ͓8,11͔, or blade cascades, ͓12,13͔. The numerical results showed a good agreement with experiments. Mainly, the code leads to a reliable simulation of the cyclic self-oscillation behavior of unsteady cavitating flows.

The complex unsteady mechanism that governs this cyclic cavitation behavior is strongly affected by the applied turbulence model, ͓9,10͔. The aim of this paper is to study the influence of different models ͑standard k-RNG; modified k-RNG; kincluding or not compressibility effects͒ on the numerical simulation of cavitating flows. Indeed, because of the barotropic state law adopted in the model, the vaporization and condensation processes correspond mainly to highly compressible flow areas. Therefore, it is of primary importance to take into account the compressibility of the vapor/liquid mixture in the turbulence model.

The paper describes the turbulence models applied in the numerical simulations and presents comparisons between different numerical results obtained in a two-dimensional Venturi-type section, whose experimental behavior has been already studied by Stutz and Rebouch ͓14,15͔. In this geometry, the flow is characterized by an unstable cavitation behavior, with almost periodical vapor cloud shedding.

In addition to the standard k-model, three other turbulence models were applied, to investigate their influence on the twophase turbulent flow. Mainly, the turbulence model proposed by Wilcox ͓16͔, which includes compressibility considerations, was applied, and promising results are obtained. Comparisons with experimental measurements of void ratio and velocities inside the cavity, ͓9͔, are presented.

The ability of these turbulence models to predict complex twophase flow is discussed, and the role of the compressibility effects is studied.

Physical Model

The present work applies a single fluid model based on previous numerical and physical studies, ͓1,8 -10͔. The fluid density varies in the computational domain according to a barotropic state law ( P) that links the density to the local static pressure ͑Fig. 1͒. When the pressure in a cell is higher than the neighborhood of the vapor pressure ( PϾ P v ϩ(⌬ P v /2)), the fluid is supposed to be purely liquid. The entire cell is occupied by liquid, and the density l is calculated by the Tait equation, ͓17͔. If the pressure is lower than the neighborhood of the vapor pressure ( PϽ P v Ϫ(⌬ P v /2)), the cell is full of vapor and the density v is given by the perfect gas law ͑isotherm approach͒. Between purely vapor and purely liquid states, the cell is occupied by a liquid/vapor mixture, which is considered as one single fluid, with the variable density . The density is directly related to the void fraction ␣ ϭ(Ϫ l )/( v Ϫ l ) corresponding to the local ratio of vapor contained in this mixture.

To model the mixture state, the barotropic law presents a smooth link in the vapor pressure neighborhood, in the interval Ϯ(⌬ P v /2). In direct relation with the range ⌬ P v , the law is characterized mainly by its maximum slope 1/A min 2 , where A min 2 ϭץP/.ץ A min can thus be interpreted as the minimum speed of sound in the mixture. Its calibration was done in previous studies, ͓10͔. The optimal value was found to be independent of the hydrodynamic conditions, and is about 2 m/s for cold water, with P v ϭ0.023 bar, and corresponding to ⌬ P v Ϸ0.06 bar ͑green chart on Fig. 1͒. That value is applied for the computations presented hereafter.

Mass fluxes resulting from vaporization and condensation processes are treated implicitly by the barotropic state law, and no supplementary assumptions are required. Concerning the momentum fluxes, the model assumes that locally velocities are the same for liquid and for vapor: In the mixture regions vapor structures are supposed to be perfectly carried along by the main flow. This hypothesis is often assessed to simulate sheet-cavity flows, in which the interface is considered to be in dynamic equilibrium, ͓3͔. The momentum transfers between the phases are thus strongly linked to the mass transfers.

Numerical Model

To solve the time-dependent Reynolds-averaged Navier-Stokes equations associated with the barotropic state law presented here above, the numerical code applies, on two-dimensional structured curvilinear-orthogonal meshes, the SIMPLE algorithm, modified to take into account the cavitation process. It uses an implicit method for the time-discretization, and the HLPA nonoscillatory second-order convection scheme proposed by Zhu ͓18͔. The numerical model is detailed in Coutier-Delgosha et al. ͓10͔: A complete validation of the method was performed, and the influence of the numerical parameters was widely investigated. The present paper describes mainly the different turbulence models applied.

Boundary Conditions. In the code, the velocity field is imposed at the computational domain inlet, and the static pressure is imposed at the outlet. Along the solid boundaries, the turbulence models are associated with laws of the wall.

Initial Transients Conditions.

To start unsteady calculations, the following numerical procedure is applied: First of all, a stationary step is carried out, with an outlet pressure high enough to avoid any vapor in the whole computational domain. Then, this pressure is lowered slowly at each new time-step, down to the value corresponding to the desired cavitation number . Vapor appears during the pressure decrease. The cavitation number is then kept constant throughout the computation.

The Geometry

Numerical simulations have been performed on a Venturi-type section whose convergent and divergent angles are, respectively, about 18 deg and 8 deg ͑Fig. 2͑a͒͒. The shape of the Venturi bottom downstream from the throat simulates an inducer blade suction side with a beveled leading edge geometry ͑Kueny et al. ͓19͔͒ and a chord length L ref ϭ224 mm.

According to experimental observations, in this geometry the flow is characterized by unsteady cavitation behavior, ͓14͔, with quasi-periodic fluctuations. Each cycle is composed of the following successive steps: The attached sheet cavity grows from the Venturi throat. A re-entrant jet is generated at the cavity closure and flows along the Venturi bottom toward the cavity upstream end. Its interaction with the cavity surface results in the cavity break off. The generated vapor cloud is then convected by the main stream, until it collapses.

For a cavitation number of about 2.4 ͑based on the timeaveraged upstream pressure͒ and an inlet velocity V ref ϭ7.2 m/s, vapor shedding frequency observed experimentally is about 50 Hz for a cavity length of 45Ϯ5 mm, ͓14͔.

The standard computational grid is composed of 160ϫ50 orthogonal cells ͑Fig. 2͑a͒͒. A special stretching of the mesh is applied in the main flow direction just after the throat, so that the two-phase flow area is efficiently simulated: about 50 grid points are used in this direction to model the 45-mm long mean cavity obtained hereafter ͑Fig. 2͑b͒͒. In the other direction, a stretching is also applied close to the walls, to obtain at the first grid point the non dimensional parameter y ϩ of the boundary layer varying between 30 and 100 and to use standard laws of the walls. The grid is finer in the bottom part of the Venturi section than in its upper part, to enhance the accuracy in the cavitation domain: Cavities obtained in the following sections contain about 30 cells across their thickness.

Turbulence Models

"a… Standard k-RNG Model. The first model applied in the numerical code is a standard k-RNG model, associated with laws of the wall, ͓9͔.

In this model, the effective viscosity applied in the Reynolds equations is defined as ϭ t ϩ l where t ϭC k 2 / is the turbulent viscosity and C ϭ0.085 ͑Yakhot et al. ͓20͔͒.

This model is originally devoted to fully incompressible fluids, and no particular correction is applied here in the case of the highly compressible two-phase mixture. Thus, the fluid compressibility is only taken into account in the turbulence equations through the mean density changes.

With this model, the unstable cavitating behavior observed experimentally is not correctly simulated: After an initial transient fluctuation of the cavity length, the numerical calculation leads to a quasi-steady behavior of the cavitation sheet, which globally stabilizes ͑Fig. 3͒.

The resulting cavity length is much too small, compared with the experimental observations ͑in this case, the error is larger than 50 percent͒. Moreover, comparisons with experimental data ob-tained by double optical probes by Stutz and Reboud ͓14,15͔ show that the numerical mean void ratio is overestimated in the main part of the cavity. Calculations give a high time-averaged void ratio in the upstream part of the cavitation sheet ͑Ͼ90 per-cent͒, abruptly falling to 0 percent in the wake, while the measured void ratio never exceeds 25 percent and decreases slowly from the cavity upstream end to its wake.

This poor agreement with the real configuration seems to be related to an overprediction of the turbulent viscosity in the rear part of the cavity. The cyclic behavior of the cloud cavitation process is strongly related to the re-entrant jet development from the cavity closure, ͓17͔. As a matter of fact, the main problem in the turbulent flow simulations consisted in the premature removal of the reverse flow along the solid wall: the re-entrant jet was stopped too early and it did not result in any cavity break off.

It is worth noting that numerical tests reported in ͓10͔ confirm that the mesh size, spatial scheme, and the time discretization applied in the model do not modify the results. Computations performed with a finer mesh (264ϫ90) and first or second-order accurate time discretization schemes still lead to the same complete stabilization of the cavitation sheet. and the vapor cloud shedding, we modified the standard k-RNG model simply by reducing the mixture turbulent viscosity, mainly in the low void ratio areas:

t ϭ f ͑ ͒C k 2 / where f ͑ ͒ϭ v ϩ ͩ v Ϫ v Ϫ l ͪ n ͑ l Ϫ v ͒ nӷ1.
Indeed, according to the experimental results ͓15͔, the reentrant jet seems to be mainly composed of liquid (␣ϳ0), and thus the reduction of the mixture turbulence viscosity leads to substantial changes in the simulation. The prediction of the unsteady re-entrant jet is now obtained, and the vapor cloud shedding is well simulated ͑See Fig. 4.͒

The transient evolution observed during this unsteady calculation is shown in Fig. 5. Fig. 5͑a͒ illustrates at a given time and for each cross section of the Venturi-type duct the value of the minimal density in the section. By comparison with Fig. 3, it gives information concerning the vapor cloud shedding process: the part of the cavity that breaks off clearly appears, and the fluctuation frequency can be evaluated. Moreover, it also supplies the maximum void ratio in each section. Curve 5b presents the transient upstream pressure evolution.

The experimental self-oscillatory behavior of the cavitation sheet is correctly simulated ͑Fig. 5͒. We obtain in this case a fluctuating cavity whose mean length is about 45 mm ͑i.e., L cav /L ref ϭ0.2). The mean fluctuation period is about 0.59 T ref , and the corresponding shedding frequency equals 55 Hz. Both length and frequency are in good agreement with measurements ͑i.e., 45 mm and 50 Hz, respectively͒. The reduced frequency based on the reference velocity and the mean cavity length equals 0.33, which only slightly overestimates the classical Strouhal number Stϭ0.3 ͑0.28 from the experimental measurements͒. The transient evolution is almost periodic, with some random disturbances affecting the oscillation regularity. This behavior is consistent with the experimental observations, which pointed out a quite regular cavitation cycle, whose frequency fluctuated around a central value, ͓14͔.

Several tests were performed to investigate the influence of the numerical parameters on the unsteady cavitating result. The whole study, including the influence of the grid size, the time-step, the order of the time discretization scheme, the ratio v / l , and the maximum slope of the barotropic law A min , is detailed in ͓10͔. Results concerning the grid resolution and the time step influence are reported in the present paper in Table 1. The effect of the mesh size appears to be small, so far it is fine enough: the cavity oscillation frequency is almost constant with the two finest grids. In contrast, the influence of the time-step cannot be completely removed: using ⌬tϭ0.005 T ref leads to a systematic slight overestimation of the cycles frequency. Nevertheless, it must be considered that the uncertainty on experimental measurements is of the same order of magnitude.

"c… k-Model: Compressibility Effects. As presented before, the proposed numerical simulations take into account a single-flow physical model to describe cavitation phenomenon. According to the adopted barotropic state law, in the vapor/liquid mixture zones the sound celerity Aϭͱץ P/ץ is very low, the fluid is locally highly compressible. Thus the flow often reaches Mach numbers larger than 5 in the vaporization or condensation areas. As a matter of fact, Birch and Eggers ͓21͔ showed that an increase of the Mach number could modify substantially the turbulence structure in some configurations, such as mixing layers. To ana-lyze and to evaluate the compressibility effects on the turbulence modeling of cavitating flows, we implemented in the numerical code the k-turbulence model proposed by Wilcox ͓16͔. Its major advantage, compared with the previous models, consists in the corrections that are proposed by Wilcox in the turbulence equations to model the effects of density fluctuations. The main features of the models are presented hereafter. All details concerning the turbulence model can be found in ͓16͔.

Basic Model Equations. Equations are similar to the governing equations resolved in the case of the k-model: All the modifications result from the use of the dissipation specific rate instead of the turbulence dissipation . The turbulent viscosity expression and the k-equation are modified, and the equation is replaced by a equation:

t ϭ k ץk ץt ϩu j ץk ץx j ϭ i j ץu i ץx j Ϫ␤ w *kϩ ץ ץx j ͫ ͑ ϩ w * t ͒ ץk ץx j ͬ ץ ץt ϩu j ץ ץx j ϭ␣ w k i j ץu i ץx j Ϫ␤ w 2 ϩ ץ ץx j ͫ ͑ ϩ w t ͒ ץ ץx j ͬ
where ␣ w ϭ5/9, ␤ w ϭ3/40, ␤ w *ϭ9/100, w ϭ0.5, and w *ϭ0.5 are the coefficients of the model.

i j ϭ2 t S i j with S i j ϭ 1 2 ͑ u i, j ϩu j,i ͒
Compressibility Effects. By analyzing the reduction of the mixing layer growth rate as a function of the Mach number, Wilcox ͓16͔ observed that the phenomenon could not be well simulated only by taking into account the mean density evolution through the mixing length, and that an explicit supplementary correction must be applied. He concludes that the Reynolds-stress transport equation is directly concerned with the compressibility effects, and introduced the influence of the Mach number in the turbulence equations.

The modifications proposed by Wilcox are based on the previous studies of Sarkar et al. ͓22͔ and Zeman ͓23͔, who aimed to take into account the compressibility effects in a k-turbulence model.

He adapted this development to the kmodel to propose the following equations giving the parameters ␤ w and ␤ w * as functions of the turbulence Mach number defined as M t 2 ϭ2k/A 2 ͑where A is the local sound celerity͒.

␤ w *ϭ␤ wi * (1ϩ*F(M t )) and ␤ w ϭ␤ wi Ϫ␤ wi * *F(M t ), where the index i indicates the values of the incompressible model given before with

F͑ M t ͒ϭ ͭ ͑ M t 2 ϪM t 0 2 ͒ if M t ϾM t 0 0 if M t рM t 0
M t 0 ϭ0.25 and *ϭ1.5.

Results. Figure 6 presents calculations performed with the k-turbulence model.

In the first case ͑Fig. 6͑a͒͒, the influence of the compressibility is not taken into account in the model. The numerical results are not satisfactory: As in the case of standard k-RNG model, the unsteady cavitation behavior is not correctly simulated and the cavity length obtained by simulations is more than 50 percent smaller than the experimental one. On the other hand, numerical results obtained by taking into account the compressibility effects are in good agreement with experimental ones. They are very similar to those presented previously in the case of the modified k-RNG model: As can be seen in Fig. 6͑b͒, the vapor cloud shedding process is correctly simulated. The mean cavity length is still correct, and the shedding frequency remains almost the same.

Comparisons are investigated with experimental data obtained by using double optical probes, ͓15͔. This intrusive technique allows measurements of the local void ratio and the velocities of the two-phase structures inside the cavitation sheet. Their timeaveraged and standard deviation values are presented for four profiles in Fig. 7.

These comparisons point out a general good agreement between experimental and numerical results. The time-averaged void ratio repartition is well predicted, and the re-entrant jet structure under the cavity appears qualitatively correct. The fluctuation levels are also in reliable close agreement with measurements. The negative velocity region at station xϭ0.065 m corresponds to the rear part of the cavity, alternatively affected by the re-entrant jet progression and by the vapor cloud shedding. Both models fail to simulate these negative velocities, because of a slight underestimation of the mean cavity length in this configuration. A second discrepancy can be observed at station 0.03 m, concerning the void ratio distribution close to the wall; nevertheless, the experimental estimation of velocities and void ratio in the re-entrant jet is associated to a high uncertainty in this area, due to important fluctuations level, ͓14͔. Thus, the comparison, along the wall, between experiments and the numerical result cannot be relevant.

Conclusion

Four turbulence models have been applied to simulate unsteady cavitating flows in a Venturi-type section.

Results obtained with the standard k-RNG model and with the k-model without taking into account compressibility effects are in poor agreement with experimental observations: the models fail to reproduce the vapor cloud shedding behavior observed experimentally. Moreover, numerical tests reported in ͓10͔ and the results presented here above indicate that the global level of dissipation induced by the numerical model seems not to be responsible for its inability to simulate the unsteady cavity behavior.

On the other hand, the modified k-RNG model and the kmodel including compressibility effects lead to a very reliable simulation of the Venturi unsteady cavitation behavior. The satisfactory results obtained with the modified k-RNG model have been confirmed by simulations performed in other geometries, such as a hydrofoil, ͓11͔, a foil cascade, ͓13͔, and another Venturitype section leading to a more stable sheet of cavitation, ͓10͔. In all cases, the general experimental behavior is correctly obtained, and oscillation frequencies are accurately predicted in unsteady configurations.

Results obtained in the present paper with kmodels pointed out the relevant role of the compressibility effects on the turbulence model. Indeed, the corrections proposed to treat compressibility effects take into account the density fluctuations: a supplementary term appears from the averaged equations, increasing the turbulent dissipation. The final effect is a reduction of the turbulent viscosity in the compressible areas, i.e., in the vapor/liquid mixture zones.

The modifications proposed to improve the k-RNG model had been also based on the reduction of the turbulent viscosity in the mixture zones, which are characterized by a very low sound celerity and large Mach number. Indeed, results obtained with the modified k-RNG model are similar to the ones calculated by the compressible kmodel.

According to the numerical calculations, the fluid compressibility has a strong effect on the turbulence structure, and must be taken into account to simulate unsteady cavitating flows.

Fig. 1

 1 Fig. 1 Barotropic state law "P…. Water 20°C.

"

  Fig. 4 Modification of the mixture viscosity "nÄ10…

Fig. 5

 5 Fig. 5 Transient evolution of unsteady cavitating flow in the Venturi-type duct. "a… Temporal evolution "in abscissa… of the cavity length "graduated in ordinate…. Instantaneous density distribution of attached and cloud cavities are drawn on the left at TÄ11T ref "velocity vectors are drawn only 1 cell over 2 in the two directions…. "b… Time evolution of the upstream pressure.

Fig. 6

 6 Fig. 6 Unsteady cavitation behavior calculated with the kmodel-"a… without and "b… with compressibility effects

Fig. 7

 7 Fig. 7 Time-averaged values and standard deviation of velocity "1,2… and void ratio "3,4…. Numerical results: "a… modified k-, "b… compressible k-"lines… and optical probe measurements "points…-V ref Ä7.2 mÕs. Mean cavity length 45 mm: dotted line Äexperimental external shape of the cavity, obtained from image processing "sigmaÄ2.4….

NomenclatureA

  min ϭ minimum speed of sound in the medium ͑m/s͒ Cp ϭ ( PϪ P ref )/(V ref 2 /2) pressure coefficient k ϭ turbulent kinetic energy per unit mass ͑m 2 s Ϫ2 ͒ L ref ϭ geometry reference length ͑m͒ P ϭ local static pressure ͑Pa͒ P ref ϭ reference pressure ͑Pa͒ P v ϭ vapor pressure ͑Pa͒ V ref ϭ reference velocity (ϭV upstream ) ͑ms Ϫ1 ͒ ␣ ϭ void ratio ϭ turbulence dissipation per unit mass ͑m 2 s Ϫ3 ͒ l , t ϭ laminar and turbulent dynamic viscosity ͑Pa•s͒ , t ϭ laminar and turbulent kinetic viscosity ͑m 2 /s͒ ⌬ P v ϭ range of mixture state law ͑Pa͒ ϭ mixture density ͑kg/m 3 ͒ l , v ϭ liquid, vapor density ͑kg/m 3 ͒ ref ϭ reference density (ϭ l ) ͑kg/m 3 ͒ ϭ ( P upstream Ϫ P v )/(V ref 2 /2) cavitation number ϭ specific dissipation rate ͑s Ϫ1 ͒

Table 1 Numerical tests performed on the Venturi-type section with the k-modified model

 1 Ä2.4, A

min Ä2 mÕs, v Õ l Ä0.01

				RESULTS
	⌬T/T ref	Mesh	Mean Vapor Volume (ϫ10 Ϫ4 m 3 )	Oscillation Frequency ͑Hz͒	Standard Deviation (ϫ10 Ϫ4 m 3 )
			Reference case
	0.005	160ϫ50	11.0	55	2.4
			Influence of the time-step
	0.002	160ϫ50	12.3	51	4.3
	0.005	160ϫ50	11.0	55	2.4
	0.01	160ϫ50	9.3	55	2.7
			Influence of the mesh
	0.002	264ϫ90	15.6	52	2.9
	0.005	264ϫ90	14.5	54	4.3
	0.005	160ϫ50	11.0	55	2.4
	0.005	110ϫ30	6.9	64	1.9

Currently at ENSTA UME/DFA, chemin de la Huniere, 91761 Palaiseau Cedex, France. O. Coutier-Delgosha, R. Fortes-Patella, LEGI-INPG, BP 53, 38041 Grenoble Cedex 9, France J. L. Reboud, ENISE-LTDS, 58, rue Jean Parot, 42023 St. Etienne, France

Acknowledgments

The authors wish to express their gratitude to the French Space Agency: Centre National d'Etudes Spatiales ͑CNES͒ for its continuous support.