Numerical Simulation of Three-Dimensional Enhanced Oil Recovery and Its Application for Horizontal Wells

In this paper, a numerical approximation is put forward to solve the simulation of three-phase (water, oil and gas) problem of permeation fluid mechanics, where the flow is driven to move by injecting ternary compound flooding in porous media. After the mechanical model formulated, then different computational algorithms are structured including a fully implicit procedure, an implicit procedure for the pressure and an implicit/explicit program respective for the pressure and the concentration. Then the pressure, the saturation and the concentrations of chemical components and petroleum acid are computed by using an upstream sequence, a difference iterative algorithm of implicit upwind fractional step and interpolation calculations of relative permeability. The software is completed and is generalized in actual industrial applications, such as in analysis and simulation of national major oil-fields, which gives the illustration with high order of accuracy, runs on a combination partition of ten-meter step and of one hundred and fifty thousand nodes and simulates on the interval of tens of years. This gives rise to outstanding economic benefits and social benefits in China. Numerical analysis and convergence argument in l norm are given for a simplified model in consideration of capillary force and an efficient tool is presented to solve an international famous problem formulated by Jr. Douglas, J.


Introduction
A mass of residual crude oil remains in the reservoir after water-driving exploiting, because the capillary force prevents the remaining oil moving, or the undesired fluidity ratio of driving phase and driven phase makes the injected water influence the flow slightly.Then a powerful tool is given to improve the displacing oil efficiency by adding some chemical agents such as the polymer, surfactant and alkali.The polymer can optimize the fluidity of driving phase, modify the ratio comparing with driven phase, and make the fluid distributing uniformly on the leading edges.Then it can weaken the motion of inner porous layers, develop the displacement efficiency and increase the gradient of pressure.Surfactant and alkali can decrease interfacial tensions of different phases, then can make the bounded oil move and gather.
The displacement of three-dimensional multi-phase (water, oil and gas), multi-component and slight compressible problem is interpreted by a nonlinear system of partial differential equations, whose solutions it is hard and complicated to find.Therefore, many modern numerical methods and techniques are discussed to approximate exact solutions.Generally speaking, the pressure is solved by an implicit scheme and the saturation is obtained implicitly or explicitly.Then numerical approximations are given and they can interpret the physical natures.During numerical simulation it is important to analyze numerical values and information and it is most helpful to describe the whole process of chemical flooding quite well.Meanwhile, it helps us understand the principles and processes of different displacements better and forecast the recovery efficiency of natural oil.We can compute the percentage of oil in producing liquid and the percentage of polymer, surfactant and alkali, then conclude the changes of different components and determine the end of injecting liquid.Finally, we can get how different parameters affect the oil recovery, and we can forecast some natures of actual applications, and optimize different recovery plans.Petroleum engineers and mathematicians pay more attention to modern new techniques and their application in oil recovery, mainly including the formulation of the mathematical model, the designing of computer software and the achievement of numerical simulation.
Yuan was engaged in theory and application of chemical flooding numerical simulation with Prof. Ewing from 1985Ewing from to 1988 in US (Ewing, 1989;Yuan, 2013).From 1991 to 1995, the research group of Yuan undertook some major national projects such as "Eighth-Five" National Key Science and Technology Tackling Program  entitled "Research on designing and application of the polymer displacement software" (Yuan, 1992(Yuan, , 1993(Yuan, , 1994))Note 1.The software was applied successfully in Daqing Oilfield such as plan design and discussion of industrial production region.From these actual applications, some conclusions about the polymer follow that the slug affects greatly, the bar is set to protect water and the dosage becomes large.These results have given rise to outstanding economic benefits and social values (Yuan, et al., 1998;Yuan, 1998Yuan, , 2000Yuan, , 2001Yuan, , 2002Yuan, , 2003;;Ewing, 1983)Note 2. Later the research group undertook a key tackling program of oil administration of Daqing Oilfield (Grant No. DQYJ-1201002-2006-JS-9565) entitled "Solving development of mathematical models and completing explain of reservoir" Note 3.This software system was also applied in numerical simulation of the polymer displacement in Zaobei fault block of Dagang Oilfield, optimization of designing plan of expanded test area of ternary compound flooding in Gudong Little Well test area of Shengli Oilfield, polymer flooding in Gudong Middle One test area, Gudong West region and feasibility of active water flooding of Gudong Eighth region, and many interesting results are obtained Note 4. In recent years the research group finished the successive key tackling project of Daqing Petroleum Administration Bureau (Grant No. DQYJ-1201002-2009-JS-1077), entitled "Research on alkali flooding principle model of chemical displacement simulator and horizontal wells model structuring and solving methods" Note 5 and presents many important theoretical results and applications.Some conclusions and extensions are discussed in this paper, mainly consisting of permeation fluid mechanical mathematical model, numerical methods, applicable software structuring and numerical simulation, and applications of horizontal wells in actual oilfields.

Permeation Fluid Mechanical Model of Chemical Flooding
The following considerations are argued in this paper for three-phase (oil, gas and water) numerical simulation.Oil phase consists of oil component and solution-gas component, water phase only has water component and gas phase only has gas component.Gas-component interchanges between oil and gas as the pressure becomes changing.The three-phase model is improved to simulate the chemical displacement as well as possible by adding chemical compound combination displacement module together with three-phase module and chemical reaction balance equations.A type of new computational software is accomplished for three-phase coupled problem with chemical compound flooding and is used successfully in numerical simulation of horizontal wells in actual oilfields.The additional chemical agents include polymer, surfactant and petroleum acid and so on, where petroleum acid exists in both water and oil and the other components only exist in water.The software is made up of three modules: three-phase solver, component equations solver and chemical balance equation solver.
Lots of subblock codes are given similarly to the program of black-oil model.The water viscosity is a constant in the black oil problem while it becomes a variable in chemical displacement due to the appearance of chemical components.On the other hand, the program structure and algorithm design for component equations and chemical reaction balance equation should be compatible with black oil model.Then additional modifications are given to finish actual computations on wedge-out district, fault and the parameter of edge-bottom water.
Two different algorithms are stated as follows.

(I)
Full-implicit algorithm (FIA).FIA, one of the most stable finite difference schemes, can obtain the values of the pressure, water saturation, gas saturation, and the ratio of solution-gas and oil implicitly.The computational time is spent mainly on a time step to complete extrapolation iterations.For a full implicit scheme is more stable than other explicit schemes invalid sometimes due to strong stability condition, so FIA adopts a large time step to decrease the total simulation cost.
(II) Implicit/explicit algorithm for the pressure and saturation (IEAFPS).This algorithm is applied to compute the pressure implicitly and get the saturations explicitly, structured on the following assumptions.The saturation does not affect distinctly the motion of the fluid.Under this hypothesis the saturation can be canceled from the discrete flow equation, and the values of the pressure are computed implicitly only at different iterations.The values of saturation are updated point by point explicitly once the pressure is given.The IEAFPS becomes unstable as the saturation changes too rapidly with from a time layer to another.While this implicit/explicit algorithm is really efficient as the change of saturation is controlled usually by 5% by taking time steps sufficiently small.
Let the variables appearing in oil-water displacement problems be denoted by the subscripts "w" and "o" respectively corresponding to the water phase and the oil phase, where the mathematical model of permeation fluid mechanics is stated as follows (see Ewing, Yuan, ＆ Li, 1989;Yuan, 2013;Yuan, Yang, ＆ Qi, 1998;Yuan, 2000;Ewing, 1983): , , ; Similarly, let "w", "o" and "g" denote the functions respectively corresponding to the water phase, the oil phase and the gas phase, and the mathematical model of three-phase displacement is followed ①②③④⑤ , ( (2.2b) (2.2c) (2.2d) In the above expressions, denotes the porosity, and means the ratio of solution-gas and oil.With respect to l-phase, the functions p l , c l , T, A l , T rl , , and denote the pressure, the concentration, the absolute permeability, the volume factor, the relative permeability, the viscosity, the density, and the source sink term (floor condition), respectively.
The water viscosity is a constant in the black oil model, while it is a function with respect to the density in three-phase chemical flooding model, , where denotes the concentration of polymer relative to water.The components move in water, and their concentration affects water-phase viscosity field, then influences the concurrent motion of three-phase fluid.Therefore, a new nonlinear system is derived naturally from the mathematical model of three-phase displacement with a convection-diffusion equation of chemical displacement.A coupled algorithm is structured to solve the equations of three-phase flow and chemical agent displacement.During an increment of time, the equation of three-phase flow is solved firstly, and the flow field is obtained.Secondly, the convection diffusion equation and the concentration field of chemical agent displacement are computed.Finally, the viscosity field of water phase is updated.The computation goes on in a similar way.
The motion laws of the components (the polymer, anions and cations) are described by convection-diffusion equations.For convenience, we don't pay more attention to the symbols of different components, so we give a mathematical description of the concentration for a certain component. (2.2e) Surface active agent is used to affect interfacial tension, and it makes capillary number increasing, then it can decrease the saturability of residual crude oil.The permeability is updated by an interpolation of the values of relative permeability of lower capillary number and those of the upper capillary number.
The primary principle of alkaline flooding is stated as follows.The alkali is injected into the fluid and reacts with petroleum acid, and surface active agents are formed.These agents can decrease interfacial tension of the fluid and make remaining oil move.Petroleum acid components exist in both water and oil, so mass transfer takes place between different phases.The computational program of three-phase chemical agent displacement runs as follows.
Step 1.At t=t 1 , the solver for the pressure and saturation runs, then the solver for the concentrations of different chemical components runs.
Step 2. The viscosity of water dependent on the concentration is modified.
Step 3. At t=t 2 , the solver for the pressure and saturation runs, then the solver for the concentrations of different chemical components run.
Step 4. The viscosity of water dependent on the concentration is corrected.

……
In the above way the pressure and saturation at t=t n+1 are obtained, the concentration of different chemical components is computed, then the viscosity of water is corrected.

……
The program ends.

Numerical Methods
Two different numerical methods are argued in this section.

Full-Implicit Algorithm of Three-Phase Flow
Cancelling spare unknown variables and formulating the equations of the pressure of oil, the concentration of water and the concentration of gas, we can present the implicit scheme.All the values on the left-hand side of the discrete scheme including the pressure, the saturation, the quantity, the relative permeability, the capillary force and other parameters are updated by the latest numerical data.The resulting implicit difference scheme is a nonlinear algebraic system, and it is solved by the method of iteration.Its computation time cost of iteration (outer iteration) is seven times that of the implicit/explicit method.The implicit scheme is easily applied to solve some difficult and complicated problems such as black-oil numerical simulation because it is unconditionally stable.In this paper let and denote two difference operators corresponding to usual time discretization A full implicit finite difference scheme for three-phase model is structured as follows (2.4) where and .Using an operator , Formulating the above expression in an expansion form, and omitting quadratic terms, then we give the remainder after the k-th iteration as follows, Continue to state the above equation with a remainder term, (2.5) The iterations are convergent as for l=w,o,g and k=1,2,….A resulting formulation is First, a linear expansion of (2.6) is given., and are chosen as unknown variables and the terms on the right-hand side are discussed later.

Notice that
, and let (2.7) then, we have Similarly, for the term on the right hand side of oil equation, and for the terms on the right-hand side of gas equation, where (2.9) The symbols and mean the derivative of volume factor and the derivative of porosity with respect to the pressure.and denote the derivatives of p c with respect to C g and C w.
Two operators are introduced to discuss the terms on the left-hand side, Then, (2.10) Only considering an expansion of M l , we have three expressions respective for water, oil and gas, and The values of conductivity coefficient in second-order difference operators are assigned by the upstream rule, and let i + and i -denote the upstream nodes between the i-th node and the (i+1)-th node and between the i-th node and the (i-1)-th node, respectively, An algebraic system is obtained by substituting the terms on both sides into initial difference equations.

Implicit/Explicit Algorithm for the Pressure and the Saturation
The implicit/explicit method is derived from a resulting equation only dependent on the pressure by combining all the flow equations, and the values of the saturation are updated explicitly by the pressure at a time level.
A discretization of (2.2) is formulated with respect to p o and the saturation, A primary assumption is given to derive the implicit/explicit scheme that the capillary pressure of the left-hand side is a constant number during a time step computation.
The coefficient S is defined by (2.12) Considering three equations of (2.11) and canceling all the terms of by multiplying water equation by A, testing gas equation by B, and adding three expressions together, we get the resulting right term The numbers A and Bare defined by (2.13) and they satisfy Therefore, the pressure equation turns into (2.14)A normal finite difference scheme is derived from a parabolic equation, whose matrix form is where K denotes a tridiagonal matrix and D is a diagonal matrix.The vector G is dependent on the gravity and capillary pressure.
If the pressure is given, and is substituted in the former two equations of (2.11), then the saturation is computed explicitly.Then the capillary pressure and are considered and used explicitly at next time level.

Numerical Method for Component Concentration Equations
The combination components including surfactant, polymer, alkali and sorts of ions are injected in water, whose physical conservation of mass is described by a convection-dominated diffusion equation.An operator division technique is applied to compute the component equation, which has great efficiency and high order of accuracy.Two simplified equations are derived from the component equation.One is a hyperbolic equation only dependent on the diffusion term and another is a parabolic equation only dependent on the dispersion term.We use an implicit upwind scheme to compute the former, where the values are obtained point by point explicitly by applying an upstream rule.The latter is solved by an alternating direction finite difference method, which can accelerate the computation speed and can improve the efficiency.For convenience, the concentration equation of components is rewritten in a typical convection-diffusion form, ).
where dispersion tensor is a diagonal matrix.When the saturation C w and the flow velocity field v w of water at t n+1 are given, the task is to compute the concentration S n+1 at the first direction, then at the second direction and the third direction alternatively.C denotes the concentration of a certain component ignoring the subscript k.A convection equation is argued by an implicit upwind scheme, (2.16) where v w =(v wx ,v wy ,v wz ) T is a velocity vector.GivenS n+1,0 , the dispersion equation is discussed alternatively in three directions.First considering the computation in x-direction, technique is introduced in actual applications to overcome some difficulties and effects of grids orientation.For example, it is considered how to keep the symmetrical computation consistently with symmetric physical problems.The values are obtained in two continuous processes, x-direction first then y-direction, and y-direction first then x-direction, then the algorithm continues in z-direction by using an average value of the above computational data.This efficient method is applied in present software.

Newton-Raphson Iterations of Chemical Reaction Balance Equations
The chemical balance equations, a nonlinear system consisting of liquid chemical agent, solid chemical agent and ions adsorbed in rocks, are solved by Newton-Raphson iterations.

Interpolation Computation of Relative Permeability
The relative permeability of low capillary number is defined by a table and the values of high capillary number are obtained by a formula, as in some curves of and .Interpolation computation is used here in the following statements.
(I) Firstly, as in Fig. 1, by geometric transformation such as stretching (or compression), the desired interval of relative permeability, is given under present capillary number by starting with (or ) of low capillary number (or high capillary number).Meanwhile, two pairs of curves of (or ) are sketched on the interval , whose values are computed by the following formulas. (2.18a) The symbols , and denote independent variables varying respectively on the intervals , and .
(II) Secondly, the graphs of relative permeability of present capillary number are illustrated in Fig. 2 by interpolation computation of above values and . (2.20a) (2.20b) In actual computations, the computation is carried out generally in a reverse order as arrow illustration in Fig. 2 and , are computed by and the formulas in (I). (2.21a) (2.21b) Then, The derivatives of relative permeability are obtained by the following expressions, (2.23a) (2.23b)

Functional Test I of Horizontal Wells
An experimental test is argued to testify the effects of horizontal wells.The computational domain is partitioned by 9×9×1, the computational time length is 5500 days and the polymer is injected at the 1460th day.The test proceeds in the following three plans.The first plan is made for five vertical wells (four injected wells and one production well), the second plan is for four injected wells and one production well (horizontal well) and the third is for three horizontal wells (two injected wells and one production well).
Numerical data of three different programs in view of streamline, saturation of the polymer and viscosity at the 3400th day are illustrated in the following figures.From these figures it is concluded that the software can simulate the process of horizontal well reservoir with high order of accuracy.Numerical simulation is well consistent with the function and the principle of horizontal wells.The values of primary physical functions distribute reasonably, some unexpected troubles such as the polymer tacking and endless loop never occur in this simulation.The numerical isograms of the streamline, saturation, phase concentration of the polymer and viscosity during three thousand and four hundred days simulation are illustrated by the figures (see Fig. ) ,

Functional Test II of Horizontal Wells (Injected Wells are Horizontal)
The displacement of polymer is illustrated on the partition 48×43×7, and computational data after 3400 days simulation are shown in the following figures (see the isograms of the pressure, the saturation, the concentration of polymer and the viscosity in Fig. 7a, Fig. 7c, Fig. 7d and Fig. 7e, the streamline in Fig. 7b).

Functional Test III of Horizontal Wells (Horizontal Wells are Production Wells)
The displacement of polymer is illustrated on the partition 48×43×7, and computational data after 3400 days simulation are shown in the following figures (see the isograms of the pressure, the saturation, the concentration of polymer and the viscosity in Fig. 8a, Fig. 8c, Fig. 8d and Fig. 8e, the streamline in Fig. 8b).

Functional Test IV of Horizontal Wells (Injected Wells are Horizontal Wells)
The displacement of polymer is illustrated on the partition 149×149×1, and streamlines of the pressure and distributions of the saturation after 3400 days simulation are shown in the following figures (see Fig. 9a and Fig. 9b).

Mechanism Analysis of Alkaline Flooding
The partitions are defined by 9×9×1 and it is tested how different concentrations of injected alkali and different naphthenic acid values affect the displacement.Three SLUGs are made and the total simulation time, 5500 days, is divided into three processes.The first process, water drive takes about 1460 days, then the polymer or the combination of the polymer and alkali flooding is injected, finally water is injected 2920 days later.
Five displacement ways are discussed here.The symbol P denotes the polymer flooding.Let ASP denote the compound combination of polymer and alkali, where the acid value is 0.0006 and the concentrations of injected Na ions and CO 3 ions are 0.3351 and 0.3929, respectively.ASP3 denotes the compound combination flooding of the polymer and alkali, where the acid value is 0.006 and the concentrations of injected Na and CO 3 are 0.3351 and 0.3929.ASP4 denotes the compound combination flooding, where the acid value is 0.0006, and the concentrations of Na and CO 3 are 3.3351 and 3.3929.The fifth case ASP5 denotes the compound combination flooding, where the acid value is 0.006, and the concentrations of injected Na and CO 3 are 3.3351 and 3.3929.
Based on numerical data of produced surface active agents, the concentration of residual oil and the generated depositing under five different cases, there four conclusions are stated as follows.
(I) From the normal concentration values at the 3000th day under different cases, the produced surface active agents becomes more and more as the acid value and the concentration are larger and larger.
(II) From the values of residual oil at the 3000th day under different cases, the residual oil quantity grows smaller as the acid value and the concentration value of injected alkali are larger.
(III) From the values of CaCO 3 under the second and the fourth cases, the values of CaCO 3 are almost independent of the concentration of alkali.
(IV) From the values of Mg(OH) 2 under the second and the fourth cases, the concentration of injected alkali gives a powerful influence to Mg(OH) 2 especial to Mg(OH) 2 of injected wells.
The numerical data of moisture content, instantaneous oil production and accumulative oil production under five different cases are compared in the following figures (see Fig. 10a-Fig.10c).

Numerical Analysis of the Model
For convenience, a simplified model is considered in this section, and its convergence analysis proceeds naturally for numerical simulation of three-dimensional problem of black-oil and chemical agent displacement in porous media.Three dimensional mathematical model of two-phase incompressible displacement problem with components in consideration of capillary force is given by a nonlinear system of partial differential equations with initial value and boundary value conditions (Ewing, Yuan & Li, 1989;Yuan, 1993Yuan, , 1993Yuan, , 1994Yuan, , 2000Yuan, , 2013;;Yuan, Yang & Qi, 1998). (3.1) where is a bounded computational domain.Let the subscripts "o" and "w" respectively denote the parameters of oil and water.The notations , , , and denote the concentration, the pressure, the relative permeability, the viscosity and the output quantity with respect to the l-phase, respectively.means the rock porosity, means the absolute permeability and denotes the component concentration.The components denote sorts of chemical agents such as the polymer, surfactant, alkali and other ions, and the number is denoted by n c .v is Darcy velocity, is diffusion coefficient, and is source sink term related with the output.The rock void space is assumed to be full of water and oil, .The capillary ( ) ], force is dependent of the concentration , , where .
The models (3.1) and (3.2) should be turned into a normal form (Ewing, 1989;Yuan, 2013).Let denote the total migration rate of two-phase fluid, and denote relative migration rate.Applying Chavent transformation (Ewing, 1989;Yuan, 2013) (3.4) The flow equation is derived from the sum of (3.1) and (3.2), (3.5) where q=q o +q w .The concentration equation is derived from the difference of (3.1) and (3.2),where q w =q and q o =0 as q is greater than 0 and equal to 0 (water injection well), and q w =λ w q and q o =λ o q as q is less than 0 (oil production well).
where denotes the outer normal unit vector.An additional condition should be introduced to determine the pressure p in consideration of no permeation case, and it is defined by The compatibility condition is The initial value condition is (3.13a) (3.13b) Douglas and other scholars present primary and foundational work (see Douglas, 1981Douglas, , 1983;;Douglas & Russell, 1982;Douglas, Russell & Wheeler, 1984;Douglas & Roberts, 1983;Yuan, 2001Yuan, , 2013) ) to discuss numerical simulation of flow displacement problems.In numerical simulation of modern oilfield exploration and development, in consideration of a large geometrical region and a long time period, the computation is greatly complicated, its scale is huge, and the program is designed on tens of thousands nodes or tens of millions nodes.Therefore, common numerical methods are unable to solve the problems of this type and a kind of alternating direction difference method is put forward by Peaceman and Douglas (Peaceman, 1980;Douglas & Gun, 1963;Marchuk, 1990;Yanenke, 1967).But there exist many natural difficulties in theoretical analysis.They use Fourier method to analyze the stability and the convergence of constant coefficient problems but cannot extend this technique for variable coefficient problems (Douglas & Gunn, 1963, 1964).Considering actual application, numerical stability and accuracy, Yuan gives fractional step schemes on upwind finite difference method and the method of characteristics and their convergence analysis for three-dimensional two-phase displacement coupled problem (Yuan, 1999(Yuan, , 2001(Yuan, , 2003;;Yuan, Cheng, Yang& Li, 2014;Yuan, Li& Sun, 2014;Yuan, Li, Sun&Liu, 2014;Yuan, Cheng, Yang, Li& Liu, 2015).Then an implicit upwind difference scheme is discussed and an interesting development appears for black oil chemical flooding in consideration of capillary force.This algorithm can overcome numerical oscillation and dispersion, and decrease the computational scale by decomposing a three-dimensional problem into three successive one-dimensional subproblems.Using the calculus of variation, energy analysis method, commutativity of the products of difference operators, decomposition of high-order difference operators and the theory of a priori estimates, the authors give second-order error estimates in l 2 -norm.So the software applied in actual oilfield exploration and development is based on mathematical and mechanical analysis (Ewing, Yuan, 1989;Yuan, 2013;Ewing, 1983;Shen, Liu & Tang, 2002;Yuan, 2015).
To give a clear illustration of theoretical analysis, we take a simple computational domain .

Let and let
The notations are defined similarly.


In the above expression, Q h denotes the cube of side-length h centered at the origin.Darcy velocity is computed by (3.15) and by similar calculations of and .
Consider the implicit upwind fractional step scheme of the saturation (3.9).Given the values c n , it is necessary to find the values c n+1 at the (n+1)-th step.Using the approximation, replacing the differential operator by the difference quotient, , then we rewrite the saturation (3.9) in the following partition form,

Conclusion and Discussion
Theoretical analysis and computational method of numerical simulation and the application of horizontal wells of three-dimensional three-phase (water, oil and gas) percolation mechanics of chemical compound combination (the polymer, surface active agents and alkali) in porous media are discussed in this paper.Summary is first stated about our project.Mathematical model of permeation fluid mechanics is presented in § 1.A full implicit numerical scheme, an implicit/explicit algorithm for the pressure and the saturation, and an interpolation computation of relative permeability are given, and an upwind difference fractional step algorithm based on upstream sequence is structured.A type of software applicable in major industries has been accomplished, mostly carried out with the spacial step of ten-meters, one hundred and fifty thousand nodes and tens of years simulation period in §2.Some experimental tests are considered for horizontal wells model, some examples occurring successfully in national major oil fields such as Daqing Oilfield, are illustrated to show the principle of alkaline flooding in §3.Then numerical analysis of simulation method is accomplished for capillary force model in this section, which gives the software a mathematical and mechanical consideration.
boundary surface of .(II)The boundary value condition for no permeation case: