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

Abstract

Multiphase flows 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 flow 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 flow and pressure fields 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 verification 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.

1. 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 used for the application of numerical research methods [10]. In turn, these numerical research methods can be divided into two main approaches: direct 3D simulation (as a rule CFD) [11,12,13] and network (pore-net) simulation [14]. The advantage of the first approach is a detailed understanding of the process, and the disadvantages are the ability to simulate only a very limited core size. However, 3D modeling allows you to consider various factors that affect the process of oil displacement. So recently, one of the methods for increasing the oil return is the use of nano-liquids and this process can be modeled in different ways. For example, by studying the physical properties of a nanofluid experimentally and setting them into a model [12], or trying to construct a simulation of particle motion by an explicitly Langrage method [13]. The second approach allows simulating a core consisting of several million pores [15]; however, it requires a complex procedure for converting the results of core tomography into a network model (based, for example, on the concept of maximal balls [16], or more complex hierarchical approaches [17]). The pore-scale network represents the space of the rock. Wide areas of space are called “pores” (hereafter, pores), which are connected by narrower areas called “throats” (hereafter, throats).

Two types of network models can be distinguished for simulating multiphase fluid flows: quasi-static and dynamic [18]. Quasi-static models are characterized by the predominance of capillary forces and the absence of viscosity gradients. Most of these models are aimed at tracking changes in the system’s state as the saturation of one phase increases. When the assumption of the dominance of capillary forces over viscous ones is no longer valid, it is necessary to use a dynamic network model. Dynamic network models take into account viscous forces explicitly: a given inflow velocity is introduced for one of the fluids, and the transient pressure characteristic as well as the corresponding position of the interface is calculated. Dynamic models of a pore network differ in three main characteristics: the network structure, the geometry of the elements, and a computational algorithm for determining the pressure field. The network structure is characterized by the location of nodes, as well as the number and orientation of links. If the nodes are located at the vertices of a regular lattice, the network is called structured, otherwise, it is unstructured. If all nodes have the same number of incident branches, the network is isotropic, and otherwise, it is anisotropic. In terms of the geometry of the elements, several main options can be distinguished. In the first option, the entire volume is concentrated in the pores but they do not have hydraulic resistance, while the throats, on the contrary, have a small volume but they resist the flow. This assumption can help save computing time and memory during simulation since there is no need to track the interface inside the throats. Moreover, the pressure drop inside porous bodies can be neglected [19]. In the second option, pores and throats have both volume and resistance. The third option concerns the case where pores and throats are combined into a single element, which is characterized by both volume and resistance. The combined elements may have different cross-sections. Connection points are not assigned any specific properties, such as volume and resistance [20]. The cross-section shape of the throats also plays an important role in forming the network model. As a rule, for the convenience of computation, a regular-shaped section, namely, round, triangular, or square is used. The fundamental difference in the cross-sections of throats with corners from those without corners (having round cross-sections) is the constant presence of the wetted fluid remaining in the corners [21] and, as a consequence, the film flow of the fluid.

The pressure field in the dynamic network model can be simulated in two ways: with a single medium pressure and pressures, specified for each phase. The first option is used under the assumption that each pore or throat section (for example, with a circular section) or network is filled with an equivalent fluid, and each throat is assigned an equivalent conductivity. In the case where the section of the throat has corners with a constant presence of wetted fluid, then for correct simulation of the film flow, models can be used in which the pressure for each phase is calculated individually.

The use of network methods for modeling the process of movement of a multiphase fluid in a porous medium continues to develop along the way, complicating models to 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 multi-scale, 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.

2. Mathematical Model of Network Simulation of Oil Displacement Processes

Before describing the mathematical model, we introduce several notions and limitations:

• 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.
• 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.
• 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.
• Fluids do not mix with each other due to a rigid phase interface.
• However, in network elements, where both components are present, the physical properties of the fluid (density and viscosity) are determined by the mixture rule.
• 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:

2.1. 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.

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).

$$D_{il}=\left\{\begin{array}{c}1, if l\in O_i\\ -1, if l\in I_i\\ 0, otherwise\end{array}\right.$$

(1)

$${\displaystyle \sum }_{l\in U_i}{D}_{il}{q}_l=Q_i, i\in N$$

(2)

$$s_l\left|q_l\right|q_l={\displaystyle \sum }_{i\in N}{D}_{il}·p_{Di}\pm P_c, l\in U$$

(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:

$$s_l=\left(\frac{{\lambda }_{\mathrm{fr}}·l}{d}+\xi \right)\frac{\rho }{2·f^2}$$

(4)

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:

$$p^k=p^{k-1}+p^{\prime }$$

(5)

where $p^k$ is the pressure value at the current iteration, $p^{k-1}$ is the pressure value at the previous iteration, $p^{\prime }$ is the pressure correction. To obtain the pressure correction equation for the network problem, we write down the difference between Equation (3) for the current and previous iterations and substitute the result into Equation (1):

$${\displaystyle \sum }_{l\in U}{D}_{jl}\left({\tau }_{\mathrm{net}}·{\displaystyle \sum }_{i\in N}{D}_{il}{p}^{\prime }\right)=q_j-{\displaystyle \sum }_{l\in U}{D}_{jl}{q}_l^{k-1}$$

(6)

here ${\tau }_{\mathrm{net}}$ are the coefficients that are obtained in the approximation of Equation (3).

The system of Equation (6) can be represented in matrix form:

$$A{p}^{\prime }=∆q$$

(7)

where A is the matrix of coefficients, $p^{\prime }$ 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⁻³, 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 ${\rho }_1$ and ${\mu }_1$ are the density and viscosity of the first fluid, ${\rho }_2$,${\mu }_2$ are, respectively, the density and viscosity of the second fluid, and φ is the volumetric fraction of the first phase:

$$\rho ={\rho }_1\phi +\left(1-\phi \right){\rho }_2,$$

(8)

$$\mu ={\mu }_1\phi +\left(1-\phi \right){\mu }_2.$$

(9)

Phase transfer along the length of the throat is modeled by solving a one-dimensional convective equation with an up-wind approximation.

$$\frac{\partial \phi }{\partial t}+\frac{\partial u\phi }{\partial z}=0,$$

(10)

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.

2.2. 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:

$${\phi }_i^{new}={\phi }_i^{old}+\frac{\left(\sum q_{in}-\sum q_{out}\right)}{V_i}·∆T,$$

(11)

where ${\phi }_i^{new}$ is the volumetric fraction at the next time step, ${\phi }_i^{old}$ is the volumetric fraction at the previous time step, $\sum q_{in}$ and $\sum q_{out}$ are the total flows bringing/carrying away the phase into/out of the element, respectively, $V$ is the volume of the element, $\Delta 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,$\sum q_{in}-\sum q_{out} \ne 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:

• If the pore is not filled with the displacing fluid, then only the displaced fluid comes out of the node (Figure 4a).
• 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.2.1. 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.2.2. 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:

$$P_c=\frac{2·\sigma ·\cos \left(\theta \right)}{r}$$

(12)

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.

2.2.3. Computational Algorithm

Considering the above-described assumptions, the following computation algorithm can be constructed, consisting of the following steps:

• Setting the initial distribution of pressure, flow rates, and phase concentrations over the network.
• Determining the list of network elements in which the phase replacement process at the first time step is possible.
• Calculating the density and pressure for all network elements by Equations (5) and (6).
• Calculate (by Equation (12)) and set the pressure for the branches corresponding to capillary pressure.
• Solving the hydraulic problem, which results in obtaining a new distribution of pressure and flow rates.
• Calculating the pressure field for the nodes.
• Determining the flow rate in the branches using the above calculated pressure field.
• Prohibiting the reverse displacement of the fluid in the network branches.
• Repeating the operation starting from point 5a until the number of blocked branches becomes zero.
• $\mathrm{Determining} ∆T_{\min }$.
• Recalculating the phase concentrations using Equation (11).
• Conducting checking for the presence of closed clusters.
• Updating the list of network elements in which the phase replacement process is possible at a time step.
• 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.

3. 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.

3.1. 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).

The physical properties of the oil were as follows: $\rho =850 \mathrm{kg}/{\mathrm{m}}^3$ and $\mu =0.02 \mathrm{Pa}$; water properties: $\rho =998.2 \mathrm{kg}/{\mathrm{m}}^3$, $\mu =0.001 \mathrm{Pa}·\mathrm{s}$. The surface tension coefficient was $\sigma =0.03 \mathrm{Pa}·\mathrm{s}$. When replacing oil with water, the average properties of the liquid change and the pressure drop along the length of the microchannel decreases. Thus, for conducting 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:

$$P_c=\frac{2·\sigma ·\cos \left(\theta \right)}{r}=\pm \frac{2·0.03·\cos \left(40\right)}{0.00025}=\pm 0.184\mathrm{kPa}$$

(13)

These analytical values of capillary pressure were used for comparison with the calculated values (see

Table 1. Total and capillary pressure drop along the length of the microchannel depending on the wetting angle at a time of 1.5 s.

Wetting Angle	Calculated pressure drop, kpa	Analytical Pressure Drop, kpa	Discrepancy of Total Pressure Drop, %	Calculated Difference with the Option with a Wetting Angle of 90°, kPa	Analytical Difference with the Option with a Wetting Angle of 90°, kPa	Discrepancy of Capillary Pressure Drop, %
40	3.437	3.437	0	0.186	0.184	0.01
90	3.251	3.253	0.06	0	0	0
140	3.07	3.069	−0.03	−0.181	−0.184	−0.02

).

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.

Figure 8 shows the pressure distribution along the length of the microchannel and the position of the waterfront at different time points (1 s, 2.5 s, and 4 s) for the variant with a wetting angle of 40°. The dotted line shows the position of the front obtained by the analytical method. The figure shows the coincidence of the calculated and analytical position of the displacement front. Figure 8 also shows a jump in capillary pressure at the displacement front.

3.2. 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\mathsf{\mu }\mathrm{L}/\mathrm{m}\mathrm{i}\mathrm{n}$. 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, 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. Parameters of a porous microfluidic chip.

Parameter	Value
Microchannel cross-section	100 × 110 µm
Cross-section of the pores	Ø85 and Ø63 µm
Inlet channel length, including bifurcations	27.7 mm
Outlet channel length, including bifurcations	99.2 mm
Total length of the porous channel	4800 mm
Inlet channel volume	0.9 µl
Outlet channel volume	3.2 µl
Porous space volume	38 µl
Roughness of channel surface	5 nm
Channel coating	hydrophilic
Microchannel cross-section	100 × 110 µm

.

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:

$$∆P=\frac{8·\xi }{{\pi }^2·\rho ·d^4}·u_l^2+\frac{128·\mu ·l}{\pi ·\rho ·d^4}·u_l^{}$$

(14)

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:

$$∆P\left(u\right)=a·u^2+b·u$$

(15)

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: ${\xi }_{\varnothing 85}=1.29$ and ${\xi }_{\varnothing 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).

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:

$$S=\frac{V_w}{V_m}·100\%$$

(16)

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%.

4. 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).