Numerical Simulation of Black Oil-three Compound Combination Flooding

Numerical methods of permeation fluid mechanics for black oil (water, oil and gas)--three compound combination flooding (polymer, surface active agent and alkali) in porous media is discussed in this paper. In view of petroleum geology, geochemistry, computational mechanics of flow and computer technology, a mechanics model of three-phase flow (water, oil and gas)--three compound combination flooding (the polymer, surface active agent and alkali) is presented firstly, then the characters and the application are stated in this paper. 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. This program runs quickly and is of high accuracy. The design of this software is given and can be applied in major industries, which is made up on ten-meters steps, hundreds of thousands 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, Shengli Oilfield and Dagang Oilfield and others, which gives rise to outstanding economic and social benefits. A precise analysis is given for a simplified model and the numerical simulation system depends on mathematics and mechanics. An idea is presented to solve this international famous problem.


Introduction
A popular numerical simulation way, water flooding, is applied in the world to keep the reservoir pressure, and the recovery efficiency is more outstanding than any other natural exploring forms.This gives more benefits and helps Chinese oil fields keep high quantity production.It continues to be more important how a strategic project works to develop the exploiting efficiency of crude oil in the way of water-flooding driving.
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 This paper concludes the former research and discusses proceeding analysis, mainly consisting of permeation fluid mechanical mathematical models of numerical simulation of black oil--three compound combination flooding in porous media, numerical methods, applicable software structuring, theoretical analysis and applications of actual oilfields.

Permeation Fluid Mechanical Model of Black Oil--three Compound Combination Displacement
The numerical simulation program of black oil 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 black-oil model is given to simulate the displacement process of three compound combinations (polymer, surface active agent and alkali) and a new simulation system is made for three compound combinations by considering the black oil module and the chemical reaction balance equation, and revising and adding three compound combinations module.
The oil, water and gas three phases 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 black oil model.The water viscosity is a variable in three 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 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 of 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.Eqs.(1a), (1b), (2e), (2f), (2g) and (2h) give a full description of the mathematical model of water and oil two-phase polymer displacement, and eqs.(2a), (2b), (2c), (2d), (2e), (2f), (2g) and (2h) represent the full explanation of the mathematical model of water-oil-gas three phases polymer displacement.
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   as l=w,o, and b V x y z     .
Then a full implicit finite difference algorithm for the black oil model is derived as follows The iterations are convergent as First, the equation is turned into a linearized expansion.Solving variables are denoted by and the right side term is considered later.
Notice that 1 ( ) and consider the right side term of the water phase equation, where Consider the right side term of oil phase equation, where .
For the right side term of the gas phase equation, where 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, Here an expansion of M l is only considered.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 coefficient C is defined by ) Consider the three equations of (11) together 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 numbers A and B are defined by the following relations.
By simple algebraic calculations, Then the pressure equation turns 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.
Given the pressure, taken in 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 considered, which 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 x-direction,   The values of S  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 F=(F 1 ,F 2 ,…,F N ) T , X=(x 1 ,x 2 ,…,x N ) T , 0=(0,0,…,0) T .Its Newton-Raphson iterations are 1 (X ) , 0,1, , (X ) 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 where denotes the partial derivative value of

Computation Program Illustration
This section illustrates three computation programs 1,2,3,5 : the program of black oil-three compound combination flooding, the data program of upstream sequence algorithm of fault-joints, 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 relative permeability and the computation program of chemical reaction balance equation (see Fig.1,Fig.2,Fig.3,Fig.4,Fig.5 and Fig. 6). 1.There are five strategies considered in the simulation and let X45ASP be the polymer flooding.X45ASP2 denotes the compound flooding combination the polymer with alkali, the acid number is 0.0006, and the concentration values of injected Na ions and CO 3 ions are 0.3351 and 0.3929 in the second way.X45ASP3 is the compound flooding combination the polymer with alkali, the acid number is 0.006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the third way.X45ASP4 is the compound flooding combination the polymer with alkali, the acid number is 0.0006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the fourth way.X45ASP5 is the compound flooding combination the polymer with alkali, the acid number is 0.006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the fifth way.The moisture content, instantaneous oil production and accumulative oil production in five different strategies are compared in the following curves (Fig. 7a, Fig. 7b and Fig. 7c).

Test II
The model grids scale in this test is defined by 119×79×9 (84609 nodes), and oil displacement efficiency is tested under different petroleum acid numbers and different concentrations of injected alkali.Three SLUGs are installed and the simulation period is 4000 days divided into three time segments: water flooding segment during the former 730 days, injected polymer and three compound combination flooding segment during the middle thirteen hundreds and seventy days and water flooding segment during the last nineteen hundreds days.There are five strategies considered in the simulation and let B1DDASP be the polymer flooding.B1DDASP2 denotes the compound flooding combination the polymer with alkali, the acid number is 0.0006, and the concentration values of injected Na ions and CO 3 ions are 0.3351 and 0.3929 in the second way.B1DDASP3 is the compound flooding combination the polymer with alkali, the acid number is 0.006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the third way.B1DDASP4 is the compound flooding combination the polymer with alkali, the acid number is 0.0006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the fourth way.B1DDASP5 is the compound flooding combination the polymer with alkali, the acid number is 0.006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the fifth way.The moisture content, instantaneous oil production and accumulative oil production in five different strategies are compared in the following curves (Fig. 8a, Fig. 8b and Fig. 8c).

Test III
The model grids scale of this test is defined by 149× 149×7 (155407 nodes), and oil displacement efficiency is tested under different petroleum acid numbers and different concentrations of injected alkali.Three SLUGs are installed and the simulation period, 2010.9.1-2015.11.1, is divided into three time segments: water flooding segment 2010.9.1-2012.11.1, injected polymer and three compound combination flooding segment 2012.11.1-2014.11.1 and water flooding segment 2014.11.1-2015.11.1.There are five strategies considered in the simulation and let P be the polymer flooding.The number 1 denotes the compound flooding combination the polymer with alkali, the acid number is 0.0006, and the concentration values of injected Na ions and CO 3 ions are 0.3351 and 0.3929 in the second way.The number 3 is the compound flooding combination the polymer with alkali, the acid number is 0.006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the third way.The number 4 is the compound flooding combination the polymer with alkali, the acid number is 0.0006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the fourth way.The number 5 is the compound flooding combination the polymer with alkali, the acid number is 0.006 and the concentration values of Na ions and CO 3 ions are 0.3351 and 0.3929 in the fifth way.The moisture content, instantaneous oil production and accumulative oil production in five different strategies are compared in the following curves (Fig. 9a, Fig. 9b and Fig. 9c).
The conclusion is derived from above figures that the more efficient the oil flooding is, the larger the acid number of petroleum acid and the concentration value of injected alkali.Error data of mass balance are shown in the following table (Table 1).
Table 1.Error results of mass balance

Theoretical Analysis of the Model
Theoretical analysis of numerical simulation is given for three dimensional three phases (oil, water and gas) -three 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, 2002, 1999, 2003, 2001 ).
where p(x,t) is the pressure of the mixture, ( , ) n c is the number of components.There are n c -1 independent components because of .Let c(x,t)=(c 1 (x,t),c 2 (x,t),…,c nc-1 (x,t)) T be the vector function of component saturations, , ( ) x  be the porosity of rock, z  be the compressible invariant of  , u be Darcy velocity of the mixture,  be the permeability, ( ) c  be the viscosity, , and let D=D(x) be the diffusion parameter.The pressure p(x,t) and the saturation vector c(x,t) are basic unknown functions to compute.
Boundary condition without permeation: where  is the outer normal vector of boundary surface  of  .
Initial conditions: (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 can be defined similarly.
The implicit finite difference fractional steps algorithm for the flow equation ( 20a) is shown as follows, Darcy velocity U=(U 1 ,U 2 ,U 3 ) T is computed by the following expression and analogous forms 1 2, The implicit upwind finite difference fractional steps algorithm for the saturation equations is ),1 , , , 1,2, , 1.
The implicit upwind fractional steps algorithm runs as follows.Given  26b) and (26c).Notice the problem is positive definite, so the solution of ( 24)-( 26) exists and is unique.
where p, c  are exact solutions and P, C  are numerical solutions.The pressure equation is considered firstly, and its equivalent difference expression is Then the error equation of the pressure is where Multiplying both sides of error equation ( 29) by Then it follows from (30) by complicated estimates, Secondly an error analysis proceeds for the saturation equations.Eq. ( 26) is rewritten in the following equivalent form From ( 20c) and (32),   Theorem Suppose that the exact solutions of ( 20)-( 22) are suitably smooth.Applying implicit upwind difference fractional steps method ( 24)-( 26) to compute layer by layer, we can conclude the following error estimates,  

Discussion
Theory, method and application of numerical simulation of three-dimensional three-phase (water, oil and gas) percolation mechanics of three compound combination flooding (the polymer, surface active agents and alkali) in porous media are discussed in this paper consisting of six sections.Summary is first stated about our project and the full process of this project.Mathematical model of permeation fluid mechanics, basic physical assumption and the characters of the model are presented in the second section.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 the third section.This program runs quickly, is of high accuracy and applied in general cases.A type of software applicable in major industries has been accomplished, mostly carried out with the spacial step of ten-meters, tens of thousands nodes and tens of years simulation period in the fourth section.Some experimental tests occurring successfully in national major oil fields such as Daqing Oilfield, Shengli Oilfield and Dagang Oilfield, are illustrated in the fifth section.This gives outstanding social benefits and economic benefits, and accelerates the development of energy sciences.Numerical analysis proceeds for the model problem and precise theoretical results are stated on mathematical and mechanical consideration in the sixth section.This research brings important theoretical values for the design and development of enhanced oil recovery simulation, the principle, method and the software.
acid dispersion tensors (including molecular diffusion and dispersion) of water phase and oil phase, respectively.The total concentration HA tot C is defined as follows .
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, Fields of Black Oil-three Compound Combination Alkali Flooding 5.1 Test I The model grids scale of experimental tests is defined by 46×83×7, and oil displacement efficiency is tested under different petroleum acid numbers and different concentrations of injected alkali.Three SLUGs are installed and the simulation period, 1970.1.1-1994.1.1, is divided into three time segments: water flooding segment 1970.1.1-1982.1.1,injected polymer and three compound combination flooding segment 1982.1.1 -1988.1.1 and water flooding segment 1988.1.1-1994.1.
and by (24c), respectively.Using (25) we can have the numerical values of 12345 , .