Non-Isothermal Compositional Model for Simulation of Multiphase Porous Media Flows

The research carried out in this work is aimed at the development of an explicit computational algorithm to simulate non-isothermal flow of multiphase multicomponent slightly compressible fluid through porous media. Such a modeling is necessary, for example, at solving practically important problems of hydrocarbon recovery via thermal methods and contaminated soil remediation via the hot steam injection. The original model and algorithm developed by the authors earlier draw the analogy with the quasi-gasdynamic set of equations. Phase continuity equations were modified to get regularizing terms, while the type of equations was changed from parabolic to hyperbolic for the approximation by an explicit scheme with a mild enough stability condition. Since an adequate description of temperature-dependent porous media processes requires an explicit consideration of the transfer of mass and energy between the phases, the present work focuses on the generalization of the proposed approach to the case of multicomponent composition of fluids. Conservation laws are currently formulated for the components in terms of the mass concentrations of components in the phases. Constants of the phase equilibrium were used to close the set of equations. The developed model and algorithm have been verified numerically via test predictions of twoand three-phase flows, and, the phase transition of water into the gas phase at heating has been reproduced correctly.


Introduction
Many industrial and environmental applications imply modeling of complex fluid flow in the subsurface: as examples one can point out problems of the hydrocarbon recovery or ecological problems concerning the contaminated soil remediation. As the corresponding technologies widely use thermal methods, in these cases non-isothermal multiphase fluid flow through porous media are to be predicted. An adequate description of temperature-dependent porous media processes requires considering the transfer of mass and energy between the phases. Therefore the governing model should account for the multicomponent composition of fluids. Besides algorithms of the model numerical implementation should allow efficient high-performance computing because the problems under consideration are large-scale and their solution is impossible without the use of modern supercomputers.
From the standpoint of parallelization, explicit difference schemes are the most attractive. In cases when calculations with the critical accuracy are really needed, the computational grid steps on space as well as on time are strongly restricted and explicit methods can gain over implicit ones [1]. Following this topic the authors proposed an approach to porous medium flow modeling based on the analogy with the quasigasdynamic (QGD) system of equations [2]: the mass balance equation for each fluid phase was modified to get a regularizing term and to change the equation type from parabolic to hyperbolic [3]. Note that hyperbolization is a modern trend in CFD [4][5][6] to increase the stability of employed explicit schemes. The QGD-based model includes also the total energy conservation equation [7,8].
The present paper reports an extension of the nonisothermal hyperbolic QGD model of porous medium flow to describe multicomponent fluids. Among primary variables to be obtained by the developed computational algorithm of the explicit type there are the mass concentrations of components in the phases.
The created approach has been verified by a number of test problems involving two-and three-phase flows, physically correct results have been observed.

Governing model of multiphase multicomponent fluid flow in a porous medium
We consider multiphase multicomponent flow through non-deformable homogeneous isotropic porous medium. By a component a single chemical compound (a pure component) or a mixture of compounds (a pseudo-component) is assumed. Let n α be the number of phases and n κ be the number of components in the system under consideration. Typically, the number of phases is no more than four (water, Non-Aqueous Phase Liquid (NAPL), gas and a solid phase), the number of components can be arbitrary. One and the same component may be present in different phases. The relative amount of the component in the phase can be expressed either in mass or in molar concentrations. The mass transfer by components between the phases occurs. We assume that components of the solid phase cannot be present in the fluid phases and vice versa. For example, the water component is a pure component, it can exist in aqueous, non-aqueous and gas phases. Air is simplified to a single pseudo-component, its composition is neglected. Further speaking of phases we mean fluid phases, subscripts w, n and g denote water, NAPL and gas phases respectively. Speaking of components (superscripts in notations) we mean components of the fluid 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 to be identical.
We choose the mass concentration of component κ in phase α as a primary variable presented as follows: where m κ α is the mass of component κ in phase α, m α is the mass of phase α.
The hyperbolic QGD model of porous medium flow proposed and substantiated by the authors in previous works [3,7,8] has been elaborated. If to pass from conservation laws for phases to conservation laws for components and to determine lacking closure relations one can formulate the next compositional model: The following notations are used: S α is the saturation (of α-phase), P α is the pressure, ρ α is the density, u α is the Darcy velocity, T is the temperature, E α is the internal energy, H α is the enthalpy, λ eff is the effective coefficient of heat conductivity, ϕ is the porosity, k is the absolute permeability, k α is the relative phase permeability, µ α is the dynamic viscosity, Q κ is the source of component κ, g is the gravity vector, ρ r = const -the rock density, c α is the sound speed in α-phase, small parameters l and τ are the minimal reference length and time respectively, Note that the derivation of kinetically-consistent finite difference (KCFD) schemes and the related QGD system is based on the principle of minimal sizes [2,9]. This principle states that at the numerical solution of continuum mechanics problems it makes no sense to consider scales smaller than some minimal reference values. In porous media the above mentioned l (the minimal reference length) is the distance at which the rock microstructure is negligible, 2 10 l : rock grain sizes. The minimal reference time τ is assumed to be the time interval for inner equilibrium establishing in the volume of linear size l. The next estimation is valid at numerical implementation [3]: / h c t : , where h is the computational grid step in space, c has the order of the sound speed in fluid.
The above model consists of • the mass balance equation for each component (2) (it accounts for the convective mass transfer and does not take into account the diffusive flux of the component yet), • the extended Darcy's law (3) (all components of the phase have one and the same velocity equals the phase Darcy velocity), • the total energy conservation law (4) (the temperature is assumed to be the same for all components), • the state equations (5) (the phase density depends on the multicomponent composition), • differences between the pressures (6) (they are defined as capillary pressures), • constitutive relations (7) and (8) (they follow from the corresponding variables' definitions), • the constansts of phase equilibrium (9) (they describe the multicomponent composition of the phases under the condition of phase equilibrium and close the system of equations).
K-values should be specified for all components to be condense, for example, if the water component according to the problem statement may evaporate and condense (phase transitions "water -gas" and "gaswater" happen) then the K-value is as follows: where ω w is the water molecule acentric factor, T c and P c are known critical values of the temperature and the pressure at which water turns to steam and they have identical properties. Some K-values can be found in [10]. Equation (4) includes the next quantities: , , , H H w n g κ α α κ α = = ∑ (12) where components' enthalpies H κ α can be found via heat capacities of the components.
The connection between the internal energy and the enthalpy is expressed by the relations: subscript r denotes the rock.
Phase densities (5) where β w is the coefficient of isothermal compressibility, η w is the coefficient of thermal expansion, the terms with subscript 0 are reference values of the corresponding variables.
It is known that the main nonlinearities in the porous medium flow equations are associated with the form of the dependence of relative phase permeabilities and capillary pressures on saturations. Generally, the capillary pressure and the relative phase permeability are temperature dependent functions but we neglect the temperature dependence and consider their dependence on saturations only. In the present research for the simulation of three-phase flow Parker's model [12] is chosen to describe capillary pressures (6), the relative phase permeability is presented by Stone's Model I [13], these models are also cited in [3,7]. For the two-phase flow van Genuchten constitutive relationships [14] are used.

Computational algorithm
A computational algorithm of the explicit type was developed for numerical implementation of the QGDbased model of fluid flow in a porous medium [3,7]. When constructing the algorithm, rectangular computational domains covered by Cartesian grids were considered. High efficiency of parallel implementation on CPU cores as well as on GPUs of a hybrid supercomputer was demonstrated while solving infiltration problems.
The algorithm is naturally generalized to predict multicomponent fluid flow. In this case the set of primary variables is enlarged due to the concentrations of components in the phases (1). Their number depends on the problem statement.
For simplicity let us formulate the algorithm for the simulation of non-isothermal two-phase two-component flow. Let the phases be water and gas, the components be water and air (superscripts w and a respectively). The water phase consists of one component -water, the gas phase consists of two components -water (as steam) and air. This composition is considered, for example, at solving the test problem on the heatpipe effect [15]: at heating water can evaporate, at cooling steam can condense.
The primary variables are P g , S w , T and C g w , the initial and boundary conditions are set for them. On each time level j = 0, 1, 2, ... the next main stages are fulfilled: • Calculation of all necessary dependent variables and coefficients including  j g S from (7),   (13) and K-value (9) (more precisely (10)). The just found values C κ (16) and E (17) are used in the right-hand sides, equality 1 w w C = is taken into account: This system is solved by Newton's method that takes only a few iterations.
• Data exchange at multiprocessor computing on inner boundaries of subdomains while the principle of geometric parallelism (i.e. data partitioning) is applied.

Test predictions
The hyperbolic QGD model of porous medium flow and algorithm of its implementation were successfully used to solve test and model problems, including those on hybrid supercomputers [16]. For single fluid flow verification was carried out using the problem of planeradial inflow to a well, the numerical solution almost coincided with the exact analytical one [3]. For verification in the two-phase case the authors used the water drainage problem studied in [14]. A rather good agreement was observed at comparison of the obtained results [8,17] with results from [14]. The physical correctness of results obtained at modeling three-phase isothermal and non-isothermal flows was qualitatively evaluated [7,8].
One of popular test problems in CFD is the problem on natural convection of air in a differently heated square cavity [17]. The similar problem has been solved at the assumption that the cavity is filled by a porous medium saturated with three-phase (water-oil-air) fluid. Earlier the problem of the phases' redistribution under the gravity was solved by the authors in the isothermal statement [3]: eventually fluids distributed over the domain according to their densities. If to apply the temperature gradient to the horizontal (top and bottom) boundaries then the infiltration process is accelerated, the corresponding results are presented in [7].
But now we are interested in the case when the temperature gradient occurs between the vertical walls,   The pressure distribution quickly becomes hydrostatic and remains almost unchanged over time.   The saturation patterns are asymmetric. One can see the mixed effect of the gravity and the temperature difference. Gradually gas accumulates at the top of the warm wall (at the upper left corner), water moves down to the cold wall (to the lower right corner), oil tends to occupy the middle part of the region (formation of the oil "tongue" is noticeable) but it strongly clings to the cold wall.
Another test problem is the heatpipe effect simulation to verify the proposed compositional model (2)-(15) and the algorithm presented in Sect. 3. We have already referred to this two-phase two-component flow problem at the algorithm formulation. A similar problem is discussed in [15], Fig.7 illustrates the problem statement. Due to the heat flux the water-air system is heated until water boiling, then steam is produced at the rigthhand boundary, thus the phase transition happens: some amount of the water component passes from the water to the gas phase. Fig. 8 and 9 show the obtained numerical results at some time moment. One can observe the described phase transition: the water phase saturation is decreased while the concentration of water steam in the gas phase is increased.  The steam flows away from the heat source, it condenses after reaching cooler regions near the lefthand boundary. Such heat transfer processes as convection, conduction and diffusion as well as capillary forces play an essential role in this experiment. One can conclude that the proposed computational algorithm ensures qualitatively correct predictions of multiphase multicomponent porous media flows.

Conclusions
The model of non-isothermal flow of multiphase multicomponent fluid in a porous medium has been formulated on the basis of the earlier developed hyperbolic QGD model of porous media flows. The algorithm of its numerical implementation has been described for the case of two-phase two-component fluid. All test predictions indicate the adequacy of the developed kinetic approach to porous medium flow simulation.
The proposed approach can be used to model promising thermal methods of oil recovery (heat carrier pumping into the stratum e.g.) aimed at increasing the oil production rate of difficult-to-recover hydrocarbon reserves. The approach is also useful for the solution of such environmental ecology problems as the thermally enhanced remediation of soils contaminated by NAPLs (halogenated organic solvents or hydrocarbon fuels e.g.).