Different-Scale Simulation of Flows in Porous Media

The paper considers the development of algorithms for an adequate description of processes of different scales in porous media. The choice of a computational technique is determined by the reference size of the problem being solved. Models of porous medium flow under Darcy's law, neglecting the medium microstructure, are used for the simulation at macro-scale. While at micro-scale, a direct description of fluid flow in porous channels with complex geometry by means of gas dynamic equations is used. In the first case the proposed model of non-isothermal multiphase multicomponent flow in a porous medium includes the mass balance and total energy conservation equations modified by analogy to the known quasi-gas dynamic equations. The model features are the introduction of minimal reference scales in space and in time and the change of the system type from parabolic to hyperbolic to increase the stability of explicit difference schemes applied for approximation. In the second case the dimensionless form of the quasi-gas dynamic system with pressure decomposition, developed by the authors earlier, is adapted to the simulation of flows in the pore space. The fictitious domain method is proposed to reproduce the core microstructure. The developed approaches have been verified by test predictions.


Introduction
The numerical simulation of fluid flows in porous media is one of urgent research topics all over the world. These studies are of great practical importance for the petroleum and metallurgical industries, since modeling the hydrocarbon recovery and minerals mining, calculations of groundwater dynamics, modeling operations in various technological units, for example, in catalytic reactors, as well as the digital analysis of small rock samples can reduce the number of expensive natural experiments. The flows under consideration are essentially multi-scale due to a large scatter of physical quantities and a variety of physical processes, such as convection, diffusion, phase transitions and reactions [1]. Therefore, models of porous media flows need further improvement for an adequate description of processes of different scales. This paper does not propose a single universal approach to modeling multiscale flows. The choice of a model and computational technique should be determined by the reference size of the applied problem being solved.
To simulate flows at macro-scale (e.g. modeling flows in the subsurface or in technological installations with reference sizes of the order of meters and kilometers) models of porous medium flow under Darcy's law are used [2][3][4].
To simulate flows at micro-scale (e.g. investigation of core samples with reference sizes of the order of micrometers and less) a direct description of fluid flow in the pore space of rocks by means of gas dynamic equations can be used [5,6].
In the both cases the present research is based on further development of the known quasi-gas dynamic (QGD) system of equations [7]. Governing models are transformed in accordance with physical peculiarities of flows under consideration. The developed approaches have been verified by a number of test predictions. We propose two approaches for modeling flows on a macro-scale within the framework of Darcy's law when the medium microstructure can be negligible. First, the original mathematical model of multiphase fluid flow has been constructed by analogy to the QGD system of equations [8,9]. The phase continuity equations were modified and got additional dissipative terms for the solution smoothing. The equations' type was transformed from parabolic to hyperbolic to increase the stability of explicit difference schemes used for the approximation.
The second approach is based on the classical model of compressible fluid flow in a porous medium [2][3][4] and the introduction of mass flux relaxation [10] into the continuity equations. Thus a hyperbolized system of equations has been also obtained [11]. The model allows for separation of the pressure and saturation variables and for the use of explicit schemes to determine the average pressure and the phase saturation. The both approaches apply logically simple computational algorithms aimed at efficient adaptation to massively parallel processing architectures.
The current version of the QGD-based hyperbolic model considers non-isothermal multiphase multicomponent flow through a non-deformable isotropic porous medium [12]. By a component a single chemical compound (a pure component) or a mixture of compounds (a pseudo-component) is assumed. Under phases we assume mobile phases, their number in the system is no more than three: these are water, Non-Aqueous Phase Liquid (NAPL) and gas, further subscripts w, n and g denote them respectively. The number of components n κ can be arbitrary, the same component may be present in different phases. Liquids are considered as slightly compressible, gas is ideal, the rock skeleton is incompressible. The temperature of all phases and the rock is considered identical. The mass and energy transfer between the phases occurs.
The mass concentration of component κ in phase α is presented as follows: where m κ α is the mass of component κ in phase α, m α is the mass of phase α.
The model consists of: • the mass balance equation for each component; • the extended Darcy's law for each phase; • the total energy conservation law; • state equations accounting for the multicomponent composition of phases; • capillary pressures as functions of saturations; • constitutive relations on the sum of saturations and the sum of concentrations; • constansts of phase equilibrium (so-called "Kvalues") specified for all components to be condensed. As the model naturally follows from the hyperbolized QGD model the second time derivatives and regularizers with small parameters τ and l are present in equations. In porous media the minimal reference length l is the distance at which the rock microstructure is negligible, 2 10 l : rock grain sizes. The minimal reference time τ is the time interval for inner equilibrium establishing in the volume of linear size l. The estimation / h c τ : is valid at numerical implementation [8] (h is the computational grid step in space, c is the sound speed in fluid).
A computational algorithm of the explicit type has been developed for numerical implementation of the model. The number of primary variables depends on the problem statement. For example, for the simulation of two-phase two-component flow the set of primary variables consists of the temperature, one phase pressure, the other phase saturation and the concentration of one component in a phase.
The algorithm specific feature is that equations for calculating the primary variables are not derived directly. The equations of conservation of mass and total energy are approximated by explicit difference schemes, and some auxiliary quantities are found from them on each time level. And then, solving the system of algebraic equations with these auxiliary quantities in the righthand sides, one can obtain the primary variables locally in each computational point.
Derivation and substantiation of the model as well as detailed description of the computational algorithm and its verification via test predictions can be found in previous works by the authors [8,9,11,12].

Test predictions
One of macro-scale test problems concerning phase transitions is the heat pipe effect simulation [13]. A horizontal column filled by a porous medium saturated with water and air is considered. This is a two-phase two component system: the water phase consists of the water component only while the gas phase consists of two components -water (as steam) and air. We choose the gas pressure P g , the water saturation S w , the temperature T and the concentration of water in the gas phase C g w (1) as primary variables. Figure 1 illustrates the problem statement.
The constant heat flux of 1000 W into the column and zero fluxes for all mass components are given at the right-hand boundary. In comparison with this test calculation from [12] the conditions are changed, in particular, the heat flux is reduced. As the water component may evaporate at heating and condense at cooling, the K-value is specified for it as follows [14]: where ω w is the water molecule acentric factor, T c and P c are critical values of the temperature and the pressure. Van Genuchten relationships [4] are used to describe relative phase permeabilities and capillary pressures.
The following process is reproduced as a result of computations: due to the heat flux the system is heated and some amount of the water component passes from the water phase to the gas phase. Figures 2-5 show the obtained numerical results at different time moments. One can observe the expected water-gas phase transition: the water saturation is decreased at the right boundary (see Fig. 4) while the concentration of water steam in the gas phase is increased respectively (see Fig. 5). As 1 w w C = during the process, w g C is completely determined by the exponent from formula (2).    The steam flows away from the high temperature region (see Fig. 2) under the gas pressure gradient (see Fig. 3). It condenses after reaching a cooler region near the left boundary. One can conclude that the proposed computational approach ensures correct predictions of multiphase multicomponent porous media flows.

Adaptation of QGD equations with decomposed pressure for modeling flow in pores
When modeling processes at micro-scale, in particular, in problems of core research, one of the main modern trends is direct numerical simulation of fluid flow in porous channels with complex geometry using gas dynamic equations (see, for example, [15]). The paper proposes to adapt for these purposes the system of equations with pressure decomposition developed earlier by the authors [16].
A well-known problem arises when simulating flows of weakly compressible fluids, which often include flows in porous media. Pressure changes are proportional to density changes times the Mach number squared. At low Mach numbers this leads to the pressure-density uncoupling. To overcome this problem, the authors developed an algorithm that treated the pressure instead of the density as the primary variable: the pressure was obtained via a system of equations, and the density was found from the equation of state [16].
In view of specific pressure behavior the pressure is decomposed into a sum of two components -the volume-averaged part ( ) p t and the dynamic part with reference values 0 P and 2 0 0 U ρ correspondingly: Then we get the dimensionless pressure: As a result, the system of QGD equations with pressure decomposition in the dimensionless form looks as follows (for brevity, the spatially one-dimensional isothermal case is given): The system does not contain singularities at 0 M → and therefore it can be used even for 0 M = . One can introduce the pressure increments (n is the time level number) 1 1 , n n n n p p p and derive finite difference equations for obtaining the intermediate values of these increments [16]. The equation of elliptic type for the intermediate value of the dynamic pressure increment has the next form: Since it is supposed to investigate flow in a porous medium, the algorithm should be adapted to computations in regions of complex shape. For this, it is proposed to use the method of fictitious domains. The entire computational domain is enclosed in a rectangle, in which a rectangular computational grid is set. On this grid, a map of the area is compiled: the indicator i map in the grid cell equals 1 if the cell belongs to the computational domain, that is, to the flow area, and 0 otherwise. Then the difference derivatives of the convective terms can be written in the form: and the pressure gradient in the form: Thus, difference equations with indicators are written uniformly throughout the domain so that the no-flow condition is automatically satisfied on a wall.

Test predictions
The developed algorithm was approved by the problem of flow around an unconnected complex-shaped body (see Fig. 6). The computational domain is a rectangular channel with rigid impermeable walls. On the left, a stream flows into the channel, the horizontal velocity has a parabolic profile with the maximum equal to 1.5. Inside the channel there is an impermeable area, consisting of two unconnected segments. The no-flow conditions are realized on all walls using the fictitious domain method described above, the grid size is 300 × 100. For the outgoing stream, the condition / 0 u x ∂ ∂ = is set. Figure 7 shows the results of calculations -the velocity fields at different Reynolds numbers. At small Re flow is more laminar, at larger Re the signs of turbulence are clearly manifested, the vorticity corresponds to the physics of the process. Currently, calculations are carried out in domains with a complex internal structure imitating rock samples.

Conclusions
The developed computational techniques have a wide range of applications. At macro-scale they include, for example: • the simulation of thermal methods of oil recovery to increase the production rate of highviscosity oil; • the solution of ecology problems on the thermally enhanced remediation of contaminated soils; • the simulation of processes in technological units at the purification of crude oil or natural gas. At micro-scale , they include the digital core analysis becoming an alternative to field and laboratory experiments in the oil and gas industry to obtain reliable data on reservoir properties and, as a result, to improve the macro-level models.