Numerical investigation of the problem of nonlinear three-phase filtration

Об’єктом дослідження є чисельне моделювання процесу фільтрації нафти, газу та води на адаптивних сітках з урахуванням деяких властивостей рідин під час їх спільної течії. Для отримання адекватного опису процесів необхідно одночасно враховувати вплив більшості з зазначених факторів на фільтрацію. Це в математичному плані призводить до вирішення систем нелінійних диференціальних рівнянь в приватних похідних, складність яких не дозволяє досить глибоко дослідити їх аналітичними методами. Експериментальне ж вивчення цих процесів пов’язано з проведенням тривалих і дорогих лабораторно-промислових експериментів. Одним з найбільш проблемних місць в теорії багатофазної фільтрації є те, що кроки по просторовій змінній повинні подрібнюватися в областях різкої зміни не тільки градієнта водонасиченості, але і градієнта газонасиченості. Це пояснюється тим, що через дуже малу в’язкість вільний газ під дією градієнта тиску обганяє інші компоненти суміші – воду та нафту. В ході дослідження використовувалися алгоритм побудови адаптивних сіток, який може бути адаптований до властивостей рішення, а також методи обчислювальної математики, у тому числі разностноітераційний метод в рухомих сітках. Були проведені чисельні експерименти для оцінювання впливу запропонованого методу на переміщення і розміри «нафтового валу» і дано порівняльний аналіз отриманих результатів на основі чисельних розрахунків. Завдяки цьому показано, що при наближенні нафтового валу до експлуатаційної свердловини з пласта виходить тільки газ, і в міру зменшення в’язкості нафти зменшується час наближення нафтового валу до експлуатаційної свердловини. А також показано, що при збільшенні в’язкості нафти довжина lv і «приріст» hv нафтового валу зменшуються, причому зменшення lv в порівнянні з hv відбувається з більшою швидкістю. A при збільшенні швидкості фільтрації і різниці тисків геометричні розміри і «приріст» нафтового валу різко збільшуються. Ключові слова: капілярні сили, трифазна фільтрація, адаптивна сітка, закон Дарсі, метод прогонки, в’язкопластична рідина. Qasimov S., Mammadov R., Karimova S.


Introduction
As is known [1], important practical problems of oil and gas production, hydraulic engineering, chemical technology, etc. are associated with the problem of multiphase filtration. This problem is the most urgent for the theory and practice of oil and gas production.
The tasks of multiphase filtering have certain specific features that do not make it possible to always use the finite difference method that has become a classic finitedifference method for their numerical solution. Therefore, there is a need to develop difference schemes on adaptive grids, taking into account the features of the solution [2,3].
Adaptive grids reduce artificial viscosity and oscillation of a numerical solution. Even if the computational grid has minimal points, they again give satisfactory results in terms of quantity and quality for the entire region (even with the inclusion of zones with solution features).
Therefore, it is relevant to study the processes of joint filtration of oil, gas and water in a porous medium, taking into account: -real geological field data (elastic properties and formation heterogeneity); -properties of filtered liquids and gas (compressibility, dependence of viscosities on pressure); -relative phase permeability and capillary forces.

The object of research and its technological audit
The object of research is the numerical simulation of the filtration process of oil, gas and water on adaptive grids, taking into account some properties of liquids during their joint flow.
The tasks of multiphase filtration have specific features that do not make it possible to always use the classical finite difference method in their numerical solution.
One of the most problematic places in the theory of multiphase filtration is that the steps of the spatial variable must be crushed in areas of a sharp change not only in the gradient of water saturation, but also in the gradient of gas saturation. This is explained by the fact that due to the very low viscosity, free gas under the action of the pressure gradient overtakes the other components of the mixture -water and oil. Therefore, there is a need to develop difference schemes on adaptive grids, which makes it possible to take into account the features of the solution.

The aim and objectives of research
The aim of research is development of efficient and sufficiently universal numerical methods for solving the problem of one-dimensional three-phase filtration. TECHNOLOGY AUDIT AND PRODUCTION RESERVES -№ 1/1(45), 2019

ISSN 2226-3780
To achieve this aim it is necessary to perform the following objectives: 1. Build a mathematical model of the process of onedimensional three-phase filtration in a porous medium.
2. Build an adaptive grid, taking into account the grinding steps in areas of abrupt change not only the gradient of water saturation, but also the gradient of gas saturation.
3. Develop an algorithm for the difference-iterative method in moving grids.
4. On the basis of numerical experiments, discuss the issues of computational implementation of the proposed method and give practical recommendations for its application.

Research of existing solutions of the problem
When conducting hydrodynamic calculations associated with the conduct and operation of oil and gas fields in the water-pressure mode, the Rapoport-Lis two-phase model is mainly used. However, due to the fact that in many fields there is a three-phase flow, it is necessary to consider a model of three-phase flow and evaluate the effect of the gas phase on the displacement process. And this, in turn, leads to the need to develop a numerical method for solving three-phase filtration problems described by a system of nonlinear partial differential equations. To study this process, numerical methods are mainly used [4,5].
In the works [6,7], the Buckley -Leverett model is used, and in [4] the effects of capillary forces are also calculated. And for the numerical solution of the problems obtained, the idea of the method is used [8]. However, this method does not take into account the characteristics of multiphase filtration, and in some cases gives qualitatively and quantitatively incorrect results [4].
Among the main directions of solving this problem, identified in the resources of the world scientific periodicals, the works [9][10][11] can be highlighted. However, in [9], the extrusion is assumed to be piston, and approximate analytical solutions are obtained. A numerical method, ideologically different from the finite difference method and finite element method, is proposed in [10]. It should be noted that modifications of these works are widely used by many authors [11].
In this paper, an effective difference-iterative method in moving grids is proposed, which has the property of adaptability to the features of problems and is distinguished by high accuracy [2].

Methods of research
5.1. Mathematical formulation of the problem. Let's suppose a circular horizontal reservoir of radius R and thickness H is considered. It is assumed that the reservoir does not let in fluid from below and above, in the center of the reservoir there is a perfect hydrodynamic production well of radius r = r c , on the contour r = R there is a battery of wells pumping compressible fluid into the reservoir.
Let's suppose that in the reservoir in question there is a three-phase compressible fluid, the phases are immiscible and at the initial moment of time t = 0 are in a state of capillary equilibrium. Based on the Rapoport-Lias model for liquids with reduced conditions, the isothermal planeradial nonlinear filtering problem is described by a system of partial differential equations in a dimensionless form:  , , , (1) The following notation is used here: λ α -phase mobility coefficient: Function: , , takes into account the viscoplastic property of oil [1,12].
Here G 1 -the initial pressure gradient.
To the system of equations (1) it is necessary to add the equations of state of liquids and gas: and ratios: . s s s Let's note that, in general, in problems of multiphase filtration, the pressure p α(α = 1,2,3) is not equal to each other and differs by the value of the corresponding capillary pressure, i. e.: where p k1 -the capillary pressure between the oil and water, and p k2 -the capillary pressure between the gas and the oil.
Knowing p k1 and p k2 , it is possible to use them to express the capillary pressure between gas and water: It should be noted that functions , , , , , f s s s Let's suppose that the functions unknown in system (1): Let's define them from the initial and boundary conditions.
Let's suppose that at time t = 0, i. e. before the start of operation, the values of the functions being defined are known, i. e.: , For the problem under consideration, the boundary conditions can be in the following form.
For a production well (i. e., at r = r c ), it is assumed that only oil is extracted from the formation, and at this point the capillary pressure jump is not taken into account: ( )= ( ) On the external contour (i. e., at r = R), it is assumed that only compressible fluid (water) is pumped into the formation and, accordingly, the first and third phases do not flow: , .
Thus, according to conditions (5)- (7), the task can be set as follows. It is necessary to find such functions: , , , so that they would satisfy the system of equations (1) and the initial boundary conditions (5)- (7). Assuming the existence and uniqueness of the solution to problem (1)-(7), let's give an approximate method for solving it.

5.2.
Numerical method for solving the problem. First of all, let's note one specific feature of multiphase filtration tasks. It has been established that during the whole study, usually a sharp change in phase pressures occurs around the wells. Therefore, when modeling multiphase nonlinear filtering processes, the boundary conditions must be approximated with high accuracy, i. e., the grid step around wells must be reduced [2,3].
Thus, let's cover the area: , : , [ ] For better convergence and accuracy of the results when solving the problem in question using an iterativedifference method, it is advisable to carry out calculations with an increment in time. Thus, the calculations begin with a very small step in t and gradually increase to a given T [2]. In accordance with this, the grid nodes ωτ n are defined as follows:  Thus, increasing the step in t occurs every k c time layer. In the calculations, k c = 100 and τ in. = 4·10 -8 .
Using an implicit conservative scheme, let's approximate the differential problem (1)-(7) on a non-uniform grid ω τ n h with the following discrete problem:   ds dp s ds dp s ds dp ds dp k k k k By approximating the derivatives ∂ ∂ p t α on the right side of system (8) as follows: Taking into account the above, the difference problem (12) can be written in the following three-diagonal canonical form:  1  11  1  1 2  2  13  3  , , , , ε would be satisfied at the same time. Here, ε 1 , ε 2 , ε 3 -the predetermined accuracy for finding the i-th component of the approximate solution.

Research results
The iterative-difference method proposed above was applied to the numerical simulation of one three-phase filtration process. The following initial data were used for the calculations: . Prior to the beginning of the process, the phase saturation was S 1 = 0.45, S 2 = 0.25, S 3 = 0.30, respectively. The dependences of relative permeabilities and capillary forces on saturations were taken following [15].
The conducted numerical experiments show that at the beginning of the impact process, an enlarged zone of oil saturation begins to form in the reservoir -an oil shaft, which dimensions increase and which, displacing the gas, moves towards the production well. The displacement of oil occurs almost after the oil shaft. Under the «width» and «increase» of «oil shaft» profile refers to the following values: , . It should be noted the following feature of the process. When the oil shaft approaches the production well, only gas escapes from the reservoir. Fig. 1 shows that as the viscosity of the oil decreases, the time for the oil shaft to approach the production well decreases.  In Fig. 2, the initial distribution of oil and gas saturation is shown by dots. Curves marked with numbers 1, 2, 3, 4, 5, 6 are the results of the calculations, respectively, at μ 1 0 01 0 02 0 03 0 1 0 2 0 4 = . ; . ; . ; . ; . ; . poise. As can be seen from Fig. 2, with an increase in oil viscosity, the length l v and the «increase» h v of the oil shaft decrease, and the decrease in l v compared to h v occurs at a higher rate. For example, from Table 1 that when increasing the value of μ 1 by 2, 10 and 40 times, l v decreases by about 41 %, 73 %, 86 %, and h v , respectively, decreases by 5 %, 36 %, 65 %.
Thus, based on the above, in the presence of a free gas phase, the geometrical dimensions of the oil shaft ISSN 2226-3780 characterize the completeness of the oil displacement, so that an increase in the shaft dimensions is the result of the oil being displaced by water. Fig. 2. The effect of oil viscosity on the size and «growth» of the oil shaft Table 1 The effect of oil viscosity on the size of the oil shaft Let's consider the effect of the displacement rate on the process of three-phase displacement ( Table 2). The results of calculations show that with an increase in the filtration rate, the geometric dimensions and the «increment» of the oil shaft increase dramatically. As can be seen from the Table 2, with an increase in the pressure difference to 0.4 kgf/cm 2 , the shaft length increases by 35.4 %, and the increase by 9.6 %.