Sage Journals: Discover world-class research

Abstract

Because of a high mobility ratio in the chemical and gas flooding for oil reservoirs, the problems of numerical dispersion and low calculation efficiency also exist in the common methods, such as IMPES and adaptive implicit methods. Therefore, the original calculation process, “one-step calculation for pressure and multistep calculation for saturation,” was improved by introducing a velocity item and using the fractional flow in a direction to calculate the saturation. Based on these developments, a new algorithm of numerical solution for “one-step calculation for pressure, one-step calculation for velocity, and multi-step calculation for fractional flow and saturation” was obtained, and the convergence condition for the calculation of saturation was also proposed. The simulation result of a typical theoretical model shows that the nonconvergence occurred for about 6,000 times in the conventional algorithm of IMPES, and a high fluctuation was observed in the calculation steps. However, the calculation step of the fast algorithm was stabilized for 0.5 d, indicating that the fast algorithm can avoid the nonconvergence caused by the saturation that was directly calculated by pressure. This has an important reference value in the numerical simulations of chemical and gas flooding for oil reservoirs.

1. Introduction

The pioneering publication by Bruce and Rice [1] opened the way for the calculations of the vadose issues in oil reservoir by numerical methods. In the past half century, the rapid development in large electronic computers has greatly promoted the application of numerical simulation methods. In the beginning of the 1960s, a multidimensional and multiphase black oil model was developed. In the beginning of the 1970s, the compositional model was studied and applied to the miscible and thermal oil recovery model. In the end of the 1970s, various oil models for chemical flooding were researched. At present, the black oil, mixed-phase, thermal, and chemical flooding models have been developed as commercial software and widely used.

The two main techniques for the numerical solution of a partial differential equation (PDE) are the finite difference method (FDM), finite element method (FEM) [2, 3], and finite volume method (FVM) [4]. The FDM is mainly applied to the problems that depend on time (parabolic and hyperbolic equations) [5]. The common numerical solution methods are implicit pressure explicit saturation (IMPES), implicit alternating solution (IMPIMS), semi-implicit, and full-implicit. Poor stability occurs when a strong nonlinear fluid flow of gas injection/steam coning/high speed gas drive is solved by the IMPES, IMPIMS, and semi-implicit methods. In the full-implicit method, the degree of implicity is high, and the stability is good [6]. However, the Newton iteration method is required to increase the workload and storing capacity in the full-implicit method [7]. In 1980, Pullman and Steger [8] proposed the adaptive implicit method for solving the contradiction between the implicit degree of an equation and the calculated amount characterized using different implicit degrees to process the different grid nodes and time steps, thus reducing the amount of computation and accelerating the computing speed under the same stability conditions. However, its calculation process is still the same as those of the IMPES, IMPIMS, and semi-implicit methods [9]. In 1985, Bertiger and Kelsey [10] proposed an approximate adaptive implicit method, in which the stability was further enhanced. Another direction in the development of the FDM for numerical solution is the multigrid (MG) and preconditioned conjugate gradient (PCG) methods [11]. The basic idea of both methods is that the obtained equation set, which is dispersed, is processed by pretreatment to reach the conjugate gradient, thus accelerating the convergence in the end [12].

At present, chemical flooding technologies such as polymer flooding, movable gel flooding, and combination flooding have been widely used in oil field for EOR. Moreover, because water injection is difficult for an oil reservoir with a low or ultra-low permeability, a series of field tests for the gas flooding such as CO₂ flooding and huff-and-puff have been conducted. Furthermore, the desired effect has also been achieved. However, still some problems exist in the numerical simulation of the chemical flooding or gas flooding with a high mobility ratio such as a large amount of computation and a severe phenomenon of numerical dispersion; therefore, the adaptive implicit algorithm for adjusting the step size is still very difficult to converge the calculation of saturation. To solve this problem, a speed item was added to the calculation process of the IMPES and adaptive implicit methods, which were improved as the “one-step calculation for pressure, one-step calculation for velocity, and multistep calculation for fractional flow and saturation” and achieved the fast algorithm of numerical solution for a strong nonlinear PDE.

2. Fast Algorithm of Numerical Solution for a Strong Nonlinear PDE

The fast algorithm of numerical solution for a strong nonlinear PDE is a sequential solution method. The calculation process is shown in Figure 1.

Figure 1:

Calculation process for the fast algorithm of numerical solution for a strong nonlinear PDE.

2.1. Discretization of Seepage Mathematical Model

The seepage mathematical equations, which describe the changes over time for the internal pressure and saturation through porous medium, are generally nonlinear [13, 14]. The general method for solving this type of equation is to discretize the equation and its definite conditions, in which the numerical method is used. At present, the numerical methods applied in engineering are the FDM and FEM, while the FDM is usually used to solve the problems that depend on time. The fast algorithm of the numerical solution proposed in the following section is specific to the PDE discretized by the FDM.

2.2. Implicit Calculations for Pressure

Using the IMPES method to calculate the pressure, the result can be expressed as follows (1) [5]:

δ p = p^{n + 1} - p^{n} .

(1)

2.3. Calculations for the Total Volume Flow Rate in the l Direction

The total volume flow rate in the l direction, v_t^{n + 1}, can be obtained according to the Darcy law using the obtained pressures, p^{n + 1}, from the calculation. Among them, for the three-dimensional problem, l ∈ {x,y,z}. Considering the gravity, the equation can be expressed as follows [15]:

\begin{matrix} v_{t}^{n + 1} = v_{o}^{n + 1} + v_{w}^{n + 1} + v_{g}^{n + 1}, \\ v_{o}^{n + 1} = - \frac{k k_{r o}^{n} A}{μ_{o}^{n}} [\nabla p_{o}^{n + 1} - (γ_{O}^{n} + γ_{g d}^{n}) \nabla D], \\ v_{w}^{n + 1} = - \frac{k k_{r w}^{n} A}{μ_{w}^{n}} (\nabla p_{w}^{n + 1} - γ_{W}^{n} \nabla D), \\ v_{g}^{n + 1} = - \frac{k k_{r g}^{n} A}{μ_{g}^{n}} (\nabla p_{g}^{n + 1} - γ_{G}^{n} \nabla D) . \end{matrix}

(2)

2.4. Calculation for Fractional Flow and Saturation

Based on the obtained total volume flow rate in the l direction, v_t^{n + 1}, the saturations of oil, water, and gas S_o^{n + 1}, S_w^{n + 1}, and S_g^{n + 1} were calculated. Because of the high mobility ratio, the changes in the saturation for each time step can be very large and the calculation is very difficult to converge. At this time, the calculation of saturation by more steps can be carried out by reducing the step size and dividing the space, Δt, of the moments, n + 1 and n, into the bits of time, m. Moreover, the saturation was calculated by velocity, which also can avoid the numerical diffusion caused by the amplitude of pressure variation, less than the amplitude of saturation variation.

The calculated formula for the mobility by more steps can be expressed as follows [15]:

\begin{matrix} λ_{o}^{n + ((i + 1) / m)} = \frac{k k_{r o}^{n + (i / m)}}{μ_{o}^{n}}, \\ λ_{w}^{n + ((i + 1) / m)} = \frac{k k_{r w}^{n + (i / m)}}{μ_{w}^{n}}, \\ λ_{g}^{n + ((i + 1) / m)} = \frac{k k_{r g}^{n + (i / m)}}{μ_{g}^{n}} . \end{matrix}

(3)

The equations of fractional flow can be obtained using the calculated mobility by the conventional methods [15]

\begin{matrix} f_{o}^{n + ((i + 1) / m)} = \frac{λ_{o}^{n + ((i + 1) / m)}}{λ_{o}^{n + ((i + 1) / m)} + λ_{w}^{n + ((i + 1) / m)} + λ_{g}^{n + ((i + 1) / m)}}, \\ f_{o}^{n + ((i + 1) / m)} = \frac{λ_{w}^{n + ((i + 1) / m)}}{λ_{o}^{n + ((i + 1) / m)} + λ_{w}^{n + ((i + 1) / m)} + λ_{g}^{n + ((i + 1) / m)}}, \\ f_{o}^{n + ((i + 1) / m)} = \frac{λ_{g}^{n + ((i + 1) / m)}}{λ_{o}^{n + ((i + 1) / m)} + λ_{w}^{n + ((i + 1) / m)} + λ_{g}^{n + ((i + 1) / m)}} . \end{matrix}

(4)

Therefore, the volume flow rate in the l direction of the oil, water, and gas phases at time n + ((i + 1)/m) can be expressed as follows:

\begin{matrix} v_{o}^{n + ((i + 1) / m)} = f_{o}^{n + ((i + 1) / m)} v_{t}^{n + 1}, \\ v_{w}^{n + ((i + 1) / m)} = f_{w}^{n + ((i + 1) / m)} v_{t}^{n + 1}, \\ v_{g}^{n + ((i + 1) / m)} = f_{g}^{n + ((i + 1) / m)} v_{t}^{n + 1} . \end{matrix}

(5)

The calculated formula for the final saturation by more steps can be expressed as follows:

\begin{matrix} S_{o}^{n + ((i + 1) / m)} = \frac{S_{o}^{n + (i / m)} V^{n + (i / m)} ρ_{o}^{n + (i / m)}}{ρ_{o}^{n + ((i + 1) / m)} V^{n + ((i + 1) / m)}} + \sum_{l} \frac{v_{o}^{n + ((i + 1) / m)} ρ_{o}^{n + ((i + 1) / m)} - v_{o}^{n + (i / m)} ρ_{o}^{n + (i / m)}}{ρ_{o}^{n + ((i + 1) / m)} V^{n + ((i + 1) / m)}}, \\ S_{w}^{n + ((i + 1) / m)} = \frac{S_{w}^{n + (i / m)} V^{n + (i / m)} ρ_{w}^{n + (i / m)}}{ρ_{w}^{n + ((i + 1) / m)} V^{n + ((i + 1) / m)}} + \sum_{l} \frac{v_{w}^{n + ((i + 1) / m)} ρ_{w}^{n + ((i + 1) / m)} - v_{w}^{n + (i / m)} ρ_{w}^{n + (i / m)}}{ρ_{w}^{n + ((i + 1) / m)} V^{n + ((i + 1) / m)}}, \\ S_{g}^{n + ((i + 1) / m)} = \frac{S_{g}^{n + (i / m)} V^{n + (i / m)} ρ_{g}^{n + (i / m)}}{ρ_{g}^{n + ((i + 1) / m)} V^{n + ((i + 1) / m)}} + \sum_{l} \frac{v_{g}^{n + ((i + 1) / m)} ρ_{g}^{n + ((i + 1) / m)} - v_{g}^{n + (i / m)} ρ_{g}^{n + (i / m)}}{ρ_{g}^{n + ((i + 1) / m)} V^{n + ((i + 1) / m)}} . \end{matrix}

(6)

Inside

\begin{matrix} V^{n + ((i + 1) / m)} = V^{n + (i / m)} + δ ϕ^{n + ((i + 1) / m)} X, \\ δ ϕ^{n + ((i + 1) / m)} = ϕ^{0} C_{R} δ p^{n + ((i + 1) / m)} = ϕ^{0} C_{R} (p^{n + ((i + 1) / m)} - p^{n + (i / m)}) . \end{matrix}

(7)

The obtained formula for Δt, in which the saturation was calculated by the step size of time, can be expressed as follows:

Δ t \leq \frac{V (1 - S_{o r} - S_{w c} - S_{g c})}{\sqrt{3} \sum_{l} v_{t}} .

(8)

2.5. Calculation of Over Bubble Point

A phase change may occur during the calculation for the numerical simulation. The certain phase can be calculated and then the inside unreasonable items of physics can be corrected. Thus, the actual cases for the old and new phases alternating under the conservation of matter conditions could be met; next, the calculation can be carried out according to the new phases. The process can be expressed as follows.

When pⁿ < p_bⁿ, the double phases of oil and water can be changed to the three phases of oil, gas, and water. Then, the parameters can be corrected.

The results before the correction can be expressed as follows:

\begin{matrix} S_{o}^{n} = 1 - S_{w}^{n} \\ S_{g}^{n} = 0 . \end{matrix}

(9)

The results after the correction can be expressed as follows [8]:

\begin{matrix} {\bar{p}}_{b}^{n} = {\bar{p}}^{n} = p^{n} + Δ p, \\ {\bar{S}}_{o}^{n} = S_{o}^{n} + Δ S_{o}, \\ {\bar{S}}_{w}^{n} = S_{w}^{n} + Δ S_{w}, \\ {\bar{S}}_{g}^{n} = Δ S_{g} = - (Δ S_{o} + Δ S_{w}) . \end{matrix}

(10)

The relationship of conservation of matter can be expressed as follows [8]:

\begin{matrix} Δ (S_{o} ρ_{o}) = 0, \\ Δ (S_{g} ρ_{g} + S_{o} ρ_{g d}) = 0, \\ Δ (S_{w} ρ_{w}) = 0 . \end{matrix}

(11)

The speed of calculation in the IMPES method has the following advantages: a simple structure, multistep calculation of flow, more stable saturation, and easier to evaluate the convergence condition. This is suitable for oil and gas reservoir chemical flooding or gas flooding simulation when routine method is in numerical nonconvergence.

3. Stability Comparison of the Algorithm

A typical model for oil and gas reservoirs was esTablished using the fast and IMPES algorithms to perform the simulation calculations. The basic parameters of the typical model for oil and gas reservoirs are listed in Table 1. The bound water is 0.35, whereas residual oil is 0.13. The injected fluid phase is gel, in which viscosity parameters are 0.0067, 3.0E − 5, and 9E − 8.

Table 1:

Basic parameters of the typical model.

Parameters	Value

Rhythm	Reverse rhythm
Permeability range (10⁻³ μm²)	0.01–3,000
K_z/K_x	0.01
Liquid production rate (PV/a)	0.025
Reservoir thickness (m)	20
Viscosity of crude oil (mPa·s)	70
Porosity range (%)	6–35
Well network	Inverted nine spots of well network
Well spacing (m)	350
Number of small layers	4

To test the calculated stability of the numerical algorithm in a high mobility ratio, the permeability of typical model was distributed in the range 0.01–3000 × 10⁻³ μm², and the porosity distribution was in the range 6–35%. Moreover, the maximum calculated step size was limited by 1 day. After the simulated calculations for 4 years (1,440 d), the comparison of the calculated results between the fast and IMPES algorithms by comparing the calculated step size indicates that the IMPES algorithm is influenced by a high mobility ratio. Furthermore, the nonconvergence occurred for about 6,000 times. As Figure 2 shows, the fluctuation in the calculated step size is slightly large. However, the calculated step size of the fast algorithm is stabilized in 0.5 days. This result shows that the simulation stability is improved remarkably because of avoiding the calculation from pressure to saturation and the fast algorithm has accuracy of calculation as shown in Figure 3.

Figure 2:

Comparison of calculated step sizes between the fast and IMPES algorithms.

Figure 3:

Comparison of water cut between the fast and IMPES algorithms.

4. Conclusions

The calculation process: “one-step calculation for pressure, one-step calculation for velocity, and multistep calculation for fractional flow and saturation” was evaluated. A fast algorithm of numerical solution for a strong nonlinear PDE was accomplished.

A typical reservoir model was used to verify that the stability of the fast algorithm was better than that of the IMPES algorithm, and the corresponding convergence criterion for the fast algorithm was also proposed.

The fast algorithm solved the problems of the dispersion for numerical value and low calculation efficiency, which usually occurred in the numerical model of a high mobility ratio for the oil and gas reservoirs of chemical and gas flooding.

For the application of the idea: “one-step calculation for pressure, one-step calculation for speed, and multistep calculation for fractional flow and saturation,” the IMPIMS, semi-implicit, full-implicit, and adaptive implicit methods can also be evaluated by the corresponding fast algorithm.

Footnotes

Nomenclature

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was supported by Major National S&T Projects (2011ZX05009, 2011ZX05054, and 2011ZX05011).

References

Bruce

G. H.

and Rice

J. D.

, “Calculations of unsteady-state gas flow through porous media,” Journal of Petroleum Technology, vol. 5, no. 3, pp. 79–92, 1953.

Karimi-Fard

and Firoozabadi

, “Numerical simulation of water injection in fractured media using the discrete-fracture model and the Galerkin method,” SPE Reservoir Evaluation and Engineering, vol. 6, no. 2, pp. 117–126, 2003.

Bramkamp

Lamby

, and Müller

, “An adaptive multiscale finite volume solver for unsteady and steady state flow computations,” journal of Computational Physics, vol. 197, no. 2, pp. 460–490, 2004.

Bourbiaux

, “Fractured reservoir simulation: a challenging and rewarding issue,” Oil & Gas Science and Technology—Revue de l'Institut Français du Pétrole, vol. 65, no. 2, pp. 227–238, 2010.

Jinfu

and Zhi

, Numerical solution of Partial Differential Equation, Tsinghua University Press, 2004.

Hickel

Adams

N. A.

, and Domaradzki

J. A.

, “An adaptive local deconvolution method for implicit LES,” Journal of Computational Physics, vol. 213, no. 1, pp. 413–436, 2006.

Lang

, “Adaptive multilevel solution of nonlinear parabolic PDE systems. Theory, algorithm, and applications,” Numerical Algorithms, vol. 53, pp. 293–307, 2010.

Pullman

T. H.

and Steger

J. L.

, “Implicit finite-difference simulations of three-dimensional compressible flow,” AIAA Journal, vol. 18, no. 2, pp. 159–167, 1980.

Kaarstad

Froyen

Bjorstad

, and Espedal

, “A massively parallel reservoir simulator,” in Proceedings of the SPE Reservoir Simulation Symposium, SPE-29139-MS, San Antonio, Tex, USA, February 1995.

10.

Bertiger

W. I.

and Kelsey

J. F.

, “Inexact adaptive Newton methods,” in Proceedings of the SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 1985.

11.

Juling

Shangjun

Yang

, “State of the art and developing tendency on reservoir numerical simulation,” PGRE, vol. 9, no. 4, pp. 69–71, 2002.

12.

Hestenes

M. R.

and Stiefel

, “Methods of conjugate gradients for solving linear systems,” Journal of Research of the National Bureau of Standards, vol. 49, pp. 409–436 (1953), 1952.

13.

Jenny

Lee

S. H.

, and Tchelepi

H. A.

, “Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media,” Journal of Computational Physics, vol. 217, no. 2, pp. 627–641, 2006.

14.

Douglas

Furtado

, and Pereira

, “On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs,” Computational Geosciences, vol. 1, no. 2, pp. 155–190, 1997.

15.

Aziz

, “Reservoir simulation grids: opportunities and problems,” Journal of Petroleum Technology, vol. 45, no. 7, pp. 658–663, 1993.