Numerical Simulation of Chemical Flooding and Its Application for Horizontal Wells

Numerical simulation of permeation fluid mechanics for three phase water, oil and gas-chemical compound combination flooding of the polymer, surface active agents and alkali in porous media is discussed in this paper. In view of petroleum geology, geochemistry, computational mechanics of flow and computer technology, a permeation fluid mechanics model of three phase chemical compound combination flooding is presented firstly, then a numerical algorithm consisting of a full implicit program, an implicit computation for the pressure and an implicit/explicit program respective for the pressure and the concentration is given by structuring an upstream sequence and an iterative algorithm of implicit fined upwind fractional step finite difference to solve the pressure equation, the saturation equation and the concentration of chemical substance components and the petroleum acid concentration equation. A type of high accuracy software applicable in major industries is made on ten-meter steps, one hundred and fifty thousand nodes and tens of years and has been carried out successfully in analysis and simulation of national major oil-fields extraction such as Daqing Oilfield, which gives rise to outstanding economic and social benefits. A precise analysis is given for a simplified model and an idea is presented to solve this international famous problem.


Introduction
There remains plenty of residual crude oil in the reservoir after water-flooding exploiting because the constraint of capillary force weakens the motion and the volume of influenced regions is small due to the disadvantageous fluidity ratio between displacement phase and driven phase.Then it is more important to develop the displacement efficiency.A popular method is considered that the mixture is injected into the underground fluid including chemical addition agents such as polymer, surface active agent and alkali.The polymer can optimize the fluidity of displacement phase, modify the ratio with respect to driven phases, balance the leading edges well, weaken the inner porous layer flow, improve the efficiency of displacement phase and increase the pressure gradient.Surface active agent and alkali can decrease interfacial tensions of different phases, then make the bounded oil move and gather.Some hypotheses are given as follows to find the mathematical models.The mixture fluid flows along isothermal curves, different phases keep equilibrium state, different components have no chemical reaction and expanded Darcy Theory holds and so on.In view of the pressure p(x,t) and the saturation c i (x,t), the flow equation and convection-diffusion equations are derived with corresponding boundary value conditions and initial conditions.nonlinear coupled partial differential equations.It is hard to solve this system since many modern numerical methods are involved in the simulation.In general speaking, based on physical meanings the pressure function is solved by an implicit scheme and the concentration values are obtained by an explicit solver or an implicit solver.The scholars try to find good ways analyzing the data and numerical results and accomplishing some research work in simulation, and the results can describe the whole process of chemistry displacements very well and help the engineers control the rules and process of displacement and forecast the recovery efficiency of natural oil and compute the oil percentage of output liquid and the percentage of the polymer, surface active agent and alkali.By numerical research the curves describing different components motion are shown, and some plans are made about the beginning and end of injected liquid and some related parameters of natural oil efficiency are derived.These conclusions, important techniques in chemistry displacements, can used in forecasting the characters of fields, choosing different optimization plans, establishing the models of chemical displacements of reservoir, completing computational software and carrying out the numerical simulation.Petroleum engineers and mathematicians pay more attention to modern new techniques of exploiting natural oil.
Yuan visited United States and accomplished some work cooperate with Prof. R. E. Ewing during 1985Ewing during to 1988, and kept a series of research in theoretical analysis and applications of numerical simulation.With the leading of Yuan, several research groups (1993,1993,1994) undertook some important projects from 1991 to 1995 such as "Eighth-Five" national key science and technology program (the Program for Tackling Key Programs) (85-203-01-087) entitled "research and application of the polymer displacement software" 1 .The software has been applied in designing plan and research work of polymer displacements in industrial production region of Daqing Oilfield.Many numerical results are concluded (1998,2000) in view of effects of fragments, fragments setting of rinsing protection, quantity of polymer, and used in actual simulations which give rise to outstanding economic and social benefits 2 .Later the authors accomplish a key tackling program of Daqing Petroleum Administration Bureau (DQYJ-1201002-2006-JS-9565)--solving development of mathematical models and completing explain of reservoir 3 .This software system is also applied in numerical simulation of the polymer displacement of Zaobei fault block of Dagang Oilfield, optimization of designing plan of expanded experimental area of three compound combination flooding of Gudong Little Well experimental region of Shengli Oilfield, polymer flooding of Gudong Middle One experimental region, Gudong West region and feasibility of active water flooding of Gudong eighth region, and many interesting results are obtained 4 .In recent years the research group finishes the successive key tackling project of Daqing Petroleum Administration Bureau (DQYJ-1201002-2009-JS-1077)-research on alkali flooding principle model of chemical displacement simulator and horizontal wells model structuring and solving methods 5 and presents many important results.This paper summarizes the former research and discusses proceeding analysis, mainly consisting of permeation fluid mechanical mathematical models of numerical simulation of chemical flooding in porous media, numerical methods, applicable software structuring and numerical simulation, applications and analysis of horizontal wells in actual oilfields.

Permeation Fluid Mechanics Model of Chemical Flooding
The numerical reservoir simulation program of three-phase model is shown based on the following factors: three phases (oil, gas and water) in the reservoir, oil-component and solution gas-component of the oil phase, water-component of the water phase, gas-component of the gas phase and the exchange of gas-component between the oil phase and the gas phase because of the changing pressure.An improved three-phase model is given to simulate the displacement process of chemical compound combinations (polymer, surface active agent and alkali) and a new simulation system is made for three-phase chemical compound combination flooding by considering three-phase module and the chemical reaction balance equation, and revising and adding chemical compound combination displacement module and is applied successfully in numerical simulation of horizontal wells in actual oilfields.
The three phase oil, water and gas are considered in the system, where the oil phase consists of oil component and solution gas component, the water phase includes water component, the chemical agent component includes polymer, surface active agent, petroleum acid and so on, and the gas phase only has gas component.Petroleum acid is contained in the water phase and oil phase and the other chemical components are only in the water phase.The exchange of gas components occurs between the oil phase and the gas phase when the pressure of neighbor environments changes.The solving system consists of three basic modules: the program for solving the three-phase flow, the program for solving the component equations and the program for solving the chemical balance equation.
The solving module of three phases inherits a partial algorithm of the three-phase model.The water viscosity is a variable in chemical compound combinations displacement because of the chemical agent component but a constant in the black oil model.Additional design structure and algorithm should be compatible with the black oil model for solving the components equation and chemical reaction balance equation, and some computational programs are modified for solving the system to satisfy engineering applications of wedge-out area, fault and edge-bottom water.
Then two different algorithms are presented in this paper.

(I)
Full implicit algorithm (FIA).It is one of the most dependable finite differences to compute the values of variables implicitly such as the pressure, the water saturation, the gas saturation and the ratio of solution gas and oil.More computational time is spent for a full implicit scheme to complete each extrapolation iterative.The full implicit scheme proceeds more stable than an implicit/explicit algorithm for the pressure and saturation (IEAFPS) which is invalid sometimes because of strong stability conditions, so a large time step is introduced to decrease the total simulation cost.
(II) Implicit/explicit algorithm for the pressure and saturation (IEAFPS).This algorithm means a coupled numerical system consisting of an implicit scheme for solving the pressure and an explicit scheme for solving the saturation dependent on an assumption: no distinct flow of the fluid occurs though the saturation varies during one time step.Under this hypothesis the unknown saturation variable of discretized convection equation can canceled and the pressure is computed coupled implicitly only at different iterations.The values of saturation is updated point by point explicitly once the pressure change rules are determined during one iteration process.While the IEAFPS is unstable as the saturation changes more large with respect to time step.This implicit/explicit algorithm is really efficient as the time step is taken sufficiently small to decrease the rangeability (the relative rangeability usually is 5%).
Let "w" and "o" denote the water phase and the oil phase of water-oil two phases problem, whose mathematical model of permeation fluid mechanics is stated as follows (see Ewing, Yuan, ＆ Li, 1989;Yuan, 2013;Yuan, Yang, ＆ Qi, 1998;Yuan, 2000): , , ; Let "w", "o" and "g" respectively denote the water phase, the oil phase and the gas phase, and the mathematical model of three phases permeation fluid mechanics is followed 12345 , , , 1 ; where means the porosity, l p is the pressure of l-phase, l S is the concentration of l-phase, K is the absolute permeability, l B is the volume factor of l-phase, rl K is the relative permeability of l-phase, l  is the viscosity of l-phase, l  is the density of l-phase, s R is the ratio of solution gas and oil, and l q is the source sink term of l-phase (floor condition).
The viscosity of water phase w  is a constant of black oil model, while the viscosity is a function with respect to the density of polymer of black oil-polymer displacement, ( ) , where pw C denotes the concentration of polymer relative to water.The compound combination flooding components move in water, the concentration affects in turn the viscosity field of the water phase, then influences the flow of the three phase fluid accompanying the motion of components.The coupled algorithm runs consistent with the nonlinear coupled system consisting a convection-diffusion equation of the polymer and the mathematical model of black oil.At a time step the motion values of three phases are computed first, the flow field is obtained, then the solution of three compound combination flooding by solving a convection diffusion equation are gotten, then the viscosity field of water phase is updated.Then the computation proceeds at the next time step.
Surface active agents can change interfacial tension, thus increase capillary number, then decrease the residual oil saturability.The curve of relative permeability is updated by an interpolation of the relative permeability curves of lower capillary number and upper capillary number.
The primary principle of alkaline flooding is stated as follows.The alkali is injected into the fluid and surface active agents arise because of a chemical reaction of the alkali and petroleum acid.These agents can decrease interfacial tension of the fluid and reduce bottom oil.Petroleum acid components exist both water phase and oil phase, so mass transfer takes place between different phases.An assumption is given that the petroleum acids in water and oil can level off instantaneously, and the equation of total concentration of petroleum acid (a type of convection) is stated.
The computation of the black oil-polymer displacement model runs in the following steps: Step 1.At the time level t 1 , solve the pressure and saturation, then solve the concentrations of different chemical components.
Step 2. Correct the viscosity of water phase dependent on the concentration.
Step 3. At the time level t 2 , solve the pressure and saturation, then solve the concentrations of different chemical components.
Step 4. Correct the viscosity of water phase dependent on the concentration.

……
In a similar way the pressure and saturation of the (n+1)th time step t n+1 are obtained, the concentration of different chemical components at t n+1 is computed, then the viscosity of water phase is corrected according to the concentration.

……
The program ends.

Numerical Methods
Two different numerical methods are presented as follows.

Full Implicit Numerical Method of the Three Phase Flow
It aims to cancel some spare unknown variables and compute only three unknown variables, generally meaning the pressure of oil phase, the concentrations of water and gas phases.All the values of the left side of the equation such as the pressure, the saturation, the quantity, the relative permeability, the capillary force and other parameters are replaced by the latest numerical values during implicit computations.The resulting implicit difference is a nonlinear algebraic equations solved by iteration and its computation scale is more rich seven times the implicit/explicit method of one iteration (outer iteration).Because the full implicit method is unconditionally stable, it is applied to compute some difficult and complicated problems such as black oil simulation.In this paper let  denote the difference of a function between the n-th time level and the (n+1)-th level and let be the difference during a time step by one iteration such as from the k-th iteration to the (k+1)-th iteration.
Applying Euler backwards finite difference method for eq.( 2), where 1   as l=g and 0 Then a full implicit finite difference algorithm for the black oil model is derived as follows The iterations are convergent as 0 k l R  for l=w,o,g and k=1,2,….A resulting formulation is 1 2 3 , , , .
For the right side term of oil phase equation, where For the right side term of the gas phase equation, where (9) In the above expressions ' l b and '  present the derivatives of volume factor and porosity with respect to the pressure, cow p ' denotes the derivative of p c with respect to S w and cgo p ' is the derivative of p c with respect to S g .Two operators are introduced to discuss the left side term, Only consider the discussion of M l .For water phase, The values of conductivity coefficient appearing in two-order difference operator are given according to 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, ), , , .
Substitute the left side and the right side terms into the initial difference equation, and get the algebraic system.

Implicit/Explicit Algorithm for the Pressure and the Saturation
The implicit/explicit method is based on an equation only involved of the pressure by combining the flow equations, and the values of the saturation are obtained explicitly as the pressure at some time level is known.
Discrete difference scheme of eq. ( 2) is written related with p o and the saturation, A basic hypothesis is assumed for deriving the implicit/explicit scheme that the capillary pressure of the flow term of the left side is a constant number during a time step computation.The values of some terms related of cow p  and at the previous time level can be computed explicitly, and .The notation p o is simplified to be replaced by p.
The parameter C is defined by ) Combine the three equations of (11) in a suitable way and cancel all the t l S  terms by multiplying the water phase equation by A, testing the gas equation by B, and adding three equations.The right term is The parameters A and B are determined by the following relations.
Therefore, the pressure equation is changed into .
It is a typical finite difference equation from a parabolic type, whose matrix form is where T denotes a tridiagonal matrix and D is a diagonal matrix.The vector G is dependent of gravity and capillary pressure.
When the pressure is solved, substituted into the former two equations of (11), the saturation 1 n l S  is obtained explicitly.Then the capillary pressure 1 n cow p  and 1 n cog p  are computed and are used explicitly in next time level.

Numerical Method for Component Concentration Equations
The components meaning sorts of chemical agents in water phase such as surface active agents, polymer, alkali and kinds of ions and so on are considered in this paper.The physical nature, conversation of mass, is described by a convection-diffusion equation but convection-dominated.It is more efficient and high order of accuracy to decompose the equation into a hyperbolic equation only related of diffusion and a parabolic equation only related of dispersion.The former is discretized by an implicit upwind scheme and an upstream rule inheriting some advantages of explicit algorithms such as solving the values point by point.The latter is solved by an alternating direction finite difference method, which can accelerate the computation speed and improve the efficiency.For simplicity, the concentration equation of components is rewritten by a typical convection diffusion type.
where dispersion tensor is a matrix of diagonal type.When the saturation S w and the flow velocity field u w of water phase at t n+1 are given, it is to compute the concentration C n+1 at the first direction, the second and then the third direction alternatively.C denotes the concentration of any component ignoring the subscript k.A convection equation is discussed by an implicit upwind scheme, where u w =(u wx ,u wy ,u wz ) T is a velocity vector.By which C n+1,0 is obtained.Then the dispersion equation is discussed alternatively in three directions.First consider the computation x-direction,   Secondly consider the diffusion in y-direction, Thirdly consider the diffusion in z-direction, then C n+1 is computed by The values of and 1 n C  are obtained, and the computation ends at present time level then continues at next time level.A technique is introduced in actual applications to overcome some difficulties and effects of grids orientation (such as the symmetrical computation kept consistently with symmetric physical problems).The values are obtained in two continuous processes, x-direction first then y-direction, then y-direction first then x-direction, then the algorithm continues in z-direction by using an average result of the above computational values.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.Considering the nonlinear equations F(X)=0, where k denotes the number of iterations, , ) , x and the convergent condition is defined by The choice of the initial values X 0 affects the convergence of Newton-Raphson iteration sequences.DF(X k ), a Jacobian matrix, is defined by ( , , , , ) 0, ( , , , , ) 0, ( , , , , ) 0, where denotes the partial derivative value of

Computation Program Illustration
This section illustrates three computation programs 1,2,3,5 : the program of three-phase chemical compound combination flooding, the implicit computational program of the water phase concentration equation, the explicit computation program of the (total) concentration equation of petroleum acid components, the computation program of chemical reaction balance equation (see Fig. 1, Fig. 2, Fig. 3 and Fig. 4).

Experimental Tests for Horizontal Wells Model
An experimental test is given for the effect of horizontal wells.The simulation domain is partitioned by 9×9×1, the simulation period is 5500 days and the polymer is injected at the 1460th day.The test proceeds in the following three programs.The first program is five vertical wells (four injected wells and one production well), the second program is four injected wells and one production well (horizontal well) and the third program is three horizontal wells (two injected wells and one production well).
The 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 (see Fig. 5, Fig. 6, Fig. 7 and Fig. 8).From these figures some inclusions are shown that the software can illustrate the process of horizontal well reservoir with high order of accuracy.The numerical simulation reflects the effect and the principle of horizontal wells correctly.The values of primary physical functions distribute reasonably and the computational accuracy is satisfactory.Some unexpected cases such as the polymer tacking and endless loop occur in this simulation.( , , , , )

Mechanism Analysis of Alkaline Flooding
The computational grids are defined by 9×9×1 and the flooding efficiency is tested under different concentrations of injected alkali and different Naphthenic acid values.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.
There five flooding plans are considered here.Let P denote the polymer flooding.Let ASP denote the compound combination of the polymer and alkali, the acid value is taken as 0.0006 and the concentrations of injected Na ions and CO 3 ions are 0.3351 and 0.3929, respectively.ASP3 denotes the third case, the compound combination flooding of the polymer and alkali, 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, 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, the acid value is 0.006, and the concentrations of injected Na and CO 3 are 3.3351 and 3.3929.Based on computational data of produced surface active agents, the concentration of residual oil and the generated depositing under five different cases, there four conclusions are given as follows.

(I)
From the normal concentration values at $3000$th day under different cases, the produced surface active agents grows more and more as the acid value and the concentration are larger and larger.

(II)
From the values of residual oil at 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 is 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. 9).

Theoretical Analysis of the Model
Theoretical analysis of numerical simulation is given for three dimensional three-phase (oil, water and gas) chemical compound combination flooding in porous media.A simplified model is discussed in theoretical analysis without loss of generality, that is to say a compressible oil water two phase displacement of three dimensional multi-components problem in porous media is discussed in a nonlinear partial differential equations with initial and boundary values (see Ewing, Yuan, ＆ Li, 1989;Yuan, 2013;Ewing, 1983;Yuan, 2002Yuan, , 1999Yuan, , 2003Yuan, , 2001 ) ). ], where  is a bounded domain, and .Let s(x,t)=(s 1 (x,t),s 2 (x,t),…,s nc-1 (x,t)) T be the vector function of component saturations.The pressure p(x,t), the saturation of water c(x,t) and the component concentration vector s(x,t) are basic unknown functions to compute.

Boundary condition without permeation
where  is the outer normal vector of boundary surface  of  .
Initial conditions: Douglas (1983) presents a fundamental paper to analyze a type of two dimensional incompressible two-phase displacement problems.Because the computation of modern reservoir exploration and development is of huge scale, its simulation region is large, its simulation time is really long and there are tens of thousands or hundreds of thousands nodes, it is impossible to solve so complicated problems of this type by using common methods.Alternating direction finite difference methods are put forward by Peaceman (1980) and Douglas (1963) and are successfully applied in two dimensional problems but there is substantial difficulty to give theoretical analysis.Using Fourier method Peaceman and Douglas discuss the stability and convergence for the problems with constant coefficients, while this method is not able to generalize in variable coefficient problems.Yanenke (1967), and Marchuk (1990) give many important results on fractional steps methods.Yuan (1999) presents a characteristic finite difference fractional steps method for compressible two-phase displacement problem and discusses convergence analysis.An implicit upwind finite difference fractional steps method is considered for black oil and polymer flooding problem and some substantial improvements are given in this paper.The three dimensional problem is turned into three one-dimensional problems and this can greatly reduce the computation and can solve actual problems.Error estimates in L 2 norm are presented by using variation, energy analysis, decomposition of high order difference operators and theory and technique of product commutativity.
For simplicity, assume computational domain  ={[0,1]} 3 and the problem is  -periodic.The nonpermeation boundary condition can be dropped.Let h=1/N, X ijk =(ih,jh,kh) T , t n =n t  and W(X ijk ,t n )= n ijk W , and let The symbols , 1/2, can be defined similarly.
The implicit finite difference fractional steps algorithm for the flow equation ( 20a) is shown as follows, and similar formations 1 2, The implicit upwind finite difference fractional steps algorithm for the saturation equations is The components concentration equation ( 20d) is computed by an upwind fractional steps algorithm implicitly in parallel, The implicit upwind fractional step algorithm runs as follows.Given .Notice the problem is positive definite, so the solution of ( 24)-( 28) exists and is unique.
Applying theories and techniques of priori estimates of differential equations, and complicated calculations, we deduce the following convergence theorem.
Theorem Suppose that the exact solutions of ( 20)-( 22) are suitably smooth.Applying implicit upwind difference fractional steps method ( 24)-( 28) to compute layer by layer, we can conclude the following error estimates,  

Conclusion and Discussion
Theory, method and application of numerical simulation of three-dimensional three-phase (water, oil and gas) percolation mechanics of chemical compound combination flooding (the polymer, surface active agents and alkali) in porous media are discussed in this paper consisting of seven sections.Summary is first stated about our project.Mathematical model of permeation fluid mechanics is presented in §2.A full implicit numerical scheme and implicit/explicit algorithm for the pressure and the saturation are given, and an upwind difference fractional steps algorithm based on upstream sequence is structured in §3.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 §4.Some tests are considered for horizontal wells model in §5.Some experimental examples occurring successfully in national major oil fields such as Daqing Oilfield, are illustrated in §6.Numerical analysis proceeds for the model problem and precise theoretical results are stated on mathematical and mechanical consideration in §7.
acid dispersion tensors (including molecular diffusion and dispersion) of water phase and oil phase, respectively.The total concentration . Rewrite the above expression as follows by using an operator , Give an expansion of the above equation and omit quadratic terms, and the remainder after the k-th iteration is Continue to express the above equation with a remainder term as, 6) First, the equation is turned into a linearized expansion.Solving variables are denoted by o p  , w S  and g S  , and the right side term is considered later.the right side term of the water phase equation, where

(
is the viscosity, and D=D(x) is the diffusion coefficient.u, p=p(x,t) and c=c(x,t) are Darcy velocity, the pressure and the saturation of water phase, respectively.( , ) s x t  is the concentration function of the -th component, where the components mean sorts of chemical agents (the polymer, surface active agent, alkali and other ions) and n c is the number of components.There are n c -1 independent components because of 1 layer by (24a) and speedup method in x 1 direction firstly.24b) and by (24c), respectively.Using (25) we can have the numerical values of . .