T Gallou Ët

R Herbin

J.-C Latch

A CONVERGENT FINITE ELEMENT-FINITE VOLUME SCHEME FOR THE COMPRESSIBLE STOKES PROBLEM PART I -THE ISOTHERMAL CASE

In this paper, we propose a discretization for the (nonlinearized) compressible Stokes problem with a linear equation of state ρ = p, based on Crouzeix-Raviart elements. The approximation of the momentum balance is obtained by usual finite element techniques. Since the pressure is piecewise constant, the discrete mass balance takes the form of a finite volume scheme, in which we introduce an upwinding of the density, together with two additional stabilization terms. We prove a priori estimates for the discrete solution, which yields its existence by a topological degree argument, and then the convergence of the scheme to a solution of the continuous problem.

Introduction

The problem addressed in this paper is the system of the so-called barotropic compressible Stokes equations, which reads:

-∆u + ∇p = f (1.1a) div(ρ u) = 0 (1.1b) ρ = ̺(p) (1.1c)
where ρ, u and p stand for the density, velocity and pressure in the flow, respectively, and f is a forcing term. The function ̺ is the equation of state used for the modelling of the particular flow at hand, which may be the actual equation of state of the fluid or may result from assumptions concerning the flow. Here, we only consider the following relation, which corresponds to an isothermal flow of a perfect gas:

(1.2) ̺(p) = A p
where A is a positive constant. Since the sound velocity is defined by c 2 = dp/dρ, A = Ma 2 /V 2 , where Ma is the Mach number and V is a characteristic velocity of the flow. This system of equations is posed over Ω, a bounded domain of R d , d ≤ 3 supposed to be polygonal (d = 2) or polyhedral (d = 3). It is supplemented by homogeneous boundary conditions for u, and by prescribing the total mass M of the fluid: where M is a positive real number.

In this paper, we study a numerical scheme for the solution of this problem, which combines low order finite element and finite volume techniques, and is very close to a scheme which was proposed for the solution of barotropic Navier-Stokes equations in [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF] and further extended to two-phase flows in [START_REF] Gastaldo | An entropy-preserving finite element-finite volume pressure correction scheme for the drift-flux model[END_REF]; the resulting code is today currently used at the French Institut de Radioprotection et de Sûreté Nucléaire (IRSN) for "real-life" studies in the nuclear safety field. Up to now, stability (in the sense of conservation of the entropy) is known for these schemes, and numerical experiments show convergence rates close to one in natural energy norms. Our goal is now to prove their convergence. This work is the first one in this direction, and we address here the probably simplest toy problem, restricting ourselves to the steady case, to creeping flows (i.e. omitting the convection term in the momentum balance equation) and to a linear equation of state. The extension to laws where ρ varies linearly with p 1/γ , where γ > 1 is a coefficient which is specific to the considered fluid, which are typically obtained for isentropic flows of perfect gases, is the object of a further paper (part II of the present one); the additional difficulty posed by this further study is to prove the strong convergence for the density, which necessitates to adapt P.L. Lions' "effective pressure trick" [START_REF] Lions | Mathematical topics in fluid mechanics -volume 2 -compressible models[END_REF] at the discrete level. Finally, for the sake of simplicity, we use here a simplified form of the diffusion term (-∆u) but it is clear from the subsequent developements that the presented theory holds for any linear elliptic operator (and in particular for the usual form of the viscous term for compressible constant viscosities flows).

The finite element -finite volume discretization which is chosen here is motivated by the fact that we wish the approximate density to be positive, as in the continuous model, in order to be compatible with the physics. Moreover, the proof of convergence of a numerical approximation to (1.1) requires estimates on both velocity and pressure or density, and the density positivity is very useful to obtain these estimates. A classical way to ensure positivity is to use a finite volume upwinding technique in the discretization of the term div(ρu). This technique is easily set up if the discrete velocities are located on the edges and densities and pressures at the cell centres, which is the reason why we choose the Crouzeix-Raviart finite elements for the velocities and cell centred finite volumes for the densities. This paper is organized as follows. The discretization is first described (section 2), and we prove an L 2 compactness result for sequences of Crouzeix-Raviart functions with bounded broken H 1 semi-norm (section 3). Then the proposed scheme is given (section 4), and the above-mentioned compactness result yields the convergence of (sub-)sequences of discrete solutions to a limit, thanks to a priori estimates which are given in section 5. Finally, this limit is shown to be a solution to the continuous problem in section 6.

To our knowledge, this convergence proof is the first one for the genuine (nonlinear) compressible Stokes problem; a linearized version of this system is adressed in previous works [START_REF] Kellog | A finite element method for the compressible Stokes equations[END_REF][START_REF] Kellog | A penalized finite-element method for a compressible Stokes system[END_REF][START_REF] Kweon | An optimal order convergence for a weak formulation of the compressible Stokes system with inflow boundary condition[END_REF][START_REF] Kweon | Optimal error estimate for a mixed finite element method for compressible Navier-Stokes system[END_REF][START_REF] Bochev | Analysis and computation of least-squares methods for a compressible Stokes problem[END_REF].

Discrete spaces and relevant lemmata

Let T be a decomposition of the domain Ω in simplices. By E(K), we denote the set of the edges (d = 2) or faces (d = 3) σ of the element K ∈ T ; for short, each edge or face will be called an edge hereafter. The set of all edges of the mesh is denoted by E; the set of edges included in the boundary of Ω is denoted by E ext and the set of internal ones (i.e. E \ E ext ) is denoted by E int . The decomposition T is supposed to be regular in the usual sense of the finite element literature (e.g. [START_REF] Ciarlet | Handbook of numerical analysis volume II : Finite elements methods -Basic error estimates for elliptic problems[END_REF]), and, in particular, T satisfies the following properties: Ω = K∈T K; if K, L ∈ T , then K ∩ L = ∅, K ∩ L = ∅ is a vertex of the mesh or K ∩ L is a common edge of K and L, which is denoted by K|L. For each internal edge of the mesh σ = K|L, n KL stands for the normal vector of σ, oriented from K to L (so that n KL = -n LK ). By |K| and |σ| we denote the measure, respectively, of the element K and of the edge σ, and h K and h σ stand for the diameter of K and σ, respectively. We measure the regularity of the mesh through the parameter θ defined by:

(2.1) θ = inf { ξ K h K ; K ∈ T } ∪ { h L h K , h K h L ; σ = K|L ∈ E int },
where ξ K stands for the diameter of the largest ball included in K. Note that the following inequality holds:

(2.2)

h σ |σ| ≤ 2 θ -d |K|, ∀K ∈ T , ∀σ ∈ E(K).
Indeed, this relation is derived by noting that h σ |σ| ≤ h d K and |K| ≥ c ξ d K with c = π/4 in 2D and c = π/6 in 3D; it will be used throughout this paper. Finally, as usual, we denote by h the quantity max K∈T h K .

The space discretization relies on the Crouzeix-Raviart element (see [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations I[END_REF] for the seminal paper and, for instance, [5, p. 199-201] for a synthetic presentation). The reference element is the unit d-simplex and the discrete functional space is the space P 1 of affine polynomials. The degrees of freedom are determined by the following set of nodal functionals:

(2.3) {F σ , σ ∈ E(K)} , F σ (v) = |σ| -1 σ v dγ.
The mapping from the reference element to the actual one is the standard affine mapping. Finally, the continuity of the average value of the discrete functions (i.e., for any function v, F σ (v)) across each face of the mesh is required, thus the discrete space V h is defined as follows:

(2.4)

V h = { v ∈ L 2 (Ω) : v| K ∈ P 1 (K), ∀K ∈ T ; F σ (v) continuous across each edge σ ∈ E int ; F σ (v) = 0, ∀σ ∈ E ext }.
The space of approximation for the velocity is the space W h of vector valued functions each component of which belongs to V h : W h = (V h ) d . The pressure is approximated by the space L h of piecewise constant functions:

L h = q ∈ L 2 (Ω) : q| K = constant, ∀K ∈ T .
Since only the continuity of the integral over each edge of the mesh is imposed, the functions of V h are discontinuous through each edge; the discretization is thus nonconforming in H 1 (Ω) d . We then define, for 1 ≤ i ≤ d and v ∈ V h , ∂ h,i v as the function of L 2 (Ω) which is equal to the (piecewise constant) derivative of v with respect to the i th space variable almost everywhere. This notation allows to define the discrete gradient, denoted by ∇ h , for both scalar and vector valued discrete functions and the discrete divergence of vector valued discrete functions, denoted by div h .

The Crouzeix-Raviart pair of approximation spaces for the velocity and the pressure is inf-sup stable, in the usual sense for "piecewise H 1 " discrete velocities, i.e. there exists c i > 0 independent of the mesh such that:

∀q ∈ L h , sup v∈W h K∈T K q divv dx ||v|| 1,b = sup v∈W h Ω q div h v dx ||v|| 1,b ≥ c i ||q -q m || L 2 (Ω) ,
where q m is the mean value of q over Ω and || • || 1,b stands for the broken Sobolev H 1 semi-norm, which is defined for any function v ∈ V h or v ∈ W h by:

||v|| 2 1,b = K∈T K |∇v| 2 dx = Ω |∇ h v| 2 dx.
This broken Sobolev semi-norm is known to control the L 2 norm by an extended Poincaré inequality [19, proposition 4.13], in the sense that for any function

v ∈ V h , ||v|| 1,b ≤ c ||v|| L 2 (Ω)
where the real number c only depends on the computational domain.

We also define a discrete semi-norm on L h , similar to the usual finite volume discrete H 1 semi-norm, weighted by a mesh-dependent coefficient:

∀q ∈ L h , |q| 2 T,β = σ∈Eint, σ=K|L (h K + h L ) β |σ| h σ (q K -q L ) 2 .
From the definition (2.3), each velocity degree of freedom can be associated to an element edge. Consequently, the velocity degrees of freedom are indexed by the number of the component and the associated edge, thus the set of velocity degrees of freedom reads: {u σ,i , σ ∈ E int , 1 ≤ i ≤ d}. We denote by φ σ the usual Crouzeix-Raviart shape function associated to σ, i.e. the scalar function of V h such that F σ (φ σ ) = 1 and F σ ′ (φ σ ) = 0, ∀σ ′ ∈ E \ {σ}.

Similarly, each degree of freedom for the pressure is associated to a cell K, and the set of pressure degrees of freedom is denoted by {p K , K ∈ T }.

Finally, we define by r h the following interpolation operator:

(2.5)

r h : H 1 0 (Ω) -→ V h v → r h v = σ∈E F σ (v) φ σ = σ∈E |σ| -1 σ v dγ φ σ .
This operator naturally extends to vector-valued functions (i.e. to perform the interpolation from H 1 0 (Ω) d to W h ), and we keep the same notation r h for both the scalar and vector case. The properties of r h are gathered in the following lemma. They are proven in [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations I[END_REF].

Lemma 2.1. Let θ 0 > 0 and let T be a triangulation of the computational domain Ω such that θ ≥ θ 0 , where θ is defined by (2.1). The interpolation operator r h enjoys the following properties:

(1) Preservation of the divergence:

(2.6) ∀v ∈ H 1 0 (Ω) d , ∀q ∈ L h , Ω q div h (r h v) dx = Ω q divv dx.
(2) Stability:

(2.7) ∀v ∈ H 1 0 (Ω), ||r h v|| 1,b ≤ c 1 (θ 0 ) |v| H 1 (Ω) .
(3) Approximation properties:

(2.8) ∀v ∈ H 2 (Ω) ∩ H 1 0 (Ω), ∀K ∈ T , ||v -r h v|| L 2 (K) + h K ||∇ h (v -r h v)|| L 2 (K) ≤ c 2 (θ 0 ) h 2 K |v| H 2 (K) .
In both above inequalities, the notation c i (θ 0 ) means that the real number c i only depends on θ 0 , and, in particular, not on the parameter h characterizing the size of the cells; this notation will be kept throughout the paper.

The following lemma is known (e.g. [5, lemma 3.32]); we give its (elementary) proof for the sake of completeness.

Lemma 2.2. Let θ 0 > 0 and let T be a triangulation of the computational domain Ω such that θ ≥ θ 0 , where θ is defined by (2.1), and V h be the space of Crouzeix-Raviart discrete functions associated to T , as defined by (2.4). Then there exists a real number c(θ 0 ) such that the following bound holds for any v ∈ V h :

σ∈E 1 h σ σ [v] 2 dγ ≤ c(θ 0 ) ||v|| 2 1,b ,
where, on any σ ∈ E int , [v] stands for the jump of v across σ and, on any

σ ∈ E ext , [v] = v.
Proof. For any control volume K of the mesh, we denote by (∇v) K the (constant) gradient of the restriction of v to K. With this notation, using the continuity of v across σ at the mass center x σ of any internal edge and the fact that v vanishes at the mass center x σ of any external edge, we get:

σ∈E 1 h σ σ [v] 2 dγ = σ∈Eint, σ=K|L 1 h σ σ (((∇v) K -(∇v) L ) • (x -x σ )) 2 dγ + σ∈Eext, σ∈E(K) 1 h σ σ ((∇v) K • (x -x σ )) 2 dγ.
We thus have:

σ∈E 1 h σ σ [v] 2 dγ ≤ 2 σ∈Eint, σ=K|L h σ |σ| (|(∇v) K | 2 + |(∇v) L | 2 ) + σ∈Eext, σ∈E(K) h σ |σ| |(∇v) K | 2 .
and the result follows by regularity of the mesh.

The proof of the following trace lemma can be found in [21, section 3].

Lemma 2.3. Let T be a given triangulation of Ω and K be a control volume of T , h K its diameter and σ one of its edges. Let v be a function of H 1 (K). Then the following inequality holds:

||v|| L 2 (σ) ≤ d |σ| |K| 1/2 ||v|| L 2 (K) + h K ||∇v|| L 2 (K) .
We will also need the following Poincaré ineqality:

(2.9)

∀K ∈ T , ∀v ∈ H 1 (K), ||v -v m,K || L 2 (K) ≤ 1 π h K ||∇v|| L 2 (K) .
where v m,K stands for the mean value of v over K. This relation is proven for any convex domain in [START_REF] Payne | An optimal Poincaré-inequality for convex domains[END_REF].

We are now in position to prove the following technical result.

Lemma 2.4. Let θ 0 > 0 and let T be a triangulation of the computational domain Ω such that θ ≥ θ 0 , where θ is defined by (2.1); let (a σ ) σ∈Eint be a family of real numbers such that ∀σ ∈ E int , a σ ≤ 1 and let v be a function of the Crouzeix-Raviart space V h associated to T . Then the following bound holds:

σ∈Eint σ a σ [v] f dγ ≤ c(θ 0 ) h ||v|| 1,b |f | H 1 (Ω) , ∀f ∈ H 1 0 (Ω).
where the real number c(θ 0 ) only depends on θ 0 and on the domain Ω.

Proof. Since the integral of the jump across any edge of the mesh of a function of V h is zero, we have, for any σ ∈ E int :

σ a σ [v] f dγ = σ a σ [v] (f -f σ ) dγ,
where f σ is any real number. Using the Cauchy-Schwarz inequality, first in L 2 (σ) then in R card(E) we thus get:

σ∈Eint σ a σ [v] f dγ ≤ σ∈Eint σ [v] 2 dγ 1/2 σ (f -f σ ) 2 dγ 1/2 ≤ σ∈Eint 1 h σ σ [v] 2 dγ 1/2 σ∈Eint h σ σ (f -f σ ) 2 dγ 1/2 T 1 .
By lemma 2.2, the first term of the latter product is bounded by c(θ 0 ) ||v|| 1,b . For the second one, choosing arbitrarily one adjacent simplex to each edge and applying the above trace lemma 2.3, we get:

T 2 1 ≤ σ∈Eint (σ∈E(K)) 2d h σ |σ| |K| ||f -f σ || 2 L 2 (K) + h 2 K ||∇f || 2 L 2 (K) .
Choosing for f σ the mean value of f on K and using (2.9), we thus get:

T 2 1 ≤ σ∈Eint (σ∈E(K)) 2d (1 + 1 π 2 ) h σ |σ| |K| h 2 K ||∇f || 2 L 2 (K) .
and the result follows by observing that the H 1 semi-norm of f on K appears at most (d + 1) times in the summation and using the regularity of the mesh.

A compactness result

The aim of this section is to state and prove a compactness result for || • || 1,b bounded sequences of discrete functions. We begin by a preliminary lemma. Lemma 3.1. Let θ 0 > 0 and let T be a triangulation of the computational domain Ω such that θ ≥ θ 0 , where θ is defined by (2.1); for σ ∈ E, let χ σ be the function defined by:

χ σ : R d × R d -→ {0, 1} (x, y) → χ σ (x, y) = 1 if [x, y] ∩ σ = ∅, χ σ (x, y) = 0 otherwise,
where x and y are two points of R d . Then there exists a family of positive real numbers (d σ ) σ∈E such that:

(1) for any σ ∈ E,

d σ = c 1 (θ 0 ) h σ , (2
) for any points x and y of R d (possibly located outside Ω), the following inequality holds:

σ∈E χ σ (x, y) d σ ≤ c 2 (θ 0 ) (|y -x| + h)
Proof. We first deal with the two-dimensional case and with quasi-uniform meshes (i.e. the bound we first prove blows up when max K∈T (h/h K ) tends to infinity).

Let T be a triangulation of a two-dimensional domain Ω, K a triangular cell of T and σ an edge of K. Without loss of generality, we suppose that σ is the segment (0, h σ ) × 0 and we denote by ξ K the diameter of the largest ball included in K and by h K the diameter of K. We denote by z σ the opposite vertex to σ; the first coordinate of z σ is necessarily lower than h K while its second coordinate is necessarily greater than ξ K (in the opposite case, no ball of diameter ξ K would be included in K). It thus follows (see figure 1):

(1) that the rectangular domain

ω σ = (h σ /3, 2h σ /3) × (0, h σ ξ K /(12h K )) is included in K, (2
) that, if the similar construction is performed for another edge σ ′ of K to obtain ω σ ′ , ω σ and ω ′ σ do not intersect. We denote by d σ the quantity d σ = h σ ξ K /(12h K ). We thus have d σ ≥ (θ/12) h σ , where θ is the parameter defined by 2.1.

We now perform this construction for each edge σ of the mesh. If σ ∈ E ext , there is only one possible choice for K (the adjacent cell to σ); if σ ∈ E int , σ = K|L, we choose either K or L. Let x and y be two points of R 2 . Let t (x,y) be the vector given by: t (x,y) = yx |y -x| and n (x,y) a normal vector to t (x,y) . We denote by S (x,y) the rectangle defined by:

S (x,y) = x + ξ 1 t (x,y) + ξ 2 n (x,y) , ξ 1 ∈ (-h, |y -x| + h), ξ 2 ∈ (-h, +h)
For each edge intersected by the segment [x, y] (i.e. for each edge σ such that χ σ (x, y) = 1), the rectangle ω σ is included in S (x,y) ; thus, since these domains ω σ and ω σ ′ are disjoint:

σ∈E χ σ (x, y) |ω σ | ≤ |S (x,y) | K ω σ ξ K h σ 12 h K ω σ ′ σ ′ z σ ( h σ 3 , 0) ( 2 h σ 3 , 0) h σ ≥ ξ K ≥ ξ K h σ 3 h K Figure 1.
Notations for the control volume K and thus:

σ∈E χ σ (x, y) 1 3 d σ h σ ≤ 2 h (|y -x| + 2h),
which concludes the proof if max K∈T (h/h K ) is supposed to be bounded.

The extension to the three-dimensional case is straightforward, since it only necessitates to adapt the construction of the domains ω σ . Finally, giving up the assumption that max K∈T (h/h K ) is bounded only needs a more careful definition of the domain S (x,y) , replacing the parameter h by a local value.

The following bound provides an estimate of the translates of a discrete function v as a function of ||v|| 1,b . Lemma 3.2. Let θ 0 > 0 and let T be a triangulation of the computational domain Ω such that θ ≥ θ 0 , where θ is defined by (2.1); let V h be the space of Crouzeix-Raviart discrete functions associated to T , as defined by (2.4). Let v be a function of V h ; we denote by ṽ the extension by zero of v to R d . Then the following estimate holds:

∀η ∈ R d , ||ṽ(• + η) -ṽ(•)|| 2 L 2 (R d ) ≤ c(θ 0 ) |η| (|η| + h) ||v|| 2 1,b .
Proof. We follow the proof of a similar result for piecewise constant functions, namely [6, Lemma 9.3, pp. 770-772]. Let η ∈ R d be given, v be a Crouzeix-Raviart discrete function and ṽ its extension by zero to R d . With the definition of the function χ σ of Lemma 3.1, the following identity holds for any x ∈ R d :

ṽ(x + η) -ṽ(x) = σ∈E χ σ (x, x + η) [v](y x,η,σ ) T 1 (x) + 1 0 ∇ h ṽ(x + sη) • η ds T 2 (x)
where y x,η,σ stands for the intersection between the line issued from x and of direction η and the hyperplane containing σ. Defining for each edge σ of the mesh a real positive number d σ such that Lemma 3.1 holds, by the Cauchy-Schwarz inequality, we get for T 1 (x):

(T 1 (x)) 2 ≤ σ∈E χ σ (x, x + η) [v](y x,η,σ ) 2 d σ σ∈E χ σ (x, x + η) d σ
Integrating now over R d , we thus obtain:

R d (T 1 (x)) 2 dx ≤ c 2 (θ 0 ) (|η| + h) σ∈E R d χ σ (x, x + η) [v](y x,η,σ ) 2 d σ dx Let Q σ,η = {x = y + sη; y ∈ σ and s ∈ [-1, 0]}. Noting that the function x → χ σ (x, x + η) is in fact the characteristic function of Q σ,η
, we get that:

R d χ σ (x, x + η) [v](y x,η,σ ) 2 dx = Qσ,η [v](y x,η,σ ) 2 dx = |n σ • η| 0 -1 σ ([v](y)) 2 dy ds,
where n σ is a unit normal vector to σ. Therefore:

R d (T 1 (x)) 2 dx ≤ c 2 (θ 0 ) (|η| + h) |η| σ∈E 1 d σ σ ([v](y)) 2 dy,
and thus, by choice of d σ :

(3.1)

R d (T 1 (x)) 2 dx ≤ c 2 (θ 0 ) c 1 (θ 0 ) (|η| + h) |η| σ∈E 1 h σ σ ([v](y)) 2 dy.
On the other hand, by the Cauchy-Schwarz inequality, we have for T 2 :

(T 2 (x)) 2 ≤ |η| 2 1 0 |∇ h ṽ(x + sη)| 2 ds,
and thus, using the Fubini theorem and remarking that ∇ h ṽ vanishes outside Ω, we get:

(3.2) R d (T 2 (x)) 2 dx ≤ |η| 2 ||v|| 2 1,b .
The result then follows thanks to the inequality |ṽ(x

+ η) -ṽ(x)| 2 ≤ 2(T 1 (x)) 2 + 2(T 2 (x)) 2
, to the bounds (3.1) and (3.2) and to Lemma 2.2.

We are now in position to state the following compactness result.

Theorem 3.3. Let (v (m)
) m∈N be a sequence of functions satisfying the following assumptions:

(1) ∀m ∈ N, there exists a triangulation of the domain

T (m) such that v (m) ∈ V (m) h
, where V (m) h

is the Crouzeix-Raviart space associated to T (m) defined by (2.4), and the parameter θ (m) characterizing the regularity of T (m) is bounded away from zero independently of m, (2) the sequence (v (m) ) m∈N is uniformly bounded with respect to the broken Sobolev H 1 semi-norm, i.e.:

∀m ∈ N, ||v (m) || 1,b ≤ C
where the real number C does not depend on m and || • || 1,b stands for the broken Sobolev H 1 semi-norm associated to T (m) (with a slight abuse of notation, namely dropping for short the index (m) pointing the dependence of the norm with respect to the mesh). Then, possibly up to the extraction of a subsequence, the sequence (v (m) ) m∈N converges strongly in L 2 (Ω) to a limit v such that v ∈ H 1 0 (Ω). Proof. The result follows from the translates estimates of lemma 3.2. The compactness in L 2 (Ω) of the sequence is a consequence of the Kolmogorov theorem (see e.g. [6, theorem 14.1, p. 833] for a statement of this result). The fact that the limit belongs to H 1 0 (Ω) follows from the particular expression for the bound of the translates and is proven in [6, theorem 14.2, pp. 833-834].

The numerical scheme

Let ρ * be the mean density, i.e. ρ * = M/|Ω| where |Ω| stands for the measure of the domain Ω. We consider the following numerical scheme for the discretization of Problem (1.1):

∀v ∈ W h , Ω ∇ h u : ∇ h v dx - Ω p div h v dx = Ω f • v dx, (4.1a) ∀K ∈ T , σ=K|L v + σ,K ̺(p K ) -v - σ,K ̺(p L ) + h α |K| (̺(p K ) -ρ * ) T stab,1 (4.1b)
+ σ=K|L (h K + h L ) β |σ| h σ |̺(p K ) + ̺(p L )| (̺(p K ) -̺(p L )) T stab,2 = 0, where v + σ,K and v - σ,K stands respectively for v + σ,K = max(v σ,K , 0) and v - σ,K = -min(v σ,K , 0) with v σ,K = |σ| u σ • n KL = v + σ,K -v - σ,K . Note that the up- winded convection term σ=K|L v + σ,K ̺(p K ) -v - σ,K ̺(p L
) may also be written:

σ=K|L v σ,K ρ σ , with (4.2) ρ σ = ρ K if v σ,K ≥ 0, ρ L otherwise.
Equation (4.1a) may be considered as the standard finite element discretization of (1.1a). Since the pressure is piecewise constant, the finite element discretization of (1.1b), i.e. the mass balance, is similar to a finite volume formulation, in which we introduce the standard first-order upwinding and two stabilizing terms. The first one, i.e. T stab,1 , guarantees that the integral of the density over the computational domain is always M (this can easily be seen by summing the second relation for K ∈ T ). The second one, i.e. T stab,2 , is useful in the convergence analysis. It may be seen as a finite volume analogue of a continuous term of the form div (|ρ|∇ρ) weighted by a mesh-dependent coefficient tending to zero as h β ; note, however, that h σ is not the distance which is usually encountered in the finite volume discretization of diffusion terms; consequently, the usual restriction for the mesh when diffusive terms are to be approximated by the two-points finite volume method (namely,the Delaunay condition) is not required here. We suppose that α ≥ 1 and the convergence analysis uses 0 < β < 2.

Remark 4.1. At first glance, leaving the weight |ρ| out, the stabilization term T stab,2 may look as a Brezzi-Pitkäranta regularisation [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF], as used in [START_REF] Eymard | On a stabilized colocated finite volume scheme for the Stokes problem[END_REF] for stabilizing the colocated approximation of the Stokes problem, which would be rather puzzling since we use here an inf-sup stable pair of approximation spaces. However, using the equation of state (1.2), we obtain:

T stab,2 = A 2 σ=K|L (h K + h L ) β |σ| h σ |p K + p L | (p K -p L )
which shows, since A 2 = Ma 4 , that this term rapidly vanishes when approaching the incompressible limit.

Existence of a solution and a priori estimates

The existence of a solution to (4.1) follows, with minor changes to cope with the diffusion stabilization term, from the theory developed in [10, section 2]. In this latter paper, it is obtained for fairly general equations of state by a topological degree argument. We only give here the obtained result, together with a proof of the a priori estimates verified by the solution, and we refer to [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF] for the proof of existence.

Theorem 5.1. Let θ 0 > 0 and let T be a triangulation of the computational domain Ω such that θ ≥ θ 0 , where θ is defined by (2.1). Problem (4.1) admits at least one solution (u, p) ∈ W h × L h ; any possible solution satisfies p K > 0, ∀K ∈ T and:

(5.1) ||u|| 1,b + ||p|| L 2 (Ω) + ||ρ|| L 2 (Ω) + |ρ| T ,β ≤ C(f , M )
where C(f , M ) ∈ R only depends on Ω, A, f , M and θ 0 .

Proof. Let (u, p) ∈ W h × L h be a solution to (4.1). Let ρ K = ̺(p K ) for any K ∈ T , and let ρ denote the vector (ρ K ) K∈T . A natural ordering of the equations and unknowns in (4.1b) leads to a linear system of the form M ρ = c, where c ∈ R N , N is the number discretization cells, c ∈ R N , c > 0, and where M is an M -matrix (in particular M -1 ≥ 0 and M -t ≥ 0). Therefore the i-th component of ρ reads

ρ i = M -1 c • e i = c • M -t e
i where e i is the i-th canonical basis vector of R N . Since M -t ≥ 0, we get M -t e i ≥ 0, and since M -t e i = 0, this proves that ρ i > 0, which, in turns, yields p K > 0, ∀K ∈ T . Let us then prove the estimate (5.1). To this end, we take v = u in (4.1a) and obtain:

(5.2)

Ω |∇ h u| 2 dx - Ω p div h u dx = Ω f • u dx.
Let us then multiply (4.1b) by A -1 [1 + log(ρ K )] (see remark 5.2 below for an explanation of this choice) and we sum over K ∈ T ; dropping the terms which vanish by conservativity, we then obtain:

T 1 + T 2 + T 3 = 0 with: T 1 = A -1 K∈T log(ρ K ) σ=K|L v + σ,K ρ K -v - σ,K ρ L , T 2 = A -1 h α K∈T |K| [1 + log(ρ K )] [ρ K -ρ * ] , T 3 = A -1 K∈T log(ρ K ) σ=K|L (h K + h L ) β |σ| h σ (ρ K + ρ L ) (ρ K -ρ L ) ,
where the term |̺(p K )+̺(p L )| has been replaced by ̺(p K )+̺(p L ) in (4.1b), thanks to the positivity of the pressure. Let us first write T 1 as:

T 1 = A -1 K∈T log(ρ K ) σ=K|L v σ,K ρ σ ,
where ρ σ is the upwind choice defined by (4.2). Adding and substracting the same quantity, T 1 equivalently reads:

T 1 = A -1 K∈T ρ K σ=K|L v σ,K + A -1 K∈T σ=K|L v σ,K (ρ σ log(ρ K ) -ρ K ).
In the first summation, we recognize Ω p div h u dx. A reordering of the second summation yields:

T 1 = Ω p div h u dx + A -1 σ∈Eint, σ=K|L v σ,K [(ρ σ log(ρ K ) -ρ K ) -(ρ σ log(ρ L ) -ρ L )] .
Let ρσ be defined by ρσ

= ρ K = ρ L if ρ K = ρ L , ρσ log(ρ K ) -ρ K = ρσ log(ρ L ) -ρ L otherwise.
With this notation, we get:

T 1 = Ω p div h u dx + A -1 σ∈Eint, σ=K|L v σ,K (ρ σ -ρσ ) (log(ρ K ) -log(ρ L )).
In the last summation, we can, without loss of generality, choose the orientation of each edge in such a way that v σ,K ≥ 0. With this convention, the term in the summation reads v σ,K (ρ Kρσ ) (log(ρ K )log(ρ L )), and is non-negative thanks to the fact that ρ σ ∈ [min(ρ K , ρ L ), max(ρ K , ρ L )] and the log function is increasing. We thus finally obtain:

(5.3) T 1 ≥ Ω p div h u dx.
Let us now turn to the estimate of T 2 . Since the function z → z log(z) is convex for positive z and its derivative is z → 1 + log(z), we simply have:

(5.4)

T 2 ≥ A -1 h α K∈T |K| [ρ K log(ρ K ) -ρ * log(ρ * )] .
Finally, reordering the sums, the term T 3 reads:

T 3 = A -1 σ∈Eint, σ=K|L (h K + h L ) β |σ| h σ (ρ K + ρ L ) (ρ K -ρ L ) (log(ρ K ) -log(ρ L )) .
By concavity of the log function, we have:

| log(ρ K ) -log(ρ L )| ≥ 1 max(ρ K , ρ L ) |ρ K -ρ L |,
and thus:

(5.5)

T 3 ≥ A -1 σ∈Eint, σ=K|L (h K + h L ) β |σ| h σ (ρ K -ρ L ) 2 .
Summing equations (5.2)-(5.5) and using Young's inequality, we obtain:

||u|| 1,b + A -1/2 |ρ| T ,β ≤ C(f , M ).
Furthermore, summing (4.1b) over K ∈ T , we obtain that the mean value of the pressure p m is given by:

p m = 1 |Ω| Ω p dx = A -1 ρ * .
Using the inf-sup stability of the discretization, we get on the other hand:

||p -p m || L 2 (Ω) ≤ 1 c i sup v∈W h 1 ||v|| 1,b Ω p div h v dx = 1 c i sup v∈W h 1 ||v|| 1,b Ω (∇ h u : ∇ h v -f • v) dx,
and the control of ||p|| L 2 (Ω) (or, equivalently, ||ρ|| L 2 (Ω) ) follows from the estimate for ||u|| 1,b .

Remark 5.2 (On the choice of (log(ρ K )) K∈T as test function). At first glance, the choice of log(ρ K ) to multiply (4.1b) in the preceding proof may seem rather puzzling. In fact, this computation is a particular case of the so-called "elastic potential identity", which is well-known in the continuous setting and is central in a priori estimates for the compressible Navier-Stokes equations [START_REF] Lions | Mathematical topics in fluid mechanics -volume 2 -compressible models[END_REF][START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF][START_REF] Feireisl | Dynamics of viscous compressible flows[END_REF]]. An analogous identity is proven at the discrete level, for the same discretization as here, in [10, theorem 2.1].

For the particular case under consideration, an elementary explanation of this choice may be given. Indeed, it is crucial in the above proof that the quantity ρσ lies in the interval [min(ρ K , ρ L ), max(ρ K , ρ L )]. Let us suppose, without loss of generality, that 0 < ρ K < ρ L and that, instead of the log function, the computation is performed with a non-specified increasing ond continuously differentiable function f ; then we get:

ρσ = ρ L -ρ K f (ρ L ) -f (ρ K )
.

The condition ρσ ≥ ρ K is equivalent to:

1 ρ K ≥ f (ρ L ) -f (ρ K ) ρ L -ρ K ,
which is verified for f = log by concavity of the latter and, letting ρ L tend to ρ K , yields f ′ (x) ≤ 1/x. Conversely, the condition ρσ ≤ ρ L yields:

1 ρ L ≤ f (ρ L ) -f (ρ K ) ρ L -ρ K ,
which, once again, is verified by the log function, and now implies f ′ (x) ≥ 1/x. This limitation for the choice of the test function is the reason for the expression of the stabilizing diffusion term.

Convergence analysis

In this section, we prove the following convergence result. Theorem 6.1. Let (T (m) ) m∈N be a sequence of triangulations of Ω such that h (m) tends to zero when m tends to +∞. Let us assume that this sequence is regular in the sense that there exists θ 0 > 0 such that θ (m) ≥ θ 0 , ∀m ∈ N, where θ (m) is defined by (2.1). For m ∈ N, we denote by W (m) h and L (m) h the discrete velocity and pressure spaces associated to T (m) and by (u

(m) , p (m) ) ∈ W (m) h × L (m) h
the corresponding solution to (4.1), with α ≥ 1 and 0 < β < 2. Then, up to a subsequence, the sequence (u (m) ) m∈N strongly converges to a limit ū in L 2 (Ω) d and (p (m) ) m∈N converges to p weakly in L 2 (Ω), where the pair (ū, p) is a solution to (1.1) in the following weak sense:

ū ∈ H 1 0 (Ω) d , p ∈ L 2 (Ω) and : Ω ∇ū : ∇ψ dx - Ω p divψ dx = Ω f • ψ dx, ∀ψ ∈ C ∞ c (Ω) d , (6.1a)
Ω p ū • ∇ψ dx = 0, ∀ψ ∈ C ∞ c (Ω), (6.1b) Ω ̺(p) = M. (6.1c)
Proof. The proof is divided in three steps: we first show the existence of the limits ū and p, then we pass to the limit in (4.1a) and (4.1b). Since the equation of state is linear, the last equation is then a straightforward consequence of the weak convergence in L 2 (Ω) of the (sub)sequence (p (m) ) m∈N to p.

Step 1: existence of a limit.

By the a priori estimates of theorem 5.1, we know that:

∀m ∈ N, ||u (m) || 1,b ≤ C(f , M ). The compactness in L 2 (Ω) d of the sequence (u (m)
) m∈N , together with the fact that the limit ū lies in H 1 0 (Ω) d , thus follows by applying theorem 3.3 to each component u

(m) i , 1 ≤ i ≤ d. Once again by theorem 5.1, we have: ∀m ∈ N, ||p (m) || L 2 (Ω) ≤ C(f , M ). which is sufficient to ensure a weak convergence in L 2 (Ω) of the sequence (p (m) ) m∈N to p ∈ L 2 (Ω).
Step 2: passing to the limit in (4.1a).

Let ψ be a function of C ∞ c (Ω) d . We denote by ψ (m) the interpolation of ψ in

W (m) h , i.e. ψ (m) = r (m)
h ψ where the operator r (m) h is defined by (2.5). Taking v = ψ (m) in (4.1a), we get:

Ω ∇ h u (m) : ∇ h ψ (m) dx - Ω p (m) div h ψ (m) dx = Ω f • ψ (m) dx, ∀m ∈ N.
Since the considered interpolation operator preserves the divergence (2.6), we have:

Ω p (m) div h ψ (m) dx = Ω p (m) divψ dx -→ Ω p divψ dx as m -→ +∞.
By the approximation properties of the interpolation operator (2.8) invoked component by component, we have:

Ω f • ψ (m) dx -→ Ω f • ψ dx as m -→ ∞.
Finally, we have:

Ω ∇ h u (m) : ∇ h ψ (m) dx = Ω ∇ h u (m) : ∇ h (ψ (m) -ψ) dx T (m) 1 + Ω ∇ h u (m) : ∇ψ dx T (m) 2
Once again by (2.8), the term T (m) 1 obeys the following estimate:

|T (m) 1 | ≤ ||u (m) || 1,b ||ψ (m) -ψ|| 1,b ≤ c(θ 0 ) h (m) ||u (m) || 1,b |ψ| H 2 (Ω) ,
and thus tends to zero as m tends to +∞. Integrating by parts over each control volume, the term T (m) 2 reads:

T (m) 2 = - Ω u (m) • ∆ψ dx + σ∈E (m) int σ [u (m) ] ∇ψ • n σ dγ,
where n σ is a normal vector to σ, with the same orientation as that of the jump through σ. Applying Lemma 2.4 for each component of ∇ψ, a σ being the relevant component of the normal vector n σ , we get:

σ∈E (m) int σ [u (m) ] ∇ψ • n σ dγ ≤ c(θ 0 ) h (m) ||u (m) || 1,b |ψ| H 2 (Ω) ,
and thus tends to zero, while the first one tends to -Ω ū • ∆ψ dx as m tends to +∞. Since ū ∈ H 1 0 (Ω) d , we may integrate by parts, and collecting the limits, we obtain (6.1a).

Step 3: passing to the limit in (4.1b).

Let ψ be a function of C ∞ c (Ω). Multiplying the second equation of (4.1) by 1/|K| K ψ(x) dx and summing over K ∈ T yields T

T (m) 3 = K∈T (m) 1 |K|   σ=K|L v (m) σ,K ρ (m) σ   K ψ dx T (m) 4 = (h (m) ) α K∈T (m) |K| ρ (m) K -ρ * ψ K T (m) 5 = K∈T (m)   σ=K|L (h K + h L ) β |σ| h σ ρ (m) K + ρ (m) L ρ (m) K -ρ (m) L   ψ K ,
where ρ (m) σ is defined by (4.2) and ψ K stands for the mean value of ψ over K. Let q (m) ∈ W h be defined as

q (m) (x) = σ∈E (m) int u (m) σ ρ (m) σ φ σ (x)
, where φ σ is the Crouzeix-Raviart basis function associated to σ. The divergence of q (m) is a piecewise constant function and reads:

∀K ∈ T (m) , divq (m) = 1 |K| σ=K|L v (m) σ,K ρ (m) σ
a.e. in K.

We thus have for T

:

T (m) 3 = K∈T (m) K ψ divq (m) dx.
Integrating by parts over each control volume, we get:

T (m) 3 = - Ω ∇ψ • q (m) dx + σ∈E (m) int σ ψ [q (m) ] • n σ dγ = - Ω ∇ψ • (ρ (m) u (m) ) dx + σ∈E (m) int σ ψ [q (m) ] • n σ dγ T (m) 6 + Ω ∇ψ • (q (m) -ρ (m) u (m) ) dx T (m) 7
.

By the respectively weak and strong convergence of (ρ (m) ) m∈N and (u (m) ) m∈N to ρ and ū in L 2 (Ω) and L 2 (Ω) d , we have:

Ω ∇ψ • (ρ (m) u (m) ) dx -→ Ω ∇ψ • (ρ ū) dx as m -→ +∞.
By the definition of q (m) , the term T (m) 6

reads:

T (m) 6 = σ∈E (m) int σ ψ [ σ ′ ∈E (m) int u (m) σ ′ ρ (m) σ ′ φ σ ′ (x)] • n σ dγ = T (m) 8 + T (m) 9
.

with:

T (m) 8 = σ∈E (m) int σ ψ ρ (m) σ [ σ ′ ∈E (m) int u (m) σ ′ φ σ ′ (x)] • n σ dγ = σ∈E (m) int ρ (m) σ σ ψ [u (m) ] • n σ dγ T (m) 9 = σ∈E (m) int σ ψ [ σ ′ ∈E (m) int \{σ} u (m) σ ′ (ρ (m) σ ′ -ρ (m) σ ) φ σ ′ (x)] • n σ dγ.
Since the integral of the jump of a Crouzeix-Raviart function over an internal edge of the mesh vanishes, the term T (m) 8

can be estimated as follows:

|T (m) 8 | ≤ c ψ h (m) σ∈E (m) int ρ (m) σ σ [u (m) ] • n σ dγ,
where c ψ only depends on ψ. Using the Cauchy-Schwarz inequality then yields:

|T (m) 8 | ≤ c ψ h (m) σ∈E (m) int |σ| 1/2 ρ (m) σ σ | [u (m) ] | 2 dγ 1/2 ≤ c ψ h (m) σ∈E (m) int h σ |σ| (ρ (m) σ ) 2 1/2 σ∈E (m) int 1 h σ σ | [u (m) ] | 2 dγ 1/2
.

By the regularity of the mesh, the first summation is bounded by ||ρ

ψ (ρ (m) σ ′ -ρ (m) σ ) [φ σ ′ (x)] u (m) σ ′ • n σ dγ ≤ c ψ h σ |σ| |ρ (m) σ ′ -ρ (m) σ | |u (m) σ ′ |,
where c ψ still only depends on ψ. Since the function φ σ ′ is non-zero over σ = K|L only when σ ′ belongs to the edges of K or L, only a limited number of terms are non-zero in T (m) 9

, in such a way that the difference ρ

(m) σ ′ -ρ (m) σ
only involves two neigbouring cells or two cells sharing the same neighbour. Splitting the difference in this last case, using the previous inequality and the regularity of the mesh (in particular the fact that the ratio of the size of two neighbouring cells is bounded) and reordering the sums, we get for T (m) 9

an estimate of the form:

|T (m) 9 | ≤ c σ∈E (m) int h d σ |u (m) σ | σ ′ ∈Nσ (σ ′ =K|L) |ρ (m) K -ρ (m) L |,
where the positive real number c only depends on ψ and the regularity of the mesh and, thanks to this regularity, the set N σ is such that a given edge K|L only appears in this sum a number of times bounded independently of m. Thus, thanks to the Cauchy-Schwarz inequality, we have:

|T (m) 9 | 2 ≤ c σ∈E (m) int h d σ |u (m) σ | 2 σ∈E (m) int (σ=K|L) h d σ ρ (m) K -ρ (m) L 2 .
By the regularity of the mesh, the first term of this product is controlled by ||u (m) || L 2 (Ω) and the second one by (h (m) ) and thus also T tend to zero as m tends to +∞, for any β < 2.

Let us then examine the term T (m) 7

:

T (m) 7 = K∈T (m) K σ=K|L (ρ (m) σ -ρ (m) K ) φ σ (x) u (m) σ • ∇ψ(x) dx
Since ∇ψ is bounded in L ∞ (Ω) d , and since the functions φ σ are bounded, we get:

|T (m) 7 | ≤ c ψ K∈T (m) |K| σ=K|L |ρ (m) σ -ρ (m) K | |u (m) σ |.
Reordering the summations and using the Cauchy-Schwarz inequality yields:

|T (m) 7 | ≤ c ψ σ∈E (m) int (σ=K|L) (|K| + |L|) |ρ (m) K -ρ (m) L | |u (m) σ | ≤ c ψ T (m) 10 1/2 T (m) 11 1/2
, with:

T (m) 10 = σ∈E (m) int (σ=K|L) (|K| + |L|) |u (m) σ | 2 T (m) 11 = σ∈E (m) int (σ=K|L) h σ (h K + h L ) (1-β) |K| + |L| |σ| (h K + h L ) (h K + h L ) β |σ| h σ (ρ (m) K -ρ (m) L ) 2 .
Once again reordering the summation, we get:

T (m) 10 = K∈T |K| σ∈EK |u (m) σ | 2 ,
and thus, the term T , we get: (σ=K|L)

T (m) 5 = σ∈E (m) int (σ=K|L) (h K + h L ) β |σ| h σ ρ (m) K + ρ (m) L ρ (m) K -ρ (m) L (ψ K -ψ L ) .
(h K + h L ) 2 |σ| h σ ρ (m) K + ρ (m) L 2      1/2 |ρ (m) | T ,β ,
which, once again since the sequence (ρ (m) ) m∈N is bounded in L 2 (Ω), tends to zero by regularity of the mesh for β > 0. The proof is thus complete.

Discussion

To our knowledge, the convergence analysis performed in this paper seems to be the first result of this kind for the compressible Stokes problem (and, of course, more widely, for the compressible Navier-Stokes equations). Beside the convergence of the scheme, it also provides an existence result for solutions of the continuous problem, which could also be derived from the continuous existence theory ingredients for the steady Navier-Stokes equations (see [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF] and references therein), but does not seem to be a direct consequence of the published literature: existence of strong solutions of the Navier-Stokes equations is known only for small data (e.g. [START_REF] Valli | On the existence of stationary solutions to compressible Navier-Stokes equations[END_REF]) and existence of weak solutions is only proven for a particular class of equations of state (typically, p = ρ γ with γ > 3/2), this limitation being due to the presence of the convection term.

A puzzling fact is that the present theory relies on two ingredients which are usually not present in actual implementations. Firstly, the stabilisation term T stab,2 is needed in our proof to ensure the convergence of the discretization of the mass convection flux div(ρu) and, to our knowledge, has never been introduced elsewhere. Secondly, the control of the pressure in L 2 (Ω) relies on the stability of the discrete gradient (i.e. the satisfaction of the so-called discrete inf-sup condition), which is not verified by colocated discretizations; note that this argument is not needed for the stability of the scheme (see the proof of a priori estimates here and [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF][START_REF] Eymard | Entropy estimate for the approximation of the compressible barotropic Navier-Stokes equations using a collocated finite volume scheme[END_REF] for stability studies of schemes for the Navier-Stokes equations). Assessing the effective relevance of these requirements for the discretization should deserve more work in the future.

An easy extension of this work consists in replacing the diffusive term -∆u in (1.1a) by its complete expression -µ ∆uµ/3 ∇(divu) with µ > 0 (i.e. the usual form of the divergence of the shear stress tensor in a constant viscosity compressible flow). Another less straightforward issue is the extension to more general state equations (for instance, p = ρ γ with γ > 1); it will be the topic of a further paper. Concerning higher order issues, let us note that the fact that the pressure is approximated by a piecewise constant function appears crucial in both stability and convergence proofs: extending this study to higher degree finite element discretizations thus certainly represents a difficult task. Finally, let us remark that the present scheme relies on the approximation of the whole velocity vector at the interfaces. A less expensive scheme would be possible with a discretization u • n at the interfaces, as in the MAC scheme which is well known for the incompressible Navier-Stokes equations. However, such a discretization does not seem straightforward on unstructured meshes.

2000

  Mathematics Subject Classification. 35Q30,65N12,65N30,76N15,76M10,76M12. Key words and phrases. Compressible Stokes equations, finite element methods, finite volume methods.

9 .

 9 (m) || L 2 (Ω) while, by Lemma 2.2, the second one is bounded by c(θ 0 ) ||u (m) || 2 1,b . Let us now turn to the study of T (m) Since, for σ ′ ∈ E (m) int \ {σ}, the integral of φ σ ′ over σ vanishes, and since the functions φ σ are bounded (namely |φ σ | ≤ 1 in 2D, |φ σ | ≤ 2 in 3D) we get:

  σ

5 .

 5 ||u (m) || L 2 (Ω) , and T (m) 11 is controlled by (h (m) ) 2-β |ρ (m) | T ,β . By the a priori estimate (5.1), T (m) 7thus tends to zero for any β < 2.We now turn to the terms T Since the sequence (ρ (m) ) m∈N is bounded in L 2 (Ω), the term T (m) 4 tends to zero for any α > 0. Reordering the summation in T (m) 5

  By regularity of ψ, |ψ Kψ L | ≤ c ψ (h K + h L ) and thus:

  2-β |ρ (m) | T ,β . Consequently, thanks

	to estimate (5.1), both T	(m) 8	and T	(m) 9

Université de Provence, France E-mail address: gallouet@cmi.univ-mrs.fr Université de Provence, France E-mail address: herbin@cmi.univ-mrs.fr Institut de Radioprotection et de Sûreté Nucléaire (IRSN) E-mail address: jean-claude.latche@irsn.fr