THEORY AND APPLICATION OF FRACTIONAL STEP CHARACTERISTIC FINITE DIFFERENCE METHOD IN NUMERICAL SIMULATION OF SECOND ORDER ENHANCED OIL PRODUCTION∗

(袁让)(程爱杰)(半丹平)(李长峰),2

1.Institute of Mathematics,Shandong University,Jinan 250100,China

2.School of Economics,Shandong University,Jinan 250100,China

E-mail:yryuan@sdu.edu.cn;aijie@sdu.edu.cn;dpyang@math.ecnu.edu.cn;cfi@sdu.edu.can

THEORY AND APPLICATION OF FRACTIONAL STEP CHARACTERISTIC FINITE DIFFERENCE METHOD IN NUMERICAL SIMULATION OF SECOND ORDER ENHANCED OIL PRODUCTION∗

Yirang YUAN(袁益让)1Aijie CHENG(程爱杰)1Danping YANG(半丹平)1Changfeng LI(李长峰)1,2

1.Institute of Mathematics,Shandong University,Jinan 250100,China

2.School of Economics,Shandong University,Jinan 250100,China

E-mail:yryuan@sdu.edu.cn;aijie@sdu.edu.cn;dpyang@math.ecnu.edu.cn;cfi@sdu.edu.can

A kind of second-order implicit fractional step characteristic fnite diference method is presented in this paper for the numerically simulation coupled system of enhanced (chemical)oil production in porous media.Some techniques,such as the calculus of variations, energy analysis method,commutativity of the products of diference operators,decomposition of high-order diference operators and the theory of a priori estimates are introduced and an optimal order error estimates in l2norm is derived.This method has been applied successfully to the numerical simulation of enhanced oil production in actual oilfelds,and the simulation results are quite interesting and satisfactory.

enhanced(chemical)oil production;three-dimensional porous coupled system; second-order implicit characteristic fractional step diferences;optimal order l2estimate;application in actual oilfelds

2010 MR Subject Classifcation65M06;65N30;76M25;76S05

1 Introduction

A mass of residual crude oil remains in the reservoir after water-fooding exploiting because of the constraint of capillary force preventing the motion and the slight infuenced regions and the fuidity ratio between displacement phase and driven phase weakening the displacement force.Then it is more important to develop the displacement efciency.A popular method is considered to add some chemical addition agents such as polymer,surfactant and alkali into the injected mixture,which can improve the fooding efciency.The polymer can optimize the fuidity of displacement phase,modify the ratio with respect to driven phases,balance the leading edges well,weaken the inner porous layer,and increase the efciency of displacementand the pressure gradient.Surfactant and alkali can decrease interfacial tensions of diferent phases,then make the bound oil move and gather[1-3]1Institute of Mathematics,Shandong University,Exploration Institute of Daqing Petroleum Administration. Research and application of the polymer fooding software(summary of“Eighth-Five”national key science and technology program(Grant No.85-203-01-08),1995.10.2Institute of Mathematics,Shandong University,Exploration and development Institute of Daqing Petroleum Limited Liability Corporation.Modifcation of solving mathematical models of the polymer and improvement of reservoir description(Grant No.DQYT-1201002-2006-JS-6765),2008.10.3Institute of Mathematics,Shandong University,Shengli Oilfeld Branch,China Petroleum&Chemical Corporation.Research on key technology of high temperature and high salinity chemical agent displacement, Chapter 4,§4.1 Numerical method,83-106.2011.3..

This paper discusses a second-order fractional step characteristic fnite diference method for numerical simulation of enhanced(chemical)oil production in porous media,and gives the theoretical analysis.Based on the former mathematical and mechanical theory,the software is accomplished,applied in national major oilfelds such as Daqing Oilfeld and Shengli Oilfeld and gives rise to important benefts and social value.

The mathematical model is a nonlinear coupled system with initial and boundary values [1-3]1,2,3:

where Ω is a bounded domain,a(c)=a(X,c)=k(X)µ(c)-1,d(c)=d(X,c),and other notations are interpreted as follows.φ(X)denotes the porosity of rock,k(X)denotes the permeability, µ(c)means the viscosity of fuid,and both D=D(X)and Kα=Kα(X)(α=1,2,···,nc)denote the difusion coefcients.uis Darcy velocity,p=p(X,t)is the pressure,c=c(X,t)means the saturation of water,and sα=sα(X,t)denotes the concentrations of diferent components. The components denote sorts of chemical agents such as the polymer,surfactant,alkali and other ions,and the number is denoted by nc.

The boundary value conditions for no permeation case:

where γ denotes the outer normal unit vector.

Initial conditions:

It is easy to compute the concentration by rewriting(1.3)as the following expression

Under an assumption of periodic condition,Douglas,Ewing,Wheeler,Russell and other scholars presented fnite diference method and fnite element method to analyze a type of two dimensional incompressible two-phase displacement problems and gave theoretical error estimates[4-7].A combination method is discussed by combining the characteristic method with normal fnite diference method or with normal fnite element method,which can refect the hyperbolic nature of one-order part of convection-difusion equations and decrease truncation error.This combination technique can also overcome numerical oscillation and dispersion, and can improve greatly the computational stability and accuracy.Douglas et al.presented mathematical model of slight compression,numerical method and theoretical analysis for twodimensional compressible displacement problem under periodic condition and began a modern numerical model research[8-11].The authors dropped the periodic condition,gave a new modifed characteristic fnite diference algorithm and fnite element algorithm,and derived optimal order error estimates in L2-norm[12-17].It is almost impossible to use traditional methods to accomplish numerical simulation of modern oilfelds exploration and production in terms of the following reasons.The computation runs on huge-scaled and three-dimensional region and the simulation consists of tens or hundreds of thousand nodes on a long period.So a novel numerical technique,fractional step method,is introduced in this paper.Douglas,Peaceman applied upwind method in incompressible two-phase(water and oil)displacements successfully [18].However it is hard naturally to proceed theoretical analysis.The stability and convergence are derived by Fourier method only for constant coefcient case and this analysis is not generalized for variable coefcient equations[19-21].Considering actual application,numerical stability and accuracy,the authors present one second-order fractional step characteristic fnite diference method for three-dimensional compressible two-phase displacement coupled problem of enhanced oil production.This algorithm can overcome numerical oscillation and dispersion, and decrease the computational scale by decomposing three-dimensional problem into three successive one-dimensional subproblems.Using the calculus of variation,energy analysis method, commutativity of the products of diference operators,decomposition of high-order diference operators and the theory of a priori estimates,the authors give the second-order convergence result of accuracy and error estimates in l2-norm,and successfully solve the international problem of Douglas and Ewing[8,11].The method discussed in this paper has been applied in numerical simulation of enhanced oil production and gives a fundamental research in energy mathematics1,2,3.

In general,the coefcients in the problem are supposed to be positive defnite,

where a∗,a∗,d∗,d∗,φ∗,φ∗,D∗,D∗K∗,K∗and A∗are positive constants.The coefcientsd(c),b(c),g(c)and Qα(c,sα)are Lipschitz continuous respect to the unknown functions nearby their ε0neighbors.

Exact solutions of(1.1)-(1.7)are assumed to be suitably smooth,

where τ denotes a vector along the characteristics.In this paper M and ε express general positive constant and general positive small constant,respectively,and they have diferent meanings in diferent areas.

2 Second-order Implicit Fractional Step Characteristic Finite Diference Method

The fractional step algorithm of fow equation(1.1)is given by

Since the fuid moves along the characteristic direction of transport,the one-order hyperbolic term of saturation equation(1.2)is treated by characteristics method.This method of characteristics has high order of accuracy,strong stability and takes large time step[4-6].Let ψ(X,u)=?φ2(X)+|u|2?1/2,∂/∂τ=φ-1{φ∂/∂t+u·∇},and apply backward diference to approximate the derivative along the characteristics,

The implicit fractional step characteristic fnite diference method of saturation equation (1.2)is considered.

An implicit fractional step characteristic fnite diference method of component concentration equation(1.6)runs in parallel,

Initial approximation:

3 Convergence Analysis

Let π=p-P,ξ=c-C,ζα=sα-Sα,where p,c and sαare exact solutions of(1.1)-(1.7), and P,C and Sαare numerical solutions of(2.3)-(2.6).Introduce the inner products and norms in discrete space l2(Ω)[22,23],

The square of weighted semi-norm in discrete space h1(Ω)corresponding to H1(Ω)=W1,2(Ω) is denoted by〈D∇hf,∇hf〉,where the weighted function D(X)is positive defnite.

The diference algorithm(2.3)-(2.6)are applied layer by layer,then error estimates hold

ProofThe fow equation is considered frst.Cancelling the transitional solutions Pn+1/3and Pn+2/3by(2.3a),(2.3b)and(2.3c),

Subtracting the diference equation(3.3)from the fow equation(1.1)at the point(xi,yj,zk, tn+1),it gives rise to the error equation of the pressure,

An induction hypothesis is given by

Continue to estimate the terms on the right-hand side of(3.6),

The third term on the right side of(3.6)is considered.The frst part is discussed here

The operators-δ¯x1(Anδx1),-δ¯x2(Anδx2),···are self-conjugate,positive defnite,bounded and the domain is a unit cube,but their products are not commutative generally.Noting that the diference operators δx1δx2=δx2δx1,δ¯x1δx2=δx2δ¯x1,δx1δ¯x2=δ¯x2δx1,δ¯x1δ¯x2=δ¯x2δ¯x1are commutative,the frst term of the right side of(3.7c)is written by

Consider the third term on the right-hand side of(3.8),

Consider the rest term on the right-hand side of(3.7c)similarly,

The estimate(3.11)can be obtained identically for the other two terms of the third right term of(3.6).

For the fourth term on the right side of(3.6),

Collecting(3.7)-(3.12),it holds for error equation(3.6)as Δt and ε are sufciently small,

Error equation of the saturation is obtained by(1.2)(t=tn+1)and(3.14)

It follows from(3.15)by using the restriction condition(3.1)and hypothesis induction(3.5),

Rewrite the above expression,

and the fact that the velocityUnis bounded because of the inequality|Un|∞≤M{1+|∇hπ|∞} and(3.5),and use the restriction(3.1),we have

Consider the fourth term on the right-hand side of(3.18),

Noting the smoothness of φ and D,we have

It is easy to obtain a result for the ffth term of the right side of(3.18),

Applying(35)-(38)to estimate(3.18),we have

Summing up the estimates of the pressure for n from 0 to L and noting that π0=0,

The frst term of the right side of(3.40)is estimated as follows

Combining(3.40)with(3.41),

Summing up(3.39)for n from 0 to L,and noting that ξ0=0 and the expression(3.42),

For π0=ξ0=0,

Collecting(3.42)and(3.43),

Using Gronwall lemma,

Error equation of component concentration is derived by(1.6)(t=tn+1)and(3.46),

In numerical analysis there exists bound water in oil reservoir,that is to say there exists a positive constant c∗such that c(X,t)≥c∗＞0.It holds as h and Δt sufciently small becauseof the convergence result(3.45),of c(X,t),

The terms on the left-hand side of(3.49)are estimated as follows,

The terms on the right-hand side of(3.49)are argued later.Using(3.45)we have

For the eighth,the ninth and the tenth terms,we have

Applying(3.50)-(3.52)on both sides of(3.49),we get

For the frst term on the right-hand side of the above expression,

Consequently,

Using Gronwall lemma,we get

The proof ends.

4 Applications

The implicit fractional step characteristic fnite diference method has been applied successfully in software design of enhanced oil production of the polymer fooding and numerical simulation and analysis of actual oil production in Daqing Oilfeld.The mathematical model is formulated as follows[1-3]1,2,3:

Let“w”and“o”denote respectively water phase and oil phase and let φ mean the porosity.The signs of the l-th phase are defned as follows.pldenotes the pressure,slmeans the saturation, Blmeans the volume factor,λldenotes the fuidity,γlis the proportion,qlis the source sink term,and pcis the pressure of capillary.

The mathematical model of the polymer and the motion of kation and anion components is described by a system of convection-difusion equations,

where cα=cα(x,t)(α=1,2,···,nc)denotes the concentration of components,and ncis the number of components.The simplifed model of(4.1)and(4.2)is turned into the system of (1.1)-(1.7)[1,2,24].

Experimental tests for Daqing Oilfeld(Xing Fourth Zone of the polymer development area)are discussed by using the polymer fooding software.The three-dimensional region of geological models are decomposed into 46×83×7 subdomains,and the efective thicknesses are distributed in Fig.1.The sketch of instantaneous oil production and the curve of water moisture are shown in Fig.2 and Fig.3.Numerical data of matter balance are analyzed in Table 1.From this it is easy to see that the numerical method and the software can keep high order of accuracy and numerical simulation can refect correctly the physical process and the principle of polymer fooding.Main physical quantities distribute reasonably,computational accuracy is high,and some results of the polymer such as accumulation,endless loop do not arise.

Table 1Analysis of numerical data of matter balance

AcknowledgementsThe authors express their deep appreciation to Prof.J.Douglas Jr,Prof.Ewing R E and Prof.Jiang L S for their many helpful suggestions in the serial of research on numerical simulation of enhanced(chemical)oil production.

[1]Ewing R E,Yuan Y R,Li G.Finite element for ehemical-fooding simulation//Proceeding of the 7th International Conference Finite Element Method in Flow Problems.The University of Alabama in Huntsville, Huntsville,Alabama:UAHDRESS,1989:1264-1271

[2]Yuan Y R.Theory and application of numerical simulation of energy sources,basis of numerical simulation of chemical production(tertiary oil recovery),Chapter 3.Beijing:Science Press,2013:257-304

[3]Yuan Y R,Yang D P,Qi L Q,et al.Research on algorithms of applied software of the polymer//Gang Qinlin.Proceedings on Chemical Flooding.Beijing:Petroleum Industry Press,1998:246-253

[4]Douglas Jr J,Russell T F.Numerical method for convection-dominated difusion problems based on combining the method of characteristics with fnite element or fnite diference procedures.SIAM J Numer Anal,1982,19(5):871-885

[5]Douglas Jr J.Simulation of miscible displacement in porous media by a modifed method of characteristic procedure//Lecture Notes in Mathematics 912.Numerical Analysis,Proceedings,Dundee,1981.Springer, 1982:64-70

[6]Douglas Jr J.Finite diference methods for two-phase incompressible fow in porous media.SIAM J Numer Anal,1983,20(4):681-696

[7]Ewing R E,Russell T F,Wheeler M F.Convergence analysis of an approximation of miscible displacement in porous media by mixed fnite elements and a modifed method of characteristics.Comp Meth Appl Mech Eng,1984,47(1/2):73-92

[8]Douglas Jr J,Roberts J E.Numerical method for a model for compressible miscible displacement in porous media.Math Comp,1983,4(164):441-459

[9]Yuan Y R.Time stepping along characteristics for the fnite element approximation of compressible miscible displacement in porous media.Math Numer Sinica,1992,14(4):385-400

[10]Yuan Y R.Finite diference methods for a compressible miscible displacement problem in porous media. Math Numer Sinica,1993,15(1):16-28

[11]Ewing R E.The Mathematics of Reservior Simulation.Philadelphia:SIAM,1983

[12]Yuan Y R.Characteristic fnite diference methods for moving boundary value problem of numerical simulation of oil deposit.Science in China,1994,37A(12):1442-1453

[13]Yuan Yi-rang,The characteristic mixed-fnite element method and analysis for three-dimensional moving boundary value problem.Sci China Math,1996,39(3):276-288

[14]Yuan Y R.Finite diference method and analysis for three-dimensional semiconductor device of heat conduction.Sci China Math,1996,39(11):1140-1151

[15]Axelsson O,Gustafasson I.A modifed upwind scheme for convective transport equations and the use of a conjugate gradient method for the solution of non-symmetric systems of equations.J Inst Maths Applics, 1979,23:321-337

[16]Ewing R E,Lazarov R D,Vassilevski A T.Finite diference shceme for parabolic problems on composite grids with refnement in time and space.SIAM J Numer Anal,1994,31(6):1605-1622

[17]Lazarov R D,Mishev I D,Vassilevski P S.Finite volume method for convection-difusion problems.SIAM J Numer Anal,1996,33(1):31-55

[18]Peaceman D W.Fundamantal of Numerical Reservoir Simulation.Amsterdam:Elsevier,1980

[19]Douglas Jr J,Gunn J E.Two order correct diference analogues for the equation of multidimensional heat fow.Math Comp,1963,17(81):71-80

[20]Douglas Jr J,Gunn J E.A general formulation of alternating direction methods,Part 1.Parabolic and hyperbolic problems.Numer Math,1964,6(5):428-453

[21]Marchuk G I.Splitting and alternating direction method//Ciarlet P G,Lions J L,eds.Handbook of Numerical Analysis.Paris:Elsevior Science Publishers BV,1990.197-460

[22]Yuan Y R.Theory and application of upwind fnite diference method for moving boundary value problem of three-dimensional percolation coupled system.Sci China Math,2010,40(2):103-126

[23]Yuan Y R.The second-order upwind fnite diference fractional step method for moving boundary value problem of nonlinear percolation coupled system.Sci China Math,2012,42(8):84-864

[24]Shen P P,Liu M X,Tang L.Mathematical problems of petroleum exploration and development:Mathematical problems of oil-gas felds development,Part III.Beijing:Science Press,2002:197-264

∗Received May 8,2014;revised February 2,2015.This work is supported by the Major State Basic Research Development Program of China(G19990328),National Tackling Key Program(2011ZX05011-004, 2011ZX05052,20050200069),National Natural Science Foundation of China(11101244,11271231,10771124, 10372052)and Doctorate Foundation of the Ministry of Education of China(20030422047).