Theory of Second Order Numerical Simulation Method of Enhanced Oil Production

A kind of second-order implicit upwind fractional steps finite difference method is presented in this paper to numerically simulate the coupled system of enhanced (chemical) oil production 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 optimal order error estimates in l2 norm are derived.


Introduction
A mass of residual crude oil stays in the reservoir after water-flooding exploiting because of the constraint of capillary force preventing the motion and the slight influence of injected water and the undesirable 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 to add 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, Yuan & Li, 1989;Institute of Mathematics, 1995, 2006, 2011;Yuan et al., 1998;Yuan, 2013;Yuan, Cheng, Yang & Li, 2014, 20142, 20143).This paper discusses a second-order upwind fractional steps difference method for numerical simulation of enhanced (chemical) oil production in porous media, and gives the theoretical analysis.Based on the previous mathematical and mechanical theory, the software is accomplished, applied in national oilfields such as Daqing Oilfield and Shengli Oilfield and give rise to important benefits and social value.
The mathematical model is a nonlinear coupled system with initial values and boundary values (Ewing, Yuan & Li, 1989;Institute of Mathematics, 1995, 2006, 2011;Yuan et al., 1998;Yuan, 2013;Yuan et al., 2014, 20142, 20143): where Ω is a bounded domain, a(c) = a(X, c) = k(X)µ(c) −1 , d(c) = d(X, c), and other notations are explained as follows.ϕ(X) denotes the porosity of rock, k(X) denotes the permeability, µ(c) means the viscosity of fluid, and both D = D(X) and ) denote the diffusion coefficients.u is Darcy velocity, p = p(X, t) is the pressure, c = c(X, t) means the saturation of water, and s α = s α (X, t) denotes the concentrations of components.The components denote sorts of chemical agents such as the polymer, surfactant, alkali and other ions, and the number is denoted by n c .
Two different boundary value conditions are considered in this paper.
(II) The boundary values condition for no permeation case: where γ denotes the outer normal unit vector.
Initial conditions: p(X, 0) = p 0 (X), X ∈ Ω, c(X, 0) = c 0 (X), X ∈ Ω, It is easy to compute the concentration by rewriting (3) as the following expression Under an assumption of periodic condition, Douglas, Ewing, Wheeler, Russell and other scholars present finite difference method and finite element method to analyze a type of two dimensional incompressible two-phase displacement problems and give theoretical error estimates (Douglas, 1981(Douglas, , 1983;;Douglas & Russell, 1982;Ewing, Russell & Wheeler, 1984).
A combination of the characteristic method with normal finite difference method or with normal finite element method is discussed, which can reflect the hyperbolic nature of one-order part of convection-diffusion equations and decrease the order of truncation error.This technique can also overcome numerical oscillation and dispersion, and can improve greatly the computational stability and accuracy.Douglas, Yuan and Ewing present mathematical model of slight compression, numerical method and theoretical analysis for two-dimensional compressible displacement problem under periodic condition and begin a new modern numerical model research (Douglas & Roberts, 19832;Ewing, 1983;Yuan, 1992Yuan, , 1993)).
The authors drop the period condition, give a new modified characteristic finite difference algorithm and finite element algorithm, and derive optimal order error estimates in L 2 -norm (Yuan, 1994(Yuan, , 1996, 19962, 19962).An interpolation computation is introduced to deal the points along the characteristics lying outside the bounded domain.The characteristics intersects the grid boundary and the corresponding values of unknown function should be computed, so the time step should vary due to the position of the grids nearby the boundary along characteristics in program design.In conclusion, the actual computation is most complicated (Yuan, 1996, 19962).
For parabolic equations, Axelsson, Ewing, and Lazarov present upwind finite differences (Axelsson & Gustafasson, 1979;Ewing, Lazarov & Vassilevski, 1994;Lazarov, Mishev & Vassilevski, 1996).It can overcome numerical oscillation and cancel extra interpolation computation of grids nearby the boundary along characteristics.Douglas and Peaceman apply upwind method successfully in incompressible two-phase (water and oil) displacements (Peaceman, 1980).While it is hard naturally to give theoretical analysis.The stability and convergence are derived by Fourier method only for constant coefficient cases and this analysis is not generalized for variable coefficient equations (Douglas & Gunn, 1963, 1964;Marchuk, 1990).Considering actual application, numerical stability and accuracy, the authors present one second-order upwind fractional steps finite difference method for three-dimensional compressible two-phase displacement coupled problem of enhanced oil production.This algorithm can overcome numerical oscillation and dispersion, and decrease the computational scale by decomposition 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 famous problem of Douglas and Ewing.
Generally, the problem is positive definite, Exact solutions of (1)∼ ( 6) are assumed to be suitably smooth, In this paper, M and ε denote a general positive constant and a general positive small constant, respectively, and they may have different meanings at different places.

Method
Let Ω h denote the partition of Ω (see Figure 1.).Let h 1 , h 2 and h 3 be three different spacial steps in x 1 -axis, x 2 -axis and x 3 -axis, respectively, and the grid points is denoted by ] /2, and the signs A n i, j+1/2,k , a n i, j+1/2,k , A n i j,k+1/2 , a n i j,k+1/2 can be defined analogously.Let The fractional steps algorithm of flow equation ( 1) is given by d(C n i jk ) d(C n i jk ) d(C n i jk ) Then the values of Darcy velocity and U n+1 2,i jk , U n+1 3,i jk are obtained similarly.The implicit upwind fractional steps method of saturation equation (2) is considered. where An implicit upwind fractional steps method, second-order accuracy, of components concentration equation ( 7) runs in parallel, D α (c) = ϕ(X)c(X, t)K α (X), where , and similarly for Ūn+1 2,i jk , Ūn+1 3,i jk .Initial approximation: The implicit program runs in the following order.Given i jk } in x 1 direction is first computed by using the method of speedup from (11a) and (11b), then the transitional solution {P n+2/3 i jk } in x 2 direction is computed according to (11c) and (11d), finally the solution {P n+1 i jk } in x 3 direction is obtained according to (11e) and (11f) similarly.Secondly, the values of Darcy velocity {U n+1 i jk } are computed by ( 12).The computation of saturation proceeds later.The transitional solution {C n+1/3 i jk } is computed by using the method of speedup in x 1 direction according to (13a) and (13b), then the transitional solution {C n+2/3 i jk } in x 2 direction is computed by ( 13c) and (13d), finally the solution {C n+1 i jk } in x 3 direction is obtained by ( 13e) and (13f) similarly.Provided {C n+1 i jk }, numerical solutions { Ũn+1 i jk } are obtained continuously.In the last process the values of concentration are computed in parallel for α,i jk } in x 1 direction is computed first by using the method of speedup according to (14a) and (14b), then {S n+2/3 α,i jk } is computed in x 2 direction by ( 14c) and (14d), finally the solution {S n+1 α,i jk } is obtained by (14e) and (14f) in x 3 direction analogously.This finite difference solution of ( 11), ( 13) and ( 14) exists, and is unique according to the positive definite condition.
Theorem 1 Assume exact solutions of ( 1)∼ ( 6) are suitably smooth: And assume that the partition parameters satisfy ∆t ≤ Mh 2 . (16) The difference algorithm ( 11)∼ ( 14) are applied layer by layer, then error estimates hold where ||g|| L∞ (J;M) = sup n∆t≤T ||g n || M , and the constants are Considering the flow equation first, cancelling the transitional solutions P n+1/3 and P n+2/3 by (11a), ( 11c) and (11e), where Eq. ( 18) subtracted from the flow equation (1) (t = t n+1 ), it gives rise to the error equation of the pressure, where An induction hypothesis is given by where Using the calculus of variation, multiplying both sides of error equation ( 19) by δ t π n i jk = d t π n i jk ∆t = π n+1 i jk − π n i jk and summing by parts, then a form in inner products holds Continue to estimate the right terms of (21), The third term of the right side of ( 21) is considered.The first part is discussed here The operators −δ x1 (A n δ x 1 ), −δ x2 (A n δ x 2 ), • • • are self-conjugate, positive definite, bounded and the domain is a unit cube, but their products are not commutative generally.Noting that the difference operators δ are commutative, the first term of the right side of ( 22c) is written by By the induction hypothesis (20) we can get that Applying the positive definiteness of A, d −1 and the decomposition of high order difference operators, we can extract high-order difference quotient term δ x 1 δ x 2 d t π n from the former two terms of the above expression, and cancel related terms by Cauchy inequality, For the third term of (23), Then, Considering the rest of (22c) similarly, The estimates (26) can be obtained analogously for the other two terms of the third right term of (21).
For the fourth term on the right hand side of ( 21), Collecting ( 22)∼ ( 27), it holds for error equation ( 21) as ∆t and ε are sufficiently small, Error estimates of the saturation equation is discussed later.Cancelling C n+1/3 i jk and C n+2/3 i jk of (13a), ( 13c) and (13e), Error equation of the saturation is obtained by (2) (t = t n+1 ) and ( 29) For Then Multiplying both sides of (30) by δ t ξ n i jk = ξ n+1 i jk − ξ n i jk = d t ξ n i jk ∆t and summing by parts, we can get the inner product expression of error equation For the second term of the left side, The other terms on the right side of (31) are considered.The velocity vector U n+1 is bounded by the induction hypothesis (20), It is to estimate the second term of the right side of (31).Note that By ε 0 −Lipschitz continuity and (20), it is derived for the third, fourth, fifth terms and the last term of the right side of (31), The sixth term of the right side of ( 31) is analyzed as follows.
For the terms of the above expression, The terms consisting of δ The first term is discussed first.It is derived that h||U n+1 || 1,∞ is bounded by the induction hypothesis and inverse theorem, and the operators Estimate the other terms of (37b) similarly.
The same result (38) can be obtained for the other parts of the sixth term.
Then the error estimates of components concentration are considered.Cancelling the transition solutions S n+1/3 α and S n+2/3 α of (14a), ( 14c) and (14e), Error equation of components concentration is derived by (7) (t = t n+1 ) and ( 48), where In numerical analysis there exists bound water in oil reservoir, that is to say there exists a positive constant c * such that c(X, t) ≥ c * > 0. It holds as h and ∆t sufficiently small because of the convergence result of c(X, t) (17a), Multiplying both sides of (49) by The left side terms of ( 51) are estimated first, by (17a), The right side terms of (51) are considered later.