An Explicit Algorithm for the Simulation of Fluid Flow through Porous Media

The work deals with the development of an original mathematical model of porous medium flow constructed by analogy with the quasigasdynamic system of equations and allowing implementation via explicit numerical methods. The model is generalized to the case of multiphase multicomponent fluid and takes into account possible heat sources. The proposed approach is verified by a number of test predictions.


Introduction
The development of mathematical methods and software for the simulation of complex fluid flow in the subsurface is one of urgent tasks of industrial mathematics. Among applications there are hydrocarbon recovery problems as well as ecological problems concerning the soil and groundwater contamination. Models and algorithms for their solution still need improvement to predict thermohydro-mechanical-chemical processes in strata with sufficient accuracy at reasonable computing time.
The work presents an original model of fluid flow through porous media constructed by analogy with the quasigasdynamic system of equations and allowing implementation via explicit numerical methods [1]. Large-scale oil recovery problems (e.g., combustion fronts, phase transitions, complicated functions of the relative phase permeability) require calculations with very small space steps what strictly constrains the time step. Then explicit schemes can surpass implicit schemes involved in standard solution methods like IMPES or SS [2]. Besides explicit algorithms are preferable for HPC.

Mathematical model and numerical implementation
For many problems of continuum mechanics the principle of minimal sizes holds: we assume the existence of minimal space and time scales which act as lower limits for description details. In porous media the minimal reference length l is a distance at which rock microstructure is negligible (l is of the order of hundred rock grain sizes), the minimal reference time τ is the time needed to reach inner equilibrium in the volume of linear size l. Starting from the basic flow equations for slightly compressible fluid [2], using the principle of minimal sizes and the differential approximation technique the next multiphase flow model has been derived (α = w (water), n (NAPL), g (gas); r denotes a rock): ∂ ∂t Here S α is the saturation, P α is the pressure, ρ α is the density, u α is the Darcy velocity, T is the temperature, E α is the internal energy, q α is the source of fluid, ϕ is the porosity, K is the absolute permeability, k α is the relative phase permeability, µ α is the dynamic viscosity, g is the gravity vector, c α is the sound speed in α-phase, β α is the coefficient of isothermal compressibility, η α is the coefficient of thermal expansion, P 0α , ρ 0α and T 0 are reference values, ρ r = const -the rock density. The phase continuity equation (1) is modified: it gets a regularizing term and a second order time derivative with small parameters. In [1,3] some estimations of these parameters are given. Particularly the value τ ∼ h/c is satisfactory (h is the space grid step). The equation type is changed from parabolic to hyperbolic, consequently a three-level explicit scheme with a rather mild stability condition can be used, convective terms are approximated by central differences. As the temperature of all phases and of the rock is identical the system involves (3) as a single equation of the total energy conservation, also modified by analogy with the quasigasdynamic system. This equation includes the effective coefficient of heat conductivity λ eff = ϕ α S α λ α + (1 − ϕ) λ r , and the enthalpy which is calculated via the heat capacity at constant pressure C P α : The connection between the internal energy and the enthalpy is expressed by the next relations: Capillary pressures are described by Parker's functions, the relative phase permeability is presented by Stone's Model I [2].
For numerical implementation of the above system an algorithm of the explicit type is proposed. P w , S w , S n and T are chosen as primary variables in the three-phase case. Main stages of the algorithm are explicit approximations of (1) and (3) as well as the solution of a system of nonlinear algebraic equations by Newton's method locally at each computational point [3].
In order to account for the complicated multicomponent structure of the fluid it is necessary to pass from conservation laws for phases to conservation laws for components. The relative quantity of j-th component in phase α can be expressed as the mass concentration: C jα = m jα /m α , α = 1, . . . , n α , j = 1, . . . , n c . Then the next equation represents the mass conservation law of each component: Equations (2) and (3) are valid for the multicomponent case but the dependences ρ α = ρ α P α , T, C jα and µ α = µ α P α , T, C jα should be noted.
The next closure relations are defined: where K jαβ are constants of the phase equilibrium.

Test predictions
For verification of the proposed model and algorithm a drainage test problem [4] has been solved. Two-phase (water/air) isothermal infiltration due to the gravity is simulated. The geometry is illustrated in the figure 1(a). Initially the given thin column is filled by a fully water-saturated porous medium. The top boundary is open to the atmosphere with a no-flow condition on water, water drains down from the column, a no-flow condition on air exists at the bottom.  ure 1(b)). A good agreement is observed that testifies to adequacy of the developed approach.
To check the full mathematical model of nonisothermal three-phase flow, the 3D problem of infiltration with a hot source has been investigated (figure 2). Initially in the cube reservoir the residual water saturation is set, oil and gas are distributed periodically, the water pressure equals the atmospheric pressure, the temperature is 285 K. A source of water (S w = 0.6, S n = 0.05) with the temperature of 320 K is set on the top. The top boundary is open to the atmosphere, the bottom is impermeable, side boundaries are permeable. Results at the moment of 2000 seconds are shown in figure 2. Water gradually forces out oil and gas and tends to the bottom, the temperature front differs from saturation fronts. Physically correct dynamics of the process is observed.

Conclusions
The created approach can be used for the simulation of perspective thermal methods of oil recovery (such as heat carrier pumping into the stratum) aimed at increasing the oil production rate of difficultto-recover hydrocarbon reserves.