MODELING OF IMPACT OF HYDRAULIC FRACTURES ON THE PROCESS OF FLUID DISPLACEMENT FROM LOW-PERMEABILITY SEDIMENTARY ROCKS

In oil and gas field development, there are a number of factors that influence the production efficiency. These are both geological characteristics and technical parameters of the reservoir [1]. Much experience in using intensive field development systems has been gained. This applies primarily to pattern waterflooding [2], in which production and injection wells are arranged in a certain way within the corresponding areas. As it is known, the use of the hydraulic fracturing technology [3, 4] is reasonable in the design of low-permeability (shale) sedimentary rocks and also due to the deterioration of reservoir properties in near-wellbore regions in the course of reservoir development. As a result, fractures extend the area of influence of production wells and form associations with high-permeability zones. Based on numerical methods for quasiconformal mappings [4–6], the algorithm for solving nonlinear boundary value problems of single-phase filtration in low-permeability sedimentary rocks in the pattern waterflooding elements considering the impact of hydraulic fractures was developed. The algorithm allows predicting the properties of the reservoir system under various impacts and studying the features of filtration in near-wellbore regions. The improvement and development of numerical methods for quasiconformal mappings for mathematical modeling of nonlinear displacement processes in oil reservoirs considering the impact of hydraulic fractures is an urgent issue. This would allow determining the time points of the displacing fluid breakthrough to production wells and complete waterflooding, the coordinates of critical “suspension” points and their quasipotential values, fluid interface position at different time points, the overall filtration rate of production wells, oil fraction dependence, the volume of the fluid displaced in the reservoir within a certain time and, accordingly, the volume of the remaining fluid and so on.


Introduction
In oil and gas field development, there are a number of factors that influence the production efficiency.These are both geological characteristics and technical parameters of the reservoir [1].Much experience in using intensive field development systems has been gained.This applies primarily to pattern waterflooding [2], in which production and injection wells are arranged in a certain way within the corresponding areas.
As it is known, the use of the hydraulic fracturing technology [3,4] is reasonable in the design of low-permeability (shale) sedimentary rocks and also due to the deterioration of reservoir properties in near-wellbore regions in the course of reservoir development.As a result, fractures extend the area of influence of production wells and form associations with high-permeability zones.
Based on numerical methods for quasiconformal mappings [4][5][6], the algorithm for solving nonlinear boundary value problems of single-phase filtration in low-permeability sedimentary rocks in the pattern waterflooding elements considering the impact of hydraulic fractures was developed.The algorithm allows predicting the properties of the reservoir system under various impacts and studying the features of filtration in near-wellbore regions.The improvement and development of numerical methods for quasiconformal mappings for mathematical modeling of nonlinear displacement processes in oil reservoirs considering the impact of hydraulic fractures is an urgent issue.This would allow determining the time points of the displacing fluid break-through to production wells and complete waterflooding, the coordinates of critical "suspension" points and their quasipotential values, fluid interface position at different time points, the overall filtration rate of production wells, oil fraction dependence, the volume of the fluid displaced in the reservoir within a certain time and, accordingly, the volume of the remaining fluid and so on.

Literature review and problem statement
Many scientists around the world to investigate the process of fluid filtering to wells in the presence of hydraulic fractures.In particular, the analytical solution of the relevant boundary value problem, where fractures are presented in the form of a section of zero thickness and finite conductivity is given in [7].Such a model does not reflect the actual filtration properties of the displacement process, unlike the case where the fracture is modeled in an ellipse form, as in [8].A more complex model is proposed in [9], which investigated deviations of fractures, depending on the pressure generated by existing microfractures in shale sedimentary rocks.The study of the impact of arrangement of several hydraulic fractures on one production well on the displacement process is proposed in [10].However, there is the problem of finding a saturation field, which would allow predicting the rate of waterflooding of production wells and identifying the features of operation of a field under the projected arrangement of wells and hydraulic fractures on them.

MODELING OF IMPACT OF HYDRAULIC FRACTURES ON THE PROCESS OF FLUID DISPLACEMENT FROM LOW-PERMEABILITY SEDIMENTARY ROCKS А . B o m b а
Doctor of technical sciences, Professor, Head of the Department* Е-mail: abomba@mail.ru

A . S i n c h u k
PhD, Lecturer* Е-mail: sinchukk@mail.ru*Department of Computer Science and Applied Mathematics Rivne State Humanitarian University S. Bandery str., 12, Rivne, Ukraine, 33028
The problems of optimization of the size, permeability coefficient and arrangement of hydraulic fractures have been examined in [13].The problems of optimization of filtration characteristics of the displacement process in the presence of hydraulic fractures have been solved in [14].However, these publications did not analyze the locations of so-called "stagnation" zones, depending on the set parameters and arrangement of fractures.
Thus, there is a need to solve a wider range of problems.In addition to the optimum arrangement of injection and production wells, it is necessary to identify efficient arrangement of hydraulic fractures in the vicinities.This would satisfy certain criteria, including by selecting the parameters of hydraulic fractures under constant parameters of the reservoir to achieve the maximum time of water breakthrough to production wells and amount of extracted oil, and the minimum water flow rate.

Research goal and objectives
The goal of the research is mathematical modeling of fluid displacement processes in oil reservoirs considering the impact of hydraulic fractures, and development of numerical methods for quasiconformal mappings for solving the relevant boundary value problems of single-phase filtration.
To achieve the goal, the following tasks were set: -to enhance the mathematical model of fluid displacement from low-permeability (shale) sedimentary rocks in the pattern waterflooding elements in the presence of hydraulic fractures; -to develop a methodology for solving boundary value problems of filtration processes of displacement from lowpermeability (shale) sedimentary rocks considering the impact of hydraulic fractures and related deformation processes in the near-wellbore region when the process under study is described by specially modified Darcy's law regarding the critical value of the pressure gradient under quasistationary filtration flow; -to develop numerical algorithms for solving relevant boundary value problems, to conduct numerical calculations and analysis of the results on this basis.

The method of comprehensive analysis of modeling of nonlinear processes of displacement in oil reservoirs considering the impact of hydraulic fractures
Let us consider the process of single-phase isothermal filtration in horizontal reservoir bed, generated by dou-bly-symmetric rectilineal rows of injection and production wells, riddled with finite-permeability hydraulic fractures (Fig. 1) without overflows between the respective rows.
Considering the symmetrical arrangement of wells in the reservoir, we have the opportunity to allocate the element ∈ z z G G , containing n * injection wells and one production well with corresponding fractures and their symmetrical parts (Fig. 2, where = + 0 * d n (r a) -the distance between the separating lines of the symmetry elements, = * n 3, 0 r -the radius of wells, a -half the distance between the injection wells, h -distance between rows).
-coefficient of absolute permeability of the soil, where the reservoir area that corresponds to the κ-fracture, κ = 1, 2, 3..., ( κ = k const).Hydraulic fractures are simulated by fragments of ellipses with semiaxes κ a , κ b and the corresponding angle of directionκ α ; L g , L -boundaries of injection and production wells, respectively, χ -the coefficient characterizing the dependence of permeability of sedimentary rocks (in complicated geological conditions of filtration, for which µ 0 k is a small size) on the pressure gradient value and is determined by the following ratio: where F -monotonically increasing function, I kr -critical value of the initial gradient.

Fig. 2. The symmetry element of the reservoir under pattern waterflooding
To construct an approximate solution of the problem, we introduce the velocity quasipotential in the form of the Laybenson function [5]: and rewrite the equation ( 1) with the corresponding boundary conditions: where Similarly to [4], by introducing the flow function ψ, complex conjugate to ϕ, the problem of constructing the hydrodynamic grid, determining filtration rate and other specific filtration parameters by the found (fixed at a given time) saturation field is reduced to the quasiconformal mapping x,y i x,y of a simply connected region G z on the corresponding area of complex quasipotential where ( ) where l=1, 2, 3..., The problem of finding the saturation field (2), according to [4], can be presented as follows: where the equation ( 6) is actually spatially one-dimensional, because the variable ψ appears as a parameter.
The difference analogue and the solution algorithm are built as in [4].At the initial stage, we find the parameter ϕ 1 H , then we consistently solve a series of intermediate problems, corresponding to Fig. 3.

Fig. 3. The region of complex quasipotential
The nodes ϕ ψ i j ( , ) of the grid area ω G are determined as follows: n n n 1, m ,n ,n N.
The equation ( 4) is approximated using the finite volume method [15] as follows: where ( ) ( ) Approximations of boundary conditions can be written as: f (x ,y ) 0, j 0,m, x 0, 2(g 1)(a r ) a y y , x 0, y y 2(g 1)(a r ) a, i 0,n , f(x ,y ) 0, f(x ,y ) 0, i 0,n, j 0,m, g 2,n .( 9) Here, as in [6], complex conjugation of harmonic functions = ϕ ψ i,j i j x x( , ), = ϕ ψ i,j i j y y( , ) is provided by the conditions of orthogonality of near-boundary normal vectors relative to corresponding tangents along the boundary of the region z G .Their difference analogues on the well boundaries are as follows: (4y 3y y )(y y ) 0, j m ,m , ( Unknown approximate values of flow rate g Q and potential ϕ g H in the flow divergence points in the process of iterations are found by the formulas: where ∆ϕ γ + ∆ϕ γ ∆ψ = γ γ 2 , and γ g l is obtained under "quasiconformal similarity in small" of respective elementary quadrangles of two regions: The equation ( 7) is approximated by the "upwind" differencing scheme [4] as follows: where τ -time step, i,j s , i,j s -saturation at the corresponding time points, υ i,j -velocity.Boundary and initial conditions for saturation in the grid area are as follows: = 0,j * s s , = i,j i,j s(x ,y ,0) i,j i,j s(x ,y ), = j 0,m, = i 1,n.By setting the step τ, the parameters of partition g 1 n , g 2 n , m g = * g 1,n , of the region ω G and the accuracy ε 1 , ε 2 of the algorithm, the initial approximations of coordinates of the boundary nodes (so as to fulfill the condition ( 8)) and the initial approximations of coordinates of internal nodes ( (0) i,j x , (0) i,j y ) according to the formulas (10), we find approximations of values γ g l .Then, refinement of coordinates of internal nodes of the hydrodynamic grid by solving (7) with respect to x i,j and y i,j is performed.After that, we correct the boundary nodes with the surrounding boundary and near-boundary nodes fixed using the orthogonality condition, and find the approximation of values Q g , ϕ g H .The conditions for completion of the construction (finding unknown filtration parameters, including the velocity field) algorithm of the hydrodynamic grid in this iterative phase are: stabilization of the flow rate g Q ( and so on.In the event of non-compliance with at least one of the conditions, the regions of violation of quasiconformality on the hydrodynamic grid are observed.According to (10), we find new distribution of saturation in the reservoir and repeat the steps of the algorithm using the velocity field and saturation field from the previous iteration step in time (considering the boundary conditions).According to the calculations, the formulas (4) and ( 5), the volume of oil in these elements of symmetry before the displacement -V=7.106, and the volume of the oil produced and the residue thereof in the reservoir during the process is V=3.965,V(18.62)=3.141(case 1), V=3.428, V(18.62)=3.678(case 2), V=3.789, V(18.62)=3.317(case 3), respectively.In view of the study (computer experiments), we see that the oil withdrawal volume in total filtration rate in each case is reduced differently after reaching the time of displacing reagent breakthrough to the production well.This is due to a significant difference in phase viscosity coefficients, relative phase permeability, parameters and locations of hydraulic fractures.In case 1, there is rapid waterflooding of the production well through one of the fractures and there is a risk of so-called of oil stagnation zones.In case 2, that risk decreases considerably, but the time for complete oil displacement increases.In case 3, we observe too slow course of the oil withdrawal process ("proximity" of withdrawals, especially for cases 1 and 3, is due to equal areas of hydraulic fractures).Also, the fact is confirmed that the "transverse direction" (with respect to injection wells) of hydraulic fractures accelerates the time of the displacing reagent breakthrough to the production well (although provides some growth of oil withdrawal values at the initial stages), and their "longitudinal" direction reduces the number of oil stagnation zones.

Numerical calculations of the model problem
At the same time, we emphasize that oil stagnation zones in these cases are close to the so-called stagnant zones (areas of reservoirs at the points of which the gradient value is less than some critical value).
The mathematical modeling of low-permeability oil field development allows predicting the rate of waterflooding of production wells and identifying the features of operation under the projected arrangement of wells and hydraulic fractures on them.In these conditions, the problem of optimizing the oil withdrawal process, depending on the given parameters of hydraulic fractures on the production well and determining the location of stagnant zones is solved.

Conclusions
The mathematical model of fluid displacement from the low-permeability oil fields in the pattern waterflooding elements in the presence of hydraulic fractures, particularly considering the coefficient that characterizes the dependence of permeability of sedimentary rocks on the pressure gradient value is improved.
Numerical methods for quasiconformal mappings are extended to solve nonlinear boundary value problems of single-phase filtration in low-permeability (shale) sedimentary rocks in the presence of hydraulic fractures.At the same time, the impact of related deformation processes in the near-wellbore region is considered.That is, under quasiconformal filtration flow, the process under study is described by specially modified Darcy's law regarding the critical value of the pressure gradient.
Numerical algorithms for the calculation of filtration characteristics: saturation field, velocity quasipotential, time of the displacing fluid breakthrough to the production well and its complete waterflooding are developed.The algorithm also allows determining the coordinates of the critical "suspension" points and their quasipotential values, fluid interface position at different time points, the overall filtration rate of the production well, the dependence of oil fraction in it.For an effective analysis of the research, calculations of the volume of the fluid displaced in the reservoir within a certain time and the volume of the remaining fluid at an arbitrary time are performed.

Fig. 1 .
Fig. 1.The diagram of reservoir development with the allocated symmetry element (( ) -injection well, ( ) -production well) When modeling this process, the law of motion and the flow continuity equation, according to [4-6], is represented as: χ υ = − µ kr k (I,I ) grad p, ρυ = div 0, flow rate of production wells.Inverse to (1) boundary value problem of quasiconformal mapping( ) ( ) ( ) = ω = ϕ ψ + ϕ ψ z zx , iy , of the region ω G on G z , and, consequently, the equation for the real parts of the characteristic flow function is written as:

Fig. 4 . 3 Fig. 6 ,
Fig. 4. The hydrodynamic grids of the symmetry elements: а -in case 1; b -in case 2; c -in case 3 Fig.6, 7 show the dependencies of the total filtration rate ( ) * n Q t and oil withdrawal values

Fig. 7 .
Fig. 7.The graphs of dependence of the oil withdrawal value on the time in the corresponding symmetry elements