Theory and Applications of Numerical Simulation of Permeation Fluid Mechanics of the Polymer-Black Oil

Theory and application of numerical simulation of permeation fluid mechanics is discussed about the polymer--black oil (water, oil and gas) in this paper. In view of petroleum geology, geochemistry, computational mechanics of flow and computer technology, the mechanics model of three-phase flow (water, oil and gas) of the polymer is presented firstly, then a numerical algorithm consisting of a full implicit program, and an implicit computation for the pressure and an explicit computation for 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 concentration equation and the saturation of chemical substance components. A type of software applicable in major industries has been made in consideration of ten-meters steps, hundreds of thousands nodes and tens of years and high accuracy and the application has been carried out successfully in analysis and simulation of 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 an idea is presented to solve this international famous problem.


Introduction
At present an effective method, water-flooding, to hold the pressure of reservoirs is popular in the world, and the recovery efficiency is more outstanding than any other natural exploring forms.It gives more benefits and helps Chinese oil fields keep high quantity production.It continues to be more important and how a strategic project works to develop the exploiting efficiency of crude oil in the way of water-flooding driving.
A mass of residual crude oil remains in the reservoir after water-flooding exploiting because of the constraint of capillary force preventing the motion and the slight influenced regions and the fluidity ratio between displacement phase and driven phase weakening the displacement force.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, increase the efficiency of displacement and the pressure gradient.Surface active agent and alkali can decrease interfacial tensions of different phases, then make the bound oil move and gather.Some hypotheses should be made 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 ( , ) (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.Most time is taken 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) 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 consists of an implicit scheme for solving the pressure and an explicit scheme for solving the saturation based on the following hypothesis: no distinct flow of the fluid occurs though the saturation changes during one time step.This hypothesis can cancel the unknown saturation of discretized convection equation.The values of saturation is updated point by point explicitly when the pressure changes keep implicit coupling and the change rules are determined during one iteration process.While the IEAFPS is unstable as the saturation changes faster with respect to time step.This implicit/explicit algorithm is very efficient as the time step is taken sufficiently small to decrease the rangeability (the relative rangeability usually be 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 polymer components motion in water, the concentration affects in turn the viscosity field of the water phase, then influence the flow of the three phase fluid accompanying the motion of components.The coupled algorithm runs consistent with the nonlinear coupled system concluding 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 convection diffusion equation are gotten, then the concentration field of polymer and the viscosity field of water phase are updated.
Then the computation proceeds at the next time step.
The following convection diffusion equation explains the motion law of the polymer, anions and cations, for convenience the different components are not distinguished and the subscript w denotes the concentration of a certain component.
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 concentration of the 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 concentration of the 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 the components at t n+1 is computed, then the viscosity of water phase is corrected according to the concentration.

……
The program ends.

Numerical Computation Method
Two different numerical methods are presented as follows.

Full Implicit Numerical Method of the Three Phases Flow
It aims to compute three unknown variables, generally the pressure of oil phase, the concentrations of water and gas phases are considered and other spare unknown numbers are cancelled.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 taken 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 one iteration (outer iteration).By that 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 values after one time step and let be the values of the (k+1)-th iteration after one 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 method for the black oil model is derived as follows where . Rewrite the operator as follows, Give an expansion of the above equation and omit quadratic terms, and the remainder after the k-th iteration is Continue to give an expression with remainder term, 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 water phase, where Consider the right side term of oil phase equation, where .  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, The difference equation is rewritten as follows, ( ) , , , .
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 adding all the flow equations, and the values of the saturation are obtained explicitly as the pressure at a certain 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 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
The components in water phase such as anions, cations, polymer molecules 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 is of high order 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 as a typical convection-diffusion expression 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, by which C n+1,0 is obtained.Then the dispersion equation is discussed alternatively in three directions.First consider x-direction,  

Computation Program Illustration
This section illustrates three computation programs (Note 1, Note2, and Note 3): the program of black oil and polymer flooding, the data program of upstream sequence algorithm connected by faults, and the implicit computational program of the concentration equation of black oil and polymer components (see Fig. 1, Fig. 2 and Figure 3).

Experimental Tests for Black Oil and Polymer Flooding
The model grids scale of experimental tests is defined by 9×9×3, and its lengths in x-direction and y-direction are DX=DY=44.5meters and its height is DZ=2.0 meters.The simulation begins at January 1, 2000 and ends at January 1, 2015.The polymer is injected from January 1, 2004 to January 1, 2008, and the concentration is 2000 mg/L.Concentration isogram of the polymer, saturation isogram of water and viscosity isogram of water phase are illustrated in nine figures (Fig. 4-Fig.6) in three different layers at the end of injection, January 1, 2008.These numerical data are shown in Fig. 7, Fig. 8 and Fig. 9 layer by layer when the moisture content is 98% and the date is August 20, 2011.During the whole injection and production, the moisture content of production wells, the curve of water yield and the curve of oil production are shown in Fig. 10 ①②③ .What numerical data conclude is that black oil and polymer flooding can precisely explain the physical process and mechanism of the displacement because the distribution curves of main physical quantities such as saturation, component concentration, viscosity of water and so on are reasonable, the order of accuracy is satisfactory, and undesired results such as polymer stacking and infinite circulation appear.Error estimates of mass balance of water, oil, gas, the polymer, anions and cations are respectively denoted by 2.71E-6, 2.70E-7, 2.70E-7, 1.58E-15, 7.84E-16, 2.00E-15.

Theoretical Analysis of the Model
Theoretical analysis of numerical simulation is given for three dimensional three phases (oil, water and gas) polymer 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 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 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: 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 .
can be defined similarly.
The implicit finite difference fractional steps algorithm for the flow equation ( 20a) is shown as follows, 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  , summation by parts, and obtaining an inner product form 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 From ( 20c) and (32), where , summation by parts, and obtaining an inner product form where Collecting (36a) and (36b), 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 polymer flooding in porous media are discussed in this paper consisting of six sections.Summary is first stated about our project.Mathematical model of permeation fluid mechanics is presented in the first 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 second section.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 third section.Some experimental tests occurring successfully in major oil fields such as Daqing Oilfield, Shengli Oilfield and Dagang Oilfield, are illustrated in the fourth section.Numerical analysis proceeds for the model problem and precise theoretical results are stated on mathematical and mechanical consideration in the fifth section.
are known, and the computation continues at the next time level.
layer by (24a) and speedup method in x 1 direction firstly.Then 2computed by (24b) and by (24c), respectively.Applying (25) we can have the numerical values of velocity 1 by (26a) and speedup method in x 1 direction. . .