NUMERICAL SOLUTION OF THE PROBLEM OF GAS FILTRATION IN POROUS MEDIA BY THE METHOD OF COORDINATE SPLITTING

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract. The article presents a three-dimensional mathematical model and numerical calculation algorithm based on the method of coordinate splitting for analysis and forecasting the process of gas filtration in porous media. The use of the method of coordinate splitting made it possible to reduce the original problem to three locally onedimensional problems and to reduce the computation time by reducing the number of computational operations. Developed mathematical apparatus can be successfully used to solve design problems or refinement of design solutions of gas fields, as well as production forecasting process.


INTRODUCTION
It is known that the basis of the real sector of the economy of any country in the world is the fuel and energy complex. The level of its development, among other things, is determined by the development index of a particular state in the world community.
The increasing volumes of consumption and consumption of fuel and energy resources can be covered not only through the development and design of new oil and gas fields, but also through the rational operation and effective management of existing fields.
Despite the fact that to date, significant theoretical and applied results have been obtained in the field of hydrodynamic studies of processes occurring in porous reservoir media, and the influence of various factors on them (filtration coefficient, porosity and thickness of the reservoir, oil and gas recovery coefficient, viscosity, etc.), this scientific direction still has many unsolved problems.
In particular, the paper [1] investigates a system of partial differential equations describing a two-phase flow in a porous medium and consisting of hyperbolic terms for advective transport, parabolic terms for dissipative effects, and relaxation nonequilibrium equations. The authors found that for several dissipative and nonequilibrium systems, the use of the stream function as a free variable instead of time divides the general system (n + 1) × (n + 1) into an auxiliary n × n system and one scalar lift equation. In many cases, when the auxiliary system admits an exact solution, the general flow problem is reduced to the numerical or semi-analytical solution of one lifting equation.
A similar technique is used in [2], where the division of the mathematical model of the enhanced oil recovery (EOR) method into thermodynamic and hydrodynamic parts is discussed.
The n × n conservation law model for two-phase flows is transformed into an (n − 1) × (n − 1) auxiliary system containing only thermodynamic variables and one lifting equation containing only hydrodynamic parameters. The authors show that splitting makes it possible to prove the independence of phase transitions that occur during a shift of phase relative permeabilities and viscosities. Consequently, the minimum miscibility pressure of the gas flooding providing insulation surfactant from produced water, can be calculated using condensed auxiliary system. Moreover, the very fact of reducing the number of equations allows the generation of new analytical model for EOR.
The authors of [3] consider the construction and analysis of algorithms for splitting operators with a large time step for the numerical simulation of multiphase flows in heterogeneous porous media. The authors present some results of solving quasilinear degenerate parabolic equations on two-dimensional and three-dimensional numerical test examples. The main conclusion drawn from numerical experiments is that the operator separation algorithms do demonstrate the property of accurate resolution of inner layers with steep gradients, give very little numerical diffusion, and, at the same time, allow the use of large time steps. In addition, these algorithms seem to cover all possible combinations of convective and diffusion forces, ranging from convection-dominated problems (including the purely hyperbolic case) to diffusion-dominated problems.
The study [4] addresses three types of two-level agreed splitting algorithms for nonstationary Navier-Stokes equations. The basic technique is to first solve the nonlinear problem in the subspace coarse level and then solve the linear equation in a subspace of a shallow level. The authors conclude that two-level methods can save time for solving problems and the number of computational operations in comparison with single-level methods. The stability and convergence of the schemes show that two-level methods can achieve optimal accuracy with the correct choice of coarse and fine mesh scales. Numerical examples show that Stokes' correction is the simplest, Newton's correction has the best accuracy, and Oseen's correction is preferable for problems with a large Reynolds number and simulation for a long time among the three methods.
The applied problems of oil and gas production are numerically investigated in [5] using mathematical models of the flow of multiphase fluids in porous media. The basic model includes continuity equations and Darcy's laws for each phase, as well as an algebraic expression for the sum of saturations. The authors highlight the main properties of the pressure problem and discuss the need for their implementation at a discrete level. The resulting elliptic problem for the pressure equation is characterized by a non-self-adjoint operator. This problem is solved by the authors approximately using iterative methods. Particular attention is paid to numerical algorithms for calculating pressure on parallel computers.
On the basis of additive schemes, effective numerical algorithms for the approximate solution of initial-boundary value problems for systems of nonstationary partial differential equations are constructed in [6]. The authors note that usually additive operator-difference schemes for systems of evolutionary equations are constructed for operators that are connected in space. Therefore, the authors focused on the study of more general problems in which there is a connection between the time derivatives for the components of the solution vector. The splitting schemes proposed in this paper are based on a triangular two-component representation of operators.
The authors of [7] describe a method based on splitting into physical processes rather than the dimensionality of the problem in the form of a nonlinear system of partial differential equations, which simulates a three-phase flow through heterogeneous porous media. In the numerical algorithm for solving the problem of convective transfer of liquid phases, a conservative central difference scheme of the second order of accuracy is used. This scheme is combined with locally conservative mixed finite elements to numerically solve parabolic and elliptic problems related to diffusion transfer of liquid phases and the pressure-velocity problem. The authors numerically investigated the existence and stability of nonclassical shock waves in two-dimensional inhomogeneous flows. The results of numerical experiments show that the proposed operator splitting technique provides computational efficiency and accuracy.

STATEMENT OF THE PROBLEM
In order to investigate and more adequately describe the process of gas filtration in porous media as well as determine the main indicators used in the development of oil and gas fields, we introduce a 3D mathematical model described by the equation and in special nodes (wells) In the equations (1) -(2) Substituting equation (3) in (2) and taking into account the variability of the capacity of the reservoir we will get [8,9,10,11,12,13,14]: Here Q -volumetric flow rate (at atmospheric pressure) in the wells, Qρ -mass flow, Ppressure; P at -atmospheric pressure, ρ -density, b -power of the stratum,b -the average power value in the grid "square", ∆x, ∆y, ∆z -steps by x, y and z coordinates; m -reservoir porosity; K, µ -respectively the filtration coefficient and viscosity of the gas, K z = f (m, g), γ ν -the set of points of the region G, in which wells may be present.
Let's assume that the gas is ideal and ρ = const · P . Then we obtain Equation (5) is valid for any law of filtration and any dependence on the density of the pressure.
If in the equation (5) all coefficients are constant, i.e., with the corresponding boundary conditions: As a result, the final form of the equation was obtained which gives us possibility to study the filtration process of any considered component in porous medium in order to determine the main parameters of the development and design of oil and gas fields.

SOLUTION OF THE PROBLEM
To solve problem (6) -(9), using the relations x = x * L; y = y * L; z = z * H; P = P * P 0 ; we bring equation (6) to a dimensionless form: To linearize equation (10) using the expression P 2 = 2P ·P − P 2 , we get Since the problem (7) -(10) is described by nonlinear partial differential equation with the corresponding internal and boundary conditions, it is complicated to obtain an analytical solution.
To solve equation (11) is necessary to set an effective approximation of the boundary condition in the well. Of great practical interest is the case when the boundary conditions are set with a known gas production at the well. To find the gas-dynamic pressure field, it is necessary to set the well radius equal to R = R c +W.
Here W --is determined from the relation ln W ∆x = const, where const = 1.57. Thus, the resulting pressure field corresponds to the operation of large wells with a radius of R o = R c + 0.208 ∆x at an appropriate flow rate.
It should be noted that the flow rate of a well with a radius is unknown, therefore, for each moment in time, it is determined in the process of solving the problem of unsteady gas filtration to an actual well with a radius of R c .
Since movement in the bottomhole formation zone occurs with violation of Darcy's filtration law, the drag law, which generalizes the linear filtration law, can be taken in the form Here ω, u, w -components of the filtration speed in the direction of the x, y, z axes; ρgas density; β * --coefficient that is introduced when describing the two-term filtration law.
When solving the problem of gas filtration in a porous medium, the main problem arises in the course of approximating the boundary conditions at wells by finite-difference relations.
To do this, select a volume with an axis of symmetry passing through a fictitious point well.
Consider the balance equation for this volume, the horizontal projection of which is a square with side ∆x = ∆y = ∆z = const. The summation of the amount of flow through the side faces in the direction of the ∆x, ∆y and ∆z axes is calculated using the expression Expression (12) for an ordinary node should be equal to where ρ and bmean values of the function in the area x i−0.5 ≤ x ≤ x i+0.5 ; y j−0.5 ≤ y ≤ y j+0.5 ; z k−0.5 ≤ z ≤ z k+0.5 .
For a special node, this value is added to the production at the well with the contour Γ κ , equal to (Q κ ρ), and Q k is determined by the formula γk Kb µ ∂ P ∂ n P P at ds = Q κ (t).
Here n --intrinsic normal to the y curve γk; ds --an arc element of a given curve.
The finite difference equation approximating the boundary condition for a special mesh node with mass flow rate Q κ ρ has the form V x (bρv) i, j,k ∆x ∆t +V y (bρv) i, j,k ∆y ∆t +V z (bρv) i, j,k ∆z ∆t = Now dividing both sides of equation (13) by b, and assuming that gas filtration in a porous medium obeys Darcy's law, we obtain 1 bV x bρ K µV x P i, j,k From equation (14) it follows that the expression Q κ ∆x ∆y∆z b -is the gas flow rate per unit volume.
For effective numerical integration of the problem on a computer, we split it in the given variables (x, y, z) and obtain three problems: The first task, where gas filtration is considered along the OX axis: where the initial calculation time P 1,0 = P 0 , P n+1 3 in time is determined at the next stages of the problem. Solving this problem, we find P 1,0 = P n+1 3 . The second task in the direction of the OY axis: The third task in the direction of the OZ axis: It should be emphasized here that the solution obtained as a result of integrating the first problem is used as an initial condition for solving the second problem, and the solution of the second problem, respectively, is used as an initial condition for solving the third problem, the obtained solution of the third problem is used as an initial condition for the first problem for integration at the next time level.
To solve problem (15) - (20), using the finite-difference method, we replace the differential operators with the finite-difference ones, we get: First task: We approximate (15) by a finite-difference scheme in Ox, we obtain [8,9,10,11,12,13,14,15,16,17,18,19,20]: We transform equation (21) to the following form: then grouping similar terms, equation (22) is present in the form of a system of algebraic equations: (23) a i, j, k (P 1 ) Replacing the boundary condition for finite difference we get To solve the system of equations (23), we use the sweep method. To do this, we present the following recurrence relation: where α i, j,k , β i, j,k -sweep coefficients.
Thus, a numerical algorithm has been obtained for solving the problem of gas filtration in porous media, with the help of which it is possible to carry out computational experiments to determine the main indicators of the development of gas fields.

RESULTS AND DISCUSSION
On the basis of the proposed mathematical apparatus, there was developed a software tool for carrying out computational experiments. Then a series of calculations was carried out using two numerical methods for different values of individual parameters of the object under study.
Based on obtained results for the coordinate splitting method based algorithm, the isobar maps were compiled for years and arbitrary filtration area (Fig. 2, 4 and 6). Then they were compared with the results of calculations using the method of alternating directions (Fig. 1, 3 and 5).
Computational experiments showed that the lower the value of the formation thickness, the faster the pressure redistribution in the filtration area occurs ( Fig. 1-6). As can be seen, the results of the alternating direction method (Fig. 1, 3 and 5) and the coordinate splitting method       Although the difference between the results of these two methods is quite small, nevertheless, the advantage of the coordinate splitting method is achieved by reducing the number of arrays and reducing the computation time by up to 25% (Fig. 7).

CONCLUSION
A three-dimensional nonlinear mathematical model of the gas filtration process in a porous medium has been developed. Also, an efficient numerical algorithm and software were developed for solving the problem of gas filtration in porous media by the method of coordinate splitting. The developed computational algorithm and software provide the possibility of reducing the computation time by 25% in comparison with other computational methods by reducing the number of cycles when calculating the gas filtration process.
The developed mathematical support of the process under consideration can be used to carry out numerical calculations on a computer for the purpose of analysis, forecasting, making managerial decisions on the development and design of oil and gas fields under various conditions of impact on the reservoir, and making specific practical recommendations depending on hydrogeological and geophysical properties porous media.

CONFLICT OF INTERESTS
The author(s) declare that there is no conflict of interests.