Numerical coupling of petroleum wellbore and reservoir models with heat transfer

c TU Delft, The Netherlands, 2006

NUMERICAL COUPLING OF PETROLEUM WELLBORE

AND RESERVOIR MODELS WITH HEAT TRANSFER

Mohamed Amara∗_{, Daniela Capatina}∗ _{and Layal Lizaik}∗,† ∗_{Laboratory of Applied Mathematics, University of Pau,}

IPRA BP 1155, 64013 Pau Cedex, France

e-mail: mohamed.amara@univ-pau.fr, daniela.capatina@univ-pau.fr †_{TOTAL, CST Jean Feger, Av. Larribau,}

64018 Pau Cedex, France e-mail: layal.lizaik@total.com

Key words: Fluid and Porous Media, Mixed Finite Elements, Thermometrics, Coupling Abstract. The talk is devoted to the modelling and the nite element approximation of the ow of a monophasic compressible uid in a petroleum reservoir and wellbore, from both a dynamic and a thermal point of view. To do that, we couple a reservoir model, including viscous dissipation and compressibility eect, with a pseudo 1D wellbore model, taking into account the particular geometry of the well and the privileged direction of the ow. Each model is validated from a numerical and a physical point of view.

1 INTRODUCTION

Due to emerging technologies such as optical ber sensors, temperature measurements are destined to play a major role in petroleum production logging interpretation. Using temperature recordings from a wellbore and a owrate history on the surface, it can be envisaged to develop new ways to predict ow repartition among each producing layer of a reservoir or to estimate virgin reservoir temperatures.

In order to solve the inverse problem, one rst needs to develop a forward model describing the ow of a monophasic compressible uid (oil or gas) in a petroleum reservoir and a well, from both a dynamic and a thermal point of view. This implies to couple a reservoir model (porous medium) and a well model (uid medium).

In this talk, we rst present the reservoir, which was thoroughly studied in 1 _{and we}

next describe a pseudo 1D wellbore model, whose rst version was introduced in 2_{. Both}

problems are written in axisymmetric form and depend on the cylindrical coordinates (r, z).

elements. Its unknowns are the specic ux, the heat ux, the pressure and the temper-ature.

The wellbore model is based on the compressible Navier-Stokes equations coupled with an energy equation. In order to take into account the privileged direction of the ow and to reduce the cost of the calculation, a pseudo 1D model is derived by constructing an explicit solution in terms of the radial coordinate r. The nonlinear time-discretized problem is then solved by means of a xed point method with respect to the density. Both proposed models were separately validated from a numerical and a physical point of view.

For the sake of clarity, we shall denote by Ω1 the domain occupied by the porous

medium, by Ω2 the domain of the uid and by Σ the interface. The boundaries of Ω1, Ω2

will be respectively denoted by Υ, respectively Γ.

We then focus on the coupling of the reservoir and the well models mentioned above. In order to take into account recorded owrates at the surface, a global resolution is envisaged. This is achieved by imposing transmission conditions at the perforations at each time-step, binding thus together the two models. This approach leads to a non classical mixed formulation which is next analyzed. Numerical tests including real study cases will be presented during the talk.

2 RESERVOIR MODEL 2.1 Physical modelling

The studied domain is a cylindrical vertical petroleum well, delimited by a casing and surrounded by a cement layer and by the reservoir. The communication between the well and the reservoir is achieved through perforations. The reservoir Ω1 is treated as a

porous medium divided into several geological layers which are characterized by their own dip and physical properties. Thus, each layer is made of a porous rock (of porosity φ), characterized by vertical and horizontal permeabilities, and saturated with both a mobile monophasic uid and a residual formation water.

The reservoir model is described by the Darcy-Forchheimer equation together with an energy equation, which notably includes the temperature eects due to the decompression of the uid (Joule-Thomson eect) and the frictional heating that occurs in the formation, cf. for instance 3_{. This energy balance thus quanties the cooling or the heating of the}

produced uid before entering the wellbore. Therefore, it stands apart from the classical models which mainly consider the wellbore thermal exchanges due to conduction and convection and assume that the produced uid enters the wellbore at the geothermal temperature.

           rφ∂ρ_∂t + div(rG) = 0 ρ−1(µK−1+ F |G| I)G + ∇p = ρg q − λ∇T = 0 r(ρc)∗∂T_∂t + rρ−1(ρc)fG · ∇T − div(rq) − rφβT∂p_∂t − rρ−1(βT − 1)G · ∇p = 0 ρ = ρ(p, T ) (1)

where G = ρv represents the specic ux and q = λ∇T the heat ux, with λ the equivalent thermal conductivity.

As usually when modelling petroleum uids, we have used the Peng-Robinson cubic state equation4 to close the system. Moreover, due to the high ltration velocity which

can arise around gas wells, we have introduced a quadratic term in the standard Darcy's equation, to take into account the kinematic energy losses; F stands for the Forchheimer's coecient. In the energy equation, (ρc)∗ characterizes the heat capacity of a virtual

medium, equivalent to the uid and the porous matrix, while (ρc)f only symbolizes the

uid properties. K stands for the permeability tensor, while β is the expansion coecient. To this system we add initial conditions for ρ and T and boundary conditions. On a given part of the boundary, we impose either the normal ux G · n or the pressure; it goes the same way for the normal heat ux and its dual variable, the temperature. We are thus able to treat teh standard cases of a closed reservoir and of a reservoir feed at constant pressure. Concerning the temperature, the geothermal gradient is imposed on the top and on the bottom of the reservoir.

2.2 Analysis and numerical approximation

A complete mathematical analysis of the reservoir model, together with many numerical tests, can be found in 1_{. We only recall here the main steps.}

In order to deal with the nonlinear system, we rst implement a time discretization based on Euler's implicit scheme and then we linearize the problem at each time step. Its unknowns are now the pressure p, the specic ux G, the temperature T and the heat ux q, the density being updated by verifying the state equation at the end of each time loop. This last step is achieved by means of a thermodynamic module.

The previous semi-discretized system was shown to be well-posed by means of the Fredholm's alternative, for smooth coecients and for a suciently small ∆t .

Concerning the space discretization, we introduce a mixed variational formulation where the convective terms are treated by means of an upwind scheme. To provide an accurate determination of the uxes, the lowest-order Raviart-Thomas nite elements are employed, whereas the pressure and the temperature are approximated by piecewise constant elements. The matrix of the discrete problem can then be written as follows :

A1 =  A1 B1 −BT 1 C1  ,

where A1 is symmetric and positive, C1 is positive but not symmetric and B1satises

an inf-sup condition, see for instance 5 _{for details. The mixed problem is shown to be}

well-posed thanks to an extension of the Babu²ka-Brezzi theory, cf. 6_.

Finally, the code has been validated from both a numerical point of view (behaviour of the solution with respect to mesh renement) and a physical one (comparison of the computed pressure with analytical pressure solutions given by well-test softwares such as PIE). A realistic example has also been treated, for an existing reservoir with high heterogeneities concerning the physical properties, and a comparison with recorded data was carried on (see 1_).

3 CONSERVATIVE PSEUDO 1D WELLBORE MODEL

Let us proceed with the petroleum wellbore model, described by the compressible Navier-Stokes equations together with an energy balance and the Peng-Robinson cubic state equation. The domain is denoted by

Ω2 = {(r, z) ; 0 ≤ r ≤ R, z ∈ I}

where I = [zmin, zmax] and where R is the radius of the well. In practice, R ' 4inch while

the length of the pipe can attend several thousands meters. Our problem is written in 2D axisymmetric form as follows :

∂

∂t(rρ) + ∇ · (rρu) = 0 (3)

 _∂

∂t(rρur) + ∇ · (rurρu) + r ∂p ∂r − ∂ ∂r(rτrr) − ∂ ∂z(rτzr) + τθθ+ rκρ|u|ur= 0 ∂

∂t(rρuz) + ∇ · (ruzρu) + r ∂p ∂z − ∂ ∂r(rτrz) − ∂ ∂z(rτzz) + rρg + rκρ|u|uz = 0 (4) ∂

∂t(rρE) + ∇ · (r(ρE + p)u) − ∇ · (rτ u) − ∇ · (rλ∇T ) + rρguz = 0 (5)

where E = cvT + |u|2

2 is the total energy, cv the specic heat, µ the viscosity, λ the thermal

conductivity and κ a positive coecient depending on the pipe's diameter. We denote ∇ = (∂

∂r, ∂ ∂z)

t _{, u = (u}

r, uz) and we recall that :

τrr = 2µ ∂ur ∂r − 2 3µ  1 r ∂ ∂r(rur) + ∂uz ∂z  , τrz = τzr = µ  ∂u_z ∂r + ∂ur ∂z  , τzz = 2µ ∂uz ∂z − 2 3µ  1 r ∂ ∂r(rur) + ∂uz ∂z  , τθθ = 2µ ur r − 2 3µ  1 r ∂ ∂r(rur) + ∂uz ∂z  . One has to add initial conditions on ρ, ρu, ρE as well as boundary conditions which will be prescribed later.

Next, in order to take into account the ow's privileged direction, the particular ge-ometry of the domain as well as the supply at the perforations, we propose a pseudo 1D modelling of the above problem (see also 2 _{for a slightly dierent approach). Thus, on}

the one hand, the calculations are lightened and one can treat pipes of several thousands meters high and on the other hand, one avoids any numerical instability due to the large aspect ratio of any 2D grid.

For this purpose, we introduce two conservative variables (the specic ux G = ρu and the heat ux q = λ∇T ) and a time discretization which leads, for a given ρ, to the following decoupled system :

div(rG) = −rρ − ρ n ∆t , (7)    div(ru) = 1 ρ(div(rG) − G · ∇ρ)

rρ_∆tu + rG · ∇u + r∇p − div(rτ ) + τθθer+ rκ|G|u = rρg + rρu

n ∆t (8)              rcv ρ_∆tT + G · ∇T
− div(rq) = rρcvT n ∆t − 1 2r 

ρ|u|2−|u_∆tn|2 + G · ∇(|u|2₎_{− div(rpu) + div(rτ u) + rg · G}

q = λ∇T.

(9)

The nonlinear system is solved by applying a xed point method with respect to ρ : more precisely, we compute G from a Petrov-Galerkin formulation of (7), we next recover (u, p) from the quasi-Stokes system (8) and nally, we compute (q, T ) from a mixed formulation of (9) and update the density by means of a thermodynamic module. We then loop until the convergence is achieved.

Let us remark that the rst equation of (8) translates the fact that div(ru) = div(r ρG),

while in the other equations we have substituted ρu by G. Therefore, solving (8) and (9) is not equivalent to solving the initial Navier-Stokes and energy equations.

Next, in order to specify the boundary conditions associated to the systems (7), (8) and (9), ∂Ω2 is divided into ve parts ∂Ω2 = Γ1∪ Γ2∪ Γ3∪ Γ4∪ Σ. It is useful to notice

that on the wellbore's symmetry axis Γ3 one has r = 0. On the perforations Σ, the normal

specic ux G · n = GΣ is given, while an impermeability condition G · n = 0 is satised

on Γ2∪ Γ3∪ Γ4.

Once G computed, we impose the owrate u · n = G·n

ρ at the top of the well, as well as

an impermeability condition u·n = 0 on Γ2∪ Γ3∪ Γ4, while on Σ we prescribe a boundary

condition of Neumann type, that is p − τn · n = PΣ. A Neumann condition τn · t = 0

is set on the top Γ1, on the bottom Γ4 and on the non-perforated lateral boundary Γ2,

whereas on Σ the tangential velocity is imposed; for the sake of simplicity, we take here u · t = 0. Concerning the energy equation, one imposes the temperature T = TΣ on the

perforations Σ and the normal heat ux q · n = 0 elsewhere. Remark 1 One may also choose to impose u · n = GΣ

ρ on Σ. In this case, the radial

velocity ur is completely determined on the whole domain Ω2, without solving the

corre-sponding momentum equation. This variant of the 1D model was developed in 2 _{and the}

corresponding code was validated numerically. However, in view of the coupling with the reservoir, we prefer here to compute both ur and uz from the Navier-Stokes equations.

A relevant issue concerns the boundary condition on the top of the wellbore Γ1. Let us

notice that, even if the owrate Q is known thanks to recorded data, one cannot impose it on the outow boundary Γ1 for the transport equation (7) since Q and GΣ are related

by a compatibility condition. However, things become dierent when dealing with the coupled problem, for which a global resolution will allow us to overcome this drawback of the sole wellbore problem.

The previous problems (7), (8) and (9) are next put under weak form and proved to have a unique solution. For the rst one, we apply the Babu²ka's theorem, while for the velocity-pressure mixed formulation the well-posedness is established thanks to the Babu²ka-Brezzi theorem, under the assumptions Gn

∈ L∞_{(Ω) and ∆t suciently}

small. One may also consult 7 _{for an exhaustive analysis and numeriacl approximation}

for the reservoir model 1 _{based on the Fredholm's alternative, one can show that the}

global problem has a unique solution, under the assumptions : TΣ ∈ Hδ(Σ) for a given

δ ∈ ]0, 1], ∇λ ∈ L∞(Ω) and ∆t suciently small. The detailed proofs cand be found in 8 (see also2 _{for a similar pseudo 1D model).}

The space discretisation is then achieved on a specic 2D grid. More precisely, we consider only one cell in the radial direction and we use a regular mesh T2

h in the z

direction. We put Ω2 = ∪K∈ThK with K rectangle of width R.

We employ classical conforming nite element spaces which are compatible with the dependence in r prescribed in (10). The pressure and the temperature are approximated by piecewise constant nite elements, of Q0 type. It goes the same way for the density

thanks to its dependency on these two previous variables. The uxes are approximated by conservative lowest-order Raviart-Thomas elements RT0, while for the velocity we employ

Q1- continuous elements (which means that ur, uz,buz are piecewise linear and continuous on I). The convective terms (in the momentum, respectively the energy equations) are treated by means of two dierent upwinding schemes.

All the discrete problems are then well-posed, thanks to the Babu²ka-Brezzi theorem for ∆t suciently small, e.g.

1 ∆t ≥ ckGn hk0,Ω2 ρhmin .

This last conditions is necessary in order to ensure the well-posedness of the discrete velocity-pressure formulation, when taking into account the convective term.

4 COUPLING DARCY-FORCHHEIMER AND NAVIER-STOKES 4.1 Transmission conditions

Our goal is to couple now the reservoir and the wellbore models previously described. For this purpose, we have to introduce transmission conditions at the perforations between the two semi-discretized problems, at each time step.

Let us rst point out the main dierences between our problem and those currently considered in the literature.

Firstly, we are dealing with a multiscale modelling, since the reservoir is described by a fully 2D model while the wellbore is described by a pseudo 1D model. Secondly, we have employed as a natural variable the specic ux G1 = rρ1u1 in the reservoir, while

the wellbore model is written by using both the conservative variable G2 = ρ2u2 and

the velocity u2 = (ur, uz). So, on the one hand, the wellbore model has an additionnal

unknown and on the other hand, we cannot work with the velocity as the main unknown in the two domains (as for instance in9,10,11_{, where the density is supposed to be constant}

unknowns related to the reservoir, respectively by 2 those related to the wellbore. All the primed letters are related to the test-functions.

An integration by parts in the reservoir model yields the following interface term : Z Σ p1G01· ndσ − Z Σ T1q01· ndσ, (11)

while for the 2D wellbore model one obtains, from the momentum and energy conservation laws : Z Σ R(p2− τ2n · n)u 0 2· ndσ − Z Σ RT2q02· ndσ − Z Σ R(τ₂n · t)u0₂· tdσ. (12) So we impose, as usually when dealing with the coupling of Stokes and Darcy equations, the mass conservation and the balance of normal forces on the interface Σ, which translate into :

G1· n = RG2· n, −p1 = −p2+ τ2n · n.

Due to the viscous context, one also has to prescribe a condition on the tangential com-ponent of the uid's velocity. Several types of conditions exist in the literature. The one which seems to be in best agreement with experimental evidence is the Beavers-Joseph-Saman law, and it states that the slip velocity along the interface is proportional to the shear stress (see 10 _{and the references therein). However, the mathematical analysis}

doesn't lose in generality if one simply takes as in 9,12_{, u}

2 · t = 0. Next, the energetic

aspects imply the continuity of the temperatures and of the normal heat uxes across Σ : T1 = T2, q1· n = Rq2· n.

Finally, we add the condition ρ2u2· n = G2· n which binds together the unknowns on Σ.

4.2 Coupled problem in weak form

In order to impose a owrate Q at the wellbore head, and thus to take into account the recorded data, we turn to a global resolution of the coupled wellbore-reservoir problem, at each time step. Similarly to10,12_{, we write a mixed weak formulation binding together}

the reservoir and wellbore time-discretized formulations. We denote the unknowns of the reservoir model by

x1 = (G1, q1, p1, T1),

belonging to the space

X1 = H(div, Ω1) × H(div, Ω1) × L2(Ω1) × L2(Ω1).

   Find x1 ∈ X1Σ∗ A1(x1, x01) = F1(x01), ∀x 0 1 ∈ X1Σ0 (13) where the ane set X∗

1Σ, respectively the subspace X1Σ0 ⊂ X1 take into account the

boundary conditions, including those on the interface.

Similarly, we denote the unknowns of the wellbore model by x2 = (G2, u2, p2, q2, T2),

which belong to X2 = W × V × M × H × M , where :

W =  w =  _r Rw¯r(z) wz(z)  ; ¯wr ∈ L2(I), wz ∈ H1(Ω2), w · n = 0 on Γ2∪ Γ3∪ Γ4} ⊂ H(div, Ω2), V =  v =  _r R¯vr(z) vz(r, z)  ; vz = r Rvz(z) + R − r R bvz(z), ¯vr, vz,bvz ∈ H 1_(I), v · n = 0 on Γ2∪ Γ4, v · t = 0 on Σ} ⊂ (H1(Ω2))2, H = {w ∈ W; w · n = 0 on ∂Ω2 \ Σ}, M = {q = q(z); q ∈ L2(I)} ⊂ L2(Ω2).

We agree to denote a generic test-function by x0

2 = (χ, u 0 2, p 0 2, q 0 2, T 0 2), belonging to Y2Σ = M × V0× M × H × M.

Taking into account the nonhomogeneous boundary conditions imposed in the wellbore leads to the following ane set

X_2Σ∗ _{= W}∗_{× V}∗_{× M × H × M.}

In order to simplify the presentation, we linearize the wellbore model with respect to G2 and we are now able to bind together the three decoupled problems describing the

pseudo 1D wellbore model, in order to obtain a global weak problem :    Find x2 ∈ X2Σ∗ A2(x2, x02) = F2(x02), ∀x02 ∈ Y2Σ. (14) Remark 2 One doesn't lose in generality due to the latter linearization with respect to G2. Indeed, thanks to the decoupling of the wellbore equations, the well-posedness of the

Next, in order to write the coupled problem, we dualize the transmission conditions on Σ by means of Lagrange multipliers and we also treat the interface terms (11) and (12) appearing in both models after integration by parts. For this purpose, we rst need to introduce the following spaces :

X = {x = (x1, x2) ∈ X1× X2; G1· n, q1 · n ∈ L2(Σ)}, Y = {x0 = (x01, x 0 2) ∈ X1× Y2Σ; G1· n, q1 · n ∈ L2(Σ)}, Y0 = {x0 ∈ Y ; G01· n = 0 on ΥG\ Σ, q01· n = 0 on Υq\ Σ, u02· n = 0 on Γ1}, X∗ = {(x1, x2) ∈ X ; G1· n = G∗ on ΥG\ Σ, q1· n = q∗ on Υq\ Σ, u2· n = Q on Γ1}

and we agree to denote by X0 the vectorial space associated to the ane set X∗_{. We also}

introduce the multipliers' spaces :

L = L2(Σ) × L2(Σ),

K = L2(Σ) × L2(Σ) × L2(Σ).

We are now able to dene the bilinear forms on L × Y, respectively K × X by : IΣ(Λ, x0) = Z Σ (G0₁· n − Ru0₂· n)θdσ − Z Σ (q0₁· n − Rq0₂· n)µdσ, JΣ(Λ0, x) = Z Σ (G1· n − Rρ2u2· n)θ0dσ + Z Σ (G1· n − RG2· n)ζ0dσ − Z Σ (q1· n − Rq2· n)µ0dσ for any Λ = (θ, µ) ∈ L, Λ0 _{= (ζ}0_{, θ}0_{, µ}0 ) ∈ K. Then, putting A(x, x0) = A1(x1, x01) + A2(x2, x02), ∀x ∈ X, ∀x 0 ∈ Y, F (x0) = F1(x01) + F2(x02), ∀x 0 ∈ Y, the coupled problem can be written as follows :

           Find x ∈ X∗ , Λ ∈ L A(x, x0_{) + I} Σ(Λ, x0) = F (x0), ∀x0 ∈ Y0 JΣ(Λ0, x) = 0, ∀Λ0 ∈ K. (15)

Remark 3 Let us note that the multiplier Λ = (θ, µ) can be interpreted as (p1, T1), or

4.3 Analysis of the continuous coupled problem

This section is devoted to the study of the mixed continuous problem (15), where the convective terms appearing in the two energy equations are neglected for the moment.

First of all, we establish (by using the Fortin's trick) the following result : Lemma 1 The next two inf-sup conditions hold :

∃b1 > 0, ∀Λ ∈ L, sup x0_∈Y0 IΣ(Λ, x0) kx0_k Y ≥ b1kΛk_0,Σ, ∃b2 > 0, ∀Λ0 ∈ K, sup x∈X0 JΣ(Λ0, x) kxk X ≥ b2kΛ0k0,Σ.

Then one knows thanks to the general theory of saddle point problems (cf. for instance

5_{) that for any x solution of :}



Find x ∈ J∗

A(x, x0_{) = F (x}0_{), ∀x}0

∈ I , (16)

there exists a unique multiplier Λ ∈ L such that the pair (x, Λ) satises the initial mixed problem (15). Therefore, in what follows we study the problem (16), where :

J∗ = {x ∈ X∗; JΣ(Λ0, x) = 0, ∀Λ0 ∈ K} ,

I = x0 ∈ Y0; IΣ(Λ, x0) = 0, ∀Λ ∈ L .

Let us note that the operator A(·, ·) is of a nonstandard mixed form. Indeed, by separating the vector functions from the scalar ones, the previous problem can be written as below :    Find (U, s) ∈ U∗ × S A(U, U0) + B(s, U0) = F1(U0), ∀U0 ∈ T0 B(s0, U) − C(s, s0) = F2(s0), ∀s0 ∈ S (17) with neither A(·, ·) nor C(·, ·) symmetric and with dierent spaces for the solution, re-spectively the test-functions. Here above, U represents (G1, q1, G2, u2, q2), U0 stands for

(G0₁, q0₁, χ, u0₂, q₂0), whereas s = (p1, T1, p2, T2). We have put :

+ a(u2, u02), B(s, U0) = − Z Ω1 p1divG01dx + Z Ω1 T1divq01dx − Z Ω2 p2div(ru02)dx + Z Ω2 T2div(rq02)dx, C(s, s0) = Z Ω1 r a ∆tp1p 0 1dx − Z Ω1 r b ∆tT1p 0 1dx + Z Ω1 r d ∆tT1T 0 1dx − Z Ω1 r f ∆tp1T 0 1dx + Z Ω2 rcvρ2 ∆t T2T 0 2dx.

In order to prove that the mixed problem (17) admits a unique solution, we combine a generalization of the Babu²ka-Brezzi theorem (cf. Nicolaides 13_{) for the uniqueness with}

a Galerkin method for the existence.

We rst state some preliminary results, which will be used for both the uniqueness and the existence. The proofs of the two following lemmas are based on those given for the separate wellbore and reservoir models.

Lemma 2 There exist two positive constants β1 and β2, independent of ∆t, such that :

∀s ∈ S, sup U∈U0 B(s, U) kUk ≥ β1ksk, ∀s ∈ S, sup U0_∈T0 B(s, U0) kU0_k ≥ β2ksk.

Lemma 3 There exists a positive constant γ, depending on ∆t, such that : ∀s ∈ S, C(s, s) ≥ γ(kp1k2_0,Ω1 + kT1k2_0,Ω1 + kT2k2_0,Ω2).

Lemma 4 For ∆t suciently small, there exists α > 0 such that, for any U ∈ U0 _one

can nd U0 ∈ T0 _{satisfying :} B(s, U) = B(s, U0), _{∀s ∈ S,} kU0k ≤ c kU0k , A(U, U0) ≥ α(kG1k20,Ω1 + kq1k 2 0,Ω1 + kG2k 2 H(div,Ω2)+ ku2k 2 1,Ω2 + kq2k 2 0,Ω2).

Proof. Let any U = (G1, q1, G2, u2, q2) ∈ U0 satisfying :

G1· n = RG2· n = Rρ2u2· n, q1· n = Rq2· n on Σ. We take U0 _{= (G}0 1, q 0 1, χ, u 0 2, q 0 2) = (G 0

1, q1, χ, u2, q2) where G01 andχ will be dened later

G0₁· n = 1 ρ2

G1· n on Σ, divG01 = divG1 in Ω1.

Then obviously U0 _{belongs to T}0 _{and B(s, U) = B(s, U}0_{)for any s ∈ S. Moreover, it comes}

that : A(U, U0) ≥ ckq1k2_0,Ω1+ ku2k2_1,Ω2 + kq2k2_0,Ω2  + Z Ω1 1 rM G1· G 0 1dx + Z Ω2 χdiv(rG2)dx

where c is a global coercivity constant, which may depend on ∆t. It is useful to note that the coercivity constant for ur, denoted in what follows by c(ur), is proportional to _∆t1 .

Let us now construct G0

1. For this purpose, we introduce the auxiliary problem

       ∆ψ = 0 in Ω1 ∂ψ ∂n = R(1 − ρ2)ur on Σ ∂ψ ∂n = 0 on ΥG\ Σ ∂ψ ∂n = k on Υp ,

where the constant k is chosen such that the compatibility condition for the previous Neumann problem be checked. So one can now put G0

1 = ∇ψ + G1 and obtain, thanks

to Young's inequality : Z Ω1 1 rM G1· G 0 1dx = Z Ω1 1 rM G1· G1dx + Z Ω1 1 rM G1· ∇ψdx ≥ a(1 − δ) kG1k2_0,Ω1 − 1 4δ|ψ| 2 1,Ω1

where a is the coercivity constant of the positive denite tensor 1

rM and δ ∈ ]0, 1[ is an

arbitrary parameter. One still has to choose χ. We recall that the term R_Ω2χdiv(rG2)dx

can also be written as

R2 2 Z I χ(∂zG2z + 2 RG2r)dz.

Then by putting χ = K(∂zG2z − _R2G2r) with K a positive constant and by using that

G2r = ρ2ur on Σ, one gets : Z Ω2 χdiv(rG2)dx ≥ KR 2  k∂zG2zk2_0,Ω2 − 4ρ R2 kurk 2 0,Ω2  ≥ c kG2zk2_1,Ω2− 2Kρ R kurk 2 0,Ω2,

thanks to the Friedrichs-Poincaré inequality for G2z. It is next possible to choose ∆t, as

well as the parameters δ ∈ ]0, 1[ and K > 0, such that : c(ur) >

C2 4δ +

which ends the proof.  The previous lemmas now allow us to prove the uniqueness of the solution of problem (16), and hence of its equivalent mixed form (17).

Theorem 1 The following statement is true : ∀x ∈ J0_{, sup} x0_∈I A(x, x0₎ kx0_k Y > 0. Therefore, problem (16) has at most one solution.

Corollary 1 The mixed problem (15) has at most one solution.

The existence of a solution of (15) will be established in the next section, by means of a Galerkin approximation based on nite elements.

4.4 Space discretization of the global problem We denote by (T1

h)h>0and (Th2)h>0the families of triangulations of Ω1and Ω2previously

introduced. In what follows, we suppose that the above meshes are matched on the perforations and we agree to denote by Ehthe set of edges (belonging to both triangulations

T1

h and Th2) situated on the interface Σ.

We also assume that

(H) ρ ≥ ρ2h(z) ≥ ρ > 0 onΣ.

uniformly with respect to h, where ρ2his a piecewise constant approximation of the density

on T2 h .

We next write a conforming approximation of problem (15) based on the nite element spaces already used for separately solving the reservoir and the wellbore models. Thus, we discretize the conservative variables G1, G2, q1 and q2 by means of the lowest-order

Raviart-Thomas nite element, the scalar unknowns p1, p2, T1 and T2 by piecewise

con-stant elements, while the uid's velocity u2 is taken Q1- continuous on Ω2. Concerning

the Lagrange multipliers on the interface, we introduce the nite dimensional space : Kh = {µ ∈ L2(Σ) ; µ ∈ P0(e), ∀e ∈ Eh}

and we put :

Lh = Kh× Kh ⊂ L, Kh = Kh× Kh× Kh ⊂ K.

where the forms Ah(·, ·)and Fh(·)are obtained after an upwinding of the convective terms

and where JΣh(·, ·) is deduced from JΣ(·, ·) by replacing ρ2 with ρ2h.

We are next interested in establishing the well-posedness of the mixed problem (18). For this purpose, we rst prove that each of the bilinear forms IΣ(·, ·)and JΣh(·, ·)satisfy

an inf-sup condition, uniformly with respect to the discretisation parameter h. Lemma 5 There exist two positive constants b∗

1 and b∗2, independent of the discretisation,

such that the next statements hold : ∀Λ ∈ Lh, sup x0_∈Y h IΣ(Λ, x0) kx0_k Y ≥ b∗₁kΛk0,Σ, ∀Λ0 _{∈ K}h, sup x∈X0 h JΣh(Λ0, x) kxk_X ≥ b ∗ 2kΛ 0_k 0,Σ.

Let us now introduce the discrete kernels of the bilinear forms : J0h = x ∈ X0h; JΣh(Λ0, x) = 0, ∀Λ0 ∈ Kh ,

Ih = {x0 ∈ Yh; IΣ(Λ, x0) = 0, ∀Λ ∈ Lh} ,

as well as the ane set :

J∗h = {x ∈ X ∗

h; JΣh(Λ0, x) = 0, ∀Λ0 ∈ Kh} .

Thanks to the previous lemma, it is now sucient to study the following discrete problem :    Find xh ∈ J∗h Ah(xh, x0) = Fh(x0), ∀x0 ∈ Ih (19) which can be equivalently written as follows :

   Find (Uh, sh) ∈ U∗h × Sh Ah(Uh, U0) + B(sh, U0) = F1h(U0), ∀U0 ∈ T0h B(s0, Uh) − Ch(sh, s0) = F2h(s0), ∀s0 ∈ Sh (20) where Ah(·, ·) takes into account the convective term of the Navier-Stokes equations in the

wellbore, while Ch(·, ·) contains the convective terms coming from the energy equation in

both the reservoir and the wellbore. All these convective terms are treated by appropriate upwinding schemes.

Lemma 6 There exist two positive constants β∗

1 and β ∗

2, independent of ∆t and the

dis-cretisation, such that :

∀s ∈ Sh, sup U∈U0 h B(s, U) kUk ≥ β ∗ 1ksk, ∀s ∈ Sh, sup U0_∈T0 h B(s, U0) kU0_k ≥ β ∗ 2ksk.

Lemma 7 There exists a positive constant γ∗_{, depending on ∆t but independent of the}

discretisation, such that :

∀s ∈ Sh, C(s, s) ≥ γ∗(kp1k2_0,Ω1 + kT1k2_0,Ω1 + kT2k2_0,Ω2).

A similar result holds for Ch(·, ·), if ∆t is now taken suciently small with respect to the

discretisation parameter.

Lemma 8 Let ∆t be suciently small. Then there exists α∗ _{> 0 independent of the}

discretisation, such that for any U ∈ U0

h, one can nd U 0 ∈ T0 h satisfying : B(s, U) = B(s, U0), _{∀s ∈ S}h, kU0k ≤ c kU0k , A(U, U0) ≥ α∗(kG1k 2 0,Ω1 + kq1k 2 0,Ω1 + kG2k 2 H(div,Ω2)+ ku2k 2 1,Ω2 + kq2k 2 0,Ω2).

A similar result holds for the bilinear form Ah(·, ·), where ∆t is now related to the

dis-cretisation parameter.

We are now able to establish :

Theorem 2 The continuous coupled problem (15), without convection in the energy laws, has at least one solution.

where in the denition of JΣh(·, ·), ρ2hstands for the piecewise constant L2(Σ)-orthogonal

projection of ρ2.

The four previous Lemmas allow us to conclude that each discrete problem (21) has a unique solution (xeh, eΛh), thanks to the nite dimensional framework.

According to Lemmas 7 and 8, the discrete solutionxeh = ( eUh,esh)satises the following estimate, uniformly with respect to h :

Ge1h _0,Ω 1 + k_eq1hk0,Ω1 + Ge2h H(div,Ω2) + ku_e2hk1,Ω2 + keq2hk0,Ω2 + kp_e1hk0,Ω1 + Te1h _0,Ω 1 + Te2h _0,Ω 2 ≤ c.

The inf-sup condition on B(·, ·) ensures (cf. Lemma 6) that pe2h is bounded in L

2_(Ω 2).

It also comes that :

kdiv eG1hk0,Ω1 + kdiveq1hk0,Ω1 + k 1 rdiv(req2h)k0,Ω2 ≤ c sup s0_∈S h F2(s0) − C(s0,esh) ks0_k ,

this last estimate resulting from the weak problem (21) with a test-function x0

∈ Ih. So

one can now conclude that divGe_1h, diveq_1h, and diveq_2hare also uniformly bounded with respect to the L2-norm.

The second equation of problem (21) gives that Ge1h · n and eq1h · n are uniformly bounded in L2_{(Σ), since :} e G1h· n|e = Rρ2h |e| Z e e urhdσ, eq1h· n|e = R(qerh)|e, ∀e ∈ Eh and uerh, qerh are both bounded in L

2_(Σ).

So the sequence (exh)his bounded in the X-norm, whereas the uniform inf-sup condition

satised by IΣ(·, ·) implies that (Λe_h)_h is bounded in the L-norm. Therefore, one can extract a subsequence, still denoted by (exh, eΛh)h, weakly convergent in the space X × L

towards (ex, eΛ). Also using the approximation properties of the nite element spaces employed, a classical passage to the limit in (21) yields that the weak limit (ex, eΛ) is in fact a solution of problem (15), which ends the theorem's proof.  Remark 4 One may equally prove that the continuous problem with convection also has a unique solution for a suciently small time-step, by using the regularity of the solution of (15) together with Fredholm's alternative (as in the analysis of the separate reservoir and wellbore models).

Theorem 3 Problem (19) has a unique solution, for ∆t satisfying : ∆t ≤ min(C1h21,min, C2h2,min)

with C1, C2 independent of the discretisation.

Proof. Since we are working in a nite dimensional framework, it is sucient to prove the uniqueness of the solution. The proof is immediate due to the previous lemmas 6, 7 and 8: the positivity of Ah(·, ·) and Ch(·, ·)gives that the solution of the homogeneous problem

satises Uh = 0, p1 = T1 = 0, T2 = 0. The inf-sup condition on B(·, ·) implies p2h = 0,

which ends the proof. 

REFERENCES

[1] M. Amara, D. Capatina, B. Denel and P. Terpolilli. Mixed nite element approxima-tion for a coupled petroleum reservoir model. M2AN, 39, No. 2, 349376 (2005). [2] B. Denel. Simulation numérique et couplage de modèles thermomécaniques

puits-milieux poreux, Ph.D. Thesis, University of Pau, (2004).

[3] F. Maubeuge, M. Didek, E. Arquis, O. Bertrand and J.P. Caltagirone. A model for interpreting thermometrics, SPE, 28588, (1994).

[4] D.Y. Peng and D.B. Robinson. A new two-constant equation of state, Ind. Eng. Chem. Fundam., 15, 5964, (1976).

[5] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods, Springer Verlag, New York, (1991).

[6] J.E. Roberts and J.-M. Thomas. Mixed and Hybrid Methods, in Handbook of Nu-merical Analysis, II, North Holland, Amsterdam, 523639, (1991).

[7] V. Girault and P.A. Raviart. Finite element method for Navier-Stokes equations, Theory and Algorithms, Springer Verlag, Berlin, (1986).

[8] M. Amara, D. Capatina and B. Denel. Numerical approximation of a conservative pseudo 1D petroleum wellbore model. In preparation, (2006).

[9] C. Bernardi, F. Hecht and O. Pironneau. Coupling Darcy and Stokes equations for porous media with cracks. M2AN, 39, No. 1, 735, (2005).

[11] B. Rivière and I. Yotov. Locally conservative coupling of Stokes and Darcy ows, SIAM J. Numer. Anal., 42, No. 5, 19591977, (2005).

[12] M. Discacciati and A. Quarteroni. Convergence analysis of a subdomain iterative method for the nite element approximation of the coupling of Stokes and Darcy equations.Comput. Visual. Sci., 6, Numbers 2-3, 93103, (2004).