Computational modelling of flow in a microporous environment of a natural reservoir with consideration of the topological structure

. The article presents the initial stage of developing a hybrid mathematical model for describing the recovery of oil from a natural sand reservoir. The hybrid model means that part of the computational domain (microchannel structure in the space between particles) is calculated by one-dimensional models is based on the hydraulic circuit theory, and another part (large cracks and caverns) is calculated by computational fluid dynamics methods (CFD). The article also describes a unique algorithm for building a network model based on preliminary CFD calculation. The last section of the article presentation compares of results of calculation CFD and network models.


Introduction
At the moment, the main methods of hydrocarbon production are the use of well. However, most of the world's proven hydrocarbon reserves (78% of gas and 58% of oil) are located in so-called terrigenous (sand) reservoirs. This is a porous medium composed of sand particles, filled with oil or gas. One of the most common methods of recovering oil from such fields is flooding, in which oil is displaced from a porous medium by water or solutions with surfactants. The average production level of oil recovery is obtained by this method for various countries and regions from 25% to 40%, which means that most of the oil remains in the reservoir and becomes even more difficult to access. Thus, increasing the oil recovery rate from the reservoir is an important task. To solve this task, it is necessary to accurately study the process of oil and gas displacement from such a reservoir.
There are two main complementary approaches to the study of oil displacement: experimental [1,2] and numerical [3][4][5][6][7][8][9][10][11]. The use of numerical methods, allows us to obtain a detailed picture of the flow of oil suspension (oil-gas-water) in a porous medium. In most cases, one-dimensional approximate models of a porous medium are used to model the process of oil displacement [5,6]. The fluid motion in this model is described by Darcy's law. Unfortunately, not all rock types can be represented as an isotropic structure. So, there are rocks that are an impenetrable medium with cracks and caverns, or with a mixed type in which the permeable medium is adjacent to cracks and caverns. Another approach is direct numerical modeling when a two-or three-dimensional model of pores space is constructed and the fluid motion is described by a system of Navier-Stokes equations. To describe a multiphase spatial flow, both the mixture model [1,3,7] and models with the calculation of the phase interface [8,9] can be used. However, the use of spatial (3D) models for methodological research for these purposes leads to unnecessarily large computational costs. Therefore, it is preferable to use methods of simplifying numerical models that will reduce computational costs without losing the required calculation accuracy.
If you look closely at the structure of a porous medium, you can see that the void space is essentially an extensive network of microchannels with characteristic dimensions from 1 to 1000 microns. Thus, to model such a structure, we can use the methods Of the theory of Hydraulic Circuits [12]. Examples of constructing such models are given in [10,11]. However, large cracks and caverns are difficult to describe as one-dimensional models and they will have to be left three-dimensional, so you need to think about the mechanism for connecting the spatial and network parts of the problem.
So, the goal is to develop a mathematical model that allows to modeling multiphase fluid system of a microporous environment, different parts of which describe different scale (1D, 3D) submodels depending on the topological properties of the plot (see Fig. 1). It will be necessary to solve three main subtasks: 1. Automatic define of topological subdomains in the computational domain, and construct 3D mesh and 1D -net.
2. Selection of mathematical models that accurately describe the process of non-stationary flow of multiphase suspension for each type of subdomain.
3. Combining different-scale subdomains of the model into a single numerical model.

A mathematical model
The main idea of the proposed method is as follows: the initial calculation domain is divided into spatial and network subdomains the interaction between which is carried out using the hybrid method based on the solution of the coupled equation for the pressure correction. This approach was developed earlier for single-phase and stationary models [13].
To model the oil recovery process, we will use a numerical method based on the liquidin-cells method, which has proven itself well for calculating macroscopic flows with a free surface [14,15]. The idea of the VOF method is that into the conservation equations, the liquid is considered as a single medium with density and viscosity depending on the value of the phase concentration in each cell. The spatial distribution of phases is set using a special function F(x, y, z, t), where F is the value that determines the volume fraction of the liquid phase. F(x, y, z, t) = 0 if the cell is filled with water and F (x, y ,z , t) = 1 if the cell is filled with oil.
To determine the position of the phase interface in space, the transport equations for the volume fraction of the water phase are solved in each cell: where v is the velocity vector of a two-phase medium obtained as a result of solving the system of hydrodynamics -Navier-Stokes equations. Which consists of the law of conservation of mass: �� �� + ∇( ⋅ ) = 0, (2) and law of momentum conservation: �� �� + ∇( × ) = −∇ + ∇( ) + , (3) where τ is the viscous stress tensor, F is the volume force vector, p is the pressure, and is the mean density of the two-phase medium.
The components of the viscous stress tensor τij are defined as: where µ is the mean dynamic viscosity of the two-phase medium, and uij are the components of the velocity vector v.
The molecular viscosity and density of a two-component water-oil medium are set depending on the volume fraction ( ) of each liquid in the cell according to the mixture rule, where � , � are the density and viscosity of water, and � , � are the density and viscosity of the oil, respectively: = � + (1 − ) � . (6) When considering two-component fluid flows in a porous medium, the surface tension force at the interface plays a significant role. The study of flows, taking into account the surface tension forces, is a complex and time-consuming task. However, for the VOF method, this problem has already been solved with the power of the Continuum surface force (CSF) algorithm [14], the essence of which is to introduce an additional volume force FS into the equations of motion. The method used is described in detail in [16]. To describe the network part of the problem, we will use the methods of the theory of Hydraulic Circuits [11]. A hydraulic network consists of a set of nodes and branches that represent a directed graph, the matrix of connections of which is represented as: Here l is the number of the branch, l ∈ Oi is the set of branches issuing from the ith node, and l ∈ Ii is the set of branches entering the ith node. By (7), the problem of the flows distribution in the network can be reduced to the combination of the mass conservation law at the node (8) and the law of resistance in the pipe (9): (9) where ql is the mass flow in the branch (kg/s), Qi is the mass source at the node (kg/s), pi is the pressure at the node (Pa), sl is the resistance coefficient, which for these flow regimes -Re~0.01, depends only on friction, determined by the following formula: where λ is the linear coefficient of friction, A is the cross-sectional area of the branch (m 2 ), and S is the area of the surface of the branch (m 2 ). To model the flow of a multiphase medium, a homogeneous flow model will be implemented at the initial stage, where the average density and viscosity will be calculated using formulas (5) and (6). Combine different topological elements into a single hydrodynamic model based on the hybrid method [13] with a common pressure field.

Building a network model based on results of spatial calculation
Building a network topology in a porous system is a special task (the stages of this process are shown in Fig. 2). An attempt to determine the "manual" network structure is not possible, so we will use the unique algorithm developed by the authors to build a network model based on the results of the preliminary spatial calculation. To do this, we will perform a spatial calculation of the studies object in a simplified formulation (a coarser grid, stationary, and single-phase flow). From this calculation, we need to know the velocity field and the cell wall distance Fig. 2a.
At the first stage, we will build streamlines using the velocity field (Fig. 2b). Then, moving along the streamlines, we will build branches with a given step (Fig. 2c). The next step is to merge the nearby nodes (Fig. 2d). Thus, we get a sparse network that has a known length of branches, and knowing the distance to the wall for each node, we can set the width of the branch (Fig. 2e).
To validation, the constructed model, a series of calculations were performed in comparison with the two-dimensional formulation. Stationary water flows through a system of microchannels were studied. A constant pressure condition was set at the upper and lower boundaries (10, 20, 30, and 40 Pa at the upper boundary and 0 Pa at the lower boundary) The dependence of the total flow rate on the pressure drop was investigated. The calculation results are shown in Fig. 3 and Fig. 4. Figure 3 shows and the pressure field. These figures show a good qualitative match between the spatial and network models. A quantitative comparison is shown in Fig. 4. the results of a series of calculations show a satisfactory match of the calculation results, which indicates the correctness of the methodology for building a network model. b Fig. 3. Pressure field of two-dimensional and network variants of water flow through a system of microchannels.

Conclusion
A mathematical model of the process of oil recovery from natural reservoirs is developed using a network and hybrid models. For its implementation, the following approaches are supposed to be used: for modelling the suspension flow in the three -dimensional part of the problem, the VOF method is supposed to be used, in particular, at the initial stage-a homogeneous model. It is planned to combine the network and spatial parts of the problem using the method of the common equation for the pressure correction. Now, the procedure for constructing the topology and geometric parameters of the network model based on the results of spatial calculation is implemented. The network model constructed using this method gave satisfactory results in comparison with the twodimensional model, but for its further use, it is necessary to research with three-dimensional models.