A semi-discrete central scheme for the approximation of two-phase flows in three space dimensions

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 81 (2011) 2296–2306

Original article

A semi-discrete central scheme for the approximation of two-phase flows in three space dimensions F. Pereira ∗ , A. Rahunanthan Department of Mathematics and School of Energy Resources, University of Wyoming, Laramie, WY 82071, United States Received 15 November 2009; received in revised form 7 June 2010; accepted 5 January 2011 Available online 17 February 2011

Abstract We present a new second-order (in space), semi-discrete, central scheme for the approximation of hyperbolic conservation laws in three space dimensions. The proposed scheme is applied to a model for two-phase, immiscible and incompressible displacement in heterogeneous porous media. Numerical simulations are presented to demonstrate its ability to approximate solutions of hyperbolic equations efficiently and accurately in petroleum reservoir simulations. © 2011 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 35L65; 65M06 Keywords: Hyperbolic conservation laws; Central schemes; Two-phase flows

1. Introduction We consider a model for two-phase, incompressible, immiscible displacement in heterogeneous porous media. In such a model, the highly nonlinear equations are of very practical importance [11,21,12]. In this paper we present a new semi-discrete central scheme for the approximation of the hyperbolic conservation law arising in two-phase, three-dimensional flows in heterogeneous porous media. The conventional theoretical description of two-phase flows in a porous medium, in the limit of vanishing capillary pressure, is via Darcy’s law coupled to the Buckley–Leverett equation. We refer the two phases, water and oil, by the subscripts w and o, respectively. We also assume that the two fluid phases saturate the pores. With no sources or sinks, and neglecting the effects of gravity, these equations become ∇ · vd = 0,

where vd = −λ(s)K(x)∇p,

(1.1)

and ∂s + ∇ · (f (s)v) = 0. (1.2) ∂t Here, vd is the Darcy velocity, and the seepage velocity v = vd /φ, where φ is the porosity which is assumed to be a constant. Furthermore, s is the water saturation, K(x) is the absolute permeability, and p is the pres∗

Corresponding author. Tel.: +1 3077665091. E-mail addresses: [email protected] (F. Pereira), [email protected] (A. Rahunanthan).

0378-4754/$36.00 © 2011 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2011.01.012

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

2297

sure. The functions λ(s) and f(s) represent the total mobility and the fractional volumetric flow of water, respectively. High-resolution, non-oscillatory central schemes are efficient numerical schemes for solving hyperbolic conservation laws arising in the simulation of multiphase flows in multidimensional heterogeneous porous media. The Lax–Friedrichs (LxF) scheme [19] is the canonical first-order central scheme, which is the forerunner of all central differencing schemes that have seen many developments in recent years. It is a non-oscillatory scheme based on piecewise constant approximate solution. However, the LxF scheme suffers from excessive numerical dissipation, which causes poor resolution at discontinuities. In 1990, Nessyahu and Tadmor [20] introduced a second-order generalization to the LxF scheme. They used a staggered form of the LxF scheme and replaced the first-order piecewise constant solution with a vanLeer’s MUSCL-type linear, second-order approximation. The central Nessyahu–Tadmor (NT) scheme offers higher resolution while retaining the simplicity of Riemann-solver-free approach. The numerical viscosity present in the NT scheme is of order O(x2r /t), where r is the formal order of the scheme. For multiphase flows in highly heterogeneous petroleum reservoirs or aquifers, this family of central schemes suffers from excessive numerical viscosity when a sufficiently small time step is enforced. In 2000, Kurganov and Tadmor [18] resolved these problems and introduced a second-order, Godunov-type central scheme. This Kurganov–Tadmor (KT) scheme has smaller dissipation of order O(x2r−1 ), and thus it can be efficiently used with time steps as small as required. It can be naturally reduced to a simple semi-discrete form. This reduced KT scheme consists of three independent building blocks – a piecewise polynomial reconstruction, a spatial flux discretization and an ODE solver. These KT central schemes can be simply extended to multidimensional systems by using a “dimension-bydimension” approach – the one-dimensional numerical flux is straightforwardly applied in all coordinate directions [18]. Over the past couple of decades, much developments have been made in the area of high-resolution central schemes to approximate solution of systems of hyperbolic conservation laws. For a review, we refer the reader to [17]. Abreu et al. in [4] analyzed central schemes for solving scalar hyperbolic conservation laws arising in the simulation of multiphase flows in heterogeneous porous media. They used the “dimension-by-dimension” extension of the KT scheme, and found out that the KT scheme is less diffusive compared to the NT scheme, particularly in the presence of high permeability field, which leads to strong restrictions on the time step selection; however, the KT scheme may produce incorrect boundary behavior. In 2001, Kurganov and Petrova [17] constructed a third-order, semi-discrete, genuinely two-dimensional central scheme for systems of conservation laws and related convection-diffusion equations. This construction is based on a multidimensional extension of the idea, introduced in [18] – the use of more precise information about the local speeds of propagation, and integration over nonuniform control volumes, which contain Riemann fans. Independently, Ribeiro [22] developed a second-order, semi-discrete, genuinely two-dimensional central scheme for two-phase flows in heterogeneous porous media to encounter the incorrect boundary behavior, which shares the same general ideas with the work of Kurganov and Petrova. Later, Balbás et al. in [9], and Balbás and Tadmor in [8] presented the extensions of the semi-discrete central schemes introduced in [17] with arbitrary order r, specially third- and fourthorder reconstructions with the possibility of additional reconstruction in the diagonal directions. More recently, it was extended for three-dimensional systems in [10]. The main contributions of this work are the derivation of a second-order (in space), semi-discrete, genuinely three-dimensional scheme based on [22], the discussion of the coupling of the hyperbolic solver with velocity fields computed using mixed finite elements, and the application of the new scheme to approximate solutions of twophase flows in petroleum reservoirs. We remark that our work is closely related to the recent work of Balbás and Qian [7], where the derivation of a semi-discrete central scheme in three space dimensions, and its application to a numerical simulation of the cloud–shock interaction modeled by the Euler equations of gas dynamics were presented. This paper is structured as follows. In Section 2, we derive the new second-order, semi-discrete, genuinely three-dimensional central scheme for hyperbolic conservation laws. In Section 3, we discuss the numerical approximation of two-phase flows using an operator splitting technique, and modify the central scheme to approximate two-phase flows in heterogeneous porous media. In Section 4, we apply the proposed central scheme for heterogeneous and homogeneous reservoirs, and present the numerical results. Section 5 contains the conclusion.

2298

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

2. Second-order semi-discrete central scheme for hyperbolic system of conservation laws In this section we outline the derivation of the second-order, semi-discrete, central scheme formulation for the hyperbolic conservation laws of the form, ut +

∂ ∂ ∂ f (u) + g(u) + h(u) = 0, ∂x ∂y ∂z

(2.1)

subject to the initial condition, u(x, y, z, t = 0) = u0 (x, y, z),

(2.2)

where u ∈ Rd , and f(u), g(u) and h(u) are nonlinear fluxes. Define xi := ix, yj := jy, zk := kz, xi±(1/2) := xi ± (x/2), yj±(1/2) := yj ± (y/2), and xk±(1/2) := zk ± (z/2). Assume that we have already computed an approximation to the solution at time level t = tn of the form,  (2.3) u˜ (x, t n ):= pni,j,k (x, y, z)χi,j,k (x, y, z), i,j,k

where pi,j,k (x, y, z) = u¯ i,j,k + (ux )i,j,k (x − xi ) + (uy )i,j,k (y − yj ) + (uz )i,j,k (z − zk ),

(2.4)

with slopes are approximated by MinMod limiter, χi,j,k is the characteristic function of the corresponding region, and u¯ i,j,k is the cell average defined by  xi+(1/2)  yj+(1/2)  zk+(1/2) 1 u¯ i,j,k := pn (x, y, z) dx dy dz. (2.5) xyz xi−(1/2) yj−(1/2) zk−(1/2) i,j,k We denote the cell interface value at the faces n uFC i,j,k :=pi,j,k (xi+(1/2) , yj , zk ),

n uBaC i,j,k :=pi,j,k (xi−(1/2) , yj , zk ),

n uRC i,j,k :=pi,j,k (xi , yj+(1/2) , zk ),

n uLC i,j,k :=pi,j,k (xi , yj−(1/2) , zk ),

n uTC i,j,k :=pi,j,k (xi , yj , zk+(1/2) ),

n uBoC i,j,k :=pi,j,k (xi , yj , zk−(1/2) ),

(2.6)

where FC, BaC, RC, LC, TC and BoC in the superscripts denote Front, Back, Right, Left, Top and Bottom Centers respectively, and the cell interface value at the vertices n uFRT i,j,k :=pi,j,k (xi+(1/2) , yj+(1/2) , zk+(1/2) ), n uFRBo i,j,k :=pi,j,k (xi+(1/2) , yj+(1/2) , zk−(1/2) ), n uFLT i,j,k :=pi,j,k (xi+(1/2) , yj−(1/2) , zk+(1/2) ), n uFLBo i,j,k :=pi,j,k (xi+(1/2) , yj−(1/2) , zk−(1/2) ),

n uBaRT i,j,k :=pi,j,k (xi−(1/2) , yj+(1/2) , zk+(1/2) ), n uBaRBo i,j,k :=pi,j,k (xi−(1/2) , yj+(1/2) , zk−(1/2) ), n uBaLT i,j,k :=pi,j,k (xi−(1/2) , yj−(1/2) , zk+(1/2) ),

(2.7)

n uBaLBo i,j,k :=pi,j,k (xi−(1/2) , yj−(1/2) , zk−(1/2) ),

where the letters F, Ba, R, L, T and Bo in the superscripts denote Front, Back, Right, Left, Top and Bottom, respectively (see Fig. 1). We compute the maximum local speeds of propagation of the discontinuities by      ∂f BaC ∂f FC x ai+(1/2),j,k (ui+1,j,k ) , ρ (ui,j,k ) := max ρ , ∂u ∂u      ∂g LC ∂g RC y ai,j+(1/2),k := max ρ (ui,j+1,k ) , ρ (ui,j,k ) , (2.8) ∂u ∂u      ∂h BoC ∂h TC z ai,j,k+(1/2) := max ρ . (ui,j,k+1 ) , ρ (u ) ∂u ∂u i,j,k As in [17,7], using these local speeds, we can define the following regions in the computational cell:

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

2299

Fig. 1. Left: The non-staggered cell [xi−(1/2) , xi+(1/2) ] × [yj−(1/2) , yj+(1/2) ] × [zk−(1/2) , zk+(1/2) ] for central differencing in three-space dimensions. Right: The reconstructing points in a magnified view of the shaded cell that is on the left.

• Regions at the cell corners where discontinuities propagate in all three directions, Vi±(1/2),j±(1/2),k±(1/2) • Regions along the edges of the cell where discontinuities propagate in two directions, Vi,j±(1/2),k±(1/2) , Vi±(1/2),j,k±(1/2) , and Vi±(1/2),j±(1/2),k • Regions across the faces of the cell where discontinuities propagate in one direction, Vi±(1/2),j,k , Vi,j±(1/2),k , and Vi,j,k±(1/2) • Interior of the cell where the solution remains smooth under an appropriate CFL condition, Vi,j,k These regions are determined by the corresponding local speeds. For instance, the region Vi+(1/2),j+(1/2),k+(1/2) is defined by Vi+ (1/2), j + (1/2), k + (1/2)  := xi+ (1/2) − Axi+ (1/2), j + (1/2), k + (1/2) t, xi+ (1/2) + Axi+(1/2),j+(1/2),k+(1/2) t  y y × yj+ (1/2) − Ai+ (1/2), j + (1/2), k + (1/2) t, yj+ (1/2) + Ai+(1/2),j+(1/2),k+(1/2) t  × zk+ (1/2) − Azi+ (1/2), j + (1/2), k + (1/2) t, zk+ (1/2) + Azi+(1/2),j+(1/2),k+(1/2) t , where



x x x x , , ai+(1/2),j+1,k , ai+(1/2),j,k+1 , ai+(1/2),j+1,k+1 Axi+(1/2),j+(1/2),k+(1/2) := max ai+(1/2),j,k 

y y y y y Ai+(1/2),j+(1/2),k+(1/2) := max ai,j+(1/2),k , ai+1,j+(1/2),k , ai,j+(1/2),k+1 , ai+1,j+(1/2),k+1 , 

z z z z Azi+(1/2),j+(1/2),k+(1/2) := max ai,j,k+(1/2) . , ai+1,j,k+(1/2) , ai,j+1,k+(1/2) , ai+1,j+1,k+(1/2)

(2.9)

(2.10)

We now compute the new cell averages in these regions, and build a new polynomial reconstruction  ˜ n+1 (x, y, z) = ˜ n+1 w w ˜ i±(1/2),j±(1/2),k±(1/2) i±(1/2),j±(1/2),k±(1/2) χ i,j,k



+

˜ n+1 ˜ n+1 w ˜ i,j±(1/2),k±(1/2) + w ˜ i±(1/2),j,k±(1/2) i,j±(1/2),k±(1/2) χ i±(1/2),j,k±(1/2) χ

i,j,k

˜ n+1 +w χ ˜ i±(1/2),j±(1/2),k i±(1/2),j±(1/2),k  n+1 n+1 ˜ n+1 ˜ ˜ + w χ ˜ + w χ ˜ + w χ ˜ i±(1/2),j,k i,j±(1/2),k i,j,k±(1/2) i±(1/2),j,k i,j±(1/2),k i,j,k±(1/2) i,j,k

+

 i,j,k

˜ n+1 w ˜ i,j,k , i,j,k χ

(2.11)

2300

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

where χ˜ i,j,k is the characteristic function of the corresponding region Vi,j,k . Then, we project this polynomial back onto the original grid to obtain the new cell averages u¯ n+1 i,j,k

1 = xyz



xi+(1/2)

xi−(1/2)



yj+(1/2)

yj−(1/2)



zk+(1/2)

˜ n+1 (x, y, z) dx dy dz. w

(2.12)

zk−(1/2)

See Fig. 1. As in [18,17], we proceed to the semi-discrete limit as t → 0, u¯ n+1 ¯ ni,j,k d i,j,k − u . u¯ i,j,k (t) = lim t→0 t dt

(2.13)

Using (2.12), we get d 1 u¯ i,j,k (t) = lim t→0 t dt



1 xyz



xi+(1/2)



yj+(1/2)



zk+(1/2)

˜ w xi−(1/2)

yj−(1/2)

n+1

zk−(1/2)

(x, y, z) dx dy dz − u¯ ni,j,k

.

(2.14)

Since the size of Vi±(1/2),j±(1/2),k±(1/2) is proportional to t3 , the size of Vi,j±(1/2),k±(1/2) , Vi±(1/2),j,k±(1/2) and Vi±(1/2),j±(1/2),k is proportional to t2 , and the size of Vi±(1/2),j,k , Vi,j±(1/2),k and Vi,j,k±(1/2) is proportional to t, we have 3 ˜ n+1 ¯ n+1 w i±(1/2),j±(1/2),k±(1/2) (x, y, z) = w i±(1/2),j±(1/2),k±(1/2) + O(t ),

(2.15)

2 ˜ n+1 ¯ n+1 w i±(1/2),j±(1/2),k (x, y, z) = w i±(1/2),j±(1/2),k + O(t ), 2 ˜ n+1 ¯ n+1 w i,j±(1/2),k±(1/2) (x, y, z) = w i,j±(1/2),k±(1/2) + O(t ),

˜ n+1 w i±(1/2),j,k±(1/2) (x, y, z)

=

¯ n+1 w i±(1/2),j,k±(1/2)

(2.16)

+ O(t 2 ),

¯ n+1 ˜ n+1 w i±(1/2),j,k (x, y, z) = w i±(1/2),j,k + O(t), ˜ n+1 ¯ n+1 w i,j±(1/2),k (x, y, z) = w i,j±(1/2),k + O(t), ˜ n+1 w i,j,k±(1/2) (x, y, z)

=

¯ n+1 w i,j,k±(1/2)

(2.17)

+ O(t),

and due to the conservation property of the reconstruction 

1 |Vi,j,k |

Vi,j,k

˜ n+1 ¯ n+1 w i,j,k (x, y, z) dx dy dz = w i,j,k .

(2.18)

We now substitute (2.15)–(2.18) in (2.14), and get   d 1 u¯ i,j,k (t) = lim t→0 txyz dt ± +

   ±

Vi,j,k±(1/2)

Vi±(1/2),j,k

¯ n+1 w i±(1/2),j,k dx dy dz +

   ±

Vi,j±(1/2),k

 1 |Vi,j,k | n+1 ¯ i,j,k − u¯ ni,j,k . w t→0 t xyz 

¯ n+1 w i,j,k±(1/2) dx dy dz + lim

Now, consider the first sum on the right hand side of the semi-discrete formulation,

¯ n+1 w i,j±(1/2),k dx dy dz

(2.19)

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

1 lim t→0 txyz

 Vi±(1/2),j,k

¯ n+1 w i±(1/2),j,k dx dy dz =

x ai±(1/2),j,k

x

2301

¯ n+1 lim w . t→0 i±(1/2),j,k

(2.20)

Given the reconstruction u˜ n , we integrate the system (2.1) over Vi+(1/2),j,k , and obtain  zk+(1/2)  yj+(1/2)

1 ¯ i+(1/2),j,k = pni+1,j,k (xi+(1/2) , y, z) + pni,j,k (xi+(1/2) , y, z) lim w t→0 2yz zk−(1/2) yj−(1/2)    ⎫ ⎤ f pni+1,j,k (xi+(1/2) , y, z) − f pni,j,k (xi+(1/2) , y, z) ⎬ dy dz⎦ . − x ⎭ ai+(1/2),j,k

(2.21)

We compute the first two integrals exactly,  zk+(1/2)  yj+(1/2) pni+1,j,k (xi+(1/2) , y, z) dy dz = yz uBaC i+1,j,k , zk−(1/2) yj−(1/2)  zk+(1/2)  yj+(1/2) pni,j,k (xi+(1/2) , y, z) dy dz = yz uFC i,j,k . zk−(1/2)

(2.22)

yj−(1/2)

For the last two integrals, we use the trapezoidal rule, and obtain  zk+(1/2)  yj+(1/2)  n  yz  f pi+1,j,k (xi+(1/2) , y, z) dy dz =

yj−(1/2) zk−(1/2) zk+(1/2)  yj+(1/2)





f pni,j,k (xi+(1/2) , y, z) dy dz = zk−(1/2)

yj−(1/2)

4



BaLT BaRBo BaLBo f (uBaRT i+1,j,k ) + f (ui+1,j,k ) + f (ui+1,j,k ) + f (ui+1,j,k ) ,

(2.23)

 yz  FLT FRBo FLBo f (uFRT i,j,k ) + f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k ) . 4

Substitute (2.22) and (2.23) in (2.21), ¯ i+(1/2),j,k = lim w

t→0

−

1 x 8ai+(1/2),j,k



 1  BaC ui+1,j,k + uFC i,j,k 2





 (2.24)

BaLT BaRBo BaLBo FRT FLT FRBo FLBo f (uBaRT i+1,j,k ) + f (ui+1,j,k ) + f (ui+1,j,k ) + f (ui+1,j,k ) − f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k )

.

¯ i,j+(1/2),k and limt→0 w ¯ i,j,k+(1/2) can be evaluated. Similarly, the terms limt→0 w For the last term in (2.19), we observe as t → 0, the region Vi,j,k becomes of the size xyz. We can now easily integrate (2.1) over Vi,j,k × [tn , tn + t]. We then use the trapezoidal rule, and get   x x ai+(1/2),j,k FC ai−(1/2),j,k BaC 1 |Vi,j,k | ¯ n+1 ¯ ni,j,k = − w ui,j,k − ui,j,k i,j,k − u t→0 t xyz x x y y z z ai,j+(1/2),k RC ai,j−(1/2),k LC ai,j,k+(1/2) ai,j,k−(1/2) − ui,j,k − − ui,j,k − uTC uBoC i,j,k i,j,k y y z z    (2.25) 1  FLT FRBo FLBo BaRT BaLT BaRBo BaLBo − f (uFRT i,j,k ) + f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k ) − f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k ) 4x   FLT  1  FRT FRBo BaRBo BaLT FLBo BaLBo − g(ui,j,k ) + g(uBaRT i,j,k ) + g(ui,j,k ) + g(ui,j,k ) − g(ui,j,k ) + g(ui,j,k ) + g(ui,j,k ) + g(ui,j,k ) 4y   FRBo  1  FRT FLT BaLT BaRBo FLBo BaLBo − h(ui,j,k ) + h(uBaRT . i,j,k ) + h(ui,j,k ) + h(ui,j,k ) − h(ui,j,k ) + h(ui,j,k ) + h(ui,j,k ) + h(ui,j,k ) 4z lim

We plug (2.24)–(2.25) in (2.19), and get y

y

z z x x − Hi−(1/2),j,k − Hi,j,k−(1/2) Hi,j+(1/2),k − Hi,j−(1/2),k Hi,j,k+(1/2) Hi+(1/2),j,k d u¯ i,j,k (t) = − − − , dt x y z

where the numerical fluxes are given by

(2.26)

2302

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306 x Hi+(1/2),j,k =

1

BaRBo BaLT f (uBaRT i+1,j,k ) + f (ui+1,j,k ) + f (ui+1,j,k ) 8

 ax i+(1/2),j,k BaC FRT FRBo FLT FLBo (ui+1,j,k − uFC + f (uBaLBo i+1,j,k ) + f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k ) + f (ui,j,k ) − i,j,k ), 2 1 FLT FLBo BaLBo FRT g(ui,j+1,k ) + g(uBaLT i,j+1,k ) + g(ui,j+1,k ) + g(ui,j,k ) + g(ui,j,k ) 8  ay i,j+(1/2),k LC FRBo BaRBo + g(uBaRT ) + g(u ) + g(u ) − (ui,j+1,k − uRC i,j,k i,j,k i,j,k i,j,k ), 2

(2.27)

y

Hi,j+(1/2),k =

1 FRBo FLBo BaLBo FRT BaRT h(ui,j,k+1 ) + h(uBaRBo i,j,k+1 ) + h(ui,j,k+1 ) + h(ui,j,k+1 ) + h(ui,j,k ) + h(ui,j,k ) 8  az i,j,k+(1/2) BoC BaLT ) + h(u ) − (ui+1,j,k − uTC + h(uFLT i,j,k i,j,k i,j,k ). 2

(2.28)

z Hi,j,k+(1/2) =

(2.29)

3. Numerical approximation of two-phase flows 3.1. Operator splitting for two-phase flow We employ an operator splitting technique to compute the solutions of the saturation (1.1) and pressure equation (1.2), which are coupled together by the seepage velocity v and the water saturation s. This technique is computationally efficient in producing accurate numerical solutions for two-phase flows [14]. Typically, for computational efficiency, larger time steps are used to compute the pressure equation (1.1). The splitting technique allows time steps used in the pressure calculation that are longer than the steps allowed under an appropriate CFL condition in the saturation calculation. We thus introduce a variable time step ts for the saturation calculation, and a longer time step tp for the pressure calculation. The pressure and thus the seepage velocity are  approximated at times tm = mtp , where m = 0, 1, . . .; the saturation is approximated at times t m,k = t m + ki=1 (ts )i for tm < tm,k ≤ tm+1 . We remark that we must specify the water saturation at t = 0. For further detail, we refer the reader to [14,2,3]. For the pressure calculation, we use a hybridized mixed finite element discretization equivalent to cell-centered finite differences [14], which effectively treats the rapidly changing permeabilities that arise from stochastic geology, and produces accurate velocity fields. The resulting algebraic system can be solved by a preconditioned conjugate gradient procedure or by a multi-grid technique. 3.2. Second-order semi-discrete central scheme for the saturation equation We use a high-resolution central scheme derived in Section 2, which is an extension of the Kurganov–Tadmor central scheme [18], for solving the scalar hyperbolic conservation law (1.2). In this section, we setup the semi-discrete form (2.26) for the saturation equation (1.2). The scalar hyperbolic conservation law (1.2) can be written as ∂s ∂ ∂ ∂ + (x vf (s)) + (y vf (s)) + (z vf (s)) = 0, ∂t ∂x ∂y ∂z

(3.1)

where x v = x v(x, y, z), y v = y v(x, y, z), and z v = z v(x, y, z) denote the x, y and z components of the seepage velocity x v, respectively. Then, the numerical flux Hi+(1/2),j,k can be written as

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

1 x FRT BaRT FRT BaRBo FRBo x FLT BaLT {f (si+1,j,k ) + f (si,j,k )} + x vFRBo [ v i,j,k {f (si+1,j,k ) + f (si,j,k )} + vi,j,k {f (si+1,j,k ) 8 i,j,k x ai+(1/2),j,k FLT BaLBo FLBo BaC FC +f (si,j,k )} + x vFLBo − si,j,k ). (si+1,j,k i,j,k {f (si+1,j,k ) + f (si,j,k )}] − 2

2303

x Hi+(1/2),j,k =

(3.2)

y

z can be modified. Similarly, the numerical fluxes Hi,j+(1/2),k and Hi,j,k+(1/2) Now, we need to have the seepage velocity at the vertices. However, the seepage velocity v, which is defined in the lowest order Raviart–Thomas space, has the components vFC , vBaC , vRC , vLC , vTC and vBoC on the faces of a rectangular parallelepiped (the components of v are pointing outward). Therefore, we use a bilinear interpolation, which preserves the null divergence necessary for the incompressible flows, to compute the seepage velocity at the vertices. For instance, to compute the seepage velocity x vFRT i,j,k , we have to use the eight cells which share the vertex FRT of the cell ijk. Thus, we define, 1  FC x FRT FC BaC FC BaC FC BaC FC vi,j,k = (vi,j,k − vBaC i,j,k ) + (vi+1,j,k − vi+1,j,k ) + (vi,j+1,k − vi,j+1,k ) + (vi,j,k+1 − vi,j,k+1 ) + (vi+1,j+1,k 16 FC BaC FC BaC FC BaC −vBaC ) + (v − v ) + (v − v ) + (v − v ) i+1,j+1,k i+1,j,k+1 i+1,j,k+1 i,j+1,k+1 i,j+1,k+1 i+1,j+1,k+1 i+1,j+1,k+1 .

(3.3)

3.3. Time marching for the semi-discrete central scheme The time integration adopted for solving (2.26) is the second-order, SSP Runge–Kutta scheme [23], s(1) = sn + tL(sn ),

(3.4)

sn+1 = (1/2)sn + (1/2)s(1) + (1/2)tL(s(1) ),

where the superscripts n and n + 1 denote time level t and t + 1, respectively, and L(s) denotes the right-hand side of the semi-discrete scheme (2.26). 4. Numerical experiments In this section, we present the results of numerical simulations of three-dimensional, two-phase flows associated with flooding problems. For the saturation calculation, we use the second-order, semi-discrete, genuinely three-dimensional central scheme (2.26). In all simulations, the reservoir contains initially 79% oil and 21% water. Water is injected uniformly into the reservoir at a constant rate of one pore-volume every five years. For the heterogeneous reservoir studies, we consider a scalar, heterogeneous absolute permeability field taken to be the logarithmic of a realization of a Gaussian random field [16,15]. The spatially variable permeability field is defined on a 32 × 16 × 8 grid with the coefficient of variation (Cv ) 1.3, where the Cv is defined as (standard deviation/mean). The water fractional volumetric flow and the total mobility are given by f (s) =

krw (s)/μw , λ(s)

where kro (s) =

 1−

λ(s) =

s (1 − sro )

krw kro + , μw μo 

2 ,

krw (s) =

(4.1)

s − srw 1 − srw

2 .

Furthermore, we use the following data in all flow simulations. Viscosity μw = 0.5 cP μo = 10 cP Porosity φ = 0.2 Residual saturations sro = 0.15 srw = 0.2

(4.2)

2304

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

Fig. 2. Water saturation plots for two-phase flow in a three-dimensional heterogeneous reservoir of 64 m × 32 m × 16 m, with Cv = 1.3, and 20 for the viscosity ratio. Top to bottom: 64 × 32 × 16, 128 × 64 × 32 and 256 × 128 × 64 grids, respectively.

For the pressure calculation, we use a hybridized mixed finite element discretization equivalent to cell-centered finite differences. The velocity field is then computed using the pressure values. For further detail, we refer the reader to [14]. We now discuss simulations in a slab geometry. We consider a three-dimensional flow in a rectangular parallelepiped, heterogeneous reservoir of size 64 m × 32 m × 16 m. The boundary conditions and the injection and production specifications for the two-phase flow equations (1.1) and (1.2) are given below. The injection is made along the face, {(x, y, z) | x = 0, 0 ≤ y ≤ 32 m and 0 ≤ z ≤ 16 m} of the reservoir, and the production is taken along the face, {(x, y, z) | x = 64 m,

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

2305

Fig. 3. Mesh refinement study indicating correct boundary behavior for the new scheme in a three-dimensional, five-spot, diagonal-grid configuration. Left to right: 32 × 32 × 32, 64 × 64 × 64 and 128 × 128 × 128 grids, respectively.

0 ≤ y ≤ 32 m and 0 ≤ z ≤ 16 m} of the reservoir; no flow is allowed along the remaining faces. We use tp = 10 days for the pressure calculation. Fig. 2 shows the mesh refinement study of the problem at 200 days. Clearly, as the mesh is refined, a better resolution of the fingering instabilities is observed. We now turn our attention to the five-spot, diagonal-grid configuration. We consider a homogeneous reservoir of size 32 m × 32 m × 32 m. The injection is made along the column, {(x, y, z) | x = 0, y = 0, 0 ≤ z ≤ 32 m} of the reservoir, and the production is taken along the column, {(x, y, z) | x = 32 m, y = 32 m, 0 ≤ z ≤ 32 m} of the reservoir; no flow is allowed along the remaining faces. Near the injection well, we use the “dimension-by-dimension” approach since the approach naturally models the injection well. Fig. 3 shows the mesh refinement study of the problem at 250 days. 5. Conclusion In this paper we presented a new second-order (in space), semi-discrete, genuinely three-dimensional central scheme for hyperbolic conservation laws. The proposed central scheme was then applied to two-phase, immiscible, incompressible, constant porosity, no gravity, no capillary pressure, three-dimensional flow in heterogeneous porous media. This scheme, which has better resolution than the “dimension-by-dimension” approach, still preserves the simplicity of the central schemes; in addition to that, it fixes the problem of incorrect boundary behavior that may result from the use of the “dimension-by-dimension” approach. We remark that the scheme presented here could also be used within operator splitting procedures to simulate miscible displacement (see, e.g. [5]), compressible flow problem (see, e.g. [13]), two-phase flow with gravity and capillary pressure (see, e.g. [6]), and three-phase flow with capillary pressure (see e.g. [1]). References [1] E. Abreu, J. Douglas, F. Furtado, D. Marchesin, F. Pereira Jr., Three-phase immiscible displacement in heterogeneous petroleum reservoirs, Mathematics and Computers in Simulation 73 (1) (2006) 2–20. [2] E. Abreu, J. Douglas, F. Furtado, F. Pereira Jr., Operator splitting based on physics for flow in porous media, International Journal of Computational Science 2 (3) (2008) 315–335. [3] E. Abreu, J. Douglas, F. Furtado, F. Pereira Jr., Operator splitting for three-phase flow in heterogeneous porous media, Communications in Computational Physics 6 (1) (2009) 72–84. [4] E. Abreu, F. Pereira, S. Ribeiro, Central schemes for porous media flows, Computational and Applied Mathematics 28 (1) (2009) 87–110. [5] C. Almeida, J. Douglas, F. Pereira Jr., A new characteristics-based numerical method for miscible displacement in heterogeneous formations, Computational and Applied Mathematics 22 (2002) 573–605. [6] J. Aquino, A.S. Francisco, F. Furtado, F. Pereira, H.P. Amaral Souto, Numerical simulation of transient water infiltration in heterogeneous soils combining central schemes and mixed finite elements, Communications in Numerical Methods in Engineering 23 (6) (2007) 491–505. [7] J. Balbás, X. Qian, Non-oscillatory central schemes for 3D hyperbolic conservation laws, in: Proceedings of Symposia in Applied Mathematics (the 12th International Conference), University of Maryland, USA, 2009. [8] J. Balbás, E. Tadmor, Non-oscillatory central schemes for one- and two-dimensional MHD equations. II. High-order semi-discrete schemes, SIAM Journal on Scientific Computation 28 (2006) 533–560.

2306

F. Pereira, A. Rahunanthan / Mathematics and Computers in Simulation 81 (2011) 2296–2306

[9] J. Balbás, E. Tadmor, C.-C. Wu, Non-oscillatory central schemes for one- and two-dimensional MHD equations, Journal of Computational Physics 201 (2004) 261–285. [10] A. Chandy, S. Frankel, Non-oscillatory central schemes for hyperbolic systems of conservation laws in three space dimensions. Available at http://www.cscamm.umd.edu/centpack/publications/, 2009. [11] G. Chavent, J. Jaffré, Mathematical methods and finite elements for reservoir simulation, in: Studies in Mathematics and its Applications, vol. 17, Elsevier, North-Holland, Amsterdam, 1986. [12] Z. Chen, G. Huan, Y. Ma, Computational Methods for Multiphase Flows in Porous Media, SIAM, Philadelphia, PA, 2006. [13] J. Douglas, D. Frias, N. Henderson, F. Pereira Jr., Simulation of single-phase multicomponent flow problems in gas reservoirs by Eulerian–Lagrangian techniques, Transport in Porous Media 50 (3) (2003) 307–342. [14] J. Douglas, F. Furtado, F. Pereira Jr., On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs, Computational Geosciences 1 (1997) 155–190. [15] F. Furtado, F. Pereira, Crossover from nonlinearity controlled to heterogeneity controlled mixing in two-phase porous media flows, Computational Geosciences 7 (2) (2003) 115–135. [16] J. Glimm, B. Lindquist, F. Pereira, Q. Zhang, A theory of macrodispersion for the scale up problem, Transport in Porous Media 13 (1993) 97–122. [17] A. Kurganov, G. Petrova, A third-order semi-discrete genuinely multidimensional central scheme for hyperbolic conservation laws and related problems, Numeriche Mathematik 88 (4) (2001) 683–729. [18] A. Kurganov, E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations, Journal of Computational Physics 160 (2000) 241–282. [19] P. Lax, Weak solutions of non-linear hyperbolic equations and their numerical computation, Communications on Pure and Applied Mathematics 7 (1954) 159–193. [20] N. Nessyahu, E. Tadmor, Non-oscillatory central differencing for hyperbolic conservation laws, Journal of Computational Physics 87 (2) (1990) 408–463. [21] D.W. Peaceman, Fundamentals of Numerical Reservoir Simulation, Elsevier, New York, 1977. [22] S. Ribeiro, New central differencing schemes for reservoir simulation, PhD thesis, IPRJ-UERJ, 2007. In Portuguese, available at http://www.labtran.iprj.uerj.br/. [23] C.-W. Shu, S. Osher, Efficient implementation of essentially nonoscillatory shock-capturing schemes II, Journal of Computational Physics 83 (1) (1989) 32–78.