Theory and Application of Characteristic Finite Di ff erence Fractional Step Method of Capillary Force Enhanced Oil Production

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


Background and Mathematical Description
Massive residual crude oil remains in the reservoir after water-flooding exploiting because the constraint of capillary force prevents the motion and the undesired fluidity ratio between displacement phase and driven phase influences the flow slightly.Then it is more important to develop the displacement efficiency.A popular method is discussed by adding some chemical addition agents such as polymer, surfactant and alkali into the injected mixture, which can improve the flooding efficiency.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, and increase the efficiency of displacement and the pressure gradient.Surfactant and alkali can decrease interfacial tensions of different phases, then make the bound oil move and gather (Ewing, et al., 1988;Yuan, 2013;Yuan, et al., 1998)1 2 3 .
In view of capillary force, and immiscible and incompressible flow, this paper discusses a second-order characteristic fractional step difference method for numerical simulation of enhanced (chemical) oil production in porous media, and gives the theoretical analysis.Based on the former mathematical and mechanical theory, the software is accomplished, applied in national oilfields such as Daqing Oilfield and Shengli Oilfield and give rise to outstanding benefits and social value.
The mathematical model is described by a nonlinear coupled system with initial-boundary values (Ewing, et al., 1988;Yuan, 2013;Yuan, et al., 1998)  1,2,3 : where Ω is a bounded computational domain.Let the subscripts "o" and "w" respectively denote the parameters of oil phase and water phase.The notations c l , p l , κ γ l (c l ), µ l and q l 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, κ(X) means the absolute permeability and s α = s α (X, t) 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 .u is Darcy velocity, K α = K α (X) is diffusion coefficient, and Q α is source sink term related with the output.The rock void space is assumed to be full of water and oil, c o + c w = 1.The the capillary force is dependent of the concentration c, defined by The models ( 1) and ( 2) should be expressed by a normal form (Ewing, et al., 1988;Yuan, 2013).Let λ(c) = κ γo (c) denote the total migration rate of two-phase fluid, and let λ l (c) = c) , l = o, w denote relative migration rate corresponding to water or oil.Applying Chavent transformation (Ewing, et al., 1988;Yuan, 2013;Douglas, 1983): The flow equation is derived from (1) and (2), where q = q o + q w .The concentration equation is derived from the difference of (1) and (2), Let u = −κ(X)λ(c)∇p, and note that the fact λ o − λ w = 2λ o − 1, we get where q w = q and q o = 0, if q ≥ 0 (water injection well), q w = λ w q and q o = λ o q, if q < 0 (oil production well).
(1) and ( 2) are rewritten as follows, It is clear to restate ( 6) and ( 7) in a normalized formula, It follows that g(X, t, c) = λ o q as q ≥ 0, and it is true that g(X, t, c) = 0 as q < 0. Using (8), we can express (3) in a computational form Two different boundary value conditions are considered in this paper.
(I) The boundary value condition for the constant pressure is given as follows p = e(X, t), X ∈ ∂Ω, t ∈ J, (11a) where ∂Ω denotes the outer boundary surface of Ω.
(II) The boundary value condition for no permeation case is defined by 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 and the initial value condition is

Relevant Development
Introducing an assumption of periodic condition, Douglas, Ewing, Wheeler, Russell and other scholars presented characteristic finite difference method and characteristic finite element method to analyze a type of two dimensional incompressible two-phase displacement problems and gave theoretical error estimates (Douglas, 1981(Douglas, , 1983;;Douglas, Russell, 1982;Ewing, Russell, Wheeler, 1984).A combination method was discussed by a whole consideration of the characteristic method and normal finite difference method or normal finite element method, which can reflect the hyperbolic nature of one-order part of convection-diffusion equation and decrease truncation error.This combination technique can also overcome numerical oscillation and dispersion, and can improve greatly the computational stability and accuracy.Douglas and other scholars presented a mathematical model of slight compression, numerical method and theoretical analysis for two-dimensional compressible displacement problem under periodic assumption and began a modern numerical model research (Douglas, Roberts, 1983;Yuan, 1992Yuan, , 1993;;Ewing, 1983).The authors dropped the period condition, gave a new modified characteristic finite difference algorithm and finite element algorithm, and derived optimal order error estimates in L 2 -norm (Yuan, 1994(Yuan, , 1996(Yuan, , 1996;;Axelsson, Gustafasson, 1979;Ewing, Lazarov, et al., 1994, 1996).In numerical simulation of modern oilfields exploration and development, the computation is of greatly huge-scale, three-dimensional region and long time interval consideration and the number of nodes maybe amount up to tens of thousands or even millions.It is impossible to solve this computation by using normal methods and a typical technique, fractional step method, is introduced in this paper (Peaceman, 1980;Douglas, Gunn, 1963, 1964;Marchuk, 1990).Though Douglas and Peaceman applied this technique into oil-water two phases displacement problem successfully (Peaceman, 1980), a type of substantial difficulty appeared in their theoretical analysis.Fourier method, considered in the proof of the stability and convergence by Douglas and Peaceman, was only applied in the constant-coefficient problems, that is to say Fourier method is not popularized in general problems with variable coefficients (Douglas, Gunn, 1963, 1964;Marchuk, 1990).Considering actual application, numerical stability and accuracy, the authors present one second-order characteristic fractional step finite difference method for three-dimensional incompressible and immiscible two-phase displacement coupled problem of enhanced oil production with consideration of capillary force.This algorithm can overcome numerical oscillation and dispersion, and can decrease the computational scale by decomposing three-dimensional problem into three successive one-dimensional subproblems.Using the calculus of variation, energy analysis method, commutativity of the products of difference operators, decomposition of high-order difference operators and the theory of a priori estimates, the authors give the second-order convergence result of accuracy and error estimates in l 2 -norm, and successfully solve the international problem formulated by Douglas and Ewing.The method discussed in this paper has been applied in numerical simulation of enhanced oil production and gives a fundamental research in energy mathematics 1,2,3 .

Hypotheses
In general, the problem is positive definite, where a * , a * , ϕ * , ϕ * , D * , D * K * , K * and A * are positive constants.b(c), g(c) and Q α (c, s α ) are Lipschitz continuous in the ε 0 neighbours of exact solutions.
Exact solutions of ( 8)∼( 14) are assumed to be suitably smooth, In this paper M and ε express a generic positive constant and a generic positive small constant, respectively, and they have different meanings at different places.

Second-order Implicit Characteristic Fractional Step Finite Difference Method
Without loss of generality, the computational domain is taken as a rectangular cube Ω = {[0, 1]} 3 and the problem is assumed to be Ω-periodic.Note that numerical simulation is considered for the flow of interior region and the boundary condition affects the flow weakly, so the periodic assumption is reasonable and related contents can be found in the references (Douglas, 1983;Douglas, Russell, 1982;Douglas, 1981;Ewing, Russell, Wheeler, 1984;Douglas, Roberts, 1983).Then The boundary value conditions of no permeation case can be dropped.Let the partition Ω h replace the domain Ω and ∂Ω h denotes the boundary.Let h = 1/N denote the uniform spacial step in three-dimensional domain and the notations are defined by and The flow equation ( 8) is discretized by where Q h denotes the cube of the side-length h centered at the origin.Darcy velocity and the other quantities U n+1 2,i jk , U n+1 3,i jk are computed analogously.The implicit characteristic fractional step algorithm is discussed to solve the saturation equation ( 9).The method of characteristics is applied in discretizing first order hyperbolic part of ( 9), which is consistent with the physical transportation of the flow along the direction of characteristics.The algorithm has high order of accuracy and great stability, and runs well in a large time step (Douglas, 1983(Douglas, , 1981;;Douglas, Russell, 1982) ] 1/2 , ∂/∂τ = ψ −1 {ϕ∂/∂t + u • ∇}, and approximate the derivative along the characteristics by backward differences, It remains to compute the values of the saturation c n+1 at the (n + 1)-th time step as c n given.Replacing the differential quotient by difference quotient, and decomposing the saturation equation ( 9) into the following form, where Then the second-order implicit characteristic finite difference program combined with fractional step is illustrated as follows, The discrete value of saturation C n (X) is defined by a piecewise threefold quadratic interpolation of the values of nodes {C n i jk }, and other symbols are defined by Ĉn Introducing the symbols φn+1 = ϕC n+1 , φn+1,−1 = ( φn+1 ) −1 , and Dn+1 α = ϕC n+1 K α , we can approximate the concentration equation ( 10) by an implicit characteristic fractional step finite difference parallel algorithm, where The implicit program runs in the following order.Given {P n i jk , C n i jk , S n α,i jk , α = 1, 2, • • • , n c }, the pressure {P n+1 i jk } is computed by elimination method or conjugate gradient method from the difference equation ( 17).Then the values of Darcy velocity {U n+1 i jk } are obtained by ( 18).The solution of transition sheaf {C n+1/3 i jk } in x 1 direction is computed by the method of speedup from (20), then {C n+2/3 i jk } proceeds in x 2 direction by ( 21), {C n+1 i jk } is obtained finally in x 3 direction by ( 22) similarly.Then the computation of the saturation goes on.The saturation {S n+1 α,i jk } is given by ( 25) after {S n+1/3 α,i jk }, {S n+2/3 α,i jk } are computed respectively in x 1 , x 2 directions by the method of speedup from ( 23) and ( 24) as shown in the above expression.It is accomplished in parallel with respect to α = 1, 2, • • • , n c .The solutions of ( 17), ( 20)∼( 22), ( 23)∼( 25) and ( 26) exist and are sole because of the positive definite condition.
Theorem 1 Assume that exact solutions of ( 8)∼( 13) are appropriately smooth, p, c The modified characteristics method is applied to compute ( 17)∼( 25) layer by layer, and the partition parameters satisfy then the following error estimates hold Proof.The error equation of the pressure is derived from the difference of (8)(t = t n+1 ) and ( 17), where Testing both sides of (29) by π n+1 , and summing by parts Then, The error of the saturation is estimated later.Canceling the transient terms C n+1/3 i jk and C n+2/3 i jk in (20), ( 21) and ( 22), introducing Ūn+1 = b(C n )U n+1 , and simplifying the dispersion term D(C n ) ≈ D(X) for convenience (this simplification does not affect the nature of numerical analysis, that is to say the following process can be applied in the theoretical analysis related of D(C n )), we can obtain an equivalent difference equation,  1.From this it is easy to see that numerical method of this type and its applicable software can keep high order of accuracy and can reflect correctly the physical process and the principle of polymer flooding.Main physical quantities distribute reasonably, computational accuracy is high, and some results of the polymer such as accumulation, endless loop don't arise.

Figure 2 .
Figure 2. Sketch of instantaneous oil production curve of the polymer in Xing Fourth Zone