Numerical simulation of two-phase filtration in the near well bore zone

1 Introduction

The western regions of Kazakhstan are rich of oil and gas. Upon that, the deposits differ from each other in the occurrence depth of structural plateaus, reservoir pressures, oil viscosity, sulphur, paraffins and other minerals content. Naturally, these parameters lead to a special relation to each well separately. In this regard, when reclaiming a particular well, it is necessary to develop a separate method of oil production taking into account real natural parameters.

A large volume of experimental and theoretical studies of relative permeability of porous medium has been done by now, including Polubarinova-Kochina P.Ya., Zhukovskyi N.E., Monakhov V.N., Bocharov O.B., Antontsev S.N., Barenblatt I.A., Entov V.M., Zubov N.V., Solovyova V.N., Konovalov A.N [1, 2, 3, 5, 6, 7].

When filtering liquids and gases in formations that are a permeable medium, it is necessary to know the nature of pressure changes at the points of the formation and at its boundaries, and in particular on the walls of the well, and also the flow of formation fluids through any confining surfaces, that is very interesting to drillers from the perspective of assessing the processes of gas-water inflow, absorption, penetration of drilling mud into reservoirs, the nearwellbore zone’s permeability deterioration, etc.

M. G. Alishayev, O. B. Bocharov, V. Ya. Bulygin, M.D.Rozepberg, V. N. Monakhov, A.E.Osokin, L. I. Rubenstein, E.V.Teslyuk, E.B.Chekalyuk’s and others’ works are devoted to mathematical modeling of immiscible liquids’ nonisothermal filtration processes [4, 8].

One of the distinctive features of maintaining and increasing the production of Kazakhstani oil is the development of high-viscosity oil fields, which include the Kenbay deposit, in particular the Moldabek Vostochnyi and Kotytras Severnyi. Experts estimate the total volume of such oil in the world over 810 billion tons. For effective development of these deposits, such reservoir stimulations as polymer water flooding, physico-chemical and microbiological types of stimulation, as well as thermal methods, including in-situ combustion, thermal gas method, steam treatment and hot water injection are used.

It is known that the main problem in the theory and practice of control in the class of objects with continuous technology is the development and creation of a mathematical model on the basis of which the structure and parameters of the control system, the laws of control, as well as the choice and justification of application of technical means of passed laws’ implementation [1].

Methods of mathematical modeling of complex filtration processes in reservoirs is currently being developed in two directions:

1. Development of well-posed models considering the liquids filtration laws in porous mediums;

2. Development of ill-posed models using simplified filtering schemes:

   1. direction leads to necessity of solving complex spatial problems of multicomponent filtration in effective formations, with subsequent implementation by numerical methods.

   2. provides the use of engineering models of calculating in-situ processes resulting from a number of simplifying assumptions, such as the scheme model of a porous medium (Figure 1).

      

      Figure 1

      Model of a porous medium.

In the modern practice of engineering knowledge, or in the so-called modern information technologies, the use of the second direction is common. Developing and building mathematical models of technological processes, including liquid filtration models in a porous medium considering well interference is an actual problem.

The article is devoted to the development and implementation of methods for mathematical and computer simulation of liquid filtration process in a porous medium, taking into account the well interference in oil and gas field development and when metals leaching in uranium deposits.

Computational algorithms have been implemented and packages of applied programs have been created for analysis, control and forecast of technological indicators of deposits in the article.

2 Solution of the filtration problem in the case of multilayer systems

3 Numerical simulation of two-phase filtration

For an axisymmetric two-phase filtration model based on equation

where Q_H – oil production rate;

μ Poisson ratio;

D- well diameter;

hϕ – height of the filter;

k – filter permeability;

δc- cumulative thickness of filter cake and colmatation zone.

A computational algorithm for calculation of the nonstationary problems of mass transfer in the near-well zone has been developed. The model takes into account formation elastic capacity, capillary forces and spatial nonuniformity of the adjacent reservoir’s hydrophysical characteristics distribution.

The equation system (16), considering the conditions of the saturation sum (s_v + s_n = 1) and the capillary jump (p_v = p_n + p_k(s_v)), can be transformed to an equivalent one consisting of a parabolic equation for the pressure (p_v) in the aqueous phase [2]

and the hyperbolic transfer equation for water saturation (s_v)

In addition, boundary and initial conditions are given:

An iterative process has been used to solve this problem, which at the differential level can be written as follows

Here n - iteration step number,

There is pressure in the aqueous phase pvn at the n-th iteration step of the first equation, and then the water saturation svn is determined from the second equation. If the conditions maxrc<r<L|svn−svn−1|<ε1,maxrc<r<L|pvn−pvn−1|<ε2 are met, then iteration process stops.

In numerical calculations by finite-difference method, implicit conservative difference schemes have been used, which have been solved by the sweep method [4]. On each time layer, the difference solution has been found using the described iterative algorithm. This iterative algorithm, when choosing sufficiently small values of time steps, converges and allows to solve the problem in a general formulation.

For an incompressible soil δ = 0, it follows from equation (17) that the total velocity

is only a function of time V = V(t), in this case the computational algorithm can be simplified [2].

The approximate value of the total velocity Vⁿ at the n-th iteration step is determined by the following formula

Here, the function,

which in the case of coefficients dependences of the equations of the following form:

can be computed explicitly.

Further, using the representation

the nonlinear parabolic equation is solved numerically

and the following iteration approximation for water saturation is determined.

The program implementation of the model has been performed on the basis of Visual C ++ 6.0, the program provides a dialog mode of entering the initial information and provides the possibility of graphical interpretation of the calculation results. The program allows solving problems in an axisymmetric and one-dimensional settings, with the conditions of the first and second kind on the right border (r = L). Hydrophysical characteristics of the reservoir can vary depending on the distance from the well, when setting them in the form of piecewise constant functions.

Local changes in the properties of the formation near the well as a result of acid treatment are given according to the formulas:

When solving the problem the results of the calculations are displayed in tables or graphs. After analyzing and modifying the data, the calculations can be continued using the counted results, interpreted as a new initial data.

On the basis of the developed program, calculations were performed that estimate the flow rate of oil entering the extraction well depending on the type of model (with and without capillary forces) and with acid treatment of the wellbore zone.

Coefficients dependences on water saturation were given in the form (18). The calculations were carried out with the following data of an axisymmetric problem.

Constant parameters of the problem: borehole radius r_c = 0.1; void structure m = 0.15; filtration coefficient k₀ = 0.2m/day; field length L = 30m; initial water saturation s_v = 0.4; the ratio of viscosities of water and oil μ₀ = 0.2; lowering the pressure in the well compared to the reservoir Δp=pv1−pv0=20(pv0=0).

The following parameters were varied: the parameter in the capillary jump formula (pk0) and the compressibility coefficient of the reservoir (δ).

The results of calculations of the untreated well flow are given in Table 1, where the numerator corresponds to the flow rate of the oil inflow, the denominator is the water flow rate.

Numerical modeling of two-phase flow while acid treatment of the well was carried out at the hydrophysical parameters of the formation δ = 0.0001, pk0 = 0, pk0 = 0.8, and the treated well bore zone a₀ = 0.05, λ = 1, R_* = 2 (Table 1).

The increase in the flowrate the treated well is observed both for the model taking into account the capillary forces (pk0 = 0.8), and without taking them into account (pk0 = 0). The last two columns of the Table 2 show the increase in the costs of acid treatment (the ratio of relevant costs with acid treatment and without it), respectively, the numerator qn0 -for oil, the denominator qv0 - for water.

Figures 2, 3 show the graphs of the distribution of oil saturation and pressure for option parameters δ = 0.0001, pk0 = 0.8 without acid treatment.

Figure 2

Distribution of oil saturation

Figure 3

Pressure distribution

It should be particularly noted the effect of compressibility of the reservoir on the two-phase flow, which qualitatively changes its character.

In this case, compared to rigid filtration mode the dynamics of the pressure change and its distribution in the reservoir significantly influence the appearance of the solution.

This is especially clearly seen in the absence of the capillary jump. The Buckley-Leverett model at δ = 0, homogeneous flooding of the reservoir at the initial time moment (s_t=0 = s(r, 0) = s₀ = const) and the boundary condition s_r=L = s₀ has a solution in the form of a constant sv(r, t) = s₀. However, the solution has a completely different form for the compressible soil: water saturation and oil saturation are non-monotonic functions that depend on time.

Figure 4 shows the calculation results for option pk0 = 0.0, δ = 0.005.

Figure 4

Oil saturation distribution

The effect of compressibility of the reservoir is also influenced by considering capillary jump, especially for large values of the compressibility coefficient of the reservoir (δ = 0.005), which may be due to the presence of the gas phase in the reservoir. In this case, when the pressure surge wave passes, a local maximum of the oil saturation distribution is formed near the borehole (the graphs are shown in Figure 5 (pk0 = 0.8) and, in addition, at the initial moments of time near it, the pressure gradients are higher (Figure 6 (pk0 = 0.8)). These factors allow to explain higher values of oil inflow into the well at the initial moment of time (Table 3) compared to the results obtained for slightly compressible soils δ = 0.0001.

Figure 5

Oil saturation distribution

Figure 6

Pressure distribution (pk0 = 0.8)

From the series of calculations performed, it is possible to draw a general conclusion.

Acid treatment of the near-well zone of the production well, under other equal conditions, can be an effective method of increasing the flow of oil. Naturally, the effectiveness of this method depends on processing technology and physicochemical parameters of reservoir formations. For example, the results presented in table 3, it is clear that a significant role in increasing the inflow of oil (≈1.5 times) plays reducing the negative impact of capillary forces.

The implementation of hydraulic fracturing Main form

The possibility of using the parallel sweep method for well systems is considered in the article. The proposed system of equations includes, firstly, a system of difference equations defined on each segment as a result of approximating differential equations, in the sense of Stefan and Verigin, and, secondly, equations defined at the vertices, which can be considered as boundary conditions.

The mathematical models under consideration were investigated in domains with complex geometry. For the comparative analysis, the following solutions to the problem of determining the radius of displacement of formation fluid were used when technological measures aimed at reducing the depth of penetration of the drilling fluid filtrate during the opening of the productive formation. We use the basic filtration equations describing the laws of fluid motion in the reservoir for the solution

We write down the equation of continuity in the case of the radially symmetric motion of the filter (liquid), assuming that the filtrate penetration process occurs without diffusion mixing; an expression for the displacement radius is obtained:

If the fluid flow is a known function, then regardless of the type of formation and the nature of the filtration, we have:

where h - bed formation thickness,

R_ϕ - radius of the displacement front of formation fluid, m - void structure,

Q(t) - flow rate at r = R_c (r - current formation radius, Rc - borehole radius),

W(T)=∫0TQ(t)dt -liquid volume entering to the reservoir in a time.

Usually (for example, during repression on a reservoir), the flow Q(t) is an unknown function, depending on the pressure drop Δp(t). In this case, it is necessary to determine the flow rate Q(t) from the solution of the corresponding non-stationary filtration problem:

where ε=khμ,χ=k2mμ,

k - permeability

μ - fluid viscosity

S - degree of growth of the surface resistance at S > 0 or its decrease at S < 0

If Δp is a variable of time values, we can use the approximate formula

where the time intervals Δt_j = t_j − t_j−1 must satisfy the condition Δt>>Rc2/4χ.

The formula s=(k2k1−1) ln R1Rc in the case ΔP = ΔP₀ = const takes the form:

where a=1.8kΔP0mχμ,τ=2.25χtR02,τ∗=τe2s,

L_i(τ_ϕ) = ∫0τϕdulnu integral logarithm, the tabulated function

At τ_ϕ >> 1, using the asymptote

we obtain a formula for finding the front radius of the computation in time:

If the flow rate is calculated using formula (20), we proceed as follows.

By value, the parameters R_c/4_χ choose a time step Δ>>Rc24χ, satisfying condition, compile the table of values ΔP_u = P(t_u) − P_m for corresponding moments of time t_u = nΔt₀ and for a given value δ we compute sequentially Q_u = Q(t_n).

Then the time interval Tis divided into n equal intervals, where it is determined from the condition n>>[4TΔERc] ([] - is the integer part of the number E). Then the function W(t) is approximately computed by the formula:

If the viscosity of the solution filtrate is not equal to the viscosity of the fluid, then in determining the flow rate Q(t) it is possible to use the quasi-steady-state method by which the filtrate flow is calculated:

where ε1=khμ,ε=khμ - Hydraulic conductivity in the zones of reservoir fluid and reservoir filters, respectively

R_k - radius

R₁ - invasion range of the filtrate into the formation.

In the formula, ∂2P∂ξ12+∂2P∂ξ22=0, assuming R₁ = R_ϕ, R_k = R_t, we write it in the form:

where μ¯=μ1μ,Rt=2.25χt,Rϕ(t)<Rt.

Substituting (21) into the right-hand side of the expression (19) we obtain the differential equation with respect to

Calculations show that with increasing parameters μ and, S, invasion range of the filtrate decreases significantly, for example, when S ≥ 10³ it is practically independent of the viscosity of the formation fluid.

Practical interest represent the following limit cases. If at the opening of the permeable formation as the drilling fluid to use a liquid containing no colloidal particles (water, oil) and its viscosity is equal to the viscosity of the fluid, then, putting the equation, we obtain the following calculation formula:

Therefore, in this case, the front radius of the computation depends, in addition to the process parameters F and T only on the properties of the formation, for example, at K = 2 · 10⁻¹³, m², m = 0.15, μ = 10⁻³PA, ΔP = 5 MPa and T = 8.64 · 10⁻⁴ s will receive R_ϕ = 13.5 m.

The equation (22) can also be used to estimate the penetration radius of a grouting mortar when fixing the absorbing horizon or acid solutions during treatment of the bottom hole zone if the colmatation zone has been destroyed. A simple calculation from the formula (22) shows that if the clay crust is removed or destroyed during the cementing of the well, the penetration radius of the filtrate of the component solution can be quite significant.

The results of numerical simulation are shown in Figure 7-12. To illustrate the results, a set of programs with real technological indicators of the above deposits was created. A series of calculations were performed to adapt mathematical models.

Figure 7

Test examples of calculating the fracture volume

Figure 8

Evaluation of the cooled thermal front of the formation

Figure 9

Changing the cooled thermal front of the formation

Figure 10

Baring of the well bore zone

Figure 11

Reservoir response

Figure 12

Baring of the well bore zone

At the same time, data on the history of field development were used. Similar packages can be used for uranium deposits i.e. their technological parameters with the parameters of the Kenbai deposit are in many respects identical.

For pressure calculation during the well bore zone bearing of the formation, the window for parameter entry is shown in Figure 10, and the solved test example is shown in Figure 11.

4 Conclusion

For the numerical implementation non-structured meshes were used for the plane-parallel flow and axisymmetric cases, in which a sufficiently thick grid is closer to the producing well. The integration step over time is determined from the generalized Courant inequality. And, numerical results are presented in the form of graphical dependencies, from which the sufficient effectiveness of the developed numerical solution algorithms follows.