在多孔介質中兩相不可壓縮不相容的流體的局部質量守恆計算法

國

立

交

通

大

學

應用數學系

碩

士

論

文

在多孔介質中兩相不可壓縮不相容的流體的

局部質量守恆計算法

A locally conservative scheme for two-phase incompressible

immiscible flows in porous media

研 究 生：楊長銘

指導教授：葉立明 教授

在多孔介質中兩相不可壓縮不相容的流體的

局部質量守恆計算法

A locally conservative scheme for two-phase incompressible

immiscible flows in porous media

研 究 生：楊長銘 Student：Chung-Ming Yang

指導教授：葉立明 Advisor：Li-Ming Yeh

國 立 交 通 大 學

應 用 數 學 系

碩 士 論 文

Submitted to Department of Applied Mathematics College of Science

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master

in

Applied Mathematics

July 2009

Hsinchu, Taiwan, Republic of China

i

在多孔介質中兩相不可壓縮不相溶的流體的

局部質量守恆計算法

學生 : 楊長銘

指導教授 : 葉立明

國 立 交 通 大 學

應 用 數 學 系

碩 士 班

摘 要

應用於本論文的水流問題的數學模型可分為兩部分。一部分就是壓力

方程式，另一部分就是 saturation 方程式。其中 saturation 方程式又分

為 transport 和 diffusion 兩部分。在此論文中我們主要著重於解 trans-

port 的部分。在此文中，我們模擬一個長兩百五十六公尺、寬兩百五十六

公尺的一個儲油槽。Locally conservative Eulerian-Lagrangian methods

(LCELM)是一個有效率的數值方法並且發展來改善在計算 transport 方程

式中水流質量守恆的部分。從數值模擬的結果，我們可以了解時間變化與

流體狀況的關係。

A locally conservative scheme for two-phase incompressible

immiscible flows in porous media

student：

Advisors：Dr.

Department of

National Chiao Tung University

ABSTRACT

The mathematical model of the waterflood problem which is applied in this paper can be divided into two sections. One is the pressure equation and the other is the saturation equation. And the saturation equation also can be pa- rtitioned into the transport stage saturation and the diffusive stage saturation. However, we will pay more attention to solve the transport stage saturation in this research. Here we construct a 256 256 meters reservoir system for simu- lation. An efficient numerical method, locally conservative Eulerian-Lagrangian methods (LCELM), is developed to compute the transport equation to improve the conservation of waterflood. From the results of the numerical simulations, we can realize the relation between temporal variation and the flow condition.

iii

Acknowledgements

在我的碩士生涯劃下完美句點之前，首先感謝我的指導教授 葉立

明教授，感謝教授的教導與鼓勵，讓我可以在完成這份論文的過程中，無

論是在理論部份或是程式部份都能給予我很大的幫助，尤其是讓我從完全

不會 C 語言變成可以獨自完成想寫的 C 語言程式，並且讓我學習了其他軟

體，讓我覺得寫程式已經不再是一個困難的問題。同時也感謝在百忙中抽

控擔任我口試委員的王夏聲老師、李華倫老師和黃聰明老師。

接著我要感謝我的父母和我的家人，因為你們的支持讓我更有自信，

同時也感謝我最好的朋友，讓我在遇到一些問題時能得到一些解決問題的

建議。

此外，要感謝我的碩士班的同學，讓我能順利完成碩士的課業，同時

也讓我學習到一些除了課業以外的東西。

Contents

中文摘要.... ... ... i

Abstract ... ii

Acknowledgements ... ii

i

Contents ... iv

Figures ... v

1 Introduction... 1

2 The Waterflood Problem... 2

3 Discretization in Temporal Domain ...4

4 Spatial Discretization ... 7

5 Algorithm ... 9

5.1 The Pressure Equation... 9

5.2 Transport... 9

5.2.1 MMOC Procedure... 9

5.2.2 LCELM Procedure... 11

5.3 Diffusive Fractional Step... 14

6 Numerical Results... 16

v

Figures

Figure 1 The point



x t

,



t

...11

Figure 2 The space-time domain

C

...13

Figure 3 The area





t

...14

Figure 4 The kind of partition (black line) is used to calculate the whole domain except the injection.

...17

Figure 5 This kind of partition (green line) is used to calculate the injection.

...18

Figure 6 The water saturation at the 40th day.

...19

Figure 7 The water saturation at the 80th day.

...19

Figure 8 The water saturation at the 120th day.

...20

Figure 9 The water saturation at the 160th day.

...20

Figure 10 The water saturation at the 200th day.

...21

1

Introduction

In this research, we develop an efficient method for the two-phase incompressible immiscible flows in porous media. We will consider the waterflooding problem in this paper. The methods in this paper should also be of value in the numerical simulation of recently developed double porosity models for swelling clay soils [1]. The numerical method combines hybridized mixed finite elements, and a new version of the modified method of characteristics, a sophisticated operator-splitting procedure for separating the transport part from the diffusion part of the saturation equation.

Our approach is to write the governing equations in terms of a global pressure [2, 3, 4]; this leads to coupled pair of equations, an elliptic equation for the global pressure and a parabolic equation for the water saturation. The primary operator splitting separates the computation of the global pressure from the saturation. This permits the use of different time steps for pressure and saturation. The second operator splitting separates the effect of transport from diffusion in the saturation calculation; this allows the use of smaller steps for the transport than for the diffusion.

In section two, we present the mathematical model of the waterflood problem. Then the temporal discretization and the spatial discretization are shown in section three and section four, respectively. And we present an overview process of numerical simulation in section three. Furthermore, the section five shows how to calculate the pressure and the saturation. And the numerical results are displayed in section six.

2

The waterflood problem

We discuss the equations for two-phase, incompressible, immiscible flow in porous media. By mass conservation, we can get the equations form [4]

Φ(x)∂S(x, t)

∂t − div(KΛΛw(S)∇xΨw(x, t)) = qext,w(x, S), (1)

− Φ(x)∂S(x, t)

∂t − div(KΛΛo(S)∇xΨo(x, t)) = qext,o(x, S), (2)

Pc(S) = Ψc(x, t) + (ρo− ρw)gz, Ψc(x, t) = Ψo(x, t) − Ψw(x, t), (3)

for x ∈ Ω, where Φ is the porosity, S is the water saturation, and K is the absolute permeability tensor of the media. The total mobility is defined as

Λ = Krw(S)

µw

+Kro(S)

µo ,

where Krw and Kro are the relative permeabilities of the water and oil phases, and µw

and µo are the viscosity of the water and oil phases. The phase mobilities are

ΛΛw(S) = Krw(S) µw , ΛΛo(S) = Kro(S) µo ,

Ψα is the α-phase potential; Pc = Pc(S) = Po− Pw is the capillary pressure ( note that

P_c0 < 0 ); ρα is the ( constant ) density of the α-phase; g is the gravitational constant;

and z is the depth. Pressures are related to potentials by the relations

Ψα = Pα− ραgz, α = w, o

The relative permeability functions Krα and the capillary pressure Pcare assumed in this

paper to be independent of x . If the global pressure [2, 3, 4] given by

P ≡ 1

2(Po+ Pw+

Z Pc

0

is introduced, equations( 1 )-( 3 ) can be written as a uniformly elliptic equation for global pressure and a convection-dominated parabolic equation for water saturation as follows :

U = −K(x)Λ(S)(∇P − (Λwρw+ Λoρo)g∇z), (5) divU = q, (6) Φ∂S ∂t + div(U Λw) + div((KΛΛwΛo)(∇Pc+ (ρw− ρo)g∇z)) = q +_{− Λ} wq−, (7)

where q = qext,w + qext,o, q+ = max(q, 0) and q− = max(−q, 0). The right-hand side of

equation( 7 ) results from assuming that only water is injected and, at a production point, the flow splits according to mobility between water and oil.

We shall assume” no flow ”boundary conditions on ∂Ω :

U · −→n |∂Ω = 0, (8)

KΛΛwΛo(∇Pc− (ρw+ ρo)g∇z) · −→n |∂Ω= 0, (9)

Where ~n is the unit outward normal vector to ∂Ω . Compatibility to the incompressibility

of the fluids requires that

Z

Ω

qdx = 0. (10)

The initial condition of the system is determined by the single relation.

S(x, 0) = Sinit(x), f or x ∈ Ω (11)

It will be convenient in the discussion of MMOC procedures to the saturation equation ( 7 ) in nondivergence form; a short calculation show that

Φ∂S ∂t + Λ 0 w(S)U · ∇S − div(KD(x, S)(∇S + (ρw − ρo)g∇z P0 c )) = (1 − Λw)q+ (12) where D(x, S) = −(ΛΛwΛoPc0)(S)

3

Discretization in temporal domain

We shall employ a time-discretization procedure based on operator splitting concepts. Assume that

∆tp = i1∆ts, ∆ts= i2∆tst

where i1 and i2 are positive integers. Let tm = m∆tp and give a function evaluated at

time tm by fm . Similarly, let tn = n∆ts and tn,k = tn+ k∆tst and let fn = f (tn) and

fn,k = f (tn,k) , respectively.

The pressure will be calculated at tm , m = 0, 1, 2 · · · . ∆ts is the time step for

the diffusive stage saturation calculation. ∆tst is the time step for the microstep for the

transport stage saturation calculation.

The time steps for the pressure and saturation variables will be allowed to be different. The numerical solution is obtained sequentially. Saturation solutions are computed for increasing values of the discretized time followed by the solution of the global pressure system after several saturation steps.

( 1 ) At the beginning, we can use the initial condition S0 _{= S}

init to determine {P0, U0}

by solving the pressure equation (in mixed form )

U0 = −K(x)Λ(S0)(∇P0− (Λw(S0)ρw+ Λo(S0)ρo)g∇z)

divU0 = q0

( 2 ) Let E1U denote the extrapolation of U given by

(E1U )(t) =          U0_{, 0 < t ≤ t}1_, t−tm−1 ∆tp U m₋ t−tm ∆tp U m−1_{, t}m _{< t ≤ t}m+1_. (13)

For tm < tn+1 ≤ tm+1 , m ≥ 0 , solve the saturation equation ( now expressed in nondivergence, mixed form )

V = −KD(x, S)(∇S + (ρw− ρo)g∇z P0 c ), (14) Φ∂S ∂t + Λ 0 w(S)(E1U ) · ∇S + divV = (1 − Λw)q+, (15) S(x, tm) = Sm(x), (16)

In ( 3.4 ), Sm(x) denotes the final values from the [tm−1, tm] . S0(x) is the initial saturation.

( 3 ) For t = tm+1 , solve for Um+1 and Pm+1

Um+1 = −K(x)Λ(Sm+1)(∇Pm+1− (Λw(Sm+1)ρw+ Λo(Sm+1)ρo)g∇z)

divUm+1 = qm+1

( 4 ) If tm+1 _{≤ T , go to ( 2 ); otherwise, stop.}

In ( 2 ) , the algorithm for this is as follows:

(i) Let tn1 and assume known {P, U, S} for t ≤ tn1 .

(ii) For n = n1, · · · , n2 = n1+ i1− 1

( a ) For k = 0, · · · , i2− 1 , compute the transport over [tn,k, tn,k+1] by solving the system

Φ∂ζn,k ∂t + Λ 0 w(ζn,k)(E1U ) · ∇ζn,k = (1 − Λw)q + , x ∈ Ω (17) (E1U ) · −→n = 0, x ∈ ∂Ω (18) ζn,k(x, tn,k) =          Sn(x), k = 0, ζn,k−1(x, tn,k), k = 1, ...i2− 1. (19)

( b ) Set Sn(x) = ζn,i2−1(x, tn,i2) = ζn,i2−1(x, tn+1) .

( c ) Compute the diffusive effects over [tn, tn+1] by solving

V = −KD(x, S)(∇S +(ρw− ρo)g∇z P0 c ), x ∈ Ω (20) Φ∂S ∂t + divV = 0, x ∈ Ω (21) V · −→n = 0, x ∈ ∂Ω (22) S(x, tn) = Sn(x), x ∈ Ω (23) (iii) Set Sm+1(x) = Sn2(x, tn2+1) = Sn2(x, t m₎

4

Discretization in spatial domain

We get the discretization from [9].

Let us return to the differential system ( 5 )-( 7 ) and restate it completely in mixed form by introducing a saturation flux variable in addition to the volumetric flow rate variable U . Then the equations take form

U = −K(x)Λ(S)(∇P − (Λw(S)ρw + Λo(S)ρo)g∇z), (24) divU = q, (25) V = −KD(x, S)(∇S + (ρw− ρo)g∇z P0 c ), (26) Φ∂S ∂t + div(ΛwU ) + divV = q +_{− Λwq}− . (27)

It should be noted that gravity terms enter into both flux variables. The no-flow boundary condition are expressed by

U · −→n = V · −→n = 0, x ∈ ∂Ω. (28)

We wish to approximate each pair, {U, P } and {V, S}, by mixed finite elements. It is com-pletely feasible and can be computationally advantageous [5, 6] to define the finite element methods for the pressure equation and saturation equation over different partitions of the domain, but in order to simplify our presentation we shall restrict our considerations to the use of the same partition for the two sets of variables. Let

Ω = [0, LX] × [0, LY ]

and set H = {HX, HY } , where HX = LX/N X and HY = LY /N Y . Then let

{Bij : i = 1, 2, · · · , N X, j = 1, 2, · · · , N Y } by Bij = [Xi−1, Xi] × [Yj−1, Yj] ; will serve for both the pressure and the saturation equations.

Let

η = η(H) = { −→ν ∈ H(div, Ω) : −→ν |Bij ∈ P1,0× P0,1 and −

→_{ν · −}→_{n = 0 on ∂Ω}}

W = W (H) = { ω : ω|Bij ∈ P0} ⊂ L

2_(Ω)

Where Pk denotes the set of polynomials of total degree k and Pk,l denotes the tensor

product of polynomials of degree k in x by those of degree l in y . Then set

M = M (H) = η × W ;

i.e., the lowest index Raviart-Thomas mixed finite element space over the partition τ .

We shall seek an approximate solution to the system ( 24 )-( 27 ) such that

(1) {Um_{, P}m_{} ∈ M, m = 0, 1, 2, · · · .}

(2) {Vn, Sn} ∈ M, n = 0, 1, 2, · · · .

In our computation, we let LX = 256, LY = 256, N X = 256 and N Y = 256. Then we can get

Ω = [0, 256] × [0, 256]

and set H = {HX, HY } , where HX = LX/N X = 1 and HY = LY /N Y = 1 . We let

Xi = i and Yj = j , and define the elements of the partition τ = τ (H) = {Bij : i =

5

Algorithm

5.1

The pressure equation

By [ 7 ], we can get the method of solving the pressure equation.

In each block Bij , Uβ, β = L, R, B, T is the out normal component of the total flux

across the edges of the block. H is the side length of block. The equations to be solved

for Pm and Um are :

U_Lm+ U_Rm+ U_Bm+ U_Tm = qmH, (1 + χgβξgβm)U m β − ξ m gβP m _{= −χ} gβξgβmUf_βm0 − ξ_gβm`fm<sub>β</sub>0 + KΛm<sub>β</sub>(Λm<sub>w</sub>ρw + Λmo ρo)βg∇z · nβ, where β = L, R, B, T , ξm gβ = 2KΛ(`msβ)/H, Λmα = Λα(`msβ), α = w, o, and `sβ = KS+ eK eS K+ eK , and g`gβ0 = KP + eKβPeβ K+ eKβ +H 2g K− eKβ K+ eKβ (Λwρw+ Λoρo)(`sβ)∇z · nβ.

In the above, ePβ and eKβ indicate the value of P and K in the element across the

β−interface, respectively.

5.2

The transport equation

The LCELM is the variational type of the MMOC. We shall explain the MMOC briefly.

5.2.1 MMOC procedure

The MMOC procedure for the waterflood problem is based on introducing a character-istic derivative for the transport part of the saturation equation written in non-divergence form. Θ(x, S, U ) =pΦ(x)2_{+ |Λ}0 w(S)U |2, (29) Θ ∂ ∂ϕ = Φ ∂ ∂t+ Λ 0 w(S)U · ∇ (30)

( Note that the characteristic direction ϕ depends on x , the saturation, and the fluid velocity.) Thus, equation( 15 ) can be written as

Θ∂S

∂ϕ + divV = (1 − Λw)q

+

and equation( 17 ) can be written as

Θ∂ζn,k

∂ϕ = (1 − Λw(ζn,k))q

+_{, x ∈ Ω, t}

n,k < t ≤ tn,k+1 (31)

We compute the transport microstep by solving equation( 31 ) with the initial values given by equation( 19 ). As seen in Figure 1, the fundamental concept in the MMOC is the discretization of the characteristic derivative by backwards differencing along the tangent

to the characteristic through the point (x, tn,k+1) back to the time level tn,k for whatever

x−point arise in the quadratures used in the finite element scheme. If we ignore for the moment sources and sink and the boundary, the transported values over a micro-step would be defined by

xn,k(x) = x −

Λ0_w(ζn,k)(E1U )(x, tn,k+1)∆tst

Φ (32)

ζn,k+1(x) = ζn,k(xn,k(x)) (33) The only reasonable criterion [ 8 ] for conservation of mass globally is that the map ( 32 ) have Jacobian identically one.

Figure 1: The point (x, tn,k+1) back to the time level tn,k.

5.2.2 LCELM procedure

Recall that the saturation equation can be written in divergence form as

∇t,x·     ΦS ΛwU     + divV = q+− Λwq− (34)

Then the transport equation can be written as

∇t,x·     ΦS ΛwU     = q+− Λwq− (35)

followed by the diffusive part given by

Φ∂S

∂t + divV = 0 (36)

Let Γ = Ω × [tn,k, tn,k+1]. Let G be a reasonable shaped, simply-connected subset of Ω ,

y(x; t) of the final value problem dy dt = ΛwU ΦS , tn,k+1 > t ≥ tn,k (37) y(x; tn,k+1) = x (38) and set xn,k(x) = y(x; tn,k) (39)

Then, let G = Gn,k be the interior of the set {xn,k(x) : x ∈ ∂G} , and let C be the tube

determined by G ,G and the integer curves ( 37 )-( 38 ); see Figure 2 and Figure 3 for

an example C in a single space variable setting. ( We let ∆tst be sufficiently small. Then

the map x → xn,k is one-to-one, so that this construction can be carried out.) Now,

denote the outward normal to ∂C by σ(x, t) and note that it is orthogonal to the vector     ΦS ΛwU    

on the lateral surface of C . Then, integrate ( 35 ) over C :

Z C ∇t,x·     ΦS ΛwU     dxdt = Z ∂C     ΦS ΛwU     · σdA = Z G ΦS(tn,k+1, x) dx − Z G ΦS(tn,k, x) dx (40) = Z C (q+− Λwq−) dxdt

Thus, mass is conserved locally in the transport step, as defined in equation( 35 ) or equations( 17 )-( 19 ) above, if Z G ΦS(tn,k+1, x) dx = Z G ΦS(tn,k, x) dx + Z C (q+− Λwq−) dxdt

The no-flow boundary condition is handled in a natural way in equation( 40 ), since the integral curves ( 32 )-( 34 ) do not exit Ω in this case. In fact, if x ∈ ∂Ω , then

the integral curve remain in ∂Ω and C has a portion of its lateral surface contained in

Γ = Ω × [tn,k, tn,k+1] . Hence, no special cares arise for subsets G close to the boundary

for these boundary conditions.

Figure 3: The area (G, tn,k+1) back to the time level tn,k.

5.3

Diffusive fractional step

We shall apply backward differencing in time over [ tn, tn+1 ]

Vn+1= −KD(Sn)(∇Sn+1+ (ρw− ρo)g∇z P0 c(Sn) ), ΦSn+1− Sn ∆ts + divVn+1 = 0,

with the no-flow boundary condition

Thus, the mixed finite element equations take the form ( 1 KD(Sn) Vn+1 , −→n1) − (Sn+1 , div−→n1) = −( 1 P0 c(Sn) (ρw− ρo)g∇z , −→n1) , −→n1 ∈ η (ΦSn+1− Sn ∆ts , n2) + (divVn+1 , n2) = 0 , n2 ∈ W.

6

Numerical results

The following data and functions are held fixed for the computational results:

Viscosity µw = 0.5 cP µo = 10 cP

Density ρw = 1 g/cm3 ρo = 0.7 g/cm3

Porosity Φ = 0.2

Residual saturations Srw = 0.2 Sro = 0.15

Absolute permeability K = 6 mdarcy

Relative permeability Krw(S) = (S−Srw) 2 (1−Srw)2 Kro(S) =  1 − S 1 − Sro 2 Capillary pressure Pc(S) = η  1 (S−Srw)2 − ζ (1−S2  ζ = S_ro2 (1 − Sro− Srw)−2 η = 3000 dynes/cm2

In this study, the geometric domain is 256 × 256 meters and the grid size of the model is 1 × 1 meters. In order to solve the transport equation, we use two kinds of partitions as shown in Figure 4 and Figure 5 to compute the injection and other sections, respectively. The red line in Figure 4 and Figure 5 represents the computational domain. In the following figures, water is injected at the lower right corner at a uniform rate and a mixture of water and oil produced at the top left corner.

Figure 4: The kind of partition (black line) is used to calculate the whole domain except the injection.

Figure 6: The water saturation at the 40th day.

Figure 8: The water saturation at the 120th day.

Figure 10: The water saturation at the 200th day.

22

References

[1] M. Murad and J. Cushman A multiscale theory of swelling porous media, II:Dial porosity models for consolidation of clay incorporating Phys- iochemical effects, Preprint #287, Center for Applied Mathematics, Purdue University, August 1996.

[2] S. N. Antontsev, On the solvability of boundary value problems for degenerate two-phase porous flow equations, Dinamika Splosnoi Sredy Vyp., 10 (1972) 28-53. In Russian.

[3] G. Chavent, A new formulation of diphasic incompressible flows in porous media, in “ Applications of Methods of Functional Analysis to Problems in Mechanics ”,Lecture Notes Mathematics. 503 (1976) 258- 270,Springer-Verlag, Berlin, New York, (P. Germain and B. Nayroles, eds.)

[4] G. Chavent and J.Jaffre, “Mathematical Models and Finite Elemenes for Reservoir Simulation” , North-Holland, Amsterdam, 1986.

[5] J. Douglas, Jr., Superconvergence in the pressure in the simulation of miscible displacement, SIAM J. Numer. Anal., 22 (1985) 962-969. [6] J. Douglas, Jr., R. E. Ewing, and M. F. Wheeler, The approximation of

the pressure by a mixed method in the simulation of miscible displacement, R.A.I.R.O., Anal. Numer., 17 (1983) 17-33.

[7] J. Douglas, F. Pereira, and L.M. Yeh. A parallelizable method for two-phase flows in naturally fractured reservoirs. Computational Geosciences, 1(3):333–368, 1997.

[8] J. Douglas Jr, C.S. Huang, and F. Pereira. The modified method of characteristics with adjusted advection. To appear in Numerische

23

Mathematik; available as Technical Report #298 , Center for Applied Mathematics, Purdue University, June 1997.

[9] J. Douglas, F. Pereira, and L.M. Yeh. A locally conservative Eulerian- Lagrangian numerical method and its application to nonlinear transport In porous media.