E A 2222

Abstract: We consider the flow of a viscous incompressible fluid in a rigid homogeneous porous medium provided with boundary conditions on the pressure around a circular well. When the boundary pressure presents high variations, the permeability of the medium depends on the pressure, so that the model is nonlinear. We propose a spectral discretization of the resulting system of equations which takes into account the axisymmetry of the domain and of the flow. We prove optimal error estimates and present some numerical experiments which confirm the interest of the discretization.

Résumé: Nous considérons l'écoulement d'un fluide visqueux incompressible dans un milieu poreux rigide lorsque la pression donnée sur la partie de la frontière correspondant à un puits circulaire présente de fortes variations. La perméabilité du milieu dépend alors de la pression, de sorte que le modèle est non linéaire. Nous proposons une discrétisation spectrale de ce modèle qui tient compte de l'axisymétrie du domaine et de l'écoulement. Nous prouvons des estimations d'erreur optimales et présentons quelques expériences numériques qui confirment l'intérêt de la discrétisation.

Introduction.

Let Ω denote the three-dimensional domain between two cylinders, namely the domain Ω = (x, y, z) ∈ R 3 ; r 0 < x 2 + y 2 < r 1 and z 1 < z < 0 ,

where r 0 and r 1 are two positive constants and z 1 a negative constant. Its boundary ∂Ω is divided into two parts Γ w = (x, y, z) ∈ R 3 ; x 2 + y 2 = r 0 and z 1 < z < 0 and Γ = ∂Ω \ Γ w , (

where the index w stands for "well". Indeed, we are interested in the following model, suggested by K.R. Rajagopal [START_REF]On a hierarchy of approximate models for flows of incompressible fluids through porous solid[END_REF],

                 α(p) u + grad p = f in Ω, div u = 0 in Ω, p = p w on Γ w , u • n = g on Γ, (1.3)
where the unknowns are the velocity u and the pressure p of the fluid. This system is an extension of Darcy's equations, which model the flow of an incompressible viscous fluid in a rigid porous medium, to the case where pressure with high variations is enforced on a well. Thus, the permeability of the medium, which is assumed to be homogeneous for simplicity, depends on p in an exponential way, which leads to the nonlinear system (1.3). Indeed, problem (1.3) is a first draft of a model for the flow in the porous medium when steam is injected in the well with very large pressure at the bottom and low pressure at the top; in this case, α is an exponential function of the pressure. We refer to [START_REF] Salgado | Finite element discretization of Darcy's equations with pressure dependent porosity[END_REF] for a first study of the finite element discretization of this problem in general geometry. In this paper, we are interested in its spectral discretization for the domain Ω defined in (1.1).

In a first step, we prove the existence of a solution for problem (1.3) in the case where Ω is a general two-or three-dimensional domain with a Lipschitz-continuous boundary. Next, we consider problem (1.3) in the simpler case of the domain Ω defined in (1.1) where the data are axisymmetric, in the sense of [START_REF] Bernardi | Azaïez -Spectral Methods for Axisymmetric Domains[END_REF]§II.3]. So, by using cylindrical coordinates, we can write a variational formulation of this problem in the meridian domain. We prove the existence of an axisymmetric solution for such a system. Next, we propose a spectral discrete problem which is constructed from this formulation by the Galerkin method with numerical integration. The analysis of this problem relies on the theorem of F. Brezzi, J. Rappaz and P.-A. Raviart [START_REF] Brezzi | Raviart -Finite dimensional approximation of nonlinear problems, Part I: Branches of nonsingular solutions[END_REF], and optimal error estimates are derived. In a final step, we propose an algorithm for solving the resulting system and present some numerical experiments. also thank V. Girault for finding an error in the first version of this paper and M. Dauge for her help in repairing this mistake.

An outline of the paper is as follows.

• In Section 2, we prove the existence of a solution to problem (1.3) in a general threedimensional domain.

• In Section 3, we write the variational formulation of problem (1.3) in the case of an axisymmetric solution and prove the existence of such a solution.

• Section 4 is devoted to the description and numerical analysis of the discrete problem.

• In Section 5, we present the algorithm that is used to solve the nonlinear system, together with some numerical experiments.

The problem in general geometry.

In this section, we assume that Ω is a bounded connected domain in R d , d = 2 or 3, with a Lipschitz-continuous boundary ∂Ω divided in two parts Γ w and Γ = ∂Ω \ Γ w such that (i) the intersection Γ w ∩ Γ is a Lipschitz-continuous submanifold of ∂Ω, (ii) Γ w has a positive (d -1)-measure in ∂Ω. We intend to study problem (1.3) when α is a continuous function from R into R and satisfies for two positive constants α 0 and α 1 ,

∀ξ ∈ R, α 0 ≤ α(ξ) ≤ α 1 . (2.1)
Even if this is not true when the function α is exponential, it does not seem restrictive to make this assumption (which is easily recovered by truncating the exponential), since in practical situations the pressure is always bounded.

We consider the full scale of Sobolev spaces H s (Ω), s ∈ R, and W m,p (Ω), m ∈ N, 1 ≤ p ≤ ∞, equipped with the standard norms and seminorms (see [2, Chap. III & VII] for instance). In order to write a variational formulation of problem (1.3), we introduce the space

H 1 w (Ω) = q ∈ H 1 (Ω); q = 0 on Γ w . (2.2)
Note that the traces of functions in H1 w (Ω) on Γ belong to H

We recall from [7, §XIII.1] that Darcy's equations even for a constant function α admit several variational formulations. We have chosen here the formulation which enables us to treat the boundary condition on p as an essential one and also seems the best adapted for handling the nonlinear term α(p) u (see [START_REF] Achdou | Coquel -A priori and a posteriori analysis of finite volume discretizations of Darcy's equations[END_REF] for more numerical reasons). So, we consider the variational problem

Find (u, p) in L 2 (Ω) d × H 1 (Ω) such that p = p w on Γ w , (2.3
)

and ∀v ∈ L 2 (Ω) d , a [p] (u, v) + b(v, p) = Ω f (x) • v(x) dx, ∀q ∈ H 1 w (Ω), b(u, q) = g, q Γ , (2.4)
where the bilinear forms a [ξ] (•, •) for any measurable function ξ on Ω and b(•, •) are defined by

a [ξ] (u, v) = Ω α ξ(x) u(x) • v(x) dx, b(v, q) = Ω v(x) • (grad q)(x) dx. (2.5)
Here, •, • Γ denotes the duality pairing between H 1 2 00 (Γ) and its dual space H

1 2 00 (Γ) ′ .
It is readily checked that the forms a

[ξ] (•, •) and b(•, •) are continuous on L 2 (Ω) d × L 2 (Ω) d and L 2 (Ω) d × H 1 (Ω)
, respectively. Thus, some density arguments yield the equivalence of this problem with system (1.3). In order to prove the existence of a solution for problem (2.3) -(2.4), we first consider the linear problem associated with a fixed measurable function ξ

Find (u ξ , p ξ ) in L 2 (Ω) d × H 1 (Ω) such that p ξ = p w on Γ w , (2.6
)

and ∀v ∈ L 2 (Ω) d , a [ξ] (u ξ , v) + b(v, p ξ ) = Ω f (x) • v(x) dx, ∀q ∈ H 1 w (Ω), b(u ξ , q) = g, q Γ , (2.7)
Indeed, the following ellipticity property follows from (2.1):

∀v ∈ L 2 (Ω) d , a [ξ] (v, v) ≥ α 0 v 2 L 2 (Ω) d , (2.8)
and the ellipticity constant α 0 is independent of ξ. The following inf-sup condition is easily derived by taking v equal to grad q:

∀q ∈ H 1 w (Ω), sup v∈L 2 (Ω) d b(v, q) v L 2 (Ω) d ≥ |q| H 1 (Ω) , (2.9)
and, since Γ w has a positive measure, the equivalence of the seminorm

, g) in L 2 (Ω) d × H 1 2 (Γ w ) × H 1 2 00 (Γ) ′ , problem (2.6) -(2.7) has a unique solution (u ξ , p ξ ) in L 2 (Ω) d × H 1 (Ω).
Moreover this solution satisfies, for a constant c independent of ξ,

u ξ L 2 (Ω) d + p ξ H 1 (Ω) ≤ c f L 2 (Ω) d + p w H 1 2 (Γ w ) + g H 1 2 00 (Γ) ′ .
(2.10) Thanks to the previous lemma, we are in a position to state the existence result. Its proof relies on Brouwer's fixed point theorem, see [12, Chap. IV, Cor. 1.1], combined with the addition of a penalization term.

Theorem 2.3. For any data (f , p w , g) in L 2 (Ω) d ×H 1 2 (Γ w )×H 1 2 00 (Γ) ′ , problem (2.3)-(2.4) admits a solution (u, p) in L 2 (Ω) d × H 1 (Ω). Moreover this solution satisfies u L 2 (Ω) d + p H 1 (Ω) ≤ c f L 2 (Ω) d + p w H 1 2 (Γ w ) + g H 1 2 00 (Γ) ′ .
(2.11) Proof: We proceed in five steps, the first four ones being devoted to the simpler case where p w = 0. 1) We introduce the space Y(Ω) = L 2 (Ω) d × H 1 w (Ω), provided with the norm

V Y(Ω) = v 2 L 2 (Ω) d + |q| 2 H 1 (Ω) 1 2 , with V = (v, q).
Let ε be a fixed positive real number. Denoting by •, • the scalar product of Y(Ω) and setting U = (u, p) and V = (v, q), we consider the mapping defined from Y(Ω) into itself by

Φ(U ), V = a [p] (u, v) + b(v + ε grad q, p) -b(u, q) - Ω f (x) • v(x) dx + g, q Γ .
It follows from the boundedness of the function α that the function Φ is continuous on Y(Ω) and moreover (see (2.8)) that

Φ(U ), U ≥ α 0 u 2 L 2 (Ω) d + ε |p| 2 H 1 (Ω) -f L 2 (Ω) d u L 2 (Ω) d -c g H 1 2 00 (Γ) ′ |p| H 1 (Ω) ,
where c denotes the norm of the trace operator from

H 1 w (Ω) into H 1 2
00 (Γ). Thus, applying a Cauchy-Schwarz inequality gives

Φ(U ), U ≥ min{α 0 , ε} U 2 Y(Ω) -f 2 L 2 (Ω) d + c 2 g 2 H 1 2 00 (Γ) ′ 1 2 U Y(Ω) .
So, Φ(U ), U is nonnegative on the sphere of Y(Ω) with radius

µ = f 2 L 2 (Ω) d + c 2 g 2 H 1 2 00 (Γ) ′ 1 2 min{α 0 , ε} .
2) Let ( V n ) n and (W n ) n be increasing sequences of finite-dimensional subspaces of L 2 (Ω) d and H 1 w (Ω), respectively, such that ∪ +∞ n=0 V n is dense in L 2 (Ω) d and ∪ +∞ n=0 W n is dense in H 1 w (Ω). We set: V n = V n + grad W n , Y n (Ω) = V n × W n
and observe that, thanks to this choice, the mapping Φ is still continuous from each Y n (Ω) onto itself and satisfies the same nonnegativity property as previously. So, it follows from Brouwer's fixed point theorem [12, Chap. IV, Cor. 1.1] that, for each n, there exists a

U n = (u n , p n ) in Y n (Ω), with U n Y(Ω) ≤ µ, such that ∀V n ∈ Y n (Ω), Φ(U n ), V n = 0.
This last equation can equivalently be written, for all m ≤ n,

∀v m ∈ V m , a [p n ] (u n , v m ) + b(v m , p n ) = Ω f (x) • v m (x) dx, ∀q m ∈ W m , b(u n , q m ) -ε b(grad q m , p n ) = g, q m Γ .
(2.12)

3) Since the sequence (u n , p n ) n is bounded by µ in L 2 (Ω) d × H 1 (Ω),
) v m ) n converges to α(p ε ) v m strongly in L 2 (Ω) d . By writing the decomposition Ω α p n (x) u n (x) • v m (x) dx = Ω u n (x) • α(p n ) -α(p ε ) (x) v m (x) dx + Ω u n (x) • α p ε (x) v m (x) dx, we thus derive lim n→+∞ Ω α p n (x) u n (x) • v m (x) dx = Ω α p ε (x) u ε (x) • v m (x) dx.
The convergence of the other terms in (2.12) (which are all linear) is readily checked. Finally, using the density of

∪ +∞ m=0 Y m (Ω) into Y(Ω), we derive that (u ε , p ε ) is a solution of the problem ∀v ∈ L 2 (Ω) d , a [p ε ] (u ε , v) + b(v, p ε ) = Ω f (x) • v(x) dx, ∀q ∈ H 1 w (Ω), b(u ε , q) -ε b(grad q, p ε ) = g, q Γ . (2.13)
4) It follows from the inf-sup condition (2.9) that

|p ε | H 1 (Ω) ≤ α 1 u ε L 2 (Ω) d + f L 2 (Ω) d ,
and from the ellipticity property (2.8) (note also that b(grad p ε , p ε ) ≥ 0) that

α 0 u ε 2 L 2 (Ω) d ≤ c f L 2 (Ω) d u ε L 2 (Ω) d + g H 1 2 00 (Γ) ′ p ε H 1 (Ω) .
Therefore, the solutions (u ε , p ε ) satisfy

u ε L 2 (Ω) d + p ε H 1 (Ω) ≤ c f L 2 (Ω) d + g H 1 2 00 (Γ) ′ ,
where the constant c is independent of ε. Thus there exists a subsequence denoted by

(u ε n , p ε n ) n which converges to a (u, p) weakly in L 2 (Ω) d × H 1 (Ω) and such that (p ε n ) n converges strongly in L 2 (Ω).
Exactly the same arguments as in part 3) of the proof yield that (u, p) is a solution of problem (2.3) -(2.4) with p w = 0. Estimate (2.11) is then easily derived from (2.8) and (2.9) as in the lines above. 5) Any datum p w in H 1 2 (Γ w ) admits a lifting p w in H 1 (Ω). Setting p 0 = p-p w , we observe that (u, p) is a solution of problem (2.3) -(2.4) if and only if (u, p 0 ) is a solution of the same problem with the right-hand side of the first line replaced by

Ω f (x) • v(x) dx -b(v, p w ),
and α replaced by the function: ξ → α(p w + ξ). Since this new function has exactly the same properties as α, the existence of (u, p 0 ) follows from the first four steps of the proof and estimate (2.11) is derived from a triangle inequality. This concludes the proof.

Unfortunately, the uniqueness result that we now prove only concerns smooth solutions, at least in dimension d = 3. Proposition 2.4. Assume that the function α is uniformly Lipschitz-continuous, with Lipschitz constant γ. Let r be a real number > d. If problem (2.3) -(2.4) admits a solution (u, p) such that u belongs to L r (Ω) d and satisfies

2α 1 α 0 c r γ u L r (Ω) d < 1, (2.14)
for an appropriate constant c r only depending on r and Ω, there is no other solution of problem (2.3) -(2.4).

Proof: Let (ũ, p) be another solution of (2.3) -(2.4). We take r * such that 1 r + 1 r * = 1 2 . The first line of (2.4) yields that, for any

v in L 2 (Ω) d : Ω α(p) ũ -α(p) u (x) • v(x) dx + Ω v(x) • grad (p -p) (x) dx = 0, or equivalently Ω α(p) (x) (ũ -u)(x) • v(x) dx + Ω α(p) -α(p) (x) u(x) • v(x) dx + Ω v(x) • grad (p -p) (x) dx = 0. (2.15)
We also observe from the second line of (2.4) that, since p and p are equal on Γ w ,

Ω (ũ -u)(x) • grad (p -p) (x) dx = 0.
Thus, taking v equal to ũu, we derive from the Lipschitz property of α and Hölder's inequalities that

α 0 ũ -u 2 L 2 (Ω) d ≤ γ p -p L r * (Ω) u L r (Ω) d ũ -u L 2 (Ω) d .
Denoting by c r the norm of the Sobolev imbedding of H 1 w (Ω) into L r * (Ω), we thus obtain

α 0 ũ -u L 2 (Ω) d ≤ c r γ u L r (Ω) d |p -p| H 1 (Ω) . (2.16)
On the other hand, we take v equal to grad (pp) in (2.15), which gives

|p -p| 2 H 1 (Ω) ≤ α 1 ũ -u L 2 (Ω) d |p -p| H 1 (Ω) + γ p -p L r * (Ω) u L r (Ω) d |p -p| H 1 (Ω) .
The same arguments as previously yield

|p -p| H 1 (Ω) ≤ α 1 ũ -u L 2 (Ω) d + c r γ u L r (Ω) d |p -p| H 1 (Ω) .
It follows from (2.14) and the inequality α 0 ≤ α 1 that

|p -p| H 1 (Ω) ≤ 2α 1 ũ -u L 2 (Ω) d .
(2.17)

Finally, inserting (2.17) into (2.16) gives

α 0 ũ -u L 2 (Ω) d ≤ 2α 1 c r γ u L r (Ω) d ũ -u L 2 (Ω) d .
Therefore, (2.14) implies that ũ is equal to u. Then (2.17) yields that p is equal to p.

We conclude with a regularity property of the solution (u, p). Its proof relies on the arguments in [START_REF] Meyers | An L p -estimate for the gradient of solutions of second order elliptic divergence equations[END_REF]. We refer to [START_REF]Dauge -Neumann and mixed problems on curvilinear polyhedra[END_REF]Def. 2.2] for the exact definition of a curvilinear polyhedron.

Proposition 2.5. If Ω is a curvilinear polygon or polyhedron, there exists a real number ρ 0 > 2 only depending on the geometry of Ω such that, for all ρ, 2 < ρ ≤ ρ 0 , and for all data

(f , p w , g) in L ρ (Ω) d × W 1-1 ρ ,ρ (Γ w ) × W -1 ρ ,ρ (Γ), any solution (u, p) of problem (2.3) -(2.4) belongs to L ρ (Ω) d × W 1,ρ (Ω).
Proof: We establish the desired property in several steps, beginning with the case of zero boundary condition. 1) Let D * denote the operator which associates with any f * in L 2 (Ω) d the part u * solution (u * , p * ) of the linear Darcy's system

     u * + grad p * = f * in Ω, div u * = 0 in Ω, p * = 0 on Γ w , u * • n = 0 on Γ, (2.18)
Following the approach in [START_REF] Meyers | An L p -estimate for the gradient of solutions of second order elliptic divergence equations[END_REF], we observe that problem (2.3) -(2.4) in the case of zero boundary conditions p w = g = 0, can equivalently be written as

u -D * (1 - α(p) α 1 ) u = D * f α 1 .
So, in order to prove the desired regularity property, it suffices to check that the operator in the left-hand side of this equation is an automorphism of the space L ρ (Ω) d .

2) Multiplying the first equation in (2.18) by a function grad q, where q vanishes on Γ w , we observe that p * is the solution of the problem:

∀q ∈ H 1 w (Ω), Ω (grad p * )(x) • (grad q)(x) dx = Ω f * (x) • (grad q)(x) dx.
Thus, we derive from [13, §2.4] and [11, §3] that there exists a real number ρ * > 2 such that, for any f * in L ρ * (Ω) d , this p * belongs to W 3) On the other hand, for any function p, the multiplication by 1 -α(p) α 1 is continuous from L r (Ω) d onto itself with norm 1 -α 0 α 1 . Combining all this gives the desired property for all ρ such that (1 - Setting u 0 = u-grad p b and p 0 = p-p b , we obsereve that the pair (u 0 , p 0 ) is a solution of Darcy's system (1.3) but now with zero boundary conditions, α replaced by the mapping α ♯ : ξ → α(p b + ξ) and f replaced by f -1 + α(p) grad p b (which obviously belongs to L ρ (Ω) d ). Since the mapping α ♯ has exactly the same properties as α and in particular still satisfies (2.1), the result in the general case follows from the first parts of the proof.

α 0 α 1 ) χ θ(ρ) < 1. 4) For any pair (p w , g) in W 1-1 ρ ,ρ (Γ w ) × W -1 ρ ,ρ (
Remark 2.6. By combining Propositions 2.4 and 2.5, we observe that, in dimension d = 2, the uniqueness of the solution is ensured for small enough data in

L ρ (Ω) d × W 1-1 ρ ,ρ (Γ w ) × W -1 ρ ,ρ (Γ), ρ > 2.
But this is no longer true in dimension d = 3.

The two-dimensional formulation.

We now consider problem (1.3) in the case of the geometry introduced in (1.1) and (1.2) and for axisymmetric data, in a sense which is made precise later on. Thus, we introduce the cylindrical coordinates r = x 2 + y 2 , θ = arccos y

x when y ≥ 0, -arccos y x when y < 0.

(3.1)

Setting ω =]r 0 , r 1 [×]z 1 , 0[, we observe that Ω = (r, θ, z); (r, z) ∈ ω and -π < θ ≤ π . (3.2)
We also introduce the two parts of the boundary of ω,

γ w = {r 0 }×]z 1 , 0[ and γ = ∂ω \ γ w . (3.3)
Next, we assume that the radial, angular and axial components of f , denoted by f r , f θ and f z , are independent of θ and that p w and g are also independent of θ. Our idea is to look for a solution (u, p) such that the three components u r , u θ , u z of u and p do not depend on θ. Thus, the pair (u, p) satisfies

                                 α(p) u r + ∂ r p = f r in ω, α(p) u θ = f θ in ω, α(p) u z + ∂ z p = f z in ω, ∂ r u r + r -1 u r + ∂ z u z = 0 in ω, p = p w on γ w , u r n r + u z n z = g on γ, (3.4)
where n = (n r , n z ) now denotes the unit outward normal vector to ω. Moreover, we observe that, when the part (u r , u z , p) of the solution is known, the component u θ is obtained by

u θ = 1 α(p) f θ . (3.5)
So from now on we only consider the reduced problem

                         α(p) u r + ∂ r p = f r in ω, α(p) u z + ∂ z p = f z in ω, ∂ r u r + r -1 u r + ∂ z u z = 0 in Ω, p = p w on γ w , u r n r + u z n z = g on γ.
(3.6)

In order to write the variational formulation of problem (3.6) and according to [5, §II.2], we introduce the weighted Sobolev spaces associated with the measure r drdz.

• The space L 2 1 (ω) is defined by

L 2 1 (ω) = v : ω → R measurable; ω v 2 (r, z) r drdz < +∞ , (3.7)
and equipped with the norm

v L 2 1 (ω) = ω v 2 (r, z) r drdz 1 2 . (3.8) • For each nonnegative integer m, H m 1 (ω) is the space of functions v in L 2 1 (ω) such that all their partial derivalives ∂ k r ∂ ℓ-k z v, 0 ≤ k ≤ ℓ ≤ m, belong to L 2 1 (ω)
and is provided with the seminorm and norm

|v| H m 1 (ω) = m k=0 ∂ k r ∂ m-k z 2 L 2 1 (ω) 1 2 , v H m 1 (ω) = m ℓ=0 |v| 2 H ℓ 1 (ω) 1 2 .
(3.9)

• For each positive real number s which is not an integer, H s 1 (ω) is defined by interpolation between the spaces H ⌊s⌋+1 1

(ω) and H ⌊s⌋ 1 (ω), where ⌊s⌋ denotes the integer part of s. It can be noted that, since r 0 is positive, the spaces H s 1 (ω) coincides with H s (ω). However working with the unweighted norms leads to the fact that the ratio r 1 r 0 appears in the estimates and, since this ratio is often very large in practical situations, we have rather avoid that.

We also introduce the space

H 1 1w (ω) = q ∈ H 1 1 (ω); q = 0 on γ w . (3.10)
For simplicity, we still denote by v the pair (v r , v z ). Indeed, it is readily checked from Proposition 2.1 (see [5, §II.2]) that problem (3.6) admits the equivalent variational formulation

Find (u, p) in L 2 1 (ω) 2 × H 1 1 (ω) such that p = p w on γ w , (3.11
)

and ∀v ∈ L 2 1 (ω) 2 , ã[p] (u, v) + b(v, p) = ω f (r, z) • v(r, z) r drdz, ∀q ∈ H 1 1w (ω), b(u, q) = γ g(τ )q(τ ) r(τ ) dτ, (3.12)
where the bilinear forms ã[ξ] (•, •) for any measurable function ξ on ω and b(•, •) are defined by

ã[ξ] (u, v) = ω α ξ(r, z) u r (r, z)v r (r, z) + u z (r, z)v z (r, z) r drdz, b(v, q) = ω v r (r, z)∂ r q(r, z) + v z (r, z)∂ z q(r, z) r drdz. (3.13)
Here and for simplicity, we assume that g belongs to L 2 1 (γ) (with obvious definition for this new space) and r(τ ) denotes the value of r at the point with tangential coordinate τ .

Note that Theorem 2.3 is not sufficient to prove the existence of a solution for problem (3.11) -(3.12), since the solution which is exhibited in this theorem can depend on θ even if the data do not. However, the ellipticity property

∀v ∈ L 2 1 (ω) 2 , ã[ξ] (v, v) ≥ α 0 v 2 L 2 1 (ω) 2 , (3.14)
and the inf-sup condition

∀q ∈ H 1 1w (ω), sup v∈L 2 1 (Ω) 2 b(v, q) v L 2 1 (ω) 2 ≥ |q| H 1 1 (ω) , (3.15)
can be derived by the same arguments as for (2.8) and (2.9). To make this last condition complete, we now prove a weighted Poincaré-Friedrichs condition which ensures the equivalence between the norms

| • | H 1 1 (ω) and • H 1 1 (ω) on H 1 1w (ω).
Lemma 3.1. The following inequality holds for all functions q in H 1 1w (ω)

q L 2 1 (ω) ≤ r 2 1 2 log( r 1 r 0 ) + r 2 0 -r 2 1 4 1 2 |q| H 1 1 (ω) . ( 3

.16)

Proof: Any function q in H 1 1w (ω) satisfies for all (r, z) in ω q(r, z) = r r 0 (∂ r q)(ρ, z) dρ, so that, thanks to a Cauchy-Schwarz inequality

|q(r, z)| ≤ r r 0 (∂ r q) 2 (ρ, z) ρ dρ 1 2 r r 0 dρ ρ 1 2 ≤ (log r r 0 ) 1 2 ( r 1 r 0 (∂ r q) 2 (ρ, z) ρ dρ) 1 2 .
Integrating the square of this inequality on ω with respect to the measure r drdz yields

q L 2 1 (ω) ≤ r 1 r 0 (log r r 0 ) r dr 1 2 ∂ r q L 2 1 (ω) . So the desired result follows from r 1 r 0 (log r r 0 ) r dr = r 2 1 2 log( r 1 r 0 ) - r 2 1 4 + r 2 0 4 .
We skip the proof of the next theorem since it relies on exactly the same arguments as for Theorem 2.3.

Theorem 3.2. For any data (f , p w , g) in L 2 1 (ω) 2 ×H 1 2 (γ w )×L 2 1 (γ), problem (3.11)-(3.12) admits a solution (u, p) in L 2 1 (ω) 2 × H 1 1 (ω). Moreover this solution satisfies u L 2 1 (ω) 2 + p H 1 1 (ω) ≤ c f L 2 1 (ω) d + p w H 1 2 (γ w ) + g L 2 1 (γ) , (3.17)
where the constant c depends on r 0 and r 1 .

Remark 3.3. Note that Theorem 3.2 would not hold if r 0 were equal to zero, since functions in H 1 1 (ω) have no traces on γ w in this case.

The discrete problem and its numerical analysis.

As standard in spectral methods, the discrete problem is constructed from (3.11) -(3.12) by the Galerkin method with numerical integration. For the choice of the discrete spaces, following the approach used in [3, §3.a] for Darcy's equations with constant coefficient, we fix an integer N and define

X N = P N (ω) 2 , M N = P N (ω), M w N = M N ∩ H 1 1w (ω), (4.1)
where, for each nonnegative integer n, P n (ω) stands for the space of restrictions to ω of polynomials with two variables r and z and degree ≤ n with respect to each of them.

On the other hand, in order to handle the non constant coefficient α and as suggested in [START_REF] Maday | Rønquist -Optimal error analysis of spectral methods with emphasis on non-constant coefficients and deformed geometries[END_REF], we have decided to use over-integration. We fix an integer M > N and recall that there exist a unique set of M + 1 nodes ξ j , 0 ≤ j ≤ M , with ξ 0 = -1 and ξ M = 1, and a unique set of M + 1 weights ρ j , 0 ≤ j ≤ M , such that the following equality holds

∀Φ ∈ P 2M -1 (-1, 1), 1 -1 Φ(ζ) dζ = M j=0 Φ(ξ j ) ρ j , (4.2)
with obvious notation for the polynomial spaces P n (-1, 1). Moreover, the ρ j are positive.

To go further, we introduce the discrete product. Let F denote the simplest mapping that sends ] -1, 1[ 2 onto ω: F is the product of a homothety and a translation in each direction. Thus, we set, for all functions u and v continuous on ω,

(u, v) M = (r 1 -r 0 )|z 1 | 4 M i=0 M j=0 u • F (ξ i , ξ j )v • F (ξ i , ξ j ) ρ i ρ j . (4.
3)

It follows from the choice of F that this product coincides with the scalar product of L 2 (ω) on all functions u and v such that uv belongs to P 2M -1 (ω). Assuming that F maps the edge {-1}×] -1, 1[ onto γ w , we also define a discrete product on γ:

(u, v) γ M = (r 1 -r 0 ) 2 M j=0 u • F (ξ j , -1)v • F (ξ j , -1) ρ j + |z 1 | 2 M j=0 u • F (1, ξ j )v • F (1, ξ j ) ρ j + (r 1 -r 0 ) 2 M j=0 u • F (ξ j , 1)v • F (ξ j , 1) ρ j . (4.4)
We now denote by ξ * j , 0 ≤ j ≤ N , the analogues of the nodes ξ j but now for M = N . Let i w N stand for the Lagrange interpolation operator at the nodes F (-1, ξ * j ), 0 ≤ j ≤ N , with values in P N (γ w ). Assuming that the data f , p w and g are continuous, the discrete problem reads

Find (u N , p N ) in X N × M N such that p N = i w N p w on γ w , (4.5)
and ∀v N ∈ X N , a

[p N ] M (u N , v N ) + b M (v N , p N ) = (f , v N r) M , ∀q N ∈ M w N , b M (u N , q N ) = (g, q N r) γ M , (4.6)
where the bilinear forms a

M (u N , v N ) = α(ξ)u N r , v N r r M + α(ξ)u N z , v N z r M , b M (v N , q N ) = v N r , ∂ r q N r M + v N z , ∂ z q N r M . (4.7)
It can be noted that the form a

[ξ]

M (•, •) only involves the values of α(ξ) at the nodes F (ξ i , ξ j ), so that computing the nonlinear term is not expensive. Moreover, the boundedness of α implies the continuity of the form a

[ξ] M (•, •).
Similar arguments as in the previous sections would yield the existence of a solution for problem (4.5) -(4.6). However, relying on the ideas in [START_REF] Brezzi | Raviart -Finite dimensional approximation of nonlinear problems, Part I: Branches of nonsingular solutions[END_REF], we have rather state and prove a more precise result. We need some further notation for that. First, we introduce the space

X (ω) = L 2 1 (ω) 2 × H 1 1 (ω). (4.8)
We set: α = α 0 +α 1 2 and introduce the operator T which associates with any data (f , p w , g)

in L 2 1 (ω) 2 × H 1 2 (γ w ) × L 2 1 (γ), the solution U * = (u * , p * ) in X (ω) of the problem p * = p w on γ w , (4.9
)

and ∀v ∈ L 2 1 (ω) 2 , a(u * , v) + b(v, p * ) = Ω f (r, z) • v(r, z) r drdz, ∀q ∈ H 1 1w (ω), b(u * , q) = γ g(τ )q(τ ) r(τ ) dτ, (4.10)
where the bilinear form a(•, •) is defined by

a(u, v) = α Ω u r (r, z)v r (r, z) + u z (r, z)v z (r, z) r drdz. ( 4.11)
The continuity and ellipticity of this new form being obvious, it follows from (3.15) that this operator is well-defined and also that the following stability property holds

T (f , p w , g) X (ω) ≤ c f L 2 1 (ω) 2 + p w H 1 2 (γ w ) + g L 2 1 (γ) . (4.12)
Thus, problem (3.11) -(3.12) can equivalently be written

U -T G(U ) = 0, with G(U ) = α -α(p) u + f , p w , g . (4.13)
Similarly, let X N (ω) stand for the space

X N (ω) = X N × M N . (4.14)
The discrete Darcy operator T N ssociates with any smooth enough data (f , p w , g) the solution

U * N = (u * N , p * N ) in X N (ω) of the problem p * N = i w N p w on γ w , (4.15)
and

∀v N ∈ X N , a M (u * N , v N ) + b M (v N , p * N ) = Ω f (r, z) • v N (r, z) r drdz, ∀q N ∈ M w N , b M (u * N , q N ) = γ g(τ )q N (τ ) r(τ ) dτ, (4.16)
where the bilinear form a M (•, •) is defined by

a M (u N , v N ) = α u N r , v N r r M + α u N z , v N z r M . (4.17)
The fact that this operator is well-defined is established in the next lemma. Finally, we observe that problem (4.5) -(4.6) can equivalently be written

U N -T N G N (U N ) = 0, (4.18)
where the three components G N 1 , G N 2 and G N 3 of the mapping G N are defined with obvious notation by

G N 1 (U N ), V N = a M (u N , v N ) -a [p M ] M (u N , v N ) + (f , v N r) M , G N 2 (U N ) = p w , G N 3 (U N ), V N = (g, q N r) γ M . (4.19)
Since M > N , the forms a(•, •) and a M (•, •) coincide on X N (ω) × X N (ω), so that the continuity and the ellipticity of the form a M (•, •) (with constants independent of M ) are obvious. The forms b(•, •) and b M (•, •) also coincide on X N × M N , whence the continuity of the form b M (•, •). Taking v N equal to grad q N leads to the following inf-sup condition

∀q ∈ M w N , sup v N ∈X N b M (v N , q N ) v N L 2 1 (ω) 2 ≥ |q N | H 1 1 (ω) , (4.20)
and, since M w N is imbedded in H 1 1w (ω), applying Lemma 3.1 makes this condition complete.

Lemma 4.1. The operator

T N is continuous from L 2 1 (ω) 2 × H τ (γ w ) × L 2 1 (γ) into X (ω) for all τ > 1 2
. Moreover, the following stability property holds

T N (f , p w , g) X (ω) ≤ c f L 2 1 (ω) 2 + i w N p w H 1 2 (γ w ) + g L 2 1 (γ) , (4.21)
for a constant c independent of N .

see [START_REF] Bernardi | Spectral Methods[END_REF]Thm 7.1] or [5, Prop. V.2.1]. Similarly, we denote by I N the Lagrange interpolation operator at all nodes F (ξ * i , ξ * j ) with values in M N , where the ξ * j , 0 ≤ j ≤ N , are the analogues of the nodes ξ j but now for M = N (see the definition of i w N ). Setting pN = I N p * , we recall from [START_REF] Bernardi | Spectral Methods[END_REF]Thm 14.2] that (note that this requires s > 1 2 )

p * -pN H 1 1 (ω) ≤ c N -s p * H s+1 1 (ω) . (4.28)
Since pN is equal to i w N p w on γ w , the pair (u * N -ũN , p * N -pN ) belongs to X N × M w N and satisfies

∀v N ∈ X N , a M (u * N -ũN , v N ) + b M (v N , p * N -pN ) = a(u * -ũN , v N ) + b(v N , p * -pN ), ∀q N ∈ M w N , b M (u * N -ũN , q N ) = b(u * -ũN , q N ).
Therefore, standard arguments give

u * -u * N L 2 1 (ω) 2 + p * -p * N H 1 1 (ω) ≤ c u * -ũN L 2 1 (ω) 2 + p * -pN H 1 1 (ω) ,
so that the desired estimate follows from (4.27) and (4.28).

Estimate (4.26) requires some regularity of T (f , p w , g) but, when combined with (4.21), it is sufficient to prove the following convergence result: For any triple (f , p w , g) in We also need a modified version of Lemma 4.3, where less regularity on T (f , p w , g) is required.

L 2 1 (ω) 2 × H τ (γ w ) × L 2 1 (γ), τ > 1 2 , lim N →+∞ (T -T N )(f , p w , g) X (ω) = 0. ( 4
Corollary 4.4. There exists a constant c such that the following estimate holds for any (f , p w , g) such that p w belongs to H τ + 1 2 (γ w ), τ > 1 2 , and T (f , p w , g) belongs to

H s 1 (ω) 2 × H s+1 1 (ω), s > 0, (T -T N )(f , p w , g) X (ω) ≤ c N -s T (f , p w , g) H s 1 (ω) 2 ×H s+1 1 (ω) + N -τ p w H τ + 1 2 (γ w ) . ( 4

.31)

Proof: Let p w be a lifting of p w in H τ +1 (ω). We set: p 0 = p *p w and take, with the same notation as above, p * N = I N p w + p 0 N , where p 0 N is the orthogonal projection of p 0 from H 1 1w (ω) onto M w N . Thus, standard arguments [6, Thm 7.2] yield that

p * -p * N H 1 1 (ω) ≤ c N -s p 0 H s+1 1 (ω) + N -τ p w H τ + 1 2 (γ w ) . itself.
2) It is proved in [START_REF] Bernardi | Blouza -Spectral discretization of a Naghdi shell model[END_REF]Lemma 3.1] for instance that there exists α ⋄ N in P N ⋄ (ω) such that

α(p ⋄ N ) -α ⋄ N L ∞ (ω) ≤ c N -s (log N ) 1 2 α(p) H s+1 1 (ω) .
This gives the first part of (4.36). Proving the second part relies on exactly the same argument.

Finally, we note that only the first components of DG(U ) and of DG N (U ⋄ N ) are not zero and that they are given by, for any ) is an isomorphism of X N (ω) and the norm of its inverse is bounded independently of N .

W N = (w N , ρ N ) in X N (ω), DG 1 (U ).W N , V = a(w N , v) -a [p] (w N , v) - Ω α ′ p(r, z) ρ N (r, z) u r (r, z)v r (r, z) + u z (r, z)v z (r, z) r drdz, DG 1N (U ⋄ N ).W N , V N = a M (w N , v N ) -a [p ⋄ N ] M (w N , v N ) -α ′ (p ⋄ N )ρ N u ⋄ N r , v N r r M -α ′ (p ⋄ N )ρ N u ⋄ N z , v N z r M . ( 4

Proof: We write

Id -T N DG N (U ⋄ N ) = Id -T DG(U ) + T -T N DG(U ) + T N DG(U ) -DG(U ⋄ N ) + T N DG(U ⋄ N ) -DG N (U ⋄ N ) .
Owing to Assumption 4.5, the desired result is obtained if the last three terms in this equality tend to zero when N tends to +∞. We now check this property, successively for the three terms. Let L X N (ω) denote the space of endormorphisms of X N (ω). 1) Let W N run through the unit sphere of X N (ω). Thus, αα(p) w N belongs to a finite-dimensional space and is bounded in L 2 1 (ω) 2 , hence belongs to a compact subset of L 2 1 (ω) 2 . On the other hand, setting 1 r = 1-s 2 (so that H s 1 (ω) is imbedded in L r (ω)) and

1 r + 1 r * = 1 2
, we use the compactness of the imbedding of H 1 1 (ω) into L r * (ω) to derive that the quantities α ′ (p)ρ N u N r and α ′ (p)ρ N u N z also belong to a compact set of L 2 1 (ω). Thus, it follows from (4.37) and (4.30) that lim

N →+∞ T -T N DG(U ) L(X N (ω)) = 0. (4.38)
2) By using the imbedding of H 1 1 (ω) into L r * (ω) for the r * introduced above and also into L 2r * (ω), we observe from (4.37) that

DG 1 (U ) -DG 1 (U ⋄ N ) .W N , V ≤ c p -p ⋄ N L ∞ (ω) + c ′ p -p ⋄ N L ∞ (ω) u H s 1 (ω) 2 + c u -u ⋄ N H s 2 1 (ω) 2 W N X (ω) v L 2 1 (ω) 2 ,
We also agree to denote by (u ℓ N (λ), p ℓ N (λ)) the solution of problem (5.3) -(5.4) with the function α replaced by α λ . Next, we choose a sample of values of λ in [0, 1], for instance the λ k = k m , 0 ≤ k ≤ m, for a fixed integer m. We also fix a positive integer L, and • for k = 0, since α 0 (ξ) is equal to α, we solve problem (4.5) -(4.6) which in this case is simply the linear Darcy system with constant permeability α;

• for 1 ≤ k ≤ m, the initial guess (u 0 N (λ k ), p 0 N (λ k )) being equal to (u L N (λ k-1 ), p L N (λ k-1
)), we iteratively solve problem (5.3) -(5.4) for 1 ≤ ℓ ≤ L. It is readily checked that, for L large enough, we obtain an initial guess (u 0 N (1), p 0 N (1)) whih satisfies (5.5). However, in practical stuations, L is chosen equal to 2 or 3.

From now on, we work in the physical case of the function α given by

α(ξ) = α * exp(ξ), (5.8)
for a constant α * > 0, see [18, §3.5] for the justification of this choice. Note that this function does not satisfy assumption (2.1). However a direct consequence of Proposition 2.5 is that the function p is bounded as a function of the data. So the previous numerical analysis can be performed by truncating the function α at appropriate values, without modifying the discrete problem (4.5) -(4.6).

In a first step, we test the convergence of the discretization on a model problem. We take α * equal to 10 2 and we work with the domain Ω defined in (1.1), with r 0 = 0.6 m, r 1 = 3 m, z 1 = -6 m.

(5.9) (5.11) Newton's algorithm is iterated until the error (u ℓ N , p ℓ N ) -(u, p) X (ω) becomes smaller than 10 -8 or at most five times. In a second step, we still work with the domain Ω defined from (5.9), but now with α * equal to 1. The datum g is given by g(r 1 , z) = g(r, z 1 ) = 0, g(r, 0) = 1.

(5.12) 4 present from left to right the isovalues of the two components of the velocity and of the pressure, obtained with N = 80 and M = 81, for g given in (5.12) and first for p w (r 0 , z) equal to 1 (Figure 3), second for p w given by p w (r 0 , z) = z 6

4

(5.13) (Figure 4). These curves are in good coherence with the physics of the problem and prove the accuracy of spectral discretizations.

  Spectral discretization of Darcy's equations with pressure dependent porosity by Mejdi Azaïez 1 , Faker Ben Belgacem 2 , Christine Bernardi 3 , and Nejmeddine Chorfi 4

  Γ) and ρ small enough, there exists a harmonic function p b in W 1,ρ (Ω) such that p b = p w on Γ w and ∂ n p b = g on Γ.

M

  (•, •) for any continuous function ξ on ω and b M (•, •) are defined by a [ξ]

  .29) This yields the further result which is useful in what follows: For any compact K of L 2 1 (ω) 2 , lim N →+∞ sup f ∈K (T -T N )(f , 0, 0) X (ω) = 0. (4.30)

  .37) We are now in a position to prove the next lemma.

Lemma 4 . 8 .

 48 If Assumptions 4.5 and 4.6 hold, there exists an integer N ♯ such that, for all N ≥ N ♯ , the operator Id -T N DG N (U ⋄ N

Figure 1 : 6 ) 4 6 ) 4

 16464 Figure 1: Convergence curves for the solution in (5.10) for µ = 1.5 and µ = 2.5

Figure 1

 1 Figure 1 presents the curves of the logarithms of the errors

Figure 2

 2 Figure 2 presents from left to right the isovalues of the radial and axial components of the velocity and of the pressure, obtained for µ = 2.5, N = 80 and M = 81. These curves are exactly the same as for the exact solution in (5.10).

Figure 2 :

 2 Figure 2: The discrete solution computed from (5.11) for µ = 2.5

Figure 3 :

 3 Figure3: The discrete solution with boundary conditions (5.12) and p w constant Figures3 and 4present from left to right the isovalues of the two components of the velocity and of the pressure, obtained with N = 80 and M = 81, for g given in (5.12) and first for p w (r 0 , z) equal to 1 (Figure3), second for p w given by

Figure 4 :

 4 Figure 4:The discrete solution with boundary conditions (5.12) and(5.13)

  Proposition 2.1. For any data (f , p w , g) in L 2 (Ω) d × H -(2.4) are equivalent, in the sense that any pair (u, p) in L 2 (Ω) d × H 1 (Ω) is a solution of system (1.3) in the distribution sense if and only if it is a solution of problem (2.3) -(2.4).

	and (2.3)	1 2 (Γ w ) × H	1 2 00 (Γ) ′ , problems (1.3)

For any measurable function ξ and for any data (f , p w

	| • | H 1 (Ω) and w (Ω) is a direct consequence of the Poincaré-Friedrichs inequality norm • H 1 (Ω) on H 1 (the equivalence constants only depending on the geometry of Ω and Γ w ). So, standard
	arguments [12, Chap. I, §4.1] lead to the following result.
	Lemma 2.2.

  × H 1 (Ω) and such that (p n ) n converges strongly in L 2 (Ω). It follows from the properties of the function α that, for any funcion v m in V m , since (p n ) n converges to p ε a.e. in Ω, (α(p n ) v m ) n also converges a.e. to α(p) v m and is upper bounded by α 1 |v m | which belongs to L 2 (Ω). So, using the Lebesgue dominated convergence theorem yields that α(p n

there exists a subsequence, still denoted by (u n , p n ) n for simplicity, which converges to a (u ε , p ε ) weakly in L 2 (Ω) d

  (Ω) d onto itself with norm 1, a simple interpolation argument, see[START_REF]Adams -Sobolev Spaces[END_REF] Thms 7.17 & 7.20], yields that it is continuous from L ρ (Ω) d onto itself for 2 ≤ ρ ≤ ρ * , with norm χ θ(ρ) where θ is a continuous increasing function from 0 to 1 onto[2, ρ

1

,ρ * (Ω). This in turn implies that D * is continuous from L ρ * (Ω) d onto itself; let χ denote its norm. Since it is continuous from

L 2 * ].

00 (Γ), see [14, Chap. 1, §11].

,(5.10)

Acknowledgements. We are deeply grateful towards K.R. Rajagopal for suggesting us to work on this model and D. Smets for finding the perfect reference in the right time. We

Proof: Since τ > 1 2 , p w is continuous on γ w , so that i w N p w is well-defined. Moreover, it follows from [START_REF] Bernardi | Rapetti -Discrétisations variationnelles de problèmes aux limites elliptiques, Collection[END_REF]Th. III.3.1] for instance that there exists a p b N in M N equal to i w N p w on γ w and such that (we use here the imbedding of H 1 (ω) into H 1 1 (ω))

On the other hand, combining the inf-sup condition (4.20) and Lemma 3.1 yields the existence of a function u b N such that

and that

Finally, it follows from the ellipticity of a M (•, •) and the inf-sup condition (4.20) that the problem: Find (u 0 N , p 0 N ) in X N × M w N such that

has a unique solution. The pair (u

) is then the unique solution of problem (4.15) -(4. [START_REF] Meyers | An L p -estimate for the gradient of solutions of second order elliptic divergence equations[END_REF], and the operator T N is well-defined. Moreover, taking v N equal to u 0 N in (4.24) yields

We need this last property in what follows.

Lemma 4.3. There exists a constant c such that the following estimate holds for any (f , p w , g) such that T (f , p w , g) belongs to We now consider a solution (u, p) of problem (3.11) -(3.12) which we wish to approximate. We are led to make the following assumption. Here, D stands for the differential operator. (ω) for a real number s > 1 4 , (ii) is such that Id -T DG(U ) is an isomorphism of X (ω).

Part (i) of this assumption seems likely but is not proved. Part (ii) is much less restrictive than the global uniqueness of the solution, see Proposition 2.4, since it only requires its local uniqueness.

From now on, we assume that there exists a real number µ, 0 < µ ≤ 1, such that M is equal to the integer part of (1 + µ)N . Thus, denoting by N ⋄ the integer part of µ N -1, standard arguments (see [START_REF] Bernardi | Spectral Methods[END_REF]Thm 7.4]) yield the existence

and

Some further approximation results which are stated in the next lemma require a new assumption (which is not at all restrictive since we work with a very smooth function α).

As usual, ⌈s⌉ denotes the smallest integer ≥ s.

Assumption 4.6. For the same real number s as in Assumption 4.5, the function α is of class C ⌈s⌉+2 on R and belongs to H ⌈s⌉+2 (R).

where the constant c(p) only depends on p H s+1 1 (ω) .

Proof: We proceed in two steps. 1) Since H s+1 (ω) is imbedded in W 1,2(s+1) (ω) for all s ≥ 0, using Assumption 4.6 and applying [8, Thm 1] yield that the function:

where c and c ′ stand for the norms of α in W 1,∞ (R) and W 2,∞ (R). Then, it follows from (4.21), the imbedding of H

(ω) into L ∞ (ω), (4.32) and (4.34) that lim

3) It follows from the exactness property (4.2) of the quadrature formula that

and also that

Adding and subtracting these equations and using severy Cauchy-Schwarz inequalities, we obtain

Combining this with (4.33), Lemma 4.7 and finally (4.25) yields lim

The lemma is then a direct consequence of (4.38) to (4.40).

Lemma 4.9. If Assumptions 4.5 and 4.6 hold, there exists a constant c depending on u such that the following property holds for any Z N in X N (ω)

We now evaluate sucessively the three lines of the right-hand side that we denote by A 1 , A 2 and A 3 for brevity.

1) The same arguments as in the proof of Lemma 4.8 give

To go further, we recall from standard inverse inequalities [START_REF]Quarteroni -Some results of Bernstein and Jackson type for polynomial approximation in L p -spaces[END_REF] combined with the fact that the norm of the imbedding of H 1 (ω) into L q (ω) behaves like √ q for large values of q, see [START_REF]Talenti -Best constant in Sobolev inequality[END_REF], that

so that taking q equal to log N gives

Combining all this yields

2) Similar arguments give

so that taking 1 r = 1-s 2 and 1 r + 1 r * = 1 2 and using appropriate Sobolev imbeddings lead to

Thus, it follows from (4.33) and (4.42) that

3) To bound the last term, we introduce the Lagrange interpolation operator I M at the nodes F (ξ i , ξ j ) with values in P M (ω). Indeed, it follows from the exactness of the quadrature formula and the boundedness of α ′ that

We recall from [6, Remark 13.5] that, for any positive integer k,

Since ρ N (u ⋄ Nz N ) belongs to P 2N (ω), applying this inequality with k = 2 for instance yields

Thus, using (4.42) once more yields

The lemma is now a consequence of (4.43) to (4.45), combined with (4.25).

Lemma 4.10. If Assumptions 4.5 and 4.6 hold and if the data (f , p w , g) belong to

, there exists a constant c depending on (u, p) such that the following property holds

(4.46)

Proof: Owing to (4.13), we have

The first term is bounded in (4.32) and (4.34), and evaluating the second one follows from Corollary 4.4. The same arguments as in the previous proofs yield that

so that the bound for the third term is derived from (4.32) and (4.34), combined with (4.25). Finally, we note that, for any polynomial αN in P N (ω), we have

where I M denotes the interpolation operator at the nodes F (ξ i , ξ j ) and i M the interpolation operator at the nodes which are contained in γ (we do not make its definition precise for simplicity), So, adding and subtracting these equations, assuming for instance that s < 1 and using standard Sobolev imbeddings and finally an extension of (4.25) give

Owing to (4.33), the final estimate follows from standard approximation and interpolation properties, see [START_REF] Bernardi | Spectral Methods[END_REF]Thms 7.3,13.4 & 14.2] (this requires the assumption s > 1 4 ).

Thanks to Lemmas 4.8 to 4.10, we now are in a position to apply the Brezzi-Rappaz-Raviart theorem [START_REF] Brezzi | Raviart -Finite dimensional approximation of nonlinear problems, Part I: Branches of nonsingular solutions[END_REF] (see also [START_REF] Crouzeix | Rappaz -On Numerical Approximation in Bifurcation Theory, Collection[END_REF]Thm 3.1] for a different version).

, there exists an integer N * and a constant c > 0 such that, for all N ≥ N * , problem (4.5) -(4.6) has a unique solution (u N , p N ) which satisfies

Moreover, the following error estimate holds between these solutions (u, p) and (u N , p N )

(4.48) Estimate (4.48) is fully optimal and the assumptions which are required on the solution (u, p) seem likely.

5. An algorithm and some numerical experiments.

Applying Newton's method to problem (4.5) -(4.6) consists in solving iteratively the equation

This of course requires that Id -T N DG N (U ℓ-1 N ) is an isomorphism of X N (ω), which seems likely when U ℓ-1 N is sufficiently close to a solution U of problem (4.13) satisfying Assumption 4.5, see the proof of Lemma 4.8. Multiplying both sides of equation (5.1) by Id -T N DG N (U ℓ-1 N ) we observe that this equation can equivalently be written as follows: Being given an initial guess (u 0

we solve the following problem, for ℓ ≥ 1,

It is readily checked that, for each value of ℓ, problem (5.3) results into a square linear system. The convergence of this method can easily be derived from [START_REF] Brezzi | Raviart -Finite dimensional approximation of nonlinear problems, Part I: Branches of nonsingular solutions[END_REF] (see also [START_REF] Girault | Raviart -Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms[END_REF]Chap. IV, Thm 6.3]) owing to Lemma 4.9. for a constant c independent of N .

As standard for Newton's method, the key point is to exhibit an initial guess (u 0 N , p 0 N ) satisfying (5.5). In order to do that, we have decided to use a continuation method. For the constant α = α 0 +α 1 2 , we define ∀λ ∈ [0, 1], α λ (ξ) = (1λ) α + λ α(ξ).

(5.7)