Gasmi and Nouri Boundary Value Problems (2015) 2015:7 DOI 10.1186/s13661-014-0256-6

0 Boundary Value Problems

a SpringerOpen Journal

RESEARCH

Open Access

Numerical simulation for two-phase flow in a porous medium

Souad Gasmi and Fatma Zohra Nouri*

"Correspondence: fz_nouri@yahoo.fr MathematicalModeling and NumericalSimulation Laboratory LAM2SIN, Badji Mokhtar University, P.O. Box 12, Annaba, 23000, Algeria

Abstract

In this paper, we introduce a numerical study of the hydrocarbon system used for petroleum reservoir simulations. This system is a simplified model which describes a two-phase flow (oil and gas) with mass transfer in a porous medium, which leads to fluid compressibility. This kind of flow is modeled by a system of parabolic degenerated non-linear convection-diffusion equations. Under certain hypotheses, such as the validity of Darcy?s law, incompressibility of the porous medium, compressibility of the fluids, mass transfer between the oil and the gas, and negligible gravity, the global pressure formulation ofChavent (Mathematical Models and Finite Elements for Reservoir Simulation: Single Phase, Multiphase and Multicomponent Flows Through Porous Media, 1986) is formulated. This formulation allows the establishment of theoretical results on the existence and uniqueness of the solution (Gasmi and Nouri in Appl. Math. Sci. 7(42):2055-2064, 2013). Furthermore, different numerical schemes have been considered by many authors, among others we can refer the reader to (Chen in Finite element methods for the black oil model in petroleum reservoirs, 1994; Chen in Reservoir Simulation: Mathematical Techniques in Oil Recovery, 2007) and (Gagneux etal. in Rev. Mat. Univ. Complut. Madr. 2(1):119-148, 1989). Here we make use of a scheme based on the finite volume method and present numerical results for this simplified system.

Keywords: compressible fluids; porous medium; multiphasic flow; finite volume method

ringer

1 Introduction

The fluids flow within porous media has an important role in various domains, such as geothermal studies, geotechnics (the mechanics of soils), chemical engineering, ground water storage, hydrocarbon exploitation (see references [1] and [2]), etc. In some cases, there are two or more fluids with different characteristics. We are concerned with the simulation of the displacement of a fluid by another one, within a porous medium, while the displacing fluid is immiscible with the fluid being displaced. Different numerical techniques, mainly based on finite elements, have been used by many authors to solve such problems, for example see [3] and [4].

In this paper we introduce a finite volume method for solving the hydrocarbon system model often used for petroleum reservoir simulations. To prevent consistency defects in the scheme, we propose to modify the mesh where the discontinuity occurs, because of the porous media. We propose to design our new scheme, called the ?diamond meshes scheme? (DMS), whose convergence is proved, and which can be used to solve the nonlin-

© 2015 Gasmi and Nouri; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the originalwork is properly credited.

ear discrete equations. Finally, numerical simulations confirm that the DMS is significantly useful for such difficult problems taking into account their physical properties.

In this case, there are only two chemical species, or components, gas and oil, where this last component may exist in both phases (gas and oil), that is to say, the heavy residual component in the oil phase and the light volatile component in the gas phase. In order to reduce confusion, we need to distinguish carefully between the ?oil component? and the ?oil phase?. This model, called a hydrocarbon system, is a simplified compositional model describing two-phase flow in a porous medium with mass interchange between them. Therefore it can predict compressibility and mass transfer effects, in the sense that it is assumed that there is mass transfer between the oil and the gas phase. In this model the ?oil component? (also called stock-tank oil) is the residual liquid at atmospheric pressure left after a differential vaporization, while the ?gas component? is the remaining fluid.

2 Mathematical model

One of the fluids wets the porous medium more than the other; we refer to this as the wetting phase fluid and we refer to the other as the non-wetting phase fluid. In an oil-gas system, oil is the wetting phase. Let us consider a bounded connex open ft of Rd (d = 2 or 3), describing the porous medium (the reservoir), with a Lipchitz boundary r, and let t be the time variable t in [0, T [, T < to. Let CGg be the mass fraction of the gas component in the gas phase, COg the mass fraction of the oil component (the light component) in the gas phase, and Coo the mass fraction of the oil component (the heavy component) in the oil phase which is equal to 1. While this distribution of the hydrocarbon components between the oil and gas phases plays an important role in a steam drive process, we cannot say that the mass of each phase is conserved because of the possibility of transfer of the oil component between the two phases. Instead, we observe that the total mass of each component must be conserved. Then the mass flux of the oil and the gas components are

Cog PgUg + PoUo, CGg pgUg.

The oil and gas mass components are

$(Cog PgSg + poSo), $CGg PgSg.

Then, for each fluid, we can write the conservation equations as d

— [$(COgPgSg + PoV] + V • (COgPgUg + PoUo) = 0

— [<P(CGgPgSg^ + V • (CGgPgUg) = 0

where (i = o, g) Si, Ui, Pi represent the saturation, the velocity, and the density of the phase i, respectively. The parameter $ is the porosity of the medium. We have

Cog + Cgg = 1,

Pg =fl(Pg, CGg), Po =f2(Po),

where Po and Pg are the oil and gaz pressure respectively.

We consider compressible fluids, with constant dynamic viscosities and where the gravity effect is neglected. Under these hypotheses, Darcy?s law combined with the mass conservation equations for each one of the component leads to the following system of partial differential equations of parabolic convection-diffusion type:

0(x)dt(PgSg + Po^'oSo) + v • (pgUg + PoaiU0) = o, (2)

0(x)-(Pou%) + v • (poO>h0Uo) = o, (3)

Uo = -K (x)—vPo, (4)

Ug = -K (x)—v Pg, (5)

where (i = o,g) Pi, —i, kri represent the pressure, the viscosity, and the relative permeability of the phase i, respectively. The parameter K is the absolute permeability of the medium and wco, c = 0,1 is the mass fraction of the component c, denoted by 0 for the heavy component and by l for the light component in the oil phase. We have

—g =f3(Pg, CGg), —o =f4(Po),

krg =f5(Sg, So), kro = f6(Sg, So).

Furthermore, we shall use the subscript S to indicate standard conditions, i.e. appropriate conditions for the temperature and the pressure of medium.

Let pOS, pGS be the density (measured at standard conditions) of the oil and the gas components, respectively. The gas formation volume factor, denoted by BG, is the ratio of the volume of free gas (all of which is gas component), measured at the reservoir conditions, to the volume of the same gas measured at standard conditions. Thus

Bg (P, r ) = Vpil.

Let WG be the weight of free gas, since VG = and VGS = W, then

Bg = —,

so that the volatility of the oil in the gas is expressed by the ratio

*V = VOS .

The mass fractions of the two components in the gas phase are

Rv PoS

C _ ' gS

By adding the last two equations and noting that CGg + COg _ 1, we obtain

(Pgs + Rv Pos)

Pg _ -R-.

The substitution of these mass fractions and densities into (1) gives, for the gas and the oil components,

д dt д dt

RV PoSSg + PoSo BG

+ V., — poSUg + PoUo) = 0,

^P^)] + v( P^) =0.

Bg )\ \ Bg We suppose that it is a saturated regime and is expressed by

So + Sg = 1.

The capillary pressure is given by

Pg - Po = Pc(So) = pc(So)pcM,

pcM = sup I Pc (So) I and 0 < pc(So) < 1.

We define the mobility of each phase by the formula

Xi = —, i = o, g,

and the total mobility X by

X = Xo + Xg.

To simplify, we set

Ph = Poыh,

P = Pg + Po,

b = Pg Xg + PoXo, d = Pg - Po.

(2) (13) (4) (15)

Let us now introduce the new unknowns, namely the reduced saturation and the global pressure in the following way: if we denote by Sim and SiM, the residual and the maximum saturations of the fluid i _ o,g, respectively; the reduced saturation is given by

S , ()

1 - S - S

1 Sg,m So,m

0 < S < 1. (17)

The ?global pressure? was first introduced by Chavent and Jaffre1[] in the following form: P _2(Pg + Po) + Y (S), (18)

S Xg — Xo

Y (S) = -l JLYJLp'M )pcMdH. (9)

Y(S) = 2/o «(t) dï, (0)

Xg (S)—Xo(S)

"(S) _ X(S) p'c(S)PCM (21)

is the capillary diffusion. Therefore, our model is given by the following simplified system:

d d $(x) - (Ph0S) + $(x)So,m - P)

- V • (K(x)p%(S)VP) + V • (K(x)Pha(S)V(S)) _ 0, (22)

d d *(x) — (pS) + <(x) — (p So,m + Pg) d t d t

- V • (K(x)b(S,P)VP) + V • (K(x)d(P)a(S)VS) _ 0, (23) X(S)VP • n _0, a(S)VS _0, on T x (0, T), (24) S(x, 0) _S0(x), P(x, 0) _ P0(x), in a, (25)

where <(x) _ ^(x)(1 - So,m - So,g).

3 Weak formulations

First we introduce the following spaces:

H(div, a)_ {v e (L2(a, (0, T)))d, div(v) e L2(a, (0, T)), d _ 2,3}, (26)

V(a) _ {v e H(div, a), v • n _ 0 on T},

W(a) _ {v e V(a), v(x, T)_ 0 in a.

The weak formulation of problem (22)-(25) is written as

($(x)p0s, - (K (x)p°0ko(S)VP, v v)ft + (K (x)p0a(S)VS, V v)ft = (fi, v), (28)

($(x)pS, ^ - (K(x)b(S, P)VP, V v)ft + (K(x)d(P)a(S)VS, V v)a = (2, v), (29)

where (•, 0ft is the inner product defined on W(ft). A theoretical study of the existence and uniqueness of the weak solution was done and the details of the results can be found in [5].

4 Finite volume approximation

In this section we study the approximation of solutions to our model in the discrete finite volume framework. This family of schemes allows very general meshes and deals with the main properties of the physical features of the treated problems. The time interval [0, T[ is divided into finite sub-intervals [tn, tn+i [ of length At„, n = 0,..., M with t0 = 0 and tM = T. The space domain (the reservoir ft) is discretized by a non-structured stitching T0.

4.1 Notations

We introduce the following notation:

- Let |K | denote the cell K surface, N(K) the set of triangles having a common side with the cell K.

- Let eK,L be the common side of the triangles K and L and "" K>L be the normal oriented from K towards L.

- "" ei is the external normal corresponding to the part of ei at the boundary r.

- Let S0 be the set of sides of the stitching T0 and S*h be the set of the interior sides.

- For a given side e, let us denote by S and N the extremities e, by W and E the two triangles where e = W n E;by xe the diamond cell associated with e given by connecting the centers of gravities of the cells W and E with the extremities N and S of e.

- ((ei)i=i,4) are the four segments forming the border of xe.

- ""s = "nn— xi, —yi) is the normal on s i outgoing of xe.

- For a given node, V(N) is the set of triangles with this node in common.

For the numerical resolution of problem (22)-(25), the first two equations will be dis-cretized separately. For more details on finite volume methods, see for example [6] and [7].

4.2 Discretization of the first equation

Let C be a cell of the stitching T0, at time tn; we integrate (28) on C to get j[<b(x)p0,n(Sn+1- Sn)] dx + ^ jc[^(x)(p0,n)' Sn (Pn+i -Pn)] dx

+ ^ jc[mSo,m(p0,n)'(P"+i-Pn)] dx

- E f [K(x)p0,n K(S)VPnVe] dx

DeN (C)J ecD

£ f K(x)phn an(S)V(Sn)ne]<

DeN (C)J eCD

_ 0. (30)

Therefore we obtain

H $C phC (SC+1- SC)

+ H ®C (phhC )'SC (PC+1- PC) + H <CSo,m (p^ohC )'(PC+1- PC) - E f [KPh,nK(S)VPnne] dx

DeN(Cr eCD

+ E f [KPh,nan(S)V(Sn)ne] dx

DeN (C) eCD

_ 0. (31)

This implies that

Hph,,CSC+1 + H(Poh,C)'(^CSC + <C Som)PC+1 _ -H ^PohCnSC - ^ (PohC)' (^C - <cSo,m)PC

E Fie " £ FD,e, (2)

DeN(C) DeN(C)

where (ph,C)' is the space derivative of (phC), and FC,e and FD,e are the numerical flows of convection and diffusion, which are approximated by

FC,e _ Kept;nK,eVePn | e|, (3)

Fno,e _ KePhh,,nanVeSn | e|, (4)

where e is a side of the stitching, "" e is the normal of e outgoing from C; Ke, p^", and kno e denote the interpolations on e of the functions K, pih, and ko, respectively, while VePn is the approximation of the gradient of the pressure on the side e. The construction of the approached gradient on e is done by the so-called Green-Gauss method. It consists of approaching the gradient by its means on a co-volume in the form of a diamond constructed around the side e. We write

phnK, _ ph(pnMsn), (5

ph^n _ ph(Pne)a(Sn), (6)

VeP = E \(PNl(z) + PN^Ielte, (7)

IXeI ee3Ze 2

where xe is the diamond cell associated with e, and Ni(s) and N2(s) are the extremities of s, one of the four segments forming dxe (the boundary of xe) and"" e is the unit normal vector. The values of P at the centers W and E are PW and PE, while the values at nodes N and S are interpolated on the boundary and denoted by PN and PS. For convenience, we omit the indication n Hence at each node N we have

Pn = E yK(N)Pk, (8)

K eV (N)

where V(N) is the set of triangles sharing the node N, PK is the value of P at the center of the cell K and yK(N) are the interpolation weights. The gradient of the saturation is given by

VeS = E i (SNi (s) + SN2(s))|s"" e, (9)

IXel sedxe 2

Sn = E yK(N)SK. (0)

KeV(N)

4.3 Discretization of the second equation

In the same way, integration of (29) implies

H $cpCSC+i + (*c (pCn)'SC + (phCn)'So,m + (phCn)')PC+i = --C PC SnC - H ($C (ph0;C )'SnC - (ph0;C )%,» + (ph0£)') PC

E Gle " E GDte, ()

DgN (C) DeN(C)

where GnCe and GD,e are the numerical flows of the convection and diffusion approximated

by C,e D,e

Gnce = KebneVePn\e\, (2)

1D,e :

GDe = KednanVeSn, (43)

bn = bn{Sne, PI), (44)

dne an = dne(Pne )a(S"e), (5)

where e is a side of the stitching that limits the cell C, VePn is given by (42) and (44), and VeSn is given by (43) and (45).

5 Numerical experiments

The problem given in (22)-(25) was said not to be typical of a hydrocarbon system simulation, but it can be designed to test the resolution capabilities of the numerical method for problems involving sharp fronts.

The numerical tests are done on a homogeneous isotropic reservoir. The physical setting in 2-D was as follows. A horizontal rectangular reservoir a _]0,300[ x ]0,60[, dis-cretized with a structured mesh made up of 3,000 cells, with an intrinsic permeability of 0.001 was initially filled with a mixture of gas and oil, their residual saturations are 0.02 and 0.03, respectively. The initial fluid (oil) distribution was taken to be uniform on the whole reservoir surface, and it had an associated initial pressure P0 _ 2,000 bar. The porosity is $ _ 0.02. The mobilities and the capillary pressure are given by Pc _ (s2 -1)/2, Xo _ 0.5, Xg _ 0.2, and ¡xw _ 1, and ¡xo _ 3, the time step is At _ 10 s. These tests were carried out using the free and open-source simulator for flow and transport processes in porous media, based on the Distributed and Unified Numerics Environment DUNE (see www.dumux.org).

5.1 Results and discussion

The numerical results shown in Figures 1-3 give the pressure for the wetting and nonwetting fluid, and the saturation for the non-wetting one. Because of the complexity of the problem, we have introduced new unknowns, namely the global pressure and the reduced saturation, in order to reduce the number of unknowns from six to two (for more details

Figure 1 Non-wetting pressure at t = 0,50,100,and 175.

Figure 2 Wetting pressure at t = 0,50,100, and 175.

see [5]). An explicit non-structured finite volume scheme has been used to solve this simplified problem with the new unknowns, and because of the non-structured meshes, we proposed a method based on ?diamond cells? to approximate the gradient. These results show that the scheme is very stable.

6 Conclusion

In this paper we proposed a new scheme based on the finite volume method for solving the displacement of a fluid by another one, within a porous medium, while the displacing fluid is immiscible with the fluid being displaced. This gives the hydrocarbon system model often used for petroleum reservoir simulations. To prevent consistency defects in the scheme, we suggested to modify the mesh where the discontinuity occurs, because of the porous media. This convergent scheme, called the ?diamond meshes scheme? (DMS) was designed, to solve the associated nonlinear discrete equations. Finally, the numerical results, which are linked to our theoretical ones in [5], confirm that the DMS is significantly useful for such difficult problems (see Figures i-3).

Competing interests

The authors declare that they have no competing interests.

Authors? contributions

All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Acknowledgements

This work was financially supported by the national research project (PNR08/23/997,2011-2013). The authors are very

grateful to the referees for their valuable and helpful comments, remarks and suggestions.

Received: 26 August 2014 Accepted: 1 December 2014 Published online: 13 January 2015

References

1. Chavent, G, Jaffre, J: Mathematical Models and Finite Elements for Reservoir Simulation: Single Phase, Multiphase and Multicomponent Flows Through Porous Media, pp. 1-40. Elsevier, Amsterdam (1986)

2. Gagneux, G, Lefevere, A-M, Madaune-Tort, M: Analyse mathématique de modèles variationnels en simulation pétrolière: le cas du modèle black-oil pseudo-compositionnel standard isotherme. Rev. Mat. Univ. Complut. Madr. 2(1), 119-148(1989)

3. Chen, Z: Finite element methods for the black oil model in petroleum reservoirs, pp. 1-28. I.M.A. Preprint Series #1238 (1994)

4. Chen, Z: Reservoir Simulation: Mathematical Techniques in Oil Recovery. SIAM, Philadelphia (2007)

5. Gasmi, S, Nouri, FZ: A study of a bi-phasic flow problem in porous media. Appl. Math. Sci. 7(42), 2055-2064 (2013)

6. Afif, M, Amaziane, B: On convergence of finite volume schemes for one-dimensional two-phase flow in porous media. J. Comput. Appl. Math. 145, 31-48 (2002)

7. Boyer, F, Hubert, F: Finite volume method for 2D linear and nonlinear elliptic problems with discontinuities. SIAM J. Numer. Anal. 43(6), 3032-3070 (2008)