Numerical Algorithms 34: 117–125, 2003.  2003 Kluwer Academic Publishers. Printed in the Netherlands.

Numerical simulation of two-phase flow through heterogeneous porous media M. Afif a and B. Amaziane b,∗ a Université Cadi Ayyad, F.S.S.M., BP 2390, Marrakech, Maroc

E-mail: [email protected] b Université de Pau, LMA–CNRS, FRE2570, Av. de l’Université, 64000 Pau, France

E-mail: [email protected]

Received 4 December 2001; accepted 4 February 2003

A mixed finite element method is combined to finite volume schemes on structured and unstructured grids for the approximation of the solution of incompressible flow in heterogeneous porous media. A series of numerical examples demonstrates the effectiveness of the methodology for a coupled system which includes an elliptic equation and a nonlinear degenerate diffusion–convection equation arising in modeling of flow and transport in porous media. Keywords: finite volume method, mixed hybrid finite element, nonlinear convection– diffusion, porous media, unstructured grids AMS subject classification: 76M12, 65M12, 35K65

1.

Introduction

Multiphase flow of fluids in porous media is physically and chemically complex. It involves heterogeneities in the porous media at many different length scales and complicated processes such as diffusion and dispersion. It is well known that the transport term in the flow equations are governed by fluid velocities. Thus accurate numerical simulation requires accurate approximations of these velocities. Often the flow properties of the porous media vary abruptly with sharp changes in lithology. Consequently, the coefficient in the flow equation is quite rough. Flow simulation in petroleum and environmental applications has been extensively studied using finite element methods in the last two decades (see, e.g., [9, and references therein]). Also, discretizations using both finite element and finite volume methods are presented in [7]. More recently, finite volume methods were developed and analyzed for immiscible two-phase flow in porous media in the case where the diffusion term is neglected (see [8, and references therein]). This approach leads to robust schemes applicable for unstructured grids and the approximate solution has various interesting ∗ Corresponding author.

118

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

properties which correspond to the properties of the physical solution. These methods have been useful for advective flow problems because they combine element by element conservation of mass with numerical stability and minimal numerical diffusion. The purpose of this paper is to discuss the applicability of the mixed finite element methods and finite volume methods to various types of flow in porous media. These problems occur namely in saltwater intrusion, in groundwater contaminant transport, nuclear waste repository studies and petroleum reservoir simulation. We focus on twophase immiscible flow in heterogeneous porous media. The equation governing these types of flow can be effectively rewritten in a fractional flow formulation; i.e., in terms of a global pressure and saturation as the primary variables; see, e.g., [6]. This formulation leads to a coupled system of partial differential equations which includes an elliptic pressure–velocity equation and a nonlinear degenerate parabolic saturation equation. The saturation equation is convection dominated and thus special care should be taken in discretization. The diffusion term is small but important and cannot be neglected. The outline of the paper is as follows. Section 2 contains a short description of the mathematical and physical model used in this study. In section 3, the numerical schemes are presented with emphasis in section 3.1 on the mixed finite element method employed for the solution of the pressure equation. Then, the vertex-centered finite volume method used for the solution of the saturation equation is illustrated. Section 4 is devoted to the presentation of the results of some selected numerical investigations for the sake of illustration. Results are given for a quarter five-spots problem with variable permeability. Additional conclusions are drawn in section 5. 2.

Mathematical and physical model

We consider the immiscible displacement of one incompressible fluid by another in a horizontal reservoir  ⊂ R2 over a time interval ]0, τ [. The governing equations can be written as follows (see [6]): Pressure equation:
q = −d(u)K(x)∇P ; div( q ) = 0 in Qτ , (1) q · n| 2 = 0; P | 3 = P0 on [0, τ [. q · n | 1 = −qd ; Saturation equation:  0
u(x, t)
1,        ∂u  + div b(u) q − div K(x)∇α(u) = 0 (x) ∂t   u| 1 = 1; K∇α(u) · n | 2 = 0; u| 3 = 0   u(x, 0) = u0 (x)

in Qτ , in Qτ , on [0, τ [, in ,

(2)

where the boundary splits up into three parts such that = 1 ∪ 2 ∪ 3 and

i ∩ j = ∅ for i = j ; 1 is the part of the boundary where the water is injected,

2 is the impervious part of the boundary and 3 is the producing part of the boundary. Here K(x) and (x) are the absolute permeability tensor and porosity of the porous

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

119

medium. α(u), b(u) and d(u) ∈ C 1 ([0, 1]) are nonlinear functions which depend on the mobilities and the capillary pressure with diffusion coefficient α  vanishing for two values of saturation: α  (0) = α  (1) = 0 (degeneration of the diffusion term) and α  (s) > 0 ∀s ∈ ]0, 1[. b is a monotone increasing function and d(s)  d0 > 0 ∀s ∈ [0, 1]. qd is a flow rate specified at 1 and P0 is a given pressure at 3 . In addition, an initial condition u0 is specified. For more details on the assumptions on the data see [4]. Also Qτ denotes ×]0, τ [. The unknowns are P (x, t) the global pressure, q (x, t) the total velocity and u(x, t) the saturation of the wetting fluid. Notice that this formulation allows the use of the same basic simulation techniques to deal with both miscible and immiscible displacement.

3.

Numerical methods

In this section, the mixed finite element method employed for the pressure–velocity equation combined to a finite volume method for the saturation equation are presented. This approach leads to robust schemes applicable for unstructured grids and approximate solutions having various interesting properties related to the physical solution. Before describing the numerical discretization of the coupled problem (1)–(2), we give some notations. Let (tn )n=0,...,Nτ be a partition of [0, τ ] with a time step t = tn+1 − tn and let h = (Ti )i=0,...,Ne be an admissible triangulation of , h = |Mi |1/2 verify supi |Mi |
δ · h2 , (Mi )i=0,...,Ns the Donald dual mesh, where h := mini the barycenter of T ∈ h is such that xT = M∩T =∅ ∂M ∈ T . We denote by and xT  xM := T ∩M=∅ ∂T ∈ M, l := ∂Mi ∩ ∂Mj ∩ T the line segment between the points xT and xij the midpoint of (xMi , xMj ) and let Lh = {l ∈ ∂M\ , for M ∈ h } (see figure 1 for an illustration). For simplicity we assume that the function K is piecewise constant, we set KT = K|T and for  ∈ L∞ (), we set M = (1/|M|) M (x) dx. 3.1. Mixed finite element method In this section, we describe a mixed finite element method for the accurate approximation of the pressure equation. This method conserves mass cell by cell and produces a direct approximation of the two variables pressure and velocity. A sequential approach is adopted to solve in time the coupled problem (1)–(2). Note that the coupling is determined by the total mobility dependent d(u). When advancing from time tn to tn+1 one can assume that d(u) := d(u(x, tn )) and solve (1) for q and P . For convenience in notation, we will present the method used for a model problem similar to the pressure–velocity equation that we have to solve at each time step . More precisely, we consider an elliptic problem with homogeneous boundary condition:
q = −A(x)∇p; div( q ) = f in , (3) p=0 on ,

120

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

where A is a symmetric positive definite tensor and f is a source term. Next we introduce the following finite element spaces (see, e.g., [5]):

 2  Wh := sh ∈ L2 () ; ∀T ∈ h , sh |T ∈ RT0 ,

 Mh := wh ∈ L2 (); ∀T ∈ h , wh |T ∈ P0 ,

 Lh := λh ∈ L2 (Eh ); λh |l ∈ P0 (l), ∀l ∈ ∂T , ∀T ∈ h and λh = 0 sur , where Eh := {l ∈ ∂T ; ∀T ∈ h } and RT0 is the lowest order Raviart–Thomas space. Our implementation for the test problems uses a triangular mesh on the domains with a piecewise constant pressure and piecewise continuous flux approximations (for other examples see [5]). A Lagrangian multiplier λ is introduced in order to get the continuity of the normal component of the velocity across the inter-element boundaries (see [5]). The mixed hybrid finite element approximation ( qh , ph , λh ) ∈ Wh × Mh × Lh of ( q , p, λ) is the solution of the following problem: 

 −1  A qh · sh dx − ph div(sh ) dx + λh sh · n T ,l dl = 0 ∀sh ∈ RT0 ,    T l  T l∈∂T    div( qh )wh dx = f wh dx ∀wh ∈ P0 ,  T T       µh qh · n T ,l dl = 0 ∀µh ∈ Lh .   l

(4)

T ∈h l∈∂T

Once the approximate formulation has been written, the resulting system of equations is very large and the next step is to choose the most efficient way to implement it. It is possible (see, e.g., [5]) to reduce the solution of the system (4) to the solution of a linear system involving only the Lagrangian multipliers vector λh , with a symmetric positive definite and sparse matrix. Once this system is solved, qh and ph are obtained through a simple post-process at the element level. We have performed test cases on both constant and variable permeability fields. The hybrid method performs well in every case and it is a robust technique for the solution of the problem related to the numerical discretization of mixed problems. 3.2. Finite volume method In this section, we describe the finite volume method used for the aproximation of the saturation equation. A Godunov scheme is used for the convection term and a special discretization of the diffusion term lead to a robust scheme which satisfies a discrete maximum principle on unstructured grids. Only a short description of the method employed in this work will be given. The interested reader is referred to [1–4] for more details.

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

121

For an initial condition u0 ∈ L∞ (), we set u0M = (1/|M|) M u0 (x) dx. Let unM be an approximation of u(xM , tn ) which will be defined precisely in the sequel. Integrating (2) over the set M × [tn , tn+1 ], where M ∈ h , we obtain the following scheme

tn+1   n+1   n b(u)( q · n M,l ) dl dt M uM − uM |M| = − l∈∂M

+

tn



l∈∂M\

l tn+1

tn



  K∇α(u) · nM,l dl dt,

(5)

l

where n M,l is the outward normal to l ∈ ∂M. Using an explicit approximation of the convection term and an implicit approximation of the diffusion term, we get t = unM − b(u)nl ql n · n M,l |l| un+1 M M |M| l∈∂M

t + ∇α(u)n+1 · KT n M,l |l|, T M |M| l∈∂M∩T \ T ∩M=∅

where b(u)nl

   b unM   = b unMl

if ql n · n M,l  0, if ql n · n M,l < 0

with Ml ∈ h is such that l ∈ Ml ∩ M. The advection term is approximated by an upwind Godunov scheme and the diffusion term is approximated in the following way using a piecewise linear approximation for ∇α:

∇α(u)n+1 · KT n M,l |l| T l∈∂M∩T \

=

  n+1    α uM − α un+1 |T |∇χMl ,T · KT ∇χM,T , Ml

l∈∂M∩T \

q = 0 we have the following where χM ∈ P1 such that χMi (xMj ) = δij . Finally, since div Semi-Implicit Finite Volume scheme (S-I.F.V.):   + t   n  = unM + b uMl − b unM − ql n · n M,l |l| un+1 M M |M| l∈∂M

     t DM,l α un+1 − α un+1 , (6) + Ml M M |M| l∈∂M\ where DM,l = −|T |∇χMl ,T · KT ∇χM,T for all l ∈ ∂M\ . Note that we can also use an explicit [respectively implicit] approximation for both the diffusion and the convection terms to obtain an Explicit (E.F.V.) (respectively Implicit (I.F.V.)) Finite Volume scheme. In [4] we proved that these schemes satisfy a discrete maximum principle under appropriate CFL conditions and some convergence

122

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

results are presented for the saturation equation with a given total flux field, i.e., the case where the pressure and the saturation equations are decoupled. We refer to [3] for more details on the comparison of these schemes for the one-dimensional problem. In this paper, we extend these ideas in two spatial dimensions on unstructured grids which are preferred when treating geometrically complex geoligical structures because of their inherent flexibility.

4.

Numerical results

In this section, we present some numerical results in 2D based on the schemes presented in this paper. The two example problems used to illustrate the methodology are (1) isotropic case and (2) anisotropic case. Reservoir simulations were performed for a standard two-phase flow problem. All simulations were performed on a quarter fivespot pattern in which the spatial and time variables have been normalized. An IMPES simulator, based on the C++ programming language, has been performed which applies the mixed finite element method for computing pressure and fluid velocity approximations and the finite volume methods for approximating the saturation. The numerical simulations were performed on a PC Intel Pentium III 256 Mo RAM. The test problem involves simulation with fixed porosity  = 0.2 and a constant pressure at the production well P0 = 0. The results obtained are satisfactory even in the convection dominated case. The computed approximate solutions satisfy the discrete maximum principle and are monotone. We also demonstrate the ability of our approach in accurately capturing the fronts typical of immiscible displacement problems. 4.1. Test problem 1: isotropic case In this example, we consider a heterogeneous isotropic reservoir with two different permeability distributions as shown in figure 1, with values ranging from 10−5 to 1: (green, K = 10−5 ) and (blue, K = 1). An injection well is placed at the lower left-hand corner of the reservoir and a production well is placed at the top right-hand corner of the reservoir. The simulation is performed with a constant flow rate at the injection well qd = 0.1. As a comparison for the three finite volume schemes, the CPU times are presented in table 1. Table 1 CPU times for test problem 1. Ne = 6034 h = 2 · 10−2 t C.P.U.

E.F.V. CFL = 1

S-I.F.V. CFL = 1

I.F.V. CFL = 1

1.6 · 10−4 20 min 40 s

1.73 · 10−2 2 min 49 s

1.73 · 10−2 5 min 25 s

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

Figure 1. Dual mesh of .

Figure 2. Total velocity.

Figure 3. Saturation contours.

Figure 4. Saturation 2D plot.

123

As expected, the E.F.V. scheme leads to large CPU time computing taking into account the restrective stability criteria. It turns out that, in general, the approximation obtained by the S-I.F.V. scheme described here is better. To get an idea of how fluid flows through the field, the total velocity (see figure 2), the water saturation contours (see figure 3) and surface plots of saturation (see figure 4) are presented. Notice that the evolution of the front behaves in a physically realistic manner. 4.2. Test problem 2: anisotropic case In this example, we consider a heterogeneous anisotropic reservoir (figure 5). The values for the permeability tensor are:   1.0 0.2 for x2 ∈ [0, 0.4[ ∪ ]0.7, 1], K(x1 , x2 ) = 0.2 1.0   1.0 10−5 for x2 ∈ [0.4, 0.7]. K(x1 , x2 ) = 10−5 1.0

124

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

Figure 5. Dual mesh of .

Figure 6. Pressure contours.

Figure 7. Saturation contours.

Figure 8. Saturation profiles.

The test problem involves simulation with constant flow rate qd = 1 at the left of the reservoir, i.e., on the boundary 1 = {0}×]0, 1[. The pressure contours are presented in figure 6, the water saturation contours in figure 7 and saturation profiles on the line x2 = 0.5 in figure 8 with a close up view of the saturation profiles near the discontinuity. The results of this simulation show that the fluid initially moves away from the left side uniformly and then is influenced by the anisotropic permeability variations in the field. 5.

Concluding remarks

A mixed finite element method was used to obtain an accurate approximation of the flow equation and a vertex-centered finite volume method for the saturation equation. Numerical simulations from 2D tests on unstructured grids show that this approach leads to a set of robust schemes. In the future, we will consider the analysis of the convergence

M. Afif, B. Amaziane / Numerical simulation of two-phase flow

125

of this coupled scheme in 2D for heterogeneous anisotropic medium. Also future work will concentrate on the extension of the approach described in this paper to 3D problems. References [1] M. Afif, Schémas volumes finis pour une équation non linéaire de diffusion–convection dégénérée, in: Actas de las VI Journadas Zaragoza–Pau de Matemática Aplicada, eds. M. Madaune-Tort et al., 1999, pp. 29–36. [2] M. Afif and B. Amaziane, Analysis of finite volume schemes for two-phase flow in porous media on unstructured grids, in: Finite Volumes for Complex Applications II – Problems and Perspectives, eds. F. Benkhaldoun and R. Vilsmeier (Hermès, Paris, 1999) pp. 387–394. [3] M. Afif and B. Amaziane, On convergence of finite volume schemes for one-dimensional two-phase flow in porous media, J. Comput. Appl. Math. 145 (2002) 31–48. [4] M. Afif and B. Amaziane, Convergence of finite volume schemes for a degenerate convection– diffusion equation arising in flow in porous media, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5265–5286. [5] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods (Springer, New York, 1991). [6] G. Chavent and J. Jaffré, Mathematical Models and Finite Elements for Reservoir Simulation (NorthHolland, Amsterdam, 1986). [7] G. Chavent, J. Jaffré and J.E. Roberts, Mixed-hybrid finite elements and Cell-centred finite volumes for two-phase flow in porous media, in: Mathematical Modeling of Flow through Porous Media, eds. A.P. Bourgeat et al. (World Scientific, London, 1995) pp. 100–114. [8] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, in: Handbook of Numerical Analysis, Vol. VII, eds. P.G. Ciarlet and J.L. Lions (North-Holland, Amsterdam, 2000) pp. 723–1020. [9] R.E. Ewing, Computational sciences in environmental applications, in: Computational Science for the 21st Century, eds. M-O. Bristeau et al. (Wiley, New York, 1997) pp. 250–259.