Development and Testing of a Mathematical Model for Dynamic Network Simulation of the Oil Displacement Process

: Multiphase ﬂows in porous media are widespread in nature and various technologies. One of the most common examples of this kind of task is the task of recovering oil from the rock. This article describes a mathematical model of the ﬂow of a two-phase (immiscible) liquid based on a new approach of network hydrodynamics for a highly branched microchannel medium (simulating a porous space in the rock). The coupling of the ﬂow and pressure ﬁelds in the network is performed using a well-proven SIMPLE algorithm in CFD problems; this approach allows us to use effective approaches to modeling 3D tasks. Phase transfer over the network is carried out by an explicit method with an adaptive time step. The article presents the results of veriﬁcation of the model, with analytical calculations and in comparison with the results of experimental studies. As an experiment, the displacement of oil from a microchip (Dolomite: 3200284) simulating a porous medium was simulated. The good qualitative and quantitative compliance with the results calculated and the results of the experiment show the correct functioning of the model. S.A.F.; software, S.A.F.; validation, S.A.F. and A.V.M.; formal analysis, S.A.F. and M.I.P.; investigation, M.I.P. and A.I.P.; writing—original draft preparation, S.A.F.; writing—review and editing, S.A.F. and M.I.P.; visual-ization, S.A.F.; supervision, A.V.M.; administration, funding acquisition,


Introduction
Multiphase fluid flows in porous media are one of the significant processes taking place in nature and technology. So, many processes of fluid filtration through the ground are multiphase in nature, covering the range from a mere transfer of moisture in the soil [1,2] to processes at waste storage sites, for example, underground carbon dioxide storage [3].
A special place in the problems of a multiphase fluid flow in porous media is held by the processes of oil displacement from natural reservoirs in which fluid flow is generally three-phase consisting of oil, water, and natural gas. In the oil extraction process, one of the most important roles is played by the formation structure which is characterized by the size and location of pores and cracks, as well as the rock material, which determines the permeability and strength of the formation. The percentage remaining after water flooding of oil is also an important indicator. When studying the process of oil recovery from a rock sample, two main lines of research can be distinguished: the experimental study of core samples, and the study of artificially created microchannel systems from regular and two-dimensional [4][5][6] to highly structured and three-dimensional [7], or the study of oil recovery from the core using X-ray microtomography in real time [8]. The first research line allows obtaining integral information about the physical properties of the core. The second approach is pore-scale simulation, which requires a detailed understanding of the physical processes occurring at the pore scale, as well as a complete description of the pore space morphology. This research line allows studying in more detail the processes occurring in microchannels, and establishing patterns expressed in mathematical form. Recently, an approach combining two noted research lines has been actively developing, namely, the geometry of the pore space is restored by core tomography [9], which is subsequently account for effects that were previously discarded. For example, taking into account the roughness of the walls of cracks [22] or the transfer of colloidal particles through the network [23] with the prediction of plugging of the necks [24]. A separate direction in the development of modeling the flow of a multiphase liquid in a porous medium is the conjugation of one-dimensional (1D) and spatial (2D, 3D) models into one hybrid model. The network model is associated with both the spatial model based on CFD models [25] and the lattice Boltzmann models [26].
This paper presents the results of the development of an original dynamic network model of oil displacement from a porous rock. The hydrodynamic part of the model is based on the theory of hydraulic circuits [27]. The peculiarity of the model is that, unlike similar models in general, it does not depend on the structure of the network. The connection of the flow and pressure fields in the network is carried out using a well-proven SIMPLE algorithm in CFD problems; this circumstance allows us to use effective approaches to modeling 3D problems. This approach allows us to implement a hybrid 1D-3D approach that allows us to obtain a continuous pressure field for the entire hydrodynamic multiscale, both in one-dimensional and three-dimensional formulations, due to the end-to-end solution of the pressure correction problem. The article is devoted to the description and testing of the proposed model.

Mathematical Model of Network Simulation of Oil Displacement Processes
Before describing the mathematical model, we introduce several notions and limitations: 1.
Two-phase motion is considered: a fluid with a wetting angle less than 90 degrees will be called wetting, and the other one is non-wetting.

2.
One mode of two-phase fluid flow is implemented, which is the piston-like displacement mode where the phase interface is perpendicular to the flow. 3.
Reverse displacement of fluid is prohibited, i.e., if the element is completely filled with the displacing fluid, then the displaced fluid cannot return back to the element.

4.
Fluids do not mix with each other due to a rigid phase interface.

5.
However, in network elements, where both components are present, the physical properties of the fluid (density and viscosity) are determined by the mixture rule. 6.
Fluids are incompressible; therefore, information about pressure changes is distributed instantly without shock waves.
The mathematical model, describing the oil displacement process from the reservoir, can be divided into two parts:

A Hydraulic Network Model for Determining the Distribution of Pressure Drops and Flow Rates throughout the Network
The concept of a hybrid 1D-3D network algorithm for simulating the flow process in networks of complex topology is described in our work [28]. To describe the network model, we introduce the appropriate notation: the set of nodes is denoted by N, and the set of branches (pipes) is U (Figure 1). Let O i be a subset of branches starting at the i-th node, and I i be a subset of branches, ending at the i-th node. Accordingly, each branch is associated with a pair of nodes with the accepted designation: N in is the initial node and N out is the end node. The direction of the branch is set from the initial node to the final one, respectively, the flow rate q l and the velocity u l on the l-th branch (l ∈ U) can take either a positive value (when the fluid flow direction coincides with the branch) or a negative one. It is also worth noting that in general, most nodes in the network are not connected: on average, each node has two to four branches. Therefore, from the entire set of branches U ij , where i, j ∈ N, it is necessary to allocate only a relatively small subset of existing branches. Uij, where i,j ∈ N, it is necessary to allocate only a relatively small subset of existing branches. This approach allows defining a matrix of connections for the entire graph in the form Equation (1). Using this expression, the problem of flow distribution in the network can be reduced to a combination of the mass conservation law in the node Equation (2) and the fluid-to-pipe resistance law Equation (3).
where ql is the carrier flow in the branch, Qi is the source of the mass concentrated in the node, pDi is the pressure in the i-th node, is the capillary pressure, whose sign depends on the distribution in the pipe, sl is the resistance coefficient, determined by the formula, known also as Darcy-Weisbach law: where λfr is the linear friction coefficient, d is the hydraulic diameter of the branch, l is the length of the branch, ρ is the density of the fluid, f is the cross-sectional area of the pipe, ξ is the local resistance coefficient. Having written these equations for each node and branch (taking into account the boundary conditions), we obtain a system of nonlinear algebraic equations describing the entire model. It is worth noting that the problem of oil recovery is characterized by very small flow velocities less than 0.1 m/s, which leads to small Reynolds numbers Re ≤ 1 corresponding to a strictly laminar flow regime.
A well-proven SIMPLE-like algorithm is used to link the pressure field and flow rates [29].
This algorithm is based on iterative obtaining of the pressure field: This approach allows defining a matrix of connections for the entire graph in the form Equation (1). Using this expression, the problem of flow distribution in the network can be reduced to a combination of the mass conservation law in the node Equation (2) and the fluid-to-pipe resistance law Equation (3).
where q l is the carrier flow in the branch, Q i is the source of the mass concentrated in the node, p Di is the pressure in the i-th node, P c is the capillary pressure, whose sign depends on the distribution in the pipe, s l is the resistance coefficient, determined by the formula, known also as Darcy-Weisbach law: where λ fr is the linear friction coefficient, d is the hydraulic diameter of the branch, l is the length of the branch, ρ is the density of the fluid, f is the cross-sectional area of the pipe, ξ is the local resistance coefficient. Having written these equations for each node and branch (taking into account the boundary conditions), we obtain a system of nonlinear algebraic equations describing the entire model. It is worth noting that the problem of oil recovery is characterized by very small flow velocities less than 0.1 m/s, which leads to small Reynolds numbers Re ≤ 1 corresponding to a strictly laminar flow regime.
A well-proven SIMPLE-like algorithm is used to link the pressure field and flow rates [29]. This algorithm is based on iterative obtaining of the pressure field: where p k is the pressure value at the current iteration, p k−1 is the pressure value at the previous iteration, p is the pressure correction. To obtain the pressure correction equation here τ net are the coefficients that are obtained in the approximation of Equation (3).
The system of Equation (6) can be represented in matrix form: where A is the matrix of coefficients, p is the vector of the required pressure corrections and ∆q is the vector of flow residuals. This matrix is a system of linear algebraic equations, which can be solved using highly efficient iterative methods. In our case, we use the Conjugate Gradient method, which has proven itself for calculating symmetric matrices and is well parallelized, which in turn makes it possible to simulate the processes of oil displacement in materials consisting of millions of pores and throats. The convergence of the iterative method is set to be 1 × 10 −3 , which is at least five orders of magnitude less than the average pressure drop in the studied systems. The density and molecular viscosity of the two-component medium under consideration are found through the volume fraction of the fluid in the cell according to the mixture rule, where ρ 1 and µ 1 are the density and viscosity of the first fluid, ρ 2 , µ 2 are, respectively, the density and viscosity of the second fluid, and ϕ is the volumetric fraction of the first phase: Phase transfer along the length of the throat is modeled by solving a one-dimensional convective equation with an up-wind approximation.

∂ϕ ∂t
where t is the time and z is the length of the throat, u is the average volume velocity on the branch. After integration with the first-order variable of the time difference scheme, a system of balance equations is obtained, described in more detail in the next paragraph. Before the calculation, it is necessary to set the initial distribution of pressure and phases over the network. The distribution data can be varied, but as a rule, at the initial moment, the network is filled with oil, except for the boundary nodes in which the source of mass Q (see Equation (2)) is greater than zero. Given the incompressible nature of the medium, the initial pressure distribution can be set to zero in the entire network, so it will recover in one iteration of the hydraulic calculation.

Transfer of the Second Phase over the Network
Let us consider the second part of the mathematical model in more detail. Assuming that in the hydraulic model of the network a node is a material point with a given pressure, and the main working volume of the network is concentrated in branches, for accounting of the accumulation/entrainment of various phases of fluids, it is necessary to assign the volume to each node. To do this, the network is divided into elements such as pores which are the part of the network where flows merge/separate, and throats, representing the part of the network connecting two neighboring pores. Both pores and throats have their volume, and the length of the throat is shorter than the length of the corresponding branch (see Figure 2). Since fluids do not mix, this substitution of one fluid with another is generally reduced to the balance equation: where is the volumetric fraction at the next time step, is the volumetric fraction at the previous time step, ∑ and ∑ are the total flows bringing/carrying away the phase into/out of the element, respectively, is the volume of the element, ∆ is the time step, i is the number of the element (both for pores and the throats).
The time step ∆T is not a constant value; it is defined as the minimum time ∆ min required to fill one of the phases of a single network element. It is logical to allocate an array of elements in which, as is expected, ∑ − ∑ ≠ 0, and work only with it, recalculating it at each time layer, thus tracking the front of the phase movement in the reservoir.
In terms of phase replacement, the throats differ from the pores only by the fact that throats have just one inlet and one outlet (see Figure 3). It is worth considering specifically the situation where two flows of different fluids inflow into the pore. Figure 4 shows an example of such a situation. Since we cannot resolve the shape of the phase interface, and therefore cannot determine which phase outflows from the pore through which throat, to maintain the mass balance, it is necessary to make the following assumptions: 1 If the pore is not filled with the displacing fluid, then only the displaced fluid comes out of the node (Figure 4a). Since fluids do not mix, this substitution of one fluid with another is generally reduced to the balance equation: where ϕ new i is the volumetric fraction at the next time step, ϕ old i is the volumetric fraction at the previous time step, ∑ q in and ∑ q out are the total flows bringing/carrying away the phase into/out of the element, respectively, V is the volume of the element, ∆T is the time step, i is the number of the element (both for pores and the throats).
The time step ∆T is not a constant value; it is defined as the minimum time ∆T min required to fill one of the phases of a single network element. It is logical to allocate an array of elements in which, as is expected, ∑ q in − ∑ q out = 0, and work only with it, recalculating it at each time layer, thus tracking the front of the phase movement in the reservoir.
In terms of phase replacement, the throats differ from the pores only by the fact that throats have just one inlet and one outlet (see Figure 3). Since fluids do not mix, this substitution of one fluid with another is generally reduced to the balance equation: where is the volumetric fraction at the next time step, is the volumetric fraction at the previous time step, ∑ and ∑ are the total flows bringing/carrying away the phase into/out of the element, respectively, is the volume of the element, ∆ is the time step, i is the number of the element (both for pores and the throats).
The time step ∆T is not a constant value; it is defined as the minimum time ∆ min required to fill one of the phases of a single network element. It is logical to allocate an array of elements in which, as is expected, ∑ − ∑ ≠ 0, and work only with it, recalculating it at each time layer, thus tracking the front of the phase movement in the reservoir.
In terms of phase replacement, the throats differ from the pores only by the fact that throats have just one inlet and one outlet (see Figure 3). It is worth considering specifically the situation where two flows of different fluids inflow into the pore. Figure 4 shows an example of such a situation. Since we cannot resolve the shape of the phase interface, and therefore cannot determine which phase outflows from the pore through which throat, to maintain the mass balance, it is necessary to make the following assumptions: 1 If the pore is not filled with the displacing fluid, then only the displaced fluid comes out of the node (Figure 4a). It is worth considering specifically the situation where two flows of different fluids inflow into the pore. Figure 4 shows an example of such a situation. Since we cannot resolve the shape of the phase interface, and therefore cannot determine which phase outflows from the pore through which throat, to maintain the mass balance, it is necessary to make the following assumptions: 1.
If the pore is not filled with the displacing fluid, then only the displaced fluid comes out of the node (Figure 4a).

2.
If the pore is filled with the displaced fluid, then an artificial ban is imposed on the non-displaced fluid flow, which is the prohibition of reverse displacement (Figure 4b).
2 If the pore is filled with the displaced fluid, then an artificial ban is imposed on non-displaced fluid flow, which is the prohibition of reverse displacement (Figu 4b).

Clustering
The condition prohibiting reverse displacement ( Figure 4b) can lead to the format of closed areas of the displaced fluid, namely, clusters, surrounded by the displaced flu ( Figure 5). Network nodes and branches that fall into such a cluster must be exclud from the hydraulic computation, since this may lead to the breakup of the numeri scheme. Therefore, it is necessary to check for the presence of such clusters at each ti step.

Capillary Pressure
When two immiscible fluids flow in the capillary, a pressure drop is formed at th interface due to surface tension, whose value for a cylindrical capillary is determined fr the Young-Laplace equation: where σ is the surface tension coefficient, θ is the contact angle between the solid wall a the fluid (Figure 6), r is the radius of the element.

Clustering
The condition prohibiting reverse displacement ( Figure 4b) can lead to the formation of closed areas of the displaced fluid, namely, clusters, surrounded by the displaced fluid ( Figure 5). Network nodes and branches that fall into such a cluster must be excluded from the hydraulic computation, since this may lead to the breakup of the numerical scheme. Therefore, it is necessary to check for the presence of such clusters at each time step.

2
If the pore is filled with the displaced fluid, then an artificial ban is imposed on the non-displaced fluid flow, which is the prohibition of reverse displacement ( Figure  4b).

Clustering
The condition prohibiting reverse displacement ( Figure 4b) can lead to the formation of closed areas of the displaced fluid, namely, clusters, surrounded by the displaced fluid ( Figure 5). Network nodes and branches that fall into such a cluster must be excluded from the hydraulic computation, since this may lead to the breakup of the numerical scheme. Therefore, it is necessary to check for the presence of such clusters at each time step.

Capillary Pressure
When two immiscible fluids flow in the capillary, a pressure drop is formed at their interface due to surface tension, whose value for a cylindrical capillary is determined from the Young-Laplace equation: where σ is the surface tension coefficient, θ is the contact angle between the solid wall and the fluid (Figure 6), r is the radius of the element.

Capillary Pressure
When two immiscible fluids flow in the capillary, a pressure drop is formed at their interface due to surface tension, whose value for a cylindrical capillary is determined from the Young-Laplace equation: where σ is the surface tension coefficient, θ is the contact angle between the solid wall and the fluid (Figure 6), r is the radius of the element. The pressure, calculated by Equation (12), is set as the pressure head in Equation (3). The direction of the pressure head is set from the non-wetting fluid toward the wetting one. If the phase interface passes through the pore, then the corresponding pressure head is set on the branch upstream from the node corresponding to the pore. When the phase interface passes through the pore, the capillary pressure does not disappear; however, it is impossible to determine its value. Therefore, we assume that the capillary force remains on the branches filled with displacing fluid. The pressure, calculated by Equation (8), is set as the pressure head in Equation (3). The direction of the pressure head is set from the non-wetting fluid toward the wetting one. If the phase interface passes through the pore, then the corresponding pressure head is set on the branch upstream from the node corresponding to the pore. When the phase interface passes through the pore, the capillary pressure does not disappear; however, it is impossible to determine its value. Therefore, we assume that the capillary force remains on the branches filled with displacing fluid.

Computational Algorithm
Considering the above-described assumptions, the following computation algorithm can be constructed, consisting of the following steps: 1. Setting the initial distribution of pressure, flow rates, and phase concentrations over the network. 2. Determining the list of network elements in which the phase replacement process at the first time step is possible. 3. Calculating the density and pressure for all network elements by Equations (5) and (6). 4. Calculate (by Equation (8)) and set the pressure for the branches corresponding to capillary pressure. 5. Solving the hydraulic problem, which results in obtaining a new distribution of pressure and flow rates. 6. Calculating the pressure field for the nodes. 7. Determining the flow rate in the branches using the above calculated pressure field. 8. Prohibiting the reverse displacement of the fluid in the network branches. 9. Repeating the operation starting from point 5a until the number of blocked branches becomes zero. 10. Determining ∆ min . 11. Recalculating the phase concentrations using Equation (7). 12. Conducting checking for the presence of closed clusters. 13. Updating the list of network elements in which the phase replacement process is possible at a time step. 14. Repeating the entire algorithm starting from the 3rd step.
The algorithm terminates when the total computation time exceeds the specified one, or there is no phase change in any element.

Validation and Testing the Developed Numerical Technique
Validation and verification of the hydraulic part of the algorithm were carried out earlier and are described in our previous papers [28]. Therefore, here we will focus in detail on testing and validation of phase transfer problems.

Oil Displacement from the Straight Microchannel
The first test concerns the oil displacement process from a straight microchannel of a rectangular cross-section. The microchannel cross-section was 500 × 500 µ m, and the length was 0.1 m. The microchannel was divided into 10 branches (Figure 7).

Computational Algorithm
Considering the above-described assumptions, the following computation algorithm can be constructed, consisting of the following steps: 1.
Setting the initial distribution of pressure, flow rates, and phase concentrations over the network.

2.
Determining the list of network elements in which the phase replacement process at the first time step is possible.

3.
Calculating the density and pressure for all network elements by Equations (5) and (6).

4.
Calculate (by Equation (12)) and set the pressure for the branches corresponding to capillary pressure.

5.
Solving the hydraulic problem, which results in obtaining a new distribution of pressure and flow rates. 6.
Calculating the pressure field for the nodes. 7.
Determining the flow rate in the branches using the above calculated pressure field. 8.
Prohibiting the reverse displacement of the fluid in the network branches. 9.
Repeating the operation starting from point 5a until the number of blocked branches becomes zero. 10. Determining ∆T min . 11. Recalculating the phase concentrations using Equation (11). 12. Conducting checking for the presence of closed clusters. 13. Updating the list of network elements in which the phase replacement process is possible at a time step. 14. Repeating the entire algorithm starting from the 3rd step.
The algorithm terminates when the total computation time exceeds the specified one, or there is no phase change in any element.

Validation and Testing the Developed Numerical Technique
Validation and verification of the hydraulic part of the algorithm were carried out earlier and are described in our previous papers [28]. Therefore, here we will focus in detail on testing and validation of phase transfer problems.

Oil Displacement from the Straight Microchannel
The first test concerns the oil displacement process from a straight microchannel of a rectangular cross-section. The microchannel cross-section was 500 × 500 µm, and the length was 0.1 m. The microchannel was divided into 10 branches (Figure 7).
These analytical values of capillary pressure were used for comparison with the calculated values (see Table 1).  the test, three computations were carried out with a wetting angle of 40, 90, and 140 • , respectively. According to Equation (12), at these conditions, the capillary pressure was: These analytical values of capillary pressure were used for comparison with the calculated values (see Table 1).  Table 1 shows the total drop and capillary pressure along the length of the microchannel at a time of 1.5 s after the start of displacement, in comparison with the pressure drop obtained by analytical calculation. The Table 1 shows that for all three options, the deviation from the analytical solution is less than 0.1%. The difference between the different variants is the value corresponding to the capillary pressure at the interface.

Oil Displacement from a Microfluidic Chip with Irregular Porosity
To conduct a more complex test, the results of the chip-based numerical simulation were compared with the results of a physical experiment [30]. The essence of this experiment is reduced to modeling the process of oil displacement from a sample with a complex porous structure. A chip with a large number of branched microchannels is used as a sample simulating a porous rock. At the initial moment time, the chip is filled with oil,

Oil Displacement from a Microfluidic Chip with Irregular Porosity
To conduct a more complex test, the results of the chip-based numerical simulation were compared with the results of a physical experiment [30]. The essence of this experiment is reduced to modeling the process of oil displacement from a sample with a complex porous structure. A chip with a large number of branched microchannels is used as a sample simulating a porous rock. At the initial moment time, the chip is filled with oil, then water with a given flow rate is supplied to the chip through the inlet openings and the dynamics of oil substitution with water is investigated. At some point in time, depending on the flow mode, the oil stops being squeezed out of the chip, and the water just continues to flow along the path that it paved. The residual amount of oil or, conversely, the maximum filling of the chip with water is the purpose of the study.
The microfluidic chip (Dolomite: 3200284- Figure 9) was used in the work, with a stable result of the experiment at water flow 10 µL/min.

Oil Displacement from a Microfluidic Chip with Irregular Porosity
To conduct a more complex test, the results of the chip-based numerical simulation were compared with the results of a physical experiment [30]. The essence of this experiment is reduced to modeling the process of oil displacement from a sample with a complex porous structure. A chip with a large number of branched microchannels is used as a sample simulating a porous rock. At the initial moment time, the chip is filled with oil, then water with a given flow rate is supplied to the chip through the inlet openings and the dynamics of oil substitution with water is investigated. At some point in time, depending on the flow mode, the oil stops being squeezed out of the chip, and the water just continues to flow along the path that it paved. The residual amount of oil or, conversely, the maximum filling of the chip with water is the purpose of the study.
The microfluidic chip (Dolomite: 3200284- Figure 9) was used in the work, with a stable result of the experiment at water flow 10 µ ⁄ . The microfluidic chip was made by etching sodium-lime glass (B270 material) with hydrogen fluoride, followed by thermal bonding. The size of the microfluidic chip was 92.5 × 15.0 mm and the thickness was 4 mm. The porous area of the chip had a size of 10 × 60 mm. The chip had a single inlet and a single outlet. The microfluidic chip was connected by a 4-sided upper interface (4 mm) (Dolomite: 3000109) and a linear 4-pin connector (Dolomite: 3000024). The porous area was formed by repeating (150 times) a square of 2 × 2 mm. The arrangement of channels in a square formed a grid of 8 × 8 channels, which had an almost circular cross-section (the channel depth and width were 100 and 110 microns, The porous area was formed by repeating (150 times) a square of 2 × 2 mm. The arrangement of channels in a square formed a grid of 8 × 8 channels, which had an almost circular cross-section (the channel depth and width were 100 and 110 microns, respectively). The channels in the grid had constrictions or pores that were randomly distributed to imitate the natural structure of the rock. The grid contained 38 pores with Ø63 µm, 40 pores with Ø85 µm, and 50 straight channels. Extra parameters are presented in Table 2. When carrying out numerical simulation, it is necessary to determine the hydrodynamic and physical parameters of the narrowing areas, namely: how much the volume of the pinched channel differs from that of the straight channel, as well as the hydraulic resistance of the pinched channel. The volume of the channel with a pore diameter of Ø85 µm amounts to 90% of the straight channel, while that of the channel with a pore diameter of Ø63 µm amounts to 86%.
To determine the hydraulic resistance of the pinched channels, the following method was used. For laminar flow in hydraulic networks, the pressure drop (4) can be represented as follows: where ζ is the local hydraulic resistance, l is the channel length, d is the channel diameter, ρ is the density, u is the average flow velocity in the channel, and µ is the viscosity. Thus, since the dependence of the pressure drop on the average velocity in the channel is approximated by an expression of the form: local resistance coefficient ξ can be expressed from the coefficient a. Three-dimensional CFD simulation was used to construct the dependence of the pressure drop on the velocity. The volumetric flow rate varied within the range from 0.176 to 1.76 m/s (which corresponds to a flow rate ranging from 0.1 to 1 mL/min). Figure 10 shows the velocity field in the central cross-section for different pore diameters.  The dependence of the pressure drop on the volumetric flow rate and the interpolation curves are shown in Figure 11. The following local hydraulic resistances have resulted from the computation: ξ ∅85 = 1.29 and ξ ∅63 = 16.4. To conduct numerical simulation, a network model was constructed (Figure 12) consisting of 9628 nodes and 18,954 branches. The distribution structure of straight and pinched channels was similar to that in the original chip ( Figure 13). To conduct numerical simulation, a network model was constructed (Figure 12) consisting of 9628 nodes and 18,954 branches. The distribution structure of straight and pinched channels was similar to that in the original chip ( Figure 13).
(a) (b) Figure 11. Pressure drop depending on the volumetric flow rate: channel with a pore diameter of Ø 85 (a) and Ø 63 (b) microns.
To conduct numerical simulation, a network model was constructed (Figure 12) consisting of 9628 nodes and 18,954 branches. The distribution structure of straight and pinched channels was similar to that in the original chip ( Figure 13).   In the process of oil displacement, the water flow front moves in the form of so-called viscous fingers [31]. Over time, the viscous fingers grow in length. Such fingers quickly break through to the exit. The water flows through separate channels. Further water flow occurs through separate channels. At the same time, approximately 40% of the initial volume of oil remains in a porous medium. The water flows through separate channels. After the main viscous fingers have reached the exit, the flow process becomes stationary. Further flooding of the microfluidic chip does not lead to an increase in the oil displacement coefficient. Figure 14 shows a qualitative comparison of the oil displacement pattern during chip flooding obtained by the computation and experiment. The analysis of the conducted test shows that in general, the computation reproduces well the main qualitative points of the flooding process, revealed in the experiment. The heterogeneous structure of the pore space leads to the winding of water flows, which is noticeable both in the computation and in the experiment. In the conditions under consideration, the water flooding front is not flat. Water moves in a porous medium, saturated with oil, in a form of separate jets. This behavior was observed both in the computation and the experiment. With this displacement, the residual oil remains in the form of separate islets, located periodically throughout the space of the microfluidic chip. So, a completely acceptable agreement of the flow pattern at compatible time points is seen in the experiment and computation using a network model.
A quantitative comparison of the computation with the experiment was carried out in terms of the time of saturation with water: where is the volume of water in the microchip, and is the volume of the microchip. The volume of the microchip is calculated as the total volume of the throats and pores, and the volume of water is calculated as the sum of the elements filled with water. If the network element is not filled, then only a part of the volume corresponding to the water phase is added to the volume of water. In the process of oil displacement, the water flow front moves in the form of so-called viscous fingers [31]. Over time, the viscous fingers grow in length. Such fingers quickly break through to the exit. The water flows through separate channels. Further water flow occurs through separate channels. At the same time, approximately 40% of the initial volume of oil remains in a porous medium. The water flows through separate channels. After the main viscous fingers have reached the exit, the flow process becomes stationary. Further flooding of the microfluidic chip does not lead to an increase in the oil displacement coefficient. Figure 14 shows a qualitative comparison of the oil displacement pattern during chip flooding obtained by the computation and experiment. The analysis of the conducted test shows that in general, the computation reproduces well the main qualitative points of the flooding process, revealed in the experiment. The heterogeneous structure of the pore space leads to the winding of water flows, which is noticeable both in the computation and in the experiment. In the conditions under consideration, the water flooding front is not flat. Water moves in a porous medium, saturated with oil, in a form of separate jets. This behavior was observed both in the computation and the experiment. With this displacement, the residual oil remains in the form of separate islets, located periodically throughout the space of the microfluidic chip. So, a completely acceptable agreement of the flow pattern at compatible time points is seen in the experiment and computation using a network model.

Conclusions
A mathematical model of the two-phase (immiscible) fluid flow in a porous medium has been developed based on the network hydrodynamics methods. The hydrodynamic part of the model is based on the well-established theory of hydrodynamic chains (graphs), which provides a significantly higher computational speed compared to the methods of computational fluid dynamics (CFD). The phase interface is taken into account by the contribution of capillary forces to the momentum conservation equation in the network branch. At that, the physical properties of the fluid in each element of the network depend on the component composition and are set according to the mixture rule. In the framework of the network approach, the passive scalar transfer technique was implemented for the first time. Phase transfer simulation was carried out by an explicit method with an adaptive time step.
The developed mathematical model of the two-phase (immiscible) fluid flow in a porous medium based on the network hydrodynamics methods was implemented in the form of a code in the C++ programming language. Detailed testing and verification of the A quantitative comparison of the computation with the experiment was carried out in terms of the time of saturation with water: where V w is the volume of water in the microchip, and V m is the volume of the microchip. The volume of the microchip is calculated as the total volume of the throats and pores, and the volume of water is calculated as the sum of the elements filled with water. If the network element is not filled, then only a part of the volume corresponding to the water phase is added to the volume of water. Figure 15 shows the time-dependent saturation process. As is seen, there is a good match between the calculated and experimental data. The deviation is mainly due to the limitations of the model, mentioned above. Quantitatively, the deviation of the saturation value in the computation from that in the experiment is ≈5%.

Conclusions
A mathematical model of the two-phase (immiscible) fluid flow in a porous medium has been developed based on the network hydrodynamics methods. The hydrodynamic

Conclusions
A mathematical model of the two-phase (immiscible) fluid flow in a porous medium has been developed based on the network hydrodynamics methods. The hydrodynamic part of the model is based on the well-established theory of hydrodynamic chains (graphs), which provides a significantly higher computational speed compared to the methods of computational fluid dynamics (CFD). The phase interface is taken into account by the contribution of capillary forces to the momentum conservation equation in the network branch. At that, the physical properties of the fluid in each element of the network depend on the component composition and are set according to the mixture rule. In the framework of the network approach, the passive scalar transfer technique was implemented for the first time. Phase transfer simulation was carried out by an explicit method with an adaptive time step.
The developed mathematical model of the two-phase (immiscible) fluid flow in a porous medium based on the network hydrodynamics methods was implemented in the form of a code in the C++ programming language. Detailed testing and verification of the developed model and the software code have been carried out for the problems of two-phase immiscible flow in straight channels and model porous media. The problems of oil displacement from the microchannel and microfluidic chips of various geometries and structures of the pore space were considered. The results of simulating using a network model are compared with the data of laboratory experiments. A good qualitative and quantitative agreement was reached between the experimental results and computations in terms of the time-dependent shape and position of the displacement front, the pressure drop, as well as the average saturation and displacement factor.
At this stage, just one of several possible modes of two-phase water-oil flow in a pore channel has been implemented considering the piston-like displacement mode. In certain cases, it is dominant, but not the only one; therefore, in the future, we plan to implement other modes of immiscible water-oil flows (film and slug flow).