Application of Regularized Hydrodynamic Equations for Direct Numerical Simulation of Micro-Scale Flows in Core Samples

The paper is devoted to numerical simulation of three-dimensional isothermal two-phase twocomponent viscous fluid flows with surface effects in the pore space of core samples. The voxel representation of the flow domain is used suitable for digital rock physics applications. The flow is described by viscous compressible Navier-Stokes-Cahn-Hilliard equations. In order to use simple and computationally efficient explicit numerical algorithms with central difference approximations of spatial terms, a quasi-hydrodynamic regularization of equations is used. Simulation results of fluid displacement in pore space of realistic core sample (sandstone) are presented. The results demonstrate the applicability and good prospects of quasi-hydrodynamic regularization technique to solve the problem.


Introduction
With the development of computed tomography methods and computer technologies, mathematical modeling is applied in the ascendant to analyze physical processes and determine different properties of materials whose microstructure is unknown in detail beforehand. In the geophysical applications, all the problems concerning the determination of the properties (including filtration characteristics) of rocks by numerical simulation techniques on the basis of their micro-tomographic image are usually called as "digital rock physics" [1].
In the last decade direct numerical simulation of pore-scale flows become an actively developing direction in mathematical modeling and numerical simulation.
Historically the first class of approaches capable to describe displacement process at such fine space resolution is based on the concept of pore network model. Essentially, the pore network is described as a geometrical graph which vertices and edges describing capacitance (pores) and transport (throats or bonds) features of the sample. The vertices of the pore-network graph correspond to separate pores, and its edges describes throats (bonds) which connect separate pores and control transport of mass, momentum and heat. Geometrical graph describing pore network geometry and topology is usually extracted from the voxel image using special geometrical and topological processing techniques (which, naturally, have a finite accuracy). The voxel image itself can be obtained from micro-CT experiments, object-based or statistical simulation methods. Contemporary pore-network type models can take into account multiphase fluid composition [2], a non-Newtonian rheology [3] and other effects [4][5][6][7][8], including the combustion processes [9].
The main drawbacks of the pore network models are (i) a need for (re)construction of the geometrical pore network which naturally can be thought as an excessive step as soon as detailed experimental data is available using CT, micro-CT or MRI techniques and (ii) development of suitable flow model for pore-network, which is a nontrivial task as soon as underlying physics is "rich" enough.
These consideration led to the development of a new approaches for direct flow simulations using CT based voxel images with none or little preprocessing. The flow itself is described by appropriate hydrodynamical model, e.g., Navier-Stokes equations. The most popular among numercial solution techniques, due to its relative simplicity, parallel scalability and quite large simulation community experience, is lattice Boltzmann method (LBM) [10,11]. Another method used for pore-scale modeling is smoothed particle hydrodynamics method [12][13][14][15]. This approach is known to be quite robust, but difficult to apply for simulation of large models. Its drawbacks include extended computational resources requirements and need for calibration of the method parameters [16]. Hybrid "pore scale" -"continuum scale" models which require more modest computational resources also have been developed [15][16][17][18].
In spite of wide application the lattice Boltzmann method, there are still a number of reasons which limit applicability of these approach. Among them we mention models and simulation techniques suitable for reliable multiphase and multicomponent flow simulations with large contrast in phases' density and viscosity, consistent treatment of reactive solid phase and changes in pore and grain geometry due to chemical reactions, dissolution and precipitation, geomechanical coupling, etc. [16,19].
Research studies conducted in Keldysh Institute of QGD approach as technique for construction of both new mathematical models of hydrodynamics and an efficient computational techniques has been developed since mid of 1980-th. Results of these research efforts were published on numerous papers and summarized in recent monographs [20][21].
In [22], the authors applied the QGD model to the direct numerical modeling of the flow of a one-phase fluid in the pore space of rock samples in order to determine the absolute permeability coefficient.
In this work, the results of simulation of the isothermal flow of a two-phase two-component fluid in the pore space of a core sample are presented. To make it possible to use relatively simple finite-difference schemes, we used the Navier-Stokes-Cahn-Hilliard equations with QHD regularization supposing that small dissipative terms are added to the original equations [20,21,23]. It should be noted that the notion of phase is not introduced explicitly in the considered model. However, due to the special form of the free energy and the terms containing gradients of order parameters, the subregions occupied by a mixture with a nearly homogeneous composition are formed in the regions of flow. This allows us to consider that these regions are occupied by a particular phase determined by the corresponding component composition (nearly homogeneous).

Mathematical model of twocomponent flow with surface effects
The basic equations of the model are the following where  , u are the mixture density and the fluid velocity; C is the mass concentration of one of the fluid components. These parameters depend on ( , ) t x ,  x ,  is a bounded region with the piecewise smooth boundary  .
Let us describe the Helmholtz free energy of the mixture by the formulae [24,25] where si c is the speed of sound of the i -component, (2)  are the free energies of components, 1 > 0 The tensor = NS    Π Π Q Π is the sum of the Navier-Stokes viscous stress tensor, the cappillary stress tensor Q , and the regularizing tensor  Π : where > 0 where n is the inner normal to  . The first two conditions in (4) ensure the absence of the mass flow on the boundary. The third condition in (4) corresponds to the absence of the mass flow of the components through the solid wall (rigid impermeable boundary). The last condition in (4) means that the angle between any level set of the function C and the solid wall at the point of their intersection is / 2  . The level set = 1/ 2 C is considered by convention as the interphase boundary.
It should be noted that in [25], the solid wall was not considered; periodic conditions were fixed on the boundaries. In this work, we consider the threedimensional problem with the presence of a solid wall, which is necessary for modeling of the fluid flow in a rock sample. For equation system (1)-(3), we applied the explicit finite-difference approximation in which all space derivatives were approximated with central differences.

Results of the numerical experiments
The numerical experiments on simulation of the considered flows were carried out using digital models of core samples which are freely accessible at the site of Imperial College London [26]. The presented results are obtained by using the sample sandstone_9. In Fig. 1, the three-dimensional binary voxel image of this sample is shown. The geometrical model of the flow region is set up by the voxel "image" of the sample on the Cartesian orthogonal grid. Fig. 2 illustrate the successive steps of the computing experiment in which the red fluid (phase) ( > 0.5 C ) displaces the blue ( 0.5 C  ) fluid (phase). They show the saturation of the sample by the red fluid S , which is defined as the ratio of the cells with concentration > 0.5 C to the total number of cells. It worth mentioning that a phenomenon of residual saturation was observed: the blue fluid is not completely "washed out" from the pore space of the sample and the saturation attains the constant value 82% S  .

Conclusions
Nowadays, the diffuse interface (phase field) methods in hydrodynamics are one of the major thermodynamically consistent ways for description of complex multiphase multicomponent flows accounting for dynamics of phase interface and multiphase effects. The model equations are high-order nonlinear differential equations. For this reason, realistic applications of the model require complex numerical algorithms.
In this paper, we propose to use the QHD regularized phase space hydrodynamics equations. The regularization consists in adding small physically based and thermodynamically consistent dissipative terms, which allows us to construct relatively simple numerical algorithms based on the central finite-difference approximations of differential operators. The results of simulation of microscale two-phase two-component flows in the pore space of core samples are presented. The results demonstrate the applicability and good prospects of quasi-hydrodynamic regularization technique, as well as corresponding algorithms and program implementation to solve the problem under consideration.