François Dubois

An introduction to finite volumes for gas dynamics

Keywords: AMS Subject Classification: 35L40, 35L60, 35L65, 35Q35, 76N15 Advection, Characteristics, Roe matrix, Van Leer method

We propose an elementary introduction to the finite volume method in the context of gas dynamics conservation laws. Our approach is founded on the advection equation, the exact integration of the associated Cauchy problem, and the so-called upwind scheme in one space dimension. It is then extended in three directions : hyperbolic linear systems and particularily the system of acoustics, gas dynamics with the help of the Roe matrix and two space dimensions by following the approach proposed by Van Leer. A special emphasis on boundary conditions is proposed all along the text.

Contents 1) Advection equation and method of characteristics

1.1 Advection equation 1.2 Initial-boundary value problems for the advection equation 1.3 Inflow and outflow for the advection equation 2) Finite volumes for linear hyperbolic systems 2.1 Linear advection 2.2 Numerical flux boundary conditions 2.3 A model system with two equations 2.4 Unidimensional linear acoustics 2.5 Characteristic variables 2.6 A family of model systems with three equations 2.7 First order upwind-centered finite volumes 3) Gas dynamics with the Roe method 3.1 Nonlinear acoustics in one space dimension 3.2 Linearization of the gas dynamics equations 3.3 Roe matrix 3.4 Roe flux 3.5 Entropy correction 3.6 Nonlinear flux boundary conditions 4) Second order and two space dimensions 4.1 Towards second order accuracy 4.2 The method of lines 4.3 The method of Van Leer 4.4 Second order accurate finite volume method for fluid problems 4.5 Explicit Runge-Kutta integration with respect to time 5) References

•

Acknowledgments. The author thanks Alexandre Gault, listener at the spring 2000 "lectures in computational acoustics" at the Conservatoire National des Arts et Métiers (Paris, France), for providing his personal manuscript notes.



1)

Advection equation and method of characteristics. 1.1 Advection equation.

•

We consider a given real number a > 0 and we wish to solve the so-called advection equation of unknown function u(x, t) :

(1.1.1) ∂u ∂t + a ∂u ∂x = 0 , t ≥ 0 , x ∈ IR .

We first look to the homogeneity coherence of the different terms of equation (1.1.1). On one hand, the ratio ∂u ∂t is homogeneous to the dimension [u] of function u(•, •) divided by the dimension [t] of the time and we have : ∂u ∂t ∼

[u]

[t] . On the other hand the expression a ∂u ∂x is homogeneous to the dimension [a] of scalar a multiplied by the ratio [u] [x] and we have a ∂u ∂x ∼ [a] [u] [x] . From equation (1.1.1), the two previous terms ∂u ∂t and a ∂u ∂x have the same dimension and we deduce from the previous formulae the equality :

1 [t] ∼ [a] [x]
. Then we have established that the constant a is homogeneous to a celerity :

(1.1.2)

[a] ∼ [x]

[t] .

•

The Cauchy problem for the model equation (1.1.1) is composed by the equation (1.1.1) itself and the following initial condition : (1.1.3) u(x, 0) = u 0 (x) , x ∈ IR , where IR ∋ x -→ u 0 (x) ∈ IR is some given function. We observe that the solution of equation (1.1.1) is constant along the characteristic (straight) lines that satisfy the differential equation (1.1.4) dx dt = a .

Proposition 1.1. The solution is constant along the characteristic lines. Let 0 ≤ λ ≤ t be some given parameter and u(•, •) a solution of equation (1.1.1). Then function u(•, •) is constant along the characteristic lines, i.e.

(1.1.5) u(xaλ, tλ) = u(x, t) , ∀ x, t, λ .

•

The proof of Proposition 1.1 is obtained as follows. We consider a fixed point (x, t) in space-time IR × [0, +∞[ and the auxiliary function [0, t] ∋ λ -→ v(λ) = u(xaλ, tλ) . We have, due to the usual chain rule for derivation of operators : dv dλ = (-a) ∂u ∂x -∂u ∂t (x-aλ, t-λ) = 0 if function u(•, •) is solution of the advection equation (1.1.1). Then v(λ) does not depend on variable λ and we have in particular v(λ) = v(0) , which exactly expresses the relation  (1.1.5). We have in particular for λ = t : u(x, t) = u(xat, 0) = u 0 (xat) as illustrated on Figure 1.1. 1.2 Initial-boundary value problems for the advection equation.

•

The second step is concerned by the so-called initial-boundary value problem considered for x > 0 and t > 0 with some given initial condition u 0 (x) for t = 0 and a boundary condition v 0 (t) for x = 0 :

(1.2.1) ∂u ∂t + a ∂u ∂x = 0 , t > 0 , x > 0 ,

(1.2.2) u(x, 0) = u 0 (x) , x > 0 , (initial condition)



(1.2.3) u(0, t) = v 0 (t) , t > 0 , (boundary condition). Proposition 1.2.

Advection in the quadrant x > 0 and t > 0 . We suppose that a > 0 . Then the solution of the advection equation (1.2.1) with the initial condition (1.2.2) and the boundary condition (1.2.3) is given by the relations (1.2.4) u(x, t) = u 0 (xat) , xat > 0 (1.2.5) u(x, t) = v 0 t -x a , xat < 0 .

The initial condition u 0 (•) is advected towards space-time point (x, t) when xat > 0 and the boundary condition v 0 (•) is activated for xat < 0 .

• Proof of Proposition 1.2. In order to solve the problem (1.2.1)-(1.2.3), we use the method of characteristics. We fix a point (x, t) of space-time domain that satisfies x > 0, t > 0 and we go upstream in time with the help of the characteristic line that goes through this point (see Figure 1.2) : (1.2.6)

x(λ) = xaλ , t(λ) = tλ .

•

First case : xat > 0 . When we take the particular value λ = t in the previous relation (1.2.6), the particular point y = x(t) = xat on the axis of abscissa is strictly positive then the initial condition u 0 (y) is well defined. The solution u(•, •) is constant on the characteristic line (see Proposition 1.1) that contains this particular point. Then relation (1.2.4) is established.

•

Second case : xat < 0 . We consider the particular value λ = x a inside the expression (1.2.6). Then the corresponding foot of the characteristic belongs to the time axis : θ = tλ = t -x a and θ > 0 due to the inequalities x < at and a > 0 . The solution is constant along the characteristic line going through this point and the relation (1.2.5) is established.

•

In the particular case where datum u 0 (x) is identically equal to zero, i.e.

(1.2.7) u 0 (x) = 0 , x > 0 , and if the boundary condition v 0 (t) is sinusoïdal for time positive to fix the ideas, (1.2.8) v 0 (t) = sin(ωt) , t > 0 , the solution of the advection equation in the domain x > 0 , t > 0 via the relations (1.2.4) and (1.2.5) can be considered with the two following view points.

(i)

We take a snap shot of the solution u(•, •) at a fixed time T > 0. We consider the partial function [0, +∞[ ∋ x -→ u(x, T ) ∈ IR and taking into account the relations (1.2.4), (1.2.5), (1.2.7) and (1.2.8), we have (1.2.9) u(x, T ) = sin ω T -x a , x < aT 0 , x > aT .

 and this function is illustrated on Figure 1.3.

(ii)

We fix a particular position X in space and we look, as time is increasing, to the solution u(•, •) at this particular point. We show on Figure 1.4 the function [0, +∞[ ∋ t -→ u(X, t) ∈ IR and taking into account the relations (1.2.4), (1.2.5), (1.2.7) and (1.2.8), we have (1.2.10) u(x, T ) = 0 , t < X a sin ω T -x a , t > X a .

x 0 u(x, T) a T 1.3 Inflow and outflow for the advection equation.

•

We still suppose that celerity a is positive and we consider the resolution of the advection (1.2.1) in the space-time domain (1.3.1) 0 < x < L , t > 0 . The relations (1.2.4) and (1.2.5) can still be applied because the proof of Proposition 1.2 remains unchanged in this particular case. As a consequence of the  previous property, we remark that no boundary condition is necessary at the particular position x = L for solving the advection problem in the space-time domain defined in relations (1.3.1). The initial condition (1.2.2) has simply to be restricted in domain ]0, L[ : (1.3.2) u(x, 0) = u 0 (x) , 0 < x < L , and the boundary condition (1.2.3) at x = 0 remains unchanged : (1.3.3) u(0, t) = v 0 (t) , t > 0 . • The difference between point x = 0 and point x = L for the resolution of the advection equation in space-time domain (1.3.1) is due to the fact that we choose an orientation of the characteristic lines xat = constant associated to an increase for the time direction. With this choice of time direction, the  characteristic lines enter inside the space-time domain (1.3.1) at x = 0 and they go outside at x = L . The boundary condition (1.3.3) is given at the input of the domain (see Figure 1.5) and at x = L , there is a free output from space time domain (1.3.1), without necessity to specify any numerical boundary condition.

•

If we change the sign of celerity a, i.e. if we suppose now (1.3.4) a < 0 , the above analysis remains unchanged, but the algebraic relations (1.2.4) and (1.2.5) have to be modified (see Figure 1.6). We still start from relation (1.1.5) that expresses that the solution of the advection equation (1.1.1) is constant along the characteristics lines. The foot of the characteristic line that contains the particular point (x, t) in space-time is either the point

(y = x -at, 0) if x -at < L , either the point L, θ = t -1 a (x -L) if x -at > L .
In the first case, we have y > 0 and θ < 0 then the initial condition (1.3.2) is advected inside the domain (1.3.1) and we have :

(1.3.5) u(x, t) = u 0 (x -at) , x -at < L . •
On the contrary, if xat > L, we have y > L and θ > 0 then the boundary condition at x = L that takes now the expression (1.3.6) u(L, t) = w L (t) , t > 0 , is advected inside the domain of study and we have :

(1.3.7) u(x, t) = w L t + L a - x a , x -at > L .
We have established the following Proposition 1.3. Advection in the domain 0 < x < L , a < 0. Under the hypothesis (1.3.4), the resolution of the advection equation (1.2.1) in the space-time domain (1.3.1) conducts to a well posed problem when we introduce the initial condition (1.3.2) on the interval ]0, L[ and the boundary condition (1.3.6) at the input region located at x = L , without any boundary condition at the output located at x = 0 . The solution of Problem (1.2.1), (1.3.2) and (1.3.6) is given by the relations (1.3.5) and (1.3.7).

2)

Finite volumes for linear hyperbolic systems. 2.1 Linear advection.

•

We still study the advection equation parameterized by some celerity a > 0 :

(2.1.1)

∂W ∂t + ∂ ∂x a W = 0 , t > 0 , x ∈ IR ,
and we search a discrete version of this mathematical model. For doing this, we introduce a space step ∆x > 0 and a space grid composed by points x j whose coordinates are multiples of this space step ∆x, id est

 (2.1.2)
x j = j ∆x , j ∈ Z Z . For a finite domain, ]0, L[ to fix the ideas, the above grid is limited to integer values j such that (2.1.3) 0 ≤ j ≤ J = L ∆x and the vertices (x j ) 0≤j≤J are usually used in the context of the finite difference method. The intervals K j+1/2 = ]x j , x j+1 [ between two vertices can be considered as finite elements (or finite volumes in our study) and they cover the entire domain ]0,

L[ : (2.1.4) [0, L] = 0≤j≤J-1 [x j , x j+1 ] ,
as proposed in the general context of meshes (see e.g. Ciarlet [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). We introduce also a time step ∆t > 0 and the discrete time values at integer multiples of the above quantum :

(2.1.5)

t n = n ∆t , n ∈ IN .
We consider now a space-time volume

V n+1/2 j+1/2
obtained by cartesian product of the two intervals ]x j , x j+1 [ and ]t n , t n+1 [ (see Figure 2.1) :

(2.1.6) V n+1/2 j+1/2 = ]x j , x j+1 [ × ]t n , t n+1 [ . •
The finite volume scheme consists simply in integrating the advection equation (2.1.1) inside the space-time domain V n+1/2 j+1/2 introduced previously :

(2.1.7) Finite volume scheme. Let IR × [0, +∞[ ∋ (x, t) -→ W (x, t) ∈ IR be a solution of the advection equation (2.1.1). We introduce the space mean value W n j+1/2 of this solution

V n+1/2 j+1/2 ∂W ∂t + ∂ ∂x a W dx dt = 0 , 0 ≤ j ≤ J , n ≥ 0 . W n+1 j+1/2 W n j+1/2 f n+1/2 j f n+1/2 j+1 x j+1 x j t n+1 t n V n+1/2 j+1/2 W n j-1/2 x j-1
W (•, •) in the cell K j+1/2 :  (2.1.8) W n j+1/2 = 1 | K j+1/2 | x j+1 x j W (x, t n ) dx
and the time mean value f n+1/2 j of the so-called flux a W (•, •) at the space position x j and between discrete times t n and t n+1 :

(2.1.9)

f n+1/2 j = 1 ∆t t n+1
t n (a W )(x j , t) dt . Then we have the following constitutive relation of finite volumes schemes :

(2.1.10)

1 ∆t W n+1 j+1/2 -W n j+1/2 + 1 ∆x f n+1/2 j+1 -f n+1/2 j = 0 .
This numerical modelling characterizes the so-called finite volume method which has been developed thanks to the work of S. Godunov [START_REF] Godunov | A Difference Method for the Numerical Computation of Discontinuous Solutions of the Equations of Fluid Dynamics[END_REF], Godunov et al [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF], Patankar [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF], Harten, Lax and Van Leer [START_REF] Harten | On Upstream Differencing and Godunov-type Schemes for Hyperbolic Conservation Laws[END_REF] or Faille, Gallouët and Herbin [START_REF] Faille | Les Mathématiciens découvrent les Volumes Finis[END_REF] among others.

•

The proof of Proposition 2.1 consists in a precise evaluation of the left hand side of equality (2.1.7). We use Fubini rule for the computation of double integrals and we begin by integrating in time for the ∂ ∂t term :

V n+1/2 j+1/2 ∂W ∂t dx dt = x j+1 x j t n+1 t n ∂W ∂t (x, t) dt dx = x j+1 x j W (x, t n+1 ) -W (x, t n ) dx = ∆x W n+1 j+1/2 -W n j+1/2
due to the definition (2.1.8). We proceed in an analogous way with the ∂ ∂x term and begin now the Fubini procedure by integrating in space ; we have

V n+1/2 j+1/2 ∂ ∂x a W dx dt = t n+1 t n x j+1 x j ∂ ∂x a W (x, t) dx dt = t n+1 t n (a W (x j+1 , t)-(a W (x j , t) dt = ∆t f n+1/2 j+1 -f n+1/2 j
according to the definition (2.1.9). We add the two previous results, use identity (2.1.7) and divide by ∆t ∆x. We obtain exactly the relation (2.1.10).

•

The relation (2.1.10) is a very general form for the evolution of the mean values W j+1/2 between two time steps. The increment

W n+1 j+1/2 -W n j+1/2
is, after correction by a multiplicative factor, equilibrated by the flux difference

f n+1/2 j+1 -f n+1/2 j
. The idea of a finite volume scheme is to consider now that the algebraic object W j+1/2 is nomore the mean value of the exact solution but an approximation of this mean value. Then the relation (2.1.10) proposes  a numerical scheme for the discrete evolution of the approximated mean values W j+1/2 , j = 0, • • • , J -1. Nevertheless, the numerical scheme is not entirely defined by the relation (2.1.10). Starting from mean values at the initial time step, i.e.

(2.1.11)

W 0 j+1/2 = 1 ∆x x j+1 x j W 0 (x) dx , j = 0, • • • , J -1 ,
we are able to increment the time step with relation (2.1.10) only if all the fluxes f n+1/2 j , j = 0, • • • , J have been a priori first determined as a functional of the previous values. In a very general way, we say that the finite volume scheme (2.1.10) is an explicit scheme if each flux f n+1/2 j is a given function Ψ j of the mean values W n k+1/2 k=1,•••, J-1 at the preceding time step number n : (2.1.12)

f n+1/2 j = Ψ j {W n k+1/2 , k = 0, • • • , J -1} , j = 0, • • • , J -1 .
The function Ψ j is called the local numerical flux function at point x j and, joined with the evolution equation (2.1.10), its choice determines the numerical scheme.

•

A natural hypothesis claims that we have translation invariance for the evaluation of the flux if we move the discrete data in the same way ; in other words, the numerical flux function Ψ j only depends on the p first neighbors of the interface x j . Then the explicit numerical flux is a given function Φ of the p first neighbors and we have : (2.1.13)

f n+1/2 j = Φ W n j+1/2-p , • • • , W n j-1/2 , W n j+1/2 , • • • , W n j+1/2+p-1 .
A very important particular case is one of a two-point scheme for the evaluation of the numerical flux. We have in this particular case : (2.1.14)

f n+1/2 j = Φ W n j-1/2 , W n j+1/2
. With this particular choice, the numerical scheme for incrementing in time of the mean values takes the form :

(2.1.15)

     1 ∆t W n+1 j+1/2 -W n j+1/2 + + 1 ∆x Φ W n j+1/2 , W n j+3/2 -Φ W n j-1/2 , W n j+1/2 = 0 .
It is also a three-point finite difference scheme. The finite volume scheme (2.1.10) (2.1.13) is said to be consistent with the advection equation (2.1.1) when the numerical flux function Φ satisfies the condition

(2.1.16) Φ W, • • • , W, W, • • • , W = a W , ∀ W ∈ IR . •
The crucial question is how to choose a numerical finite volume scheme. The simplest choice consists in a two point explicit scheme such that the finite difference scheme is identical to the upstream-centered scheme (see e.g. Richtmyer-Morton [START_REF] Richtmyer | Difference Methods for Initial-Value Problems[END_REF]). It takes the following expressions :

(2.1.17)

1 ∆t W n+1 j+1/2 -W n j+1/2 + a W n j+1/2 -W n j-1/2 = 0 , a > 0 (2.1.18) 1 ∆t W n+1 j+1/2 -W n j+1/2 + a W n j+3/2 -W n j+1/2 = 0 , a < 0 .
The corresponding flux function is called the first order upstream-centered flux, is simply given by the following relations :

(2.1.19) Φ(W l , W r ) = a W l , a > 0 a W r , a < 0 . When this flux function acts at a given point x j of the mesh, we have :

(2.1.20) f n+1/2 j = Φ W n j-1/2 , W n j+1/2 = a W n j-1/2 , a > 0 a W n j+1/2
, a < 0 . If a > 0, the exact solution of the advection equation propagates the information from the left to the right ; the flux at the interface x j is issued from the cell at the left of the interface and this cell at the number j -1/2 . If a < 0, the propagation of the information with the advection equation is from right to left ; the interface flux at the abscissa x j is due to the control volume on the right, i.e. with number j +1/2 as depicted on Figure 2.2.

x j W n j-1/2 a > 0 W n j+1/2 a < 0
x = 0 1 x ... j-1 x j x ... J-1 x J x = L W 1/2 f 0 f J W J-1/2 W j-1/2 j-1/2 x

•

Recall that practical use of the upwind finite volume scheme like (2.1.17) when a > 0 or (2.1.18) if a < 0 is restricted to the usual Courant-Friedrichs-Lewy stability condition :

( J ∆x = L . Note that the j th cell is exactly the interval ]x j-1 , x j [ and it is centered at point x j-1/2 as shown on are not a priori defined because states W n -1/2 or W n J+1/2 does not exist. The situation is more complex with numerical fluxes that use four points or more as proposed in (2.1.13) and will not be detailed in this section. Even if the formula giving the numerical flux at the boundaries has to be specifically studied, the finite volume scheme remains defined by the relation (2.1.10) and we have for the two cells encountering the boundary :

(2.2.3) 1 ∆t W n+1 1/2 -W n 1/2 + 1 ∆x f n+1/2 1 -f n+1/2 0 = 0 , (2.2.4) 1 ∆t W n+1 J-1/2 -W n J-1/2 + 1 ∆x f n+1/2 J -f n+1/2 J-1 = 0 .
• The question is now to adapt the relation (1.2.14) in order to determine the two boundary fluxes f n+1/2 0 at the left of the domain and f n+1/2 J at the right. For the advection equation with celerity a > 0, we have observed in the first section that some boundary condition v 0 (t) has to be assigned at x = 0 and it is not the case for x = L. It is therefore natural to take into account this information at the input of the domain and to set :

(2.2.5)

f n+1/2 0 = 1 ∆t t n+1 t n a v 0 (t) dt  or simply (2.2.6) f n+1/2 0 = a v 0 (n + 1 2 )∆t , a > 0 , if function t -→ v 0 (t)
has a slow time variation at the scale defined by the time step. At the output x = L, no numerical datum has to be assigned to set correctly the continuous mathematical problem. We must maintain this property if we wish the numerical method to follow the mathematical physics as efficiently as possible. A simple boundary flux is associated with the previous numerical upwind scheme. For x = x J = L and a > 0, we observe that the upwind scheme (2.1.20) is simply written as : (2.2.7)

f n+1/2 J = a W n J-1/2
, a > 0 , and this relation (2.2.7) defines a first order extrapolated boundary flux.

•

The roles are reversed when a < 0. The abscissa x = 0 corresponds to an output for the advection equation and the right boundary x = L is an input where a time field t -→ w L (t) is given. In the first case, the upwind scheme (2.1.20) can be applied without modification : (2.2.8)

f n+1/2 0 = a W n 1/2
, a < 0 , and it corresponds to a first order extrapolation of the internal data

W n j-1/2 , j = 1, • • • , J
at the boundary at time step n∆t. For x = L, the boundary flux f n+1/2 J uses the given information between the two time steps : (2.2.9)

f n+1/2 J = a w L (n + 1 2 )∆t , a < 0 . Proposition 2.2.
Flux boundary conditions for the advection equation. When we approach the advection equation (2.1.1) with the finite volume method, the numerical boundary conditions induces a choice for the two boundary fluxes f n+1/2 0 and f n+1/2 J . When a > 0, the boundary condition v 0 (t) at the input can be introduced into the boundary with the relation (2.2.6) and the free output at the right can be treated with an extrapolation of the type (2.2.7). When a < 0, the free output at the left of the domain can be taken into account with the help of relation (2.2.8) whereas the input condition w L (t) at the right can be introduced thanks to relation (2.2.9).

A model system with two equations

•

Let a > 0 and b > 0 be two positive real number. We study in this section a model problem that is composed by the juxtaposition of an advection equation with celerity a and an advection with celerity -b. We explicit the associated algebra :

(2.3.1) ∂u ∂t + a ∂u ∂x = 0 , t > 0 , x ∈ IR , (2.3.2) ∂v ∂t -b ∂v ∂x = 0 , t > 0 , x ∈ IR .
We associate the two equations (2.3.1) and (2.3.2) and consider a unique problem with a vector field as unknown. We set :

(2.3.3) ϕ = u v and the set of equations (2.3.1)-(2.3.2) can naturally be written as a system :

(2.3.4) ∂ϕ ∂t + a 0 0 -b ∂ϕ ∂x = 0 .
By introducing the flux function F (ϕ) according to the relation

(2.3.5) F (ϕ) = a u -b v the system (2.3.4) takes the general conservative form : (2.3.6) ∂ϕ ∂t + ∂ ∂x F (ϕ) = 0 .
• The approximation of system (2.3.6) with a grid parameterized by a space step ∆x and a time step ∆t is conducted exactly as in the case of the advection equation. The following property is a straightforward generalization of Proposition 2.1. We left the proof to the reader. Proposition 2.3.

Finite volume scheme. Let IR × [0, +∞[ ∋ (x, t) -→ ϕ(x, t) ∈ IR × IR be a solution of the linear conservation law (2.3.6). We define the space mean value ϕ n j+1/2 of this solution

ϕ(•, •) in the cell K j+1/2 : (2.3.7) ϕ n j+1/2 = 1 | K j+1/2 | x j+1 x j ϕ(x, t n ) dx
and the time mean value f n+1/2 j of the flux function introduced in (2.3.5) at the space position x j between discrete times t n and t n+1 :

(2.3.8) f n+1/2 j = 1 ∆t t n+1 t n F ϕ(x j , t) dt .
We have the following relation that characterizes the finite volumes schemes :

(2.3.9)

1 ∆t ϕ n+1 j+1/2 -ϕ n j+1/2 + 1 ∆x f n+1/2 j+1 -f n+1/2 j = 0 .

•

We have now to propose a precise numerical flux function analogous to the relation (2.1.12) to transform the conservation property (2.3.9) into a finite volume numerical scheme able to propagate the discrete values ϕ n j+1/2 up to the discrete time t n+1 . For internal interfaces x j , j = 1, • • • , J -1 , it is natural to apply the upwinding scheme (2.1.20) with a left upwinding for the first equation and a right upwinding for the equation (2.3.2). Figure 2.4 illustrates the associated algebra : • At the left boundary x = 0 , we have an input for the variable u and we suppose given the associated datum

(2.3.10) f n+1/2 j = Φ ϕ n j-1/2 , ϕ n j+1/2 = a u n j-1/2 -b v n j+1/2 , j = 1, • • • , J -1 . u n j+1/2 v n j-1/2 x j -b < 0 a > 0 u n j-1/2 v n j+1/2
[0, +∞[ ∋ t -→ u 0 (t) ∈ IR : (2.3.11) u(0, t) = u 0 (t) ,
t > 0 whereas it is an output for the v variable. By association of relations (2.2.6) and (2.2.8), we obtain (2.3.12)

f n+1/2 0 = a u 0 (n + 1 2 )∆t -b v n 1/2 .
At the other boundary of the interval ]0, L[ , we have an output for the first variable u and an input for the second one, and an associated boundary condition

[0, +∞[ ∋ t -→ v L (t) ∈ IR is supposed to have been given : (2.3.13) v(L, t) = v L (t) , t > 0  as illustrated on Figure 2.5.
The numerical flux at the right is evaluated by association of the relations (2.2.7) and (2.2.9) :

(2.3.14)

f n+1/2 L = a u n J-1/2 -b v L (n + 1 2 )∆t
.

Unidimensional linear acoustics

• We consider a gas in a pipe of uniform section at normal conditions of temperature and pressure. The reference density is denoted by ρ 0 and the reference pressure is named p 0 . The sound celerity c 0 of this gas satisfies the relation (2.4.1) c 0 = γp 0 ρ 0 with γ = 1.4 as proved e.g. in the book of Landau and Lifchitz [START_REF] Landau | Fluid Mechanics[END_REF]. A sound wave is a small perturbation of this reference state. The differences of density, pressure and velocity fields are denoted respectively by ρ, p and u. The hypothesis of a small perturbation implies that the entropy of the reference state is maintained for all the time evolution and in consequence, it is easy to establish the following relation between the perturbations of density and pressure : (2.4.2) p = c 2 0 ρ .

•

The conservation of mass leads to a first order linear conservation law : (2.4.3) ∂ρ ∂t + ρ 0 ∂u ∂x = 0 and the conservation of momentum links the time evolution of velocity with the spatial gradient of pressure :

(2.4.4)

ρ 0 ∂u ∂t + ∂p ∂x = 0 .
We introduce the vector W = p u of unknowns. Then the equations (2.4.3) and (2.4.4) can be written as a linear hyperbolic system of conservation laws :

(2.4.5)

∂W ∂t + A ∂W ∂x = 0 with (2.4.6) A = 0 ρ 0 c 2 0 1 ρ 0 0 .
• When we consider the eigenvalues and eigenvectors of matrix A, it is natural to introduce the characteristic variables defined respectively by (2.4.7)

ϕ + = p + ρ 0 c 0 u (2.4.8)
ϕ -= pρ 0 c 0 u and the quantity ρ 0 c 0 is named the acoustic impedance. We have from the relations (2.4.3) and (2.4.4) :

∂ϕ + ∂t + c 0 ∂ϕ + ∂x = ∂p ∂t + ρ 0 c 0 ∂u ∂t + c 0 ∂p ∂x + ρ 0 c 2 0 ∂u ∂x = c 2 0 ∂ρ ∂t + ρ 0 ∂u ∂x + c 0 ρ 0 ∂u ∂t + ∂p ∂x = 0 , ∂ϕ - ∂t -c 0 ∂ϕ - ∂x = ∂p ∂t -ρ 0 c 0 ∂u ∂t -c 0 ∂p ∂x -ρ 0 c 0 ∂u ∂x = c 2 0 ∂ρ ∂t + ρ 0 ∂u ∂x -c 0 ρ 0 ∂u ∂t + ∂p ∂x = 0 ,
and we recover a system of the type (2.3.4) studied previously :

(2.4.9) ∂ ∂t

ϕ - ϕ + + -c 0 0 0 c 0 ∂ ∂t ϕ - ϕ + = 0 .
• A typically physical problem is the following : a given acoustic pressure wave [0, +∞[ ∋ t -→ Π(t) > 0 is injected at the left x = 0 of the pipe and the waves go away freely at the right boundary x = L . At t = 0, the velocity and pressure of the fluid are given : (2.4.10)

u(x, 0) = u 0 (x) , 0 < x < L (2.4.11) p(x, 0) = p 0 (x) , 0 < x < L .
From a mathematical viewpoint, the boundary conditions have to respect the dynamics of this system of acoustic equations written in diagonal form (2.4.9) : the variable ϕ + must be given at x = 0 and the variable ϕ -at the abscissa x = L. From (2.4.7) and (2.4.8), we determine the pressure as a function of the two characteristics variables ϕ + and ϕ -:

(2.4.12) p = 1 2 ϕ + + ϕ - and if the pressure is imposed at x = 0, the relation (2.4.12) can be written under the form :

(2.4.13) ϕ + (0, t) =ϕ -(0, t) + 2 Π(t) , x = 0 , t > 0 , that makes in evidence a reflection operator : the input variable ϕ + is a given affine function of the output variable ϕ -. At the other boundary x = L , the notion of free output expresses that the waves that go outside of the domain of study have no reflection at the boundary. When x = L , the characteristic variable ϕ + is going outside and there is no boundary condition for this variable. We have to express also that this wave has no influence on the characteristic ϕ - that wish to go inside the domain ]0, L[. In other terms, the input value ϕ -is independent of the variable ϕ + and also of time. We have in consequence

(2.4.14) ∂ ∂t ϕ -(L, t) = 0 .
We have established  Proposition 2.4. Boundary conditions for acoustic problem. The mathematical boundary conditions associated with the datum of a given acoustic pressure wave [0, +∞[ ∋ t -→ Π(t) > 0 at the left of the domain ]0, L[ admits the expression (2.4.13) and a condition of free output of the waves at the right boundary x = L can be expressed by the relation (2.4.14).

0 L x L c 0 p = Π(t) t ϕ (L , t) = ϕ (L , 0) - - Figure 2
.6. Solution of the acoustic equations in one space dimension for a model problem with two equations

•

The above acoustic problem associated with the first order partial differential equations (2.4.3) (2.4.4), the initial conditions (2.4.10) (2.4.11) and the boundary conditions (2.4.13) (2.4.14) is illustrated on Figure 2.6. The initial conditions are active in the beginning of the evolution in time (t ≤ L c 0 ) and have a trace for higher times due to the boundary conditon (2.4.13), that gives, due to (2.4.8) and (2.4.13) :

(2.4.15) ϕ -(x, t) ≡ p(x, t) -ρ 0 c 0 u(x, t) = p 0 (L) -ρ 0 c 0 u 0 (L) , t ≥ L c 0 .
On the other hand, the inflow boundary condition (2.4.12) and the second row of matrix equation (2.4.9) implies :

(2.4.16)

ϕ + (x, t) ≡ p(x, t) + ρ 0 c 0 u(x, t) = = 2 Π t - x c 0 -ϕ -0, t - x c 0 , t ≥ L c 0 .
We deduce from the relations (2.4.15) (2.4.16) joined with the definitions (2.4.7) and (2.4.8) :

(2.4.17

) p(x, t) = Π t -x c 0 , 0 ≤ x ≤ L , t ≥ L c 0 (2.4.18) u(x, t) = u 0 (L) + 1 ρ 0 c 0 Π t-x c 0 -p 0 (L) , 0 ≤ x ≤ L , t ≥ L c 0 . •
We turn now to the numerical finite volume scheme. We have to determine the internal fluxes f

F (W ) = ρ 0 c 2 0 u 1 ρ 0 p with W = p u .
Proposition 2.5. Upwind scheme for computational acoustics. The extension of the upwind finite volume scheme (2.3.10), (2.3.12) and (2.3.14) is determined by the following relations :

(2.4.20)

f n+1/2 j = ρ 0 c 2 0 2 u n j-1/2 + u n j+1/2 -c 0 2 p n j+1/2 -p n j-1/2 1 2 ρ 0 p n j-1/2 + p n j+1/2 -c 0 2 u n j+1/2 -u n j-1/2
for the internal fluxes, i.e. for indexes j that satisfy 1 ≤ j ≤ J -1 . The two boundary fluxes follow the following relations :

(2.4.21)

f n+1/2 0 = ρ 0 c 2 0 u n 1/2 + c 0 Π (n + 1 2 )∆t -p n 1/2 1 ρ 0 Π (n + 1 2 )∆t (2.4.22) f n+1/2 J = ρ 0 c 2 0 2 u n J-1/2 + u 0 J-1/2 -c 0 2 p 0 J-1/2 -p n J-1/2 1 2 ρ 0 p n J-1/2 + p 0 J-1/2 -c 0 2 u 0 J-1/2 -u n J-1/2
.

•

The internal fluxes are determined with the scheme (2.3.10) applied with the diagonal form of relation (2.4.9). We have (2.4.23)

ϕ n+1/2 +, j = ϕ n +, j-1/2 ≡ p n j-1/2 + ρ 0 c 0 u n j-1/2 (2.4.24) ϕ n+1/2 -, j = ϕ n -, j+1/2 ≡ p n j+1/2 -ρ 0 c 0 u n j+1/
= 2Π (n + 1 2 )∆t -ϕ n+1/2 -, 0
and use a first order extrapolation of the outgoing characteristic variable :

(2.4.26)

ϕ n+1/2 -, 0 = ϕ n -, 1/2 ≡ p n 1/2 -ρ 0 c 0 u n 1/2 .
ϕ n+1/2 -, J = ϕ 0 -, J ≡ p 0 (L) -ρ 0 c 0 u 0 (L) ≈ p 0 J-1/2 -ρ 0 c 0 u 0 J-1/2
and the output characteristic variable is extrapolated from the interior of the domain :

(2.4.28)

ϕ n+1/2 +, J = ϕ n +, J-1/2 ≡ p n J-1/2 + ρ 0 c 0 u n J-1/2
. The relation (2.4.22) follows after two steps of elementary algebra.

•

We remark that both relations (2.4.20) and (2.4.22) are identical, except that the boundary state W 0 (L) ≈ W 0 J-1/2 has replaced the right state W n j+1/2 . Moreover the flux boundary condition (2.4.21) that involves the pressure is a natural discretization of the exact characteristic solution (2.4.17) (2.4.18) at x = 0 .

Characteristic variables.

•

We suppose now to fix the ideas that the unknown vector

W (•, •) (2.5.1) [0, L] × [0, +∞[ ∋ (x, t) -→ W (x, t)
∈ IR 3 has three real components w 1 , w 2 and w 3 . We suppose also that the function

W (•, •) is solution of a conservation law of the type (2.5.2) ∂W ∂t + ∂ ∂x F (W ) = 0
where the flux F (W ) is a linear function of vector W :

(2.5.3)

F (W ) = A • W and A is a 3 by 3 diagonalizable real matrix. •
We first detail the fact that matrix A is a diagonalizable matrix. There exists three non null real vectors r 1 , r 2 , r 3 and three real scalars λ 1 , λ 2 , λ 3 in such a way that (2.5.4) A • r j = λ j r j , j = 1, 2, 3. ¿From a matricial viewpoint, we denote by R k j the k 0 component of the eigenvector r j , i.e.

(2.5.5)

r j =   R 1 j R 2 j R 3 j   ≡   r j 1 r j 2 r j 3  
and we introduce the 3 by 3 matrix R composed by the scalars R k j . The vector r j is the k 0 column of matrix R. The relation (2.5.4) can also be written as (2.5.6) andΛ is the diagonal matrix whose diagonal terms are equal to the eigenvalues

A • R = R • Λ ,
λ j : (2.5.7) Λ =   λ 1 0 0 0 λ 2 0 0 0 λ 3   .

•

We consider now two distinct bases for linear space IR 3 : on one hand the canonical basis e j j=1, 2, 3 defined by (2.5.8)

e 1 =   1 0 0   , e 2 =   0 1 0   , e 3 =   0 0 1   
where the vector W admits the natural decomposition introduced above :

(2.5.9) W = k=3 k=1 w k e k , and on the other hand the basis of IR 3 composed by the eigenvectors (r j ) j=1, 2, 3 . In the latter, the vector W can be decomposed with a formula of the type (2.5.10) W = j=3 j=1 ϕ j r j and the scalar ϕ j define the characteristic variables associated with the system (2.5.2) (2.5.3). The link between the relations (2.5.9) and (2.5.10) is classical : we consider the components R k j of vector r j inside the canonical basis and we get from the relation (2.5.5) :

(2.5.11)

w k = j=3 j=1 ϕ j R k j .
Then the relation (2.5.11) can be re-written under a matricial form :

(2.5.12) W = R • ϕ .

•

The relation (2.5.12) proposes to change the unknown function, i.e. to replace the research of W (x, t) ∈ IR 3 by the equivalent research of the characteristic vector ϕ(x, t) ∈ IR 3 and defined by :

(2.5.13)

ϕ = R -1 • W .
Proposition 2.6. Characteristic variables satisfy advection equations.

The vector [0, L]×[0, +∞[ ∋ (x, t) -→ ϕ(x, t) ∈ IR 3 of characteristic variables satisfy the matrix equation

(2.5.14) ∂ϕ ∂t + Λ • ∂ϕ ∂x = 0
that takes also the equivalent scalar form :

(2.5.15) ∂ϕ j ∂t + λ j ∂ϕ j ∂x = 0 , j = 1, 2, 3 .

•

We have from (2.5.2), (2.5.3), (2.5.6) and (2.5.12) :

∂W ∂t + A ∂W ∂x = R • ∂ϕ ∂t + A • R • ∂ϕ ∂x = R • ∂ϕ ∂t + R -1 • A • R • ∂ϕ ∂x = R • ∂ϕ ∂t + Λ • ∂ϕ ∂x = 0 ,
and since the matrix R is invertible, we deduce from the previous calculus the relation (2.5.14). The relation (2.5.15) is an immediate consequence of (2.5.14) and (2.5.7).

λ 2 λ 3 λ 1 λ 3 λ 2 λ 1 x = 0 x = L t x
Figure 2.7. Linear hyperbolic system with three equations and eigenvalues satisfying λ 1 < 0 < λ 2 < λ 3 .

•

To fix the ideas, we suppose that the eigenvalues λ j of matrix A are distinct, enumerated with an increasing order and with distinct signs as illustrated on Figure 2.7 :

(2.5.16)

λ 1 < 0 < λ 2 < λ 3 .
The propagation of the first variable ϕ 1 goes from right to left (because λ 1 < 0 ) with celerity | λ 1 |, the second characteristic variable ϕ 2 from left to right with celerity λ 2 and the same property holds for variable ϕ 3 with eigenvalue λ 3 .

• A set of well posed boundary conditions is a consequence of the diagonal form (2.5.15) of the equations and of the particular choice (2.5.16) for the signs. The directions associated with eigenvalues λ 2 and λ 3 are ingoing at x = 0 and we have to give some boundary condition for ϕ 2 and ϕ 3 at this point :

(2.5.17)

ϕ 2 (x = 0, t) = β 0 (t) (2.5.18) ϕ 3 (x = 0, t) = γ 0 (t) .
The direction associated with the eigenvalue λ 1 is ingoing at the abscisssa x = L, and this condition imposes to have some datum concerning ϕ 1 at this particular point :

(2.5.19)

ϕ 1 (x = L, t) = α L (t) .
The previous boundary conditions (2.5.17) to (2.5.19) define a well posed problem. Nevertheless, the introduction of physically relevant boundary conditions (as a pressure condition as seen in the previous section) requires a more general formulation of the boundary condition. In the linear case, the stability study developed by Kreiss [START_REF] Kreiss | Initial Boundary Value Problems for Hyperbolic Systems[END_REF] shows that the ingoing characteristic can be an affine function of the outgoing characteristic through a reflection operator at the boundary. We can explicit the former with the above example. • At x = 0 , the first characteristic is outgoing and the two last ones are going inside the domain of study. Then we can replace the conditions (2.5.17) and (2.5.18) by the following ones :

(2.5.20) ϕ 2 (x = 0, t) = β 0 (t) + p(t) ϕ 1 (x = 0, t) (2.5.21) ϕ 3 (x = 0, t) = γ 0 (t) + q(t) ϕ 1 (x = 0, t) , where t -→ p(t) and t -→ q(t) are given fixed real functions of time. The conditions (2.5.20) and (2.5.21) are illustrated on Figure 2.8. We can also write them (2.5.22)

ϕ in (x, t) = g(t) + S(t) • ϕ out (x, t) , x point on the boundary, with ϕ in = ϕ 2 ϕ 3 , g(t) = β 0 (t) γ 0 (t) , S(t) = p(t) q(t) , ϕ out = ϕ 1 when x = 0 . • When x = L, the relation (2.5.19
) is replaced by a more general one (2.5.23) 2.9 and including an affine component of the outgoing characteristic variables. The boundary condition (2.5.23) takes again a form of  the type (2.5.22) with this time the following relations :

ϕ 1 (x = L, t) = α L (t) + θ(t) ϕ 2 (x = L, t) + σ(t) ϕ 3 (x = L, t) illustrated on Figure
ϕ in = ϕ 1 , g(t) = α L (t) , S(t) = θ(t) σ(t) , ϕ out = ϕ 2 ϕ 3 when x = L .
2.6 A family of model systems with three equations

•

We still study a 3 by 3 linear hyperbolic system of the type (2.5.2) (2.5.3) with the condition (2.5.16) to fix a particular example. We suggest in this section to explicit a way for evaluation of the numerical flux f n+1/2 j that is the key point for the discrete evolution in time of the mean values W j+1/2 :

(2.6.1)

1 ∆t W n+1 j+1/2 -W n j+1/2 + 1 ∆x f n+1/2 j+1 -f n+1/2 j = 0 .

•

We change the notations and wish to determine the numerical flux Φ(W l , W r ) for W l = W left and W r = W right given respectively at the left and at the right of the interface (see Figure 2.10). When we consider the advection equation (and in that case the variables W l and W r are real numbers) the relation (2.1.19) gives the result : Φ(W l , W r ) = a W l when a > 0 and Φ(W l , W r ) = a W r when a < 0 . We have to generalize this study when the field W (•, •) is three-dimensional. We first decompose the vector Φ(W l , W r ) with the basis r j of eigenvectors and introduce its (scalar) components

ψ j (W l , W r ) : (2.6.3) Φ(W l , W r ) = j=3 j=1 ψ j (W l , W r ) r j i.e. (2.6.4) Φ k (W l , W r ) = j=3 j=1 R k j ψ j (W l , W r
) . For j = 1, we have λ 1 < 0 then the numerical scheme has to be upwinded in the right direction :

(2.6.5) ψ 1 (W l , W r ) = λ 1 ϕ 1, r whereas for j = 2 or j = 3, we have λ 2 > 0 and λ 3 > 0 and the scheme must be upwinded to the left. It comes (2.6.6)

ψ 2 (W l , W r ) = λ 2 ϕ 2, l , ψ 3 (W l , W r ) = λ 3 ϕ 3, l .
In consequence of the relations (2.6.3) to (2.6.6), the numerical flux function Φ(•, •) can be written globally : (2.6.7) Φ(W l , W r ) = λ 1 ϕ 1, r r 1 + λ 2 ϕ 2, l r 2 + λ 3 ϕ 3, l r 3 , or in an equivalent way with introducing the Cartesian components :

(2.6.8) Φ k (W l , W r ) = λ 1 ϕ 1, r R k 1 + λ 2 ϕ 2, l R k 2 + λ 3 ϕ 3, l R k 3 , k = 1, 2, 3 .

•

We can also re-write the relation (2.6.8) for the particular interface x j : (2.6.9)

W l = W left = W n j-1/2 , W r = W right = W n j+1/2
. We first decompose the vector W on the eigenvectors of matrix A as in (2.5.11) :

(2.6.10)

W n j+1/2 k = i=3 i=1 ϕ n i, j+1/2 R k i , k = 1, 2, 3 , j = 1, • • • , J-1 , then we introduce the component number k of the flux f n+1/2 j , i.e. (f n+1/2 j ) k = Φ k (W n j-1/2 , W n j+1/2 ) at the interface x j : (2.6.11) f n+1/2 j k = λ 1 ϕ n 1, j+1/2 R k 1 + λ 2 ϕ n 2, j-1/2 R k 2 + λ 3 ϕ n 3, j-1/2 R k 3 . •
We detail in this sub-section the determination of the numerical flux f n+1/2 0 at the boundary x = 0. We first recall that the continuous boundary conditions at this point take the form given in (2.5.20) (2.5.21). The idea is to try to apply the upwind scheme (2.6.11) at the particular vertex j = 0 :

f n+1/2 0 = λ 1 ϕ n
1, 1/2 r 1 + + λ 2 ϕ n 2, -1/2 r 2 + λ 3 ϕ n 3, -1/2 r 3 and then to replace the characteristic values ϕ n 2, -1/2 and ϕ n 3, -1/2 (that are not defined on the mesh) by their values evaluated after a rough discretization of relations (2.5.20) and (2.5.21) :

ϕ n 2, -1/2 = β n+1/2 0 + + p n+1/2 ϕ n 1, 1/2 , ϕ n 3, -1/2 = γ n+1/2 0 + q n+1/2 ϕ n 1, 1/2 .
We obtain in consequence the following expression for the boundary flux at x = 0 :

(2.6.12)

f n+1/2 0 = λ 1 ϕ n 1, 1/2 r 1 + λ 2 β n+1/2 0 + p n+1/2 ϕ n 1, 1/2 r 2 + + λ 3 γ n+1/2 0 + q n+1/2 ϕ n
1, 1/2 r 3 or in an equivalent way :

(2.6.13)

f n+1/2 0 = ϕ n 1, 1/2 λ 1 r 1 + λ 2 p n+1/2 r 2 + λ 3 q n+1/2 r 3 + + λ 2 β n+1/2 0 r 2 + λ 3 γ n+1/2 0 r 3 .

•

The determination of the boundary flux f n+1/2 J can be conducted in the same way. Starting from the expression of the upwind scheme (2.6.11) when j = J , i.e. formally f

n+1/2 J = λ 1 ϕ n 1, J+1/2 r 1 + λ 2 ϕ n 2, J-1/2 r 2 +
λ 3 ϕ n 3, J-1/2 r 3 , we replace the first characteristic variable that appears external of the domain by its value given by the boundary condition (2.5.23) :

ϕ n 1, J+1/2 = α n+1/2 L + θ n+1/2 ϕ n 2, J-1/2 + + σ n+1/2 ϕ n 3, J-1/2
. We deduce :

(2.6.14)

f n+1/2 J = λ 1 α n+1/2 L + θ n+1/2 ϕ n 2, J-1/2 + σ n+1/2 ϕ n 3, J-1/2 r 1 + λ 2 ϕ n
2, J-1/2 r 2 + λ 3 ϕ n 3, J-1/2 r 3 or in an equivalent manner :

(2.6.15)

f n+1/2 J = λ 1 α n+1/2 L r 1 + ϕ n 2, J-1/2 λ 1 θ n+1/2 r 1 + λ 2 r 2 + + ϕ n 3, J-1/2 λ 1 σ n+1/2 r 1 + λ 3 r 3 .

First order upwind-centered finite volumes

•

We consider now a general system of conservation laws (2.7.1)

∂W ∂t + ∂ ∂x F (W ) = 0 with an unknown vector W (•, •) that belongs to linear space IR m : (2.7.2) [0, L] × [0, +∞[ ∋ (x, t) -→ W (x, t) ∈ IR m and a linear flux function F (•) (2.7.3) F (W ) = A • W
associated with a diagonalizable matrix A with eigenvalues λ j and eigenvectors r j (2.7.4)

A • r j = λ j r j , j = 1, 2, • • • , m .
Introducing the m × m matrix R as in relation (2.5.5) and the diagonal matrix Λ of eigenvalues as in relation (2.5.7), we have : (2.7.5)

A • R = R • Λ .

•

We propose here to determine a first order upwind flux Φ(W l , W r ) between the two states W left = W l and W right = W r that generalizes the relation (2.6.7) when we have not done any hypothesis of the type (2.5.16) concerning the sign of the eigenvalues λ j . We decompose any state W on the basis of space IR m characterized by the eigenvectors r j : (2.7.6) W = j=m j=1 ϕ j r j , W l = j=m j=1 ϕ j, l r j , W r = j=m j=1 ϕ j, r r j , and due to the structure introduced at Proposition 2.6, we obtain an advection equation for the j o characteristic variable ϕ j :

(2.7.7)

∂ϕ j ∂t + λ j ∂ϕ j ∂x = 0 , j = 1, 2, • • • , m .
Therefore it is natural to introduce the components ψ j (W l , W r ) of the numerical flux on the basis of the eigenvectors :

(2.7.8) Φ(W l , W r ) = j=3 j=1 ψ j (W l , W r ) r j Proposition 2.7. Three expressions of the upwind first order scheme. Let Φ(W l , W r ) the upwind flux defined by the relations (2.7.8) and (2.7.9). Then we have the three equivalent expressions : (2.7.15)

Φ(W l , W r ) = F (W l ) + j=m j=1 λ - j ϕ j, r -ϕ j, l r j (2.7.16) Φ(W l , W r ) = F (W r ) -j=m j=1 λ + j ϕ j, r -ϕ j, l r j (2.7.17) Φ(W l , W r ) = 1 2 F (W l ) + F (W r ) -1 2 | A | • (W r -W l . •
We write the relation (2.7.9) under the form : (2.7.18) ψ j (W l , W r ) = λ + j ϕ j, l + λ - j ϕ j, r and we have :

Φ(W l , W r ) = j=m j=1 λ + j ϕ j, l + λ - j ϕ j, r r j = j=m j=1
(λ jλ - j ) ϕ j, l + λ - j ϕ j, r r j due to (2.7.11) = j=m j=1 λ j ϕ j, l r j + j=m j=1 λ - j ϕ j, rϕ j, l r j and the relation (2.7.15) is established. In an analogous way, we have : Φ(W l , W r ) = j=m j=1 λ + j ϕ j, l + λ - j ϕ j, r r j = j=m j=1 λ + j ϕ j, l + (λ jλ + j ) ϕ j, r r j due to (2.7.11) Φ(W l , W r ) = j=m j=1 λ j ϕ j, r r j -j=m j=1 λ + j ϕ j, rϕ j, l r j and the relation (2.7.16) holds. We remark that

 | A | • (W r -W l ) = R • | Λ | • R -1 • R • (ϕ r -ϕ l ) due to (2.7.14) and (2.5.12) = R • | Λ | • (ϕ r -ϕ l ) = k=m k=1 j=m j=1 R k j | λ j | (ϕ j, r -ϕ j, l ) e k then (2.7.19) | A | • (W r -W l ) = j=m j=1 | λ j | (ϕ j, r -ϕ j, l ) r j .
We add the previous results (2.7.15) with (2.5.16), and we divide by two. We obtain :

Φ(W l , W r ) = 1 2 F (W l ) + F (W r ) -1 2 j=m j=1 λ + j -λ - j ϕ j, r -ϕ j, l r j = 1 2 F (W l ) + F (W r ) -1 2 j=m j=1 | λ j | ϕ j, r -ϕ j, l r j due to (2.7.12) = 1 2 F (W l ) + F (W r ) -1 2 | A | • (W r -W l ) due

3)

Gas dynamics with the Roe method. 3.1 Nonlinear acoustics in one space dimension.

•

We propose here to describe quickly a physical problem that comes from the theoretical modelling of trombone, detailed for instance in the work of Hirschberg et al [START_REF] Hirschberg | Shock waves in trombones[END_REF] or in our study [START_REF] Msallam | Mathematical model for coupling a quasi-unidimensional perfect flow with an acoustic boundary layer[END_REF] with R. Msallam. In a first approximation, the duct of a trombone is a long cylinder with a constant section and the acoustic waves propagate only in the longitudinal direction. We can use a one-dimensional description of the geometry (see Figure 3.1) and in what follows, the trombone is modelled by a real space variable x that ranges from x = 0 at the input to x = L at the output.

•

At the input x = 0, a given non-stationary pressure wave t -→ Π(t) is emitted ; this wave is a perturbation of the ambiant pressure p 0 of the air : (3.1.1)

| Π(t)p 0 | << p 0 , t > 0 . At the output x = L, the waves go outside without any reflection due to the presence of a pavilion and the boundary condition is a "free output" and a nonreflecting boundary condition has to be used. At the initial time t = 0, we can consider that the air satisfies the usual conditions of pressure p(x, 0) ≡ p 0 , temperature T (x, 0) ≡ T 0 and density ρ(x, 0) ≡ ρ 0 . We study in this section a finite volume method able to treat nonlinearities in the acoustic modelling and based on the characteristic decompositions developed in the previous section. 3.2 Linearization of the gas dynamics equations.

•

We study a perfect gas subjected to a motion with variable velocity in space and time. We have noticed that the primitive unknowns of this problem are the scalar fields that characterize the thermodynamics of the gas, i.e. density ρ, internal energy e, temperature T, and pressure p. In what follows, we suppose that the gas is a polytropic perfect gas ; it has constant specific heats at constant volume C v and at constant pressure C p . These two quantities do not depend on any thermodynamic variable like temperature or pressure ; we denote by γ their ratio :

(3.2.1) γ = C p C v (= constant) .
We suppose that the gas satisfies the law of perfect gas that can be written with the following form :

(3.2.2) p = (γ -1) ρ e . As usual, internal energy and temperature are linked together by the Joule-Thomson relation : (3.2.3) e = C v T .

•

In the formalism proposed by Euler during the 18 th century, the motion is described with the help of an unknown vector field u which is a function of space x and time t : (3.2.4) u = u(x, t) . In the following, we will suppose that space x has only one dimemsion (x ∈ IR). We have four unknown functions (density, velocity, pressure and internal energy) linked together by the state law (3.2.2). In consequence, we need three complementary equations in order to define a unique solution of the problem. The general laws of Physics assume that mass, momentum and total energy are conserved quantities, at least in the context of classical physics associated to the paradigm of invariance for the Galileo group of space-time transformations (see e.g. Landau and Lifchitz [START_REF] Landau | Fluid Mechanics[END_REF]). When we write the conservation of mass, momentum and energy inside an infinitesimal volume dx advected with celerity u(x, t), which is exactly the mean velocity of particules that compose the gas, it is classical [START_REF] Landau | Fluid Mechanics[END_REF] to write the fundamental conservation laws of Physics with the help of divergence operators :

(3.2.5) ∂ρ ∂t + ∂ ∂x ρ u = 0 (3.2.6) ∂ ∂t ρ u + ∂ ∂x ρ u 2 + p = 0 (3.2.7) ∂ ∂t 1 2 ρ u 2 + ρ e + ∂ ∂x 1 2 ρ u 2 + p γ -1 u + p u = 0 .
• We introduce the specific total energy E by unity of volume (3.2.8) E = 1 2 u 2 + e , the sound celerity c following the classical expression :

(3.2.9) c = γ p ρ , and total enthalpy H defined according to (3.2.10)

H ≡ E + p ρ = 1 2 u 2 + 1 γ-1 c 2 .
The vector W is therefore composed by the "conservative variables" or more precisely by the "conserved variables" :

(3.2.11) W = ρ , ρ u , ρ E t ≡ ρ , q , ǫ t .
The conservation laws (3.2.5)-(3.2.7) take the following general form of a system of conservation laws : (3.2.12)

∂W ∂t + ∂ ∂x F (W ) = 0
where the flux vector W -→ F (W ) satisfies the following algebraic expression :

(3.2.13) F (W ) = ρ u , ρ u 2 + p , ρ u H t that can be explicited as a true function of state vector W, on one hand with the pressure law P (W ) computed with (3.2.2), (3.2.8) and (3.2.11) :

(3.2.14) P (W ) = (γ -1) ǫ -q 2 2 ρ and on the other hand with an explicit use of the conserved variables ρ, q and ǫ. We obtain :

(3.2.15)

F (W ) = q , q 2 ρ + P (W ) , q ǫ ρ + P (W ) q ρ .
Proposition 3.1. Jacobian matrix of gas dynamics. • The Jacobian matrix dF (W ) of the flux function W -→ F (W ) for the Euler equations of the gas dynamics admits the following expression :

(3.2.16) dF (W ) =   0 1 0 (γ -1) H -u 2 -c 2 (3 -γ) u γ -1 (γ -2) u H -u c 2 H -(γ -1) u 2 γu   .
• The matrix dF (W ) is diagonalizable ; the eigenvalues λ j (W ) satisfy the relations (3.2.17) λ 1 (W ) ≡ uc < λ 2 (W ) ≡ u < λ 3 (W ) ≡ u + c . and the associated eigenvectors r j (W ) are proportional to the following ones :

(3.2.18) r 1 (W ) =   1 u -c H -u c   , r 2 (W ) =   1 u 1 2 u 2   , r 3 (W ) =   1 u + c H + u c   .

•

We first differentiate the pressure law W -→ P (W ) given in (3.2.14) :

(3.2.19) ∂P ∂ρ = γ -1 2 u 2 = (γ-1) H -c 2 , ∂P ∂q = -(γ-1) u , ∂P ∂ǫ = (γ-1)
and the second row of the matrix (3.2.16) is a direct consequence of the relations ∂ ∂ρ

q 2 ρ = -u 2 and ∂ ∂q q 2 ρ = 2 u .

•

The calculus of the third row of matrix in (3.2.16) demands first evaluation of the gradient of ρ u E = u ǫ relatively to the state W. We get

(3.2.20) ∂ ∂ρ q ǫ ρ = -u E , ∂ ∂q q ǫ ρ = E , ∂ ∂ǫ q ǫ ρ = u .
We have also

∂ ∂W (P u) = ∂P ∂W u + p ∂ ∂W q ρ
then we deduce from (3.2.19) and the following expressions for the gradient of velocity ∂ ∂ρ q ρ = -u ρ and

∂ ∂q q ρ = 1 ρ : (3.2.21)      ∂ ∂ρ P q ρ = γ -1 2 u 3 - u p ρ , ∂ ∂q P q ρ = -(γ -1) u 2 + p ρ , ∂ ∂ǫ P q ρ = (γ -1) u .
We add the relations (3.2.20) and (3.2.21) ; then the third row of matrix (3.2.16) admits the following expression :

γ-1 2 u 3u H , H -(γ-1) u 2 , γ u and this result is exactly the third row of the right hand side of (3.2.16) when we take into account the relation (3.2.10) between H, u 2 and c 2 . The relations (3.2.17) and (3.2.18) are elementary to satisfy ; they express simply the three relations : (3.2.22) dF (W ) • r j (W ) = λ j (W ) r j (W ) , j = 1, 2, 3 and Proposition 3.1 is established.

•

We keep into memory the following expression of the Jacobian matrix dF (W ) :

(3.2.23) dF (W ) =   0 1 0 γ-3 2 u 2 (3 -γ) u γ -1 γ-1 2 u 3 -u H H -(γ -1) u 2 γu  
that needs only the datum of velocity u and total enthalpy H of the state W.

Roe matrix.

• We consider two states W left ≡ W l and W right ≡ W r relatively to the gas dynamics, i.e. they both belong to space IR 3 and have an expression of the form (3.2.11). By definition, a Roe matrix A(W l , W r ) between these two states is a 3 by 3 matrix that satisfy the three following properties : (3.3.1)

A(W l , W r ) is a diagonalizable matrix on the field IR of real numbers (3.3.2)

A(W, W ) = dF (W ) (3.3.3) F (W r ) -F (W l ) = A(W l , W r ) • (W r -W l )
. In his original article, P. Roe [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors and Difference Schemes[END_REF] has proposed a very simple algebraic way to construct a Roe matrix for the dynamics of polytropic gas. We propose it in the following Proposition. Proposition 3.2.

Algebraic construction of a Roe matrix [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors and Difference Schemes[END_REF]. Let W l and W r be two states for gas dynamics, defined by their densities ρ l and ρ r , their velocities u l and u r and their total enthalpies H l and H r . We introduce an intermediate state W * (W l , W r ) by its density ρ * , its velocity u * and its total enthalpy H * according to the following relations :

(3.3.4) ρ * = √ ρ l ρ r (3.3.5) u * = √ ρ l u l + √ ρ r u r √ ρ l + √ ρ r (3.3.6) H * = √ ρ l H l + √ ρ r H r √ ρ l + √ ρ r .
Then the matrix A(W l , W r ) defined as the Jacobian matrix of the flux for the intermediate state W * (W l , W r ), i.e.

(3.3.7)

A(W l , W r ) = dF W * (W l , W r ) is a Roe matrix. •
Due to the expression (3.2.23) of the Jacobian matrix of gas dynamics, we remark that the formula (3.3.4) giving the density ρ * is not necessary for the determination of the matrix dF (W * (W l , W r )) and an entire family of states W * (W l , W r ) define a Roe matrix according to the relations (3.3.5), (3.3.6) and (3.3.7). Nevertheless, we keep this definition of density ρ * by convenience and simplicity for future algebraic expressions. The proof of Proposition 3.2 needs some algebraic developments. We begin by the following technical lemma. Proposition 3.3. Under the hypotheses of Proposition 3.2, we have the following relations :

(3.3.8) (u * ) 2 (ρ r -ρ l ) -2 u * (ρ r u r -ρ l u l ) + (ρ r u 2 r -ρ l u 2 l ) = 0 (3.3.9) -u * H * (ρ r -ρ l ) + H * (ρ r u r -ρ l u l ) + u * (ρ r H r -ρ l H l ) = = ρ r u r H r -ρ l u l H l .
• We first evaluate the left hand side of relation (3.3.8) :

(u * ) 2 (ρ r -ρ l ) -2 u * (ρ r u r -ρ l u l ) + (ρ r u 2 r -ρ l u 2 l ) = = u * ( √ ρ r - √ ρ l ) ( √ ρ r u r + √ ρ l u l ) -2 u * (ρ r u r -ρ l u l ) + (ρ r u 2 r -ρ l u 2 l ) = u * √ ρ l ( √ ρ l + √ ρ r ) u l - √ ρ r ( √ ρ l + √ ρ r ) u r + (ρ r u 2 r -ρ l u 2 l ) = ( √ ρ l u l + √ ρ r u r ) ( √ ρ l u l - √ ρ r u r ) + (ρ r u 2 r -ρ l u 2 l ) = 0
and the relation (3.3.8) is established.

•

We work on the left hand side of (3.3.9) as follows :

-u * H * (ρ r -ρ l ) + H * (ρ r u r -ρ l u l ) + u * (ρ r H r -ρ l H l ) = = -u * ( √ ρ r - √ ρ l ) ( √ ρ l H l + √ ρ r H r ) + H * (ρ r u r -ρ l u l ) + u * (ρ r H r -ρ l H l ) = √ ρ l ρ r u * (H r -H l ) + H * (ρ r u r -ρ l u l ) = √ ρ l √ ρ r ( √ ρ l u l + √ ρ r u r ) (-H l + H r ) + (-ρ l u l + ρ r u r ) ( √ ρ l H l + √ ρ r H r ) √ ρ l + √ ρ r = 1 √ ρ l + √ ρ r -ρ l ( √ ρ l + √ ρ r ) u l H l + ρ r ( √ ρ l + √ ρ r ) u r H r = ρ r u r H r -ρ l u l H l
and the proposition 3.3 is established.

•

The proof of Proposition 3.2 consists in satisfying the three hypotheses that define a Roe matrix. First, due to the fact that the relation (3.3.7) defines the matrix A(W l , W r ) as a Jacobian of some state, this matrix is diagonalizable with real elements due to the result of Proposition 3.1 and the first property (3.3.1) is satisfied. The second property (3.3.2) is a simple consequence of the fact that if W l = W r = W, then we have from the relations (3.3.4) to (3.3.6) : W * (W l , W r ) = W and the property results from (3.3.7).

•

The third property (3.3.3) needs more work. We remark that the first row of this matricial relation is clear. For the second row, we have :

Second row of matrix A(W l , W r ) • (W r -W l ) = = γ -3 2 (u * ) 2 (ρ r -ρ l ) + (3-γ) u * (ρ r u r -ρ l u l ) + (γ -1) (ρ r E r -ρ l E l ) = γ -3 2 (u * ) 2 (ρ r -ρ l ) -2 u * (ρ r u r -ρ l u l ) + γ -1 2 (ρ r u 2 r -ρ l u 2 l ) + (p r -p l ) = (ρ r u 2 r -ρ l u 2 l ) + (p r -p l ) due to (3.3.8) = second row of the flux difference F (W r ) -F (W l ) . •
We have also, in consequence of (3.2.23), Third row of matrix

A(W l , W r ) • (W r -W l ) = = u * γ -1 2 (u * ) 2 -H * (ρ r -ρ l ) + (H * -(γ -1) (u * ) 2 ) (ρ r u r -ρ l u l ) + + γ u * (ρ r E r -ρ l E l ) = γ -1 2 u * (u * ) 2 (ρ r -ρ l ) -2 u * (ρ r u r -ρ l u l ) +  + -H * u * (ρ r -ρ l ) + H * (ρ r u r -ρ l u l ) + γ u * (ρ r E r -ρ l E l ) = γ -1 2 u * (u * ) 2 (ρ r -ρ l ) -2 u * (ρ r u r -ρ l u l ) -u * (ρ r H r -ρ l H l ) + + (ρ r u r H r -ρ l u l H l ) + γ u * (ρ r E r -ρ l E l ) due to (3.3.9) = γ -1 2 u * (u * ) 2 (ρ r -ρ l ) -2 u * (ρ r u r -ρ l u l ) + ρ r u 2 r -ρ l u 2 l + + u * (-γ ρ r e r + γ ρ l e l + γ ρ r e r -γ ρ l e l ) + ρ r u r H r -ρ l u l H l = ρ r u r H r -ρ l u l H l due to (3.3.8) = third row of the flux difference F (W r ) -F (W l )
in the view of relation (3.2.13). The proposition 3.2 is established.

Roe flux.

•

The principal interest of the Roe matrix is to be able to use all what has been developed for linear hyperbolic systems in Section 2. In particular, the following linear hyperbolic system defined with a given Roe matrix

A(W l , W r ) (3.4.1) ∂W ∂t + A(W l , W r ) • ∂W ∂x
= 0 can be treated with the upwind scheme defined at proposition 2.7. We obtain by doing this the following Proposition 3.4.

Three formulae for a flux. • Let W l and W r be two fluid states and W * the intermediate state defined by the relations (3.3.4) to (3.3.6). The sound celerity c * of state W * is defined with the help of relation (3.2.10), i.e.

(3.4.2)

c * = (γ -1) H * - (u * ) 2 2 ,
and the eigenvalues λ * j of the Roe matrix A(W l , W r ) ≡ dF (W * (W l , W r )) are given by a relation analogous to (3.2.17).

(3.4.3)

λ * 1 ≡ u * -c * < λ * 2 ≡ u * < λ * 3 ≡ u * + c * .
The associated eigenvectors r * j ≡ r j (W * ) are proportional to the following ones :

(3.4.4) r * 1 =   1 u * -c * H * -u * c *   , r * 2 =   1 u * 1 2 (u * ) 2   , r * 3 =   1 u * + c * H * + u * c *   .
• We introduce the decomposition of vector W r -W l in the basis r * j : (3.4.5) W r -W l = j=3 j=1 α j r * j . The three following relations define a unique numerical flux Φ(W l , W r ) named the Roe flux between the two states W l and W r :

(3.4.6) Φ(W l , W r ) = F (W l ) + j=3 j=1 (λ * j ) -α j r * j (3.4.7) Φ(W l , W r ) = F (W r ) - j=3 j=1 (λ * j ) + α j r * j  (3.4.8) Φ(W l , W r ) = 1 2 F (W l ) + F (W r ) -1 2 | A(W l , W r ) | • (W r -W l . •
The first non-obvious point is to verify that the relation (3.4.2) defines a real number c * . We have

H * - (u * ) 2 2 = √ ρ l H l + √ ρ r H r √ ρ l + √ ρ r - 1 2 √ ρ l u l + √ ρ r u r √ ρ l + √ ρ r 2 = = 1 ( √ ρ l + √ ρ r ) 2 ( √ ρ l + √ ρ r ) √ ρ l 1 2 u 2 l + 1 γ -1 c 2 l + √ ρ r 1 2 u 2 r + 1 γ -1 c 2 r - 1 2 ρ l u 2 l + 2 ρ * u l u r + ρ r u 2 r = 1 ( √ ρ l + √ ρ r ) 2 1 2 ρ * u 2 l -2 ρ * u l u r + ρ * u 2 r + ρ l + ρ * γ -1 c 2 l + ρ * + ρ r γ -1 c 2 r = 1 ( √ ρ l + √ ρ r ) 2 1 2 ρ * (u r -u l ) 2 + ρ l + ρ * γ -1 c 2 l + ρ * + ρ r γ -1 c 2 r > 0 and (3.4.9) c * = γ-1 2 ρ * (u r -u l ) 2 + (ρ l + ρ * ) c 2 l + (ρ * + ρ r ) c 2 r √ ρ l + √ ρ r .
• We make the difference between the right hand sides of (3.4.6) and (3.4.7). We get :

F (W r ) -F (W l ) - j=3 j=1 (λ * j ) + + (λ * j ) -α j r * j = = A(W l , W r ) • (W r -W l ) -

•

We make explicit the parameters α j introduced in relation (3.4.5) in order to be complete for the implementation of the above formulae on a computer.

Proposition 3.5.

New acoustic impedance. With the notations introduced at Proposition 3.4, and denoting by p l and p r the respective pressures of states W l and W r , we have the following relations for the scalar components α j of the state difference W r -W l in relation (3.4.5) :

(3.4.10)

α 1 = 1 2 (c * ) 2 (p r -ρ * c * u r ) -(p l -ρ * c * u l ) (3.4.11) α 2 = - 1 (c * ) 2 (p r -(c * ) 2 ρ r ) -(p l -(c * ) 2 ρ l )  (3.4.12) α 3 = 1 2 (c * ) 2 (p r + ρ * c * u r ) -(p l + ρ * c * u l )
with acoustic impedance ρ * c * that is nomore the one ρ 0 c 0 of a reference state as in traditional acoustics but an impedance associated with the Roe intermediate state W * (W l , W r ) of relations (3.3.4) to (3.3.6).

•

We have just to explicit the three components of the relation (3.4.5). It comes :

  ρ r -ρ l ρ r u r -ρ l u l ρ r E r -ρ l E l   =   α 1 + α 2 + α 3 α 1 (u * -c * ) + α 2 u * + α 3 (u * + c * ) α 1 (H * -u * c * ) + α 2 (u * ) 2 2 + α 3 (H * + u * c * )  
(3.4.13) α 1 + α 2 + α 3 = ρ rρ l , and we deduce after multiplying the equation (3.4.13) by -u * and adding to the second equation of the above matrix equality :

c * (α 3 -α 1 ) = ρ r u r -ρ l u l -u * (ρ r -ρ l ) = ρ r u r -ρ l u l -( √ ρ r - √ ρ l ) ( √ ρ l u l + √ ρ r u r ) , then (3.4.14) c * (-α 1 + α 3 ) = ρ * (u r -u l ) . •
We deduce from the third equation of relation (3.4.5) :

(c * ) 2 γ -1 (α 1 + α 3 ) = ρ r E r -ρ l E l - 1 2 (u * ) 2 (ρ r -ρ l ) -u * c * (α 3 -α 1 ) = 1 γ -1 (p r -p l ) + 1 2 (ρ r u 2 r -ρ l u 2 l ) + 1 2 (u * ) 2 (ρ r -ρ l ) -u * (ρ r u r -ρ l u l ) = 1 γ -1 (p r -p l )
due to (3.3.8). Then we have :

(3.4.15) (c * ) 2 (α 1 + α 3 ) = p r -p l . •
The solution of the 2 by 2 linear system with unknowns α 1 and α 3 defined by the relations (3.4.14) and (3.4.15) directly gives the relations (3.4.10) and (3.4.12). The expression (3.4.11) of variable α 2 is a direct consequence of the relations (3.4.10), (3.4.12) and (3.4.13) and Proposition 3.5 is proven. • We introduce as above two states W l ≡ W 0 and W r ≡ W 3 and the Roe matrix A(W l , W r ) described in the preceding sub-sections. We have in particular the relation (3.5.1) W r -W l ≡ W 3 -W 0 = j=3 j=1 α j r * j and we do not make the confusion between the eigenvalues λ j (W 0 ) of the left state, λ j (W 3 ) of the right state, and λ * j of the Roe matrix. We introduce the following two intermediate states W 1 and W 2 according to (3.5.2)

W 1 = W 0 +α 1 r * 1 , W 2 = W 0 +α 1 r * 1 + α 2 r * 2 = W 3 -α 3 r * 3
and illustrated on Figure 3.2. Note that λ j (W k ) is well defined for j = 1, 2, 3 and k = 0, 1, 2, 3 : it is the j 0 eigenvalue of the k 0 intermediate state W k . We define now the set S of sonic indices by the condition that the sign of the j 0 eigenvalue is increasing from negative to positive values accross some j-wave :

(3.5.3) S = j ∈ {1, 2, , 3}, λ j (W j-1 ) < 0 < λ j (W j ) . The modification of the Roe flux is active only for the sonic indices and we introduce a polynomial p j of degree 3 by the classical Hermite interpolation conditions (3.5.4)

p j (0) = 0 , p ′ j (0) = λ j (W j-1 ) , p j (α j ) = λ * j α j , p ′ j (α j ) = λ j (W j ) , j ∈ S
that defines explicitely the polynomial p j (•) by the algebraic relation (3.5.5)

       p j (ξ) = λ j (W j ) + λ j (W j-1 ) -2 λ * j (α j ) 2 ξ 3 + + 3 λ * j -2 λ j (W j-1
)λ j (W j ) α j ξ 2 + λ j (W j-1 ) ξ . • With the hypothesis that j ∈ S, it is not difficult to see [START_REF] Dubois | A non-parameterized entropy correction for Roe's approximate Riemann solver[END_REF] that the polynomial p j (•) has a unique minimum inside the interval (0, α j ). The argument ξ * j of this point of minimum is given according to :

(3.5.6) ξ * j = -λ j (W j-1 ) α j 3 λ * j -2 λ j (W j-1 )λ j (W j ) + + 3 λ * jλ j (W j )λ j (W j-1 ) 2λ j (W j-1 ) λ j (W j )

.

Since p j (ξ * j ) is the unique minimum of the polynomial p j (•) on the interval (0, α j ), we have 3.6 Nonlinear flux boundary conditions.

• At the two extremities x = 0 and x = L of the pipe, we have to express on one hand the datum of a given nonstationary pressure Π(t) at x = 0 and on the other hand a free output of the waves at x = L.

•

For the numerical boundary condition for pressure, we follow a general approach founded on the so-called partial Riemann problem [START_REF] Dubois | Partial Riemann problem, Boundary conditions and gas dynamics[END_REF] that generalizes to nonlinear hyperbolic systems the reflection operator of relation (2.5.22). For a given discrete time t n = n ∆t, and a given state W n 1/2 ≡ W r in the first cell of the unidimensional mesh, we construct a boundary state W n 0 ≡ W l that satisfies the boundary constraint (3.6.1) p(W n 0 ) = Π n+1/2 , n ≥ 0 , and moreover, we impose that the state W n 1/2 present in the first cell is issued from the boundary state W n 0 with an ingoing 3-wave, i.e. we impose the relation (3.6.2)

W n 1/2 -W n 0 = α 3 r *

•

The proof of Proposition 3.7 is a consequence precisely of the preceding subsections about the Roe flux. We first remark that the relations (3.4.5) and (3.6.4) are absolutly identical. Then we deduce that necessarily α 1 = α 2 = 0 and according to the relations (3.4.11) and (3.4.12), we get (3.6

.7) Π -ρ * c * u l = p r -ρ * c * u r (3.6.8) Π -(c * ) 2 ρ l = p r -(c * ) 2 ρ r .
We deduce simply u ru l = p r -Π ρ * c * = c * (ρ rρ l ) ρ * and due to the relation (3.4.9), we get :

(c * ) 2 = 1 ( √ ρ r + √ ρ l ) 2 γ -1 2 ρ * c * (ρ r -ρ l ) ρ * 2 + γ ρ l + ρ * ρ l Π+ ρ * + ρ r ρ r p r .
Then after multiplication by (ρ rρ l ), we obtain with the help of (3.6.8) :

0 = (p r -Π) - γ -1 2 (p r -Π) ( √ ρ r - √ ρ l ) 2 ρ * -γ √ ρ r - √ ρ l √ ρ r + √ ρ l ρ l + ρ * ρ l Π + ρ * + ρ r ρ r p r = γ (p r -Π) - γ -1 2 (p r -Π) √ ρ r √ ρ l + √ ρ l √ ρ r -γ √ ρ r √ ρ l -1 Π -γ 1 - √ ρ l √ ρ r p r .



Franc ¸ois Dubois

We multiply the previous equality by ρ l ρ r and we get :

- γ -1 2 (p r -Π) ( 1 + ρ l ρ r ) -γ Π + γ ρ l ρ r p r = 0 , id est - γ -1 2 (p r -Π) + γ p r ρ l ρ r = γ Π + γ -1 2 (p r -Π)
and the relation (3.6.5) is established.

• We deduce from the previous relation :

ρ * = (γ +1) Π + (γ -1) p r (γ -1) Π + (γ +1) p r ρ r , ρ r -ρ l = 2 (p r -Π) (γ -1) Π + (γ +1) p r ρ r
and due to the relation (3.6.8) :

(c * ) 2 = (γ -1) Π + (γ +1) p r 2 ρ r , ρ * c * = (γ -1) Π + (γ +1) p r 2 √ ρ r
and the relation (3.6.6) is an easy consequence of the last equality joined with (3.6.7). The proposition 3.7 is established.

•

The determination of a nonlinear nonreflecting boundary condition at x = L is still an open mathematical problem. We recommand for deriving a flux boundary condition for such a situation to impose that no wave are present at the interaction for the last interface j = J. We just write (3.6.9) f n+1/2 J = F (W n J-1/2 ) , n ≥ 0 which is equivalent of introducing a right boundary state W n J according to the simple relation W n J = W n J-1/2 and then making these two states interacting with the Roe flux :

f n+1/2 J = Φ(W n J-1/2 , W n J )
. This last definition is equivalent to the one proposed in (3.6.9) due to the property (3.3.2) of the Roe matrix.

4)

Second order and two space dimensions. 4.1 Towards second order accuracy.

•

The finite volume method described in the previous sections is a natural method for the discretization of systems of m conservation laws. It conducts to an explicit scheme in time : the evaluation of the field W n+1 at time step (n+1) ∆t needs only the knowledge of the field W n j+1/2 for j = 0, • • • , J -1 at the preceding time step n ∆t. This evaluation needs a certain number of auxiliary computations without the resolution of any linear system involving the new field. The method is parameterized by the choice of a numerical flux and a great flexibility can be adopted at this level. We have proposed two fluxes for nonlinear problems related to nonlinear acoustics and gas dynamics, the Roe flux Φ(•, ) is first order accurate :

(4.1.4) 1 ∆x f n+1/2 j+1 -f n+1/2 j = ∂F (W ) ∂x n+1/2 j+1/2 + O ∆t + ∆x .
In a similar way, the use of an explicit scheme in time conducts to

(4.1.5) 1 ∆t (W n+1 j+1/2 -W n j+1/2 ) + 1 ∆x f n+1/2 j+1 -f n+1/2 j = 0
and maintains this first order accuracy for the finite volume scheme.

• We develop in this section the fact that it is possible to improve the method, i.e. to define a method with a relation of the type (4.1.5), and that conduct to a troncation error of second order : (4.1.6)

       1 ∆t (W n+1 j+1/2 -W n j+1/2 ) + 1 ∆x f n+1/2 j+1 -f n+1/2 j = = ∂W ∂t + ∂ ∂x F (W ) n+1/2 j+1/2 + O ∆t 2 + ∆x 2 .
The price to pay is to develop flux formulae much more complicated than the simple relation (4.1.1). When the second order precision (4.1.6) is achieved with a stable scheme, the precision is sufficient to develop predictive computations in acoustics and aerodynamics, whereas that is not the case with the initial scheme (4.1.1) (4.1.5).

4.2

The method of lines.

•

The simplest way to extend the first order finite volume scheme is first to develop a new vision of the method with emphasis more on abstraction. We have presented a method founded on the integration of the conservation law (4.1.3)

W(x, t ) n x x j+1 x j x j-1 W j+1/2 W j-1/2 W j - W j +

•

The simplest case is illustrated on Figure 4.1. It imposes simply the reconstructed function W (x) to be constant in each interval : (4.3.2) W (x) ≡ W j+1/2 , x j < x < x j+1 . The two limit values on each side of the interface x j are the following ones : W - j = W j-1/2 , W + j = W j+1/2 and the explicit version of this finite volume scheme is the standard first order numerical flux (4.1.1) as seen at Proposition 4.1. In the context of the method of lines, we obtain : (4.3.3) f j = Φ(W j-1/2 , W j+1/2 ) .

•

The second order accurate Muscl method consists first in restricting the methodology to a scalar field W (•) and to construct an affine function in each interval K j+1/2 instead of a constant function as in (4.3.2). We set (4.3.4) W (x) ≡ W j+1/2 + p j+1/2 xx j+1/2 , x j < x < x j+1 . The simplest choice for a slope is the one of the centered scheme : • The choice of a numerical flux given according to the relations (4.3.8) f j = Φ W j-1/2 + 1 4 (W j+1/2 -W j-3/2 ) , W j+1/2 -1 4 (W j+3/2 -W j-1/2 ) lead to an unstable scheme when we consider the particular case of the advection equation with the first order explicit scheme in time.

Proposition 4.2.

Linear Muscl scheme is unstable. We apply the linear Muscl approach for the advection equation. Then the numerical scheme obtained by association of (4.3.8) and the upwind scheme (2.1.19) conducts to the following explicit first order scheme : This scheme is unstable for each ∆t > 0.

• Due to the expression (2.1.19) of the upwind scheme, the discrete first order in time advection equation can be written :

W n+1 j+1/2 -W n j+1/2 + a ∆t ∆x W -, n j+1 -W -, n j = 0
and the expression (4.3.9) is a consequence of the left extrapolation (4.3.6). For the study of stability, we introduce a profile of the type W n j+1/2 = e (i k (j+1/2) ∆x) with a wave number k. The scheme (4.3.9) can be written as

W n+1 j+1/2 = g k ∆x, a ∆t ∆x W n j+1/2
with an amplification coefficient g(ξ, σ) (ξ = k ∆t, σ = a ∆t ∆x ) given simply by the expression g(ξ, σ) = 1σ 1e -i ξσ 4 e i ξ -1e -i ξ + e -2 i ξ , then For ξ arbitrarily small, we deduce from (4.3.10) : | g(ξ, σ) | 2 = 1 + σ 2 ξ 2 + O ξ 4 which establishes the instability for all σ = 0.

• The above remark motivates the introduction of so-called slope limiters, intensively studied during the period 1980-90 after the pioneering work of Van Leer [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme IV. A New Approach to Numerical Convection[END_REF]. The idea is to search an interpolation W - j of the field W (•) at the left of the point x j from the neighbouring mean values W j-3/2 , W j-1/2 and W j+1/2 and by left-right invariance of the procedure, to construct an interpolated value W + j from the first right neighbours W j-1/2 , W j+1/2 and W j+3/2 as suggested on Figure 4.2. We replace the relations (4.3.6) and (4. We remark that the relations (4.3.6) and (4.3.7) are a particular case of the general nonlinear relations (4.3.11) and (4.3.12) with the particular choice ϕ(r) = Among all the possible choices, we have adopted for fluid mechanics [START_REF] Dubois | Solution of the Euler Equations Around a Double Ellipsoïdal Shape Using Unstructured Meshes and Including Real Gas Effects[END_REF] the so-called STS-limiter defined by the relations The control volumes are exactly the elements of the mesh.

•

We focus now on the use of unstructured meshes for the extension to second order accuracy of the finite volume method. As in the one-dimensional case, the domain of study is decomposed into finite elements (or control volumes) K ∈ E T than can be structured in a Cartesian way (Figure 4.4) or with a cellular complex as in Figure 4.5. In both cases, the intersection of two finite elements define an interface f ∈ F T . We denote by n f the normal at the interface f that separates a left control volume K l (f ) and a right control volume K r (f ). The ordinary differential equation (4.2.6) is replaced by a multidimensional version : (4.4.5)

| K | dW K dt + f ⊂∂K | f | Φ W K , n f , W K r (f ) = 0 , K ∈ E T .
For internal interfaces, the function Φ • , n f , •) is equal e.g. to the Roe flux between states W K l (f ) and W K r (f ) in the one-dimensional direction along normal n f in order to take into account the invariance by rotation of the equations of gas dynamics (see [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]).

•

We consider now a finite element K internal to the domain. The extension to second order accuracy of the finite volume scheme consists in replacing the arguments W K l (f ) and W K r (f ) in relation (4.4.5) by nonlinear extrapolations W K l (f ), f and W K r (f ), f on each side of the boundary of state data and evaluated as described in what follows. We first introduce the set N (K) of neighbouring cells of given finite element K ∈ E T , as illustrated on Figure 4.6 : (4.4.6) N (K) = L ∈ E T , ∃ f ∈ F T , f ⊂ ∂K ∩ ∂L . For L ∈ N (K), we suppose by convention that the normal n f to the face f ⊂ ∂K ∩ ∂L is external to the element K id est K r (f ) = K, K l (f ) = L. We introduce also the point y K, f on the interface f ⊂ ∂K that links the barycenters x K and x K r (f ) :

(4.4.7)

y K, f ≡ (1 -θ K, f ) x K + θ K, f x K r (f ) , y K, f ∈ f , f ⊂ ∂K ,
K finite element internal to mesh T . Then, following Pollet [START_REF] Pollet | Méthodes de calcul relatives aux interfaces missiles-propulseurs[END_REF], for z equal to one scalar variable of the family : (4.4.8) z ∈ { ρ , ρ u , ρ v , p } we evaluate a mean value z K, f on the interface f : (4.4.9) z K, f = (1θ K, f ) z K + θ K, f z K r (f ) and the gradient ∇z(K) of field z(•) in volume K with a Green formula : Three neighbouring cells are necessary to determine the gradient in triangle K and to limit eventually its variation.

•

The numerical integration of such kind of system is done with a Runge-Kutta scheme as presented in [START_REF] Chargy | Méthodes numériques pour le calcul d'écoulements compressibles, applications industrielles[END_REF]. We have used with success in [START_REF] Dubois | Solution of the Euler Equations Around a Double Ellipsoïdal Shape Using Unstructured Meshes and Including Real Gas Effects[END_REF] the Heun scheme of second order accuracy for discrete integration of (4.5.1) between time steps n ∆t and (n+1) ∆t : (4.5.2)

| K | ∆t W K -W n K + f ⊂∂K | f | Φ W n K, f , n f , W n K r (f ), f = 0 , K ∈ E T (4.5.3) | K | ∆t W K -W K + f ⊂∂K | f | Φ W K, f , n f , W K r (f ), f = 0 , K ∈ E T (4.5.4) W n+1 K = 1 2 W K + W n K , K ∈ E T .

5)

References.

Figure 1 . 1 .Figure 1

 111 Figure 1.1. The solution u(x, t) of the advection equation is constant along the characteristic lines.

Figure 1 . 3 .Figure 1 . 4 .

 1314 Figure 1.3. Snap shot of the solution of the advection equation at time t = T.

Figure 1 . 5 .Figure 1 . 6 .

 1516 Figure 1.5. Initial-boundary value problem for the advection equation with a > 0 in the domain 0 < x < L and t > 0.

Figure 2 . 1 .

 21 Figure 2.1. Space-time grid for the finite volume method.

Figure 2 . 2 .

 22 Figure 2.2. Upwinding of the information for the advection equation.

  0

Figure 2 . 3 .

 23 Figure 2.3. Notations for the one-dimensional finite volume method.

  Figure 2.3. • At time step n ∆t , the discrete field is entirely known and is composed of all the values W n j-1/2 for j = 1, • • • , J . With a flux function Φ(•, •) as proposed at relation (2.1.14), we observe that the two boundary fluxes f

Figure 2 . 4 .Figure 2 . 5 .

 2425 Figure 2.4. Interface upwind numerical flux for a model problem with two equations.t

.

  n+1/2 j , j = 1, • • • , J -1 and the boundary fluxes  f Recall first that the physical flux F (W ) function for the acoustic equation (2.4.5) is equal to (2.4.19)

Figure 2 . 8 .Figure 2 . 9 .

 2829 Figure 2.8. Reflection operator at x = 0.

Figure 2 . 10 .

 210 Figure 2.10. Discontinuity at the interface between two cells.

Figure 3 . 1 .

 31 Figure 3.1. Long unidimensional pipe for the modelling of a trombone.

  A(W l , W r ) • r * j = λ * j r * j for each integer j. The proof of relation (3.4.8) is obtained by taking the half sum of (3.4.6) and (3.4.7). It is analogous to the one done for Proposition 2.7. The proof of Proposition 3.4 is completed.

  Proposition 3.6. An algorithm for the Roe flux. Let W l and W r be two compressible fluid states. The computation of the Roe flux Φ(W l , W r ) of relations (3.4.6)-(3.4.8) between these two states is summarized by the following points : • Evaluation of density ρ * , velocity u * and total enthalpy H * of the intermediate state W * with the relations (3.3.4) to (3.3.6), • Determination of the sound celerity c * of the intermediate state W * from the previous data with the relation (3.4.2), • Eigenvectors r * j of the Roe matrix with the relations (3.4.4),

Figure 3 . 2 .

 32 Figure 3.2. Intermediate states for the entropy correction of the Roe upwind scheme.

.

  Then the modified flux Φ modif (W l , W r ) is defined from the Roe flux Φ(W l , W r ) by the relation (3.5.7) Φ modif (W l , W r ) = Φ(W l , W r ) + added numerical viscosity explicit.

3ΠFigure 3 . 3 .

 33 Figure 3.3. Nonlinear boundary condition for given pressure at x = 0 with the Roe upwind scheme.

Figure 4 . 1 .Figure 4 . 2 .

 4142 Figure 4.1. First order and second order interpolation at the interface x j . W(x, t ) n

  /2 -W j-1/2 ) , and due to (4.3.4) and (4.3.5), the extrapolated values W - j and W + j on each side of the interface located at the position x j are the following ones : /2 -W j-1/2 ) .

   g(ξ, σ) = 1σ ( 1cos ξ ) -σ 4 ( -1 + cos 2ξ )i σ sin ξ + cos ξ ) 2i σ 2 sin ξ ( 3cos ξ ) .

(

  3.7) by a nonlinear interpolation parameterized by a slope limiter IR∋ r -→ ϕ(r) ∈ IR : 1/2 -W j-3/2 W j+1/2 -W j-1/2 W j+1/2 -W j-1/2 2 -W j+1/2 W j+1/2 -W j-1/2 W j+1/2 -W j-1/2 .

12

  (1 + r) . The limiter function satisfies very often the functional relation (4.3.13) ϕ(r) ≡ r ϕ 1 r , r > 0 .

  Figure 4.4. Structured Cartesian mesh. The control volumes are exactly the elements of mesh T .

Figure 4 . 5 .

 45 Figure 4.5. Unstructured mesh composed by triangular elements.The control volumes are exactly the elements of the mesh.

Figure 4 . 6 .

 46 Figure 4.6. Cellular complex mesh with triangles and quadrangles.Three neighbouring cells are necessary to determine the gradient in triangle K and to limit eventually its variation.

  Then we solve the system (2.4.25) (2.4.26) and find finally the relation (2.4.21).

	The process is analogous for the right boundary. The input datum is imposed
	according to the relation (2.4.14) :
	(2.4.27)

  to the relation (2.7.19). Then the relation (2.7.17) is established and the proposition 2.7 is proven.

and the first order upwind finite volume scheme is defined by the way we evaluate the coefficient ψ j (W l , W r ) with the upwind scheme associated with the advection equation (2.7.7) :

(2.7.9) ψ j (W l , W r ) = λ j ϕ j, l if λ j > 0 λ j ϕ j, r if λ j < 0 .

• For any real number µ , we introduce the positive part µ + and the negative part µ -by the relations (2.7.10)

We remark that we have (2.7.11) µ ≡ µ + + µ -, ∀ µ ∈ IR (2.7.12)

| µ | ≡ µ +µ -, ∀ µ ∈ IR . We introduce also the absolute value | Λ | of the diagonal matrix Λ by the condition :

(2.7.13)

and due to the relation (2.7.5), the absolute value | A | of the matrix A is defined by : (2.7.14)

• Computation of the characteristic variables α j in (3.4.5) for the difference W r -W l with the relations (3.4.10) to (3.4.12),

• Final computation of the Roe flux Φ(W l , W r ) with the minimum of work :

(3.4.16)

The proof of the relations (3.4.16)-(3.4.19) is obtained by starting from the expression of the Roe flux given in (3.4.6). We know that ) is established. We remark also that the algebraic expression (3.4.11) for α 2 is not necessary for the implementation of the algorithm.

Entropy correction.

•

The Roe flux replaces the nonlinear waves of the gas dynamics, i.e. the rarefactions and the shock waves by linear waves that are the contact discontinuities. If sufficiently weak shock waves occur for a given discontinuity between two states W left and W right , the Roe flux presented above is a good approximation, but if a rarefaction containing a sonic point is present among the nonlinear waves that solves the discontinuity problem between W left and W right , it has been early remarked that for this very particular situation, the Roe flux does not satisfy the entropy condition (see e.g. Godlewski and Raviart [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]).

•

A popular response has been proposed by Harten [START_REF] Harten | High Resolution Schemes for Hyperbolic Conservation Laws[END_REF] with a tuning parameter that plays in fact the role of an artificial viscosity and P. Roe himself [START_REF] Roe | Some contributions to the Modelling of Discontinuous Flows[END_REF] has proposed a nonparameterized entropy correction for his flux. With G. Mehlman, we have treated the same subject by the introduction of hyperbolic nonlinear models with nonconvex flux functions and have proved a discrete entropy inequality if sufficiently weak nonlinear waves are present in the problem [START_REF] Dubois | A non-parameterized entropy correction for Roe's approximate Riemann solver[END_REF]. We detail here the modification of the algorithm that we have proposed and tested numerically for various gas dynamics problems. inside the space-time domain V n+1/2 j+1/2 = ]x j , x j+1 [ × ]t n , t n+1 [ as suggested in (2.1.6). With the method of lines, we just integrate the conservation (4.1.3) in space in each control volume K j+1/2 = ]x j , x j+1 [. It is straightforward to introduce the mean value W j+1/2 (t) in this finite element :

then we integrate the conservation law (4.1.3) in space in the cell K j+1/2 and taking into account the relation

As usual with the finite volume method, a numerical scheme can be obtained from the relations (4.2.2) (4.2.3) by replacing the relation (4.2.3) by some explicit function over the set of all discrete variables introduced for the relation (4.2.1). To fix the ideas, we introduce a dynamic state vector Z(t) composed by all the dynamic variables on the finite mesh : (4.2.4)

The discretization in space is achieved if we are able to determine the numerical flux f j (t) with the help of both the dynamic state vector Z(•) and the boundary conditions, that is the input pressure Π(t) in the example considered in the last section. As in relation (2.1.12), we introduce a local numerical flux function Ψ j (•, •) relative to the vertex x j : (4.2.5) f j (t) = Ψ j Π(t), Z(t) . We replace the relation (4.2.3) by the numerical approximation (4.2.5) inside the equation (4.2.2) of dynamic evolution of the state variable W j+1/2 (•). We obtain the following ordinary differential equation

• The method of lines is a semi-discrete version of the finite volume method. It is obtained by integration in space of the conservation law without integration in time. The result is not a numerical scheme but just an ordinary differential equation for the dynamic state vector Z(•) described component by component with the equation (4.2.6). The method is parameterized by the local numerical flux functions Ψ j •, •) and take the general form of a dynamical system parameterized by the pressure function t -→ Π(t) :

The discrete dynamic function

and it is defined from the (J +1) local numerical fluxes (Ψ j ) j=0, •••, J with the very simple algebra relative to the finite volume method :

(4.2.9)

Proposition 4.1. Explicit Euler scheme. With the choice of the first order scheme in space, that is (4.2.10)

, and the first order explicit forward Euler scheme for the ordinary differential equation (4.2.7), id est

we recover the previous first order finite volume scheme (4.2.12)

• We write the relation (4.2.10) for the particular control volume K j+1/2 and we get :

) c.f. (4.2.10) and the relation (4.2.12) is established.

4.3

The method of Van Leer.

•

We turn now to the construction of a second order accurate version of the finite volume method as proposed initially with the "Multidimensional Upwindcentered Scheme for Conservation Laws" of B. Van Leer [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme V. A Second Order Sequel to Godunov's Method[END_REF]. The fundamental idea of this scheme is the reconstruction of a function IR ∋ x -→ W (x) ∈ IR from his mean values W j+1/2 in each cell K j+1/2 . The reconstructed function is regular inside each control volume K j+1/2 and is discontinuous at the interfaces x l between two control volumes. The application to the finite volume method replaces the scheme (4.1.1) by the same Roe flux interaction Φ(•, •) considered for the two extrapolated data W - j and W + j on each side of the boundary : (4.3.1) 4.4 Second order accurate finite volume method for fluid problems.

•

We detail in this section a generalization for unstructured meshes of the Muscl scheme proposed by Van Leer [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme V. A Second Order Sequel to Godunov's Method[END_REF]. At one space dimension on a uniform mesh, it is classical to consider a scalar field z among the primitive variables, id est

. This nonlinear interpolation is done with a slope limiter ϕ(•) that operates on each variable proposed in (4.4.1) and we have typically when a left-right invariance is assumed [START_REF] Dubois | Nonlinear Interpolation and Total Variation Diminishing Schemes[END_REF] :

• An ideal extrapolation of field z(•) at the interface f would be : (4.4.11) z K, f = z K + ∇z(K) • y K, fx K but the corresponding scheme is unstable as explicited at Proposition 4.2. When the variation ∇z(K) • y K, fx K is very important, it has to be "limited" as first suggested by Van Leer [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme IV. A New Approach to Numerical Convection[END_REF]. For doing this in a very general way, we introduce the minimum m K (z) and the maximum M K (z) of field z(•) in the neighbouring cells : (4.4.12)

We introduce a nonlinear extrapolation of the field z(•) between center x K and boundary face

is chosen as large as possible and less than or equal to 1 in order to satisfy the constraints (4.4.16) :

(4.4.17)

•

In the one dimensional case with a regular mesh, it is an exercice to re-write the extrapolation (4.4.15) under the usual form (4.4.3) in the context of finite differences. In this particular case, some limiter functions r -→ ϕ k (r) associated with particular parameters k are shown on Figure 4.3. For k = 1, we recover the initial limiter proposed by Van Leer in the fourth paper of the family "Towards the ultimate finite difference scheme..." [START_REF] Van Leer | Towards the Ultimate Conservative Difference Scheme IV. A New Approach to Numerical Convection[END_REF] ; for this reason, we have named it the "Towards 4" limiter (see Figure 4.3). When k = 1 2 we obtain the "min-mod" limiter proposed by Harten [START_REF] Harten | High Resolution Schemes for Hyperbolic Conservation Laws[END_REF]. The intermediate value k = 3 4 is a good compromise between the "nearly unstable" choice k = 1 and the "too compressive" min-mod choice. We have named it STS (see also (4.3.14)) and it has been chosen for our Euler computations in [START_REF] Dubois | Solution of the Euler Equations Around a Double Ellipsoïdal Shape Using Unstructured Meshes and Including Real Gas Effects[END_REF].

•

We explain now the way the preceding scheme is adapted near the boundary. We first consider a fluid boundary. When K is a finite element with some face g ⊂ ∂K lying on the boundary, we still define the set N (K) of neighbouring cells by the relation (4.4.6) as shown on Figure 4.7. The number of neighbouring cells is just less important in this case. Then points y K, f are introduced by relation (4.4.7) if face f does not lie on the boundary and by taking the barycenter of face g if it is lying on the boundary. The only difference is the way the values z K, g are extrapolated for the face g that is on the boundary ; we set (4.4.18) z K, g = z K , g ⊂ ∂K , g face lying on the boundary of the domain. When values z K, f are determined for all the faces f ⊂ ∂K, the gradient  Franc ¸ois Dubois ∇z(K), the minimal m K (z) and maximal M K (z) values among the neighbouring cells are still determined with the relations (4.4.10), (4.4.12) and (4.4.13) respectively. The constraints (4.4.16) remain unchanged except that no limitation process is due to the faces lying on the boundary. In a precise way, we set :

Then the interpolated values z K, f for all the faces f ⊂ ∂K are again predicted with the help of relation (4.4.15).

• For a rigid wall, the limitation process is a little modified, as presented at Figure 4.8. We first introduce the limit face g inside the set of neighbouring cells :

(4.4.20)

∂K, g on the boundary . For the face(s) g ⊂ ∂K lying on the solid boundary, we determine preliminary values z K, g by taking in consideration at this level the impenetrability boundary condition u • n g = 0. We introduce the two components n x g and n y g of the normal n g at the boundary and we set, in coherence with variables (4.4.8) :

(4.4.21)

We consider also these values for the limitation algorithm. We define "external values" z L for L = g and face g lying on the boundary as equal to the ones defined in relation (4.4.21) : (4.4.22) z g ≡ z K, g , z(•) field defined in (4.4.21), g ⊂ ∂K on the boundary. Then the extrapolation algorithm that conducts to relation (4.4.15) for extrapolated values z K, f is used as in the internal case.

4.5 Explicit Runge-Kutta integration with respect to time.

•

When all values z K, f are known for all control volumes K ∈ E T , all faces f ⊂ K and all fields z(•) defined at relation (4.4.8), extrapolated states W K, f are naturally defined by going back to the conservative variables. Then we introduce these states as arguments of the flux function Φ • , n f , •) and obtain by this way a new system of ordinary differential equations : (4.5.1)