Simulation of a gas-condensate mixture passing through a porous medium in depletion mode

One of important tasks in a development of gas-condensate fields is to minimize hydrocarbons loss arising from the gas condensation in pores of the gas-bearing layer. The search for the optimal gas production regime is carried out both on the basis of laboratory experiments and on the base of computer simulation. In this regard, the relevant is the verification of the constructed mathematical models by means of comparison of numerical results with experimental data obtained on the laboratory models of a hydrocarbon reservoirs. Within the classical approach on the basis of the Darcy law and the law continuity for flows, the model is formulated that describes the passing a multicomponent gas-condensate mixture through a porous medium in the depletion mode. The numerical solution of the corresponding system of nonlinear partial differential equations is implemented on the basis of the combined use of the C++ programming language and the Maple software. Shown that the approach used provides an agreement of results of numerical simulations with experimental data on the dynamics of hydrocarbon recoverability depending on the pressure obtained at VNIIGAZ, Ukhta.


Introduction
Recently, due to the decrease in easily recoverable natural gas reserves in traditional fields, the development of areas with unconventional hydrocarbon reserves, including gas-condensate fields, has attracted attention [1]- [3]. Gas production in such fields is difficult due to the presence of the condensate [4]- [6]. Methods for increasing the gas production based on reducing a condensation is reviewed in the recent paper [7]. The behavior of a multi-component gascondensate mixture in a porous medium can be difficult to predict. Indeed, the phase states of the components can change. Gas has a much lower viscosity than liquid, thus a boiling of any component gives it an advantage in the speed of passage. The movement of the mixture is caused by the pressure drop at the point of its extraction. The pressure drop leads, in turn, to the possibility of boiling, that is, the phase transition of any component from a liquid to a gaseous state. A paradoxical result may be the formation of a condensate plug if the gas thus "erode".
In this regard, the development of realistic mathematical models can help predict the physical and technological parameters that ensure optimization of the gas production regime in the gas-condensate reservoirs. One of actual problems in this direction is the verification of mathematical models developed on the basis of the comparison of numerical results with the data of practical measurements, including the data of laboratory experiments. Traditionally, an approach based on the classical Darcy law and conservation laws, described in detail in the literature, is used to model the processes of passage of multicomponent gas-condensate mixtures through a porous medium (see, for example, [8]- [13] and references therein). In this framework, one formulates the system of equations with appropriate initial and boundary conditions tailored to the specific modeling process, including the geometry of the system, physical-chemical parameters of the studied process and other factors. In [14], basing on the numerical solution of a stationary system of equations we obtained an agreement of numerical results with experimental data from [15] on stabilization of a two-component gas-condensate mixture passing through a porous medium. In [16], the experimental results of [17] on the dynamics of extraction of heavy components of a multicomponent hydrocarbon mixture in the depletion mode are numerically reproduced in the assumption of a homogeneous spatial distribution of the pressure and the density of hydrocarbons.
This work is a continuation of our studies on the numerical analysis of measurements at VNIIGAZ, Ukhta [17]. As previously, our consideration is based on the mathematical model described by the Darcy law and the law of continuity of flows. The approach developed in [16] is generalized to the case of taking into account the coordinate dependence of physical characteristics of the system modeled. Based on a numerical solution of the formulated initialboundary value problem for a system of nonlinear partial differential equations, an adequate agreement is obtained between the numerical results and the experimental data from [17] on the yield of both heavy and light hydrocarbons in the laboratory model of the gas-bearing formation at a temperature of 25 ∘ C. In the paper, the mathematical formulation of the problem is given, the computational scheme implemented as a package of C++ and Maple programs is described, and the results of computer simulations are presented.

General statement of the problem
In general, the dynamic process of the passing the -component gascondensate mixture of hydrocarbons through a porous medium is described by a system of equations Here, the first equation (1) corresponds to the Darcy law, and the the second one (2) -to the flows continuity law, -interfacial transition rates. The process involves components: = 1 … , located at pressure in the -phase, where = corresponds to the liquid phase, and = to the gaseous phase.
For components in the phase, the is the linear flow velocity, and are, respectively, the parameters of permeability and viscosity of the -phase, is the molar density of the -th component in -phase (defined as the local volume average).
It is assumed that the gas and liquid in the pores occupy a separate volumes divided by the interfacial surface, on which there are surface tension forces. Concentrations of components in gas and liquids are different, depending on time and coordinate. Inside the gas and liquid, the pressure is described by the Peng-Robins equation: Here is the gas constant, is the temperature, is the molar volume taking into account for a porosity and the gas saturation The constants and are defined as follows [18], [19]: Here is the acentric factor of the -th component, and are the critical pressure and the temperature. The gas saturation is defined as the maximal root of the equation The molar fractions of the components in the liquid and gas phases in the formulas for and are calculated as Interfacial transition rates satisfy the condition and are determined from the following relation: where , are the Gibbs chemical potentials of the -th component, is the interfacial coefficient transition, depending on many factors, including the structure of the rock.
Potential formulae are derived in [19] from equation (3): Note that in case capillar = 0, taking into account for (10), the gas saturation can be calculated from the equation Here are the roots of the equation (following from (3)) such that In the experiment [17], where a gas-condensate mixture passes through a long pipe filled with a porous substance, it is natural to consider the spatially one-dimensional case with a single coordinate along the pipe. In this case, the differentiation operator will have the form ∇ = . Note that the problem has a simple generalization to the case of a thin gas-bearing formation with circular symmetry. In this case, when passing to the polar coordinate system, the differentiation operator takes the form ∇ = 1 ( ).

Laboratory experiment
In the experiment [17], the laboratory model (LM) of the gas-bearing formation is a thermostatic cylinder with a narrow pipe (core holder) of 3 16 H 34 in molar concentrations, respectively, of 87.01%; 7.00%; 1.11%; 0.70%; 0.86%; 1.19%; 0.94%; 1.02%; 0.17%, and both ends of the pipe are closed. After the achievement of the phase equilibrium in LM, one end of the tube is opened that provides a gradual release of the substance through this end while maintaining a constant rate of its consumption by means of a regulating the pressure reduction not more than 0.2 MPa/h. The yield of heavy ( 5+ ) and light ( 2−4 ) hydrocarbons in both liquid and gas phases depending on the pressure were measured. Aggregate C 5+ includes components nC 6  = 25 ∘ C and = 60 ∘ C. Here, we present numerical results only in case = 25 ∘ C.

Computational scheme and implementation
Taking into account the conditions the laboratory experiment, we consider here a spatially one-dimensional case. Since the time and the coordinate do not enter explicitly into the system (1), (2), we pass to arbitrary units, assuming, in particular, ∈ [0, 1]. The end of = 1 is always closed; the mixture is extracted through the point = 0.
To numerically solve the system of equations (1), (2), a discrete mesh in coordinate is introduced with the step ℎ = 1/ , with main nodes = ℎ /2 + ( − 1) ⋅ ℎ and intermediate nodes ±1/2 = ± ℎ /2, where = 1, … , , is a number of the discrete mesh nodes. The difference equations at the nodes take the form: are, respectively, the density of -components at the boundaries between cells and the linear velocity of the -components across the boundaries. Note that the scheme (14) allows a simple generalization on the spatially two-dimensional case for a plane gas-bearing layer with circular symmetry.
The boundary condition for at the closed end = 1 has the form | =1 = 0 or, in terms of the difference scheme, ( +1/2) = 0. As for the end = 0, it is also closed until = 0, i.e. | =0, ⩽0 = 0 or (1/2) = 0. For > 0, taking into account the phase equilibrium achieved under the experimental conditions and the constant yield of the hydrocarbon substance, the boundary condition is formulated as ( )| =0 = | =0 , where | =0 is the given constant. The initial conditions for the functions are chosen so that the pressure at the initial time is equal to 35 MPa.
To the numerical solution of the Cauchy problem for the system (14), the fourth-order Runge-Kutta method is implemented as a combined Maple/C++ code, where The Maple part is responsible for the input data, saving and visualization of the results, while the numerical solution of the system (14) is written in C++, which, as practice has shown, provides a faster calculation compared to the previously developed pure Maple implementation.
As in the previous work [16], the numerical simulation continued until the pressure as a result of MP depletion was reduced to 1 atm. The calculations were performed with the coordinate step ℎ = 0.1 and the time step of 0.0001. The values of the parameters determining the physicochemical properties of components in the hydrocarbon mixture are taken from [20]. In calculations we assumed = . The interphase transition coefficient was chosen in the form = 0 , where is the molar weight of -th component, 0 and are constants to be varied in order to adequately reproduce the experimental data. The ratios of viscosity to permeability parameters in the and phases in equation (1) ( = and = ) were also adjusted to the experimental data.
The best agreement with the set of experimental data from [17] for the case = 25 ∘ was obtained at 0 = 0.001713, = 1.370349, = 0.004631 × 10 −7 , = 2.073255 ⋅ 10 −7 . Figure 1(a) and Figure 1(b) show the pressure dynamics and gas saturation at the points = 0.05, = 0.45, and = 0.85 which are, respectively, near the exit from the LM (solid curves), in the central part of the LM (dashed), and near the closed end (dotted). As in the laboratory experiment, the calculated decreases during the LM depletion. It can be seen that in different parts of the LM is different: closer to the exit region of LM, is less and faster tends to zero. The saturation values in different regions of the LM also differ. When drops in process of depletion, the saturation slightly decreases, and at the final of the depletion process, at low , it noticeably increases in the region of the LM exit.

Results of computer simulations
The phase equilibrium reached according to the experimental conditions by the time = 0, begins to be disturbed when the transport factor is turned on, especially at the of the depletion process. This is confirmed by the increase in the chemical potential difference (J/mol) shown in Figure 1(c) for the components nC 9 H 20 (solid), C 3 H 8 (dashed) and CH 4 (dotted) at the point = 0.05 corresponding to the exit region of LM. Figure 2 shows the evolution of the molar densities of the gas components nC 16 H 26 (dash-dotted line), nC 10 H 22 (solid), nC 9 H 20 (dashed), nC 6 H 14 (dotted line) at the points = 0.05, = 0.35, and = 0.85 which correspond to the exit region from LM (Figure 2(a)), the internal region of LM (Figure 2(b)) and the closed end region (Figure 2(c)). During the depletion process, the density of each hydrocarbon component decreases. The "bursts" of increasing concentrations in the graphs can be explained by the transition of hydrocarbon components from the condensed phase to the gaseous state (boiling) due to the pressure drop. It can be seen that the boiling begins with the heaviest hydrocarbons and is observed first in the exit region of LM (Figure 2(a)), where the pressure is the lowest. This effect requires the further studies.  Figure 3 shows the results of calculating the yields of heavy ( 5+ ) and light ( 2−4 ) hydrocarbons depending on the pressure in comparison with experimental data from [17] at temperature 25 ∘ C. It is seen that the developed approach provides adequate reproduction of the laboratory experiments. In particular, computer simulation confirms some increase in the recoverability of hydrocarbons at low pressure, which occurs, as can be seen from Figure 2, due to the transition of part of the condensate into a gaseous state.  Figure 3. (a) Dependence of the yield (g/m3) of heavy hydrocarbons C 5+ on the pressure at the outlet of the MP in the depletion mode at a temperature of 25 C in comparison with the experimental data from [17]. (b) Dependence of the molar fraction (in %) of heavy hydrocarbons C 5+ on the pressure at the outlet of the MP in the depletion mode at a temperature of 25 C in comparison with the experimental data from [17]. (c) Dependence of the molar fraction (in %) of light hydrocarbons C 2−4 on the pressure at the outlet of the MP in the depletion mode at a temperature of 25 ∘ C in comparison with the experimental data from [17]

Conclusions
The mathematical formulation of the problem and the Maple/C++ implementation has been developed for numerical simulation of the process of extraction of multi-component hydrocarbon gas-condensate composition in the depletion mode. Accounting for the hydrocarbons distribution along the length of the LM allowed one to well reproduce the experimental data on the yield of heavy and light hydrocarbons as a function of pressure obtained at laboratory model of the reservoir at the temperature of 25 ∘ C (VNIIGAZ, Ukhta) [17]. Currently, similar calculations are carried out to analyze experimental data obtained at the temperature of 60 ∘ C. In conclusion, we note that the calculations presented were carried out with a relatively small number of nodes of a discrete mesh at . Increasing the calculation accuracy, as well as modeling two-dimensional systems based on the proposed difference scheme, requires a significant increase in computer time, which makes the parallel implementation of the presented computational scheme using parallel programming technique relevant.