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

**Authors:**Volokhova A.V.^{1}, Zemlyanaya E.V.^{1}^{,2}, Kachalov V.V.^{3}, Rikhvitsky V.S.^{1}, Sokotushchenko V.N.^{2}^{,3}**Affiliations:**- Joint Institute for Nuclear Research
- State University “Dubna”
- Joint Institute for High Temperatures of Russian Academy of Sciences

**Issue:**Vol 27, No 3 (2019)**Pages:**205-216**Section:**Mathematical Modeling**URL:**http://journals.rudn.ru/miph/article/view/22700**DOI:**https://doi.org/10.22363/2658-4670-2019-27-3-205-216- Cite item

# Abstract

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.

# Full Text

1. 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. 2. 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 = - ∇, (1) + ∇( ) = . (2) 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: = () = - - ( + ) + ( - . (3) ) Here is the gas constant, is the temperature, is the molar volume taking into account for a porosity and the gas saturation = , = ∑ . (4) =1 The constants and are defined as follows [18], [19]: = 2 2 , = , = ∑ (1 - )( )1/2 , = ∑ , ,=1 =1 2 1/2 2 = 0.4572355 [(0.37464 + 1.5422 - 0.02699 )(1 - ( ) ) + 1] , = 0.077796074 / , = / , = / . 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 = ( ) = ( (1 - ) ) + capillar, (5) + = 1, 0 ⩽ ⩽ 1. The molar fractions of the components in the liquid and gas phases in the formulas for and are calculated as = . (6) Interfacial transition rates satisfy the condition + = 0 (7) and are determined from the following relation: = (̄ - ), where ̄ = , ̄ = , (8) 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): = ln( ) + ( - 1) - ln( - )+ 1/2 ⎛⎜ ∑(1 - )( ) ⎞ + √ ⎜⎜ =1 2 ⎜ - 2 ⎟ ⎟⎟ × ⎟ ⎝ √⎠ × ln ( + (1 - √2) ) , (9) where + (1 + 2) = . (10) Note that in case capillar = 0, taking into account for (10), the gas saturation can be calculated from the equation (11) (1 - ) = . Here are the roots of the equation (following from (3)) () = 3 + ( - 1)2 + ( - 32 - 2 )+ + (- + 2 + ) = 0, (12) such that = max{ ∶ () = 0}, = min{ ∶ () = 0}. (13) 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 gas-b has a simple generalization to the case of a thin earing formation with circular symmetry. In this case, when passing to the polar coordinate system, the differentiation operator takes the form ∇ = 1 ( ). 3. 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 cm in diameter and 93.27 cm long inside. A pipe filled with a terrigenous filler with porosity of = 0.1377, has been saturating under the pressure 8 of 35 MPa with a nine-component hydrocarbon mixture of CH4, C3H , nC H 4 10 , nC6 H14 7 , nC H 16 , nC9 H20 , nC10 H , nC 22 12 H26 16 , nC 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 C5+ includes components nC6H14, nC7H16, nC9H20, 10 nC H 22 , nC12 H26 , nC16 H34 of the hydrocarbon mixture, and the aggregate C 2-4 consists of components CH4 3 , C H 8 , nC4 H . The experiment was carried 10 out at two temperature values: = 25∘C and = 60∘C. Here, we present numerical results only in case = 25∘C. 4. 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: (-1) () (-1/2) = - - , (14a) -1 - () (+1/2) (+1/2) (-1/2) (-1/2) - + +1/2 - -1/2 = , () (14b) (-1/2) 1 () (-1) 1 (-1) () (+1) = 2 ( + ) - 8 ( - 2 + ). (14c) () () Here = (, ), = (, ), is the number of a coordinate node. The difference equations (14) can be interpreted as the modeling the original () () (one-dimensional) object as a set of cells ∈ [-1/2, +1/2], where , and () are, respectively, the average density of -component, the pressure and the saturation of the -phase inside each cell. Values (±1/2) and (±1/2) 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,⩽0 = 0, it is also closed until = 0, i.e. | = 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. 5. Results of computer simulations 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. 17 0.8 1500 13 P SG×10 6 0.79 1000 (c) 9 (a) 5 0.78 (b) L G φi -φi 500 0 4 5 6 7 t 0.77 4 5 6 7 t -500 4 5 6 7 t Figure 1. Evolution of the pressure = = in MPa (a) and of the gas saturation 1. in the region of exit from the LM (solid curve), in the central region of the LM (dashed curve) and in the closed end region (dotted curve). 2. Evolution of the chemical potential difference (J/mol) for the components nC9H20 3 (solid), C H 8 (dashed) and CH4 (dotted) at the point = 0.05 in the exit region of the LM 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. 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 nC9H20 (solid), C3H8 (dashed) and CH4 (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 16 nC H 26 10 (dash-dotted line), nC H 22 9 (solid), nC H 20 6 (dashed), nC 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. 101 101 101 ρG 100 ρG 100 ρG 100 10-1 (a) 10-1 (b) 10-1 (c) 10-2 4 4.5 5 5.5 6 6.5 7 t 10-2 4 4.5 5 5.5 6 6.5 7 t 10-2 4 4.5 5 5.5 6 6.5 7 t Figure 2. Evolution of molar densities of gas components nC16H26 (dash-dotted line), 9 nC10H22 (solid), nC H20 (dashed), nC6 H14 (dotted line) in the exit region of the LM (a), in the inner region of the LM (b) and in the closed end region (c) 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. 50 (a) 40 2 (b) 20 (c) 15 C , g/m3 C , mol % C , mol % 30 1 5+ 5+ 2-4 10 20 5 10 0 0 0 5 10 15 P, MPa 0 5 10 15 P, MPa 0 0 5 10 15 P, MPa Figure 3. (a) Dependence of the yield (g/m3) of heavy hydrocarbons C5+ on the pressure at the outlet of the MP in the depletion mode at a temperature of 25C in comparison with the experimental data from [17]. 3. Dependence of the molar fraction (in %) of heavy hydrocarbons C5+ on the pressure at the outlet of the MP in the depletion mode at a temperature of 25C in comparison with the experimental data from [17]. 4. 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] 6. 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.

# About the authors

### Alina V. Volokhova

Joint Institute for Nuclear Research
**Author for correspondence.**

Email: bskr@yandex.ru

6 Joliot-Curie St., Dubna, Moscow Region 141980, Russian Federation

Junior Researcher of Joint Inststute of Nuclear Research

### Elena V. Zemlyanaya

Joint Institute for Nuclear Research; State University “Dubna”
Email: elena@jinr.ru

6 Joliot-Curie St., Dubna, Moscow Region 141980, Russian Federation; 19 Universitetskaya St., Dubna, Moscow Region 141980, Russian Federation

Doctor of Physical and Mathematical Sciences, Head of sector of Joint Inststute of Nuclear Research, Professor of State University “Dubna”

### Vladimir V. Kachalov

Joint Institute for High Temperatures of Russian Academy of Sciences
Email: ongk@mail.ru

13 Bd.2 Izhorskaya St., Moscow 125412, Russian Federation

Candidate of Technical Sciences, Senior Researcher of Joint Institute for High Temperatures of Russian Academy of Sciences

### Victor S. Rikhvitsky

Joint Institute for Nuclear Research
Email: rqvtsk@mail.ru

6 Joliot-Curie St., Dubna, Moscow Region 141980, Russian Federation

Leading Programmer of Joint Inststute of Nuclear Research

### Vadim N. Sokotushchenko

State University “Dubna”; Joint Institute for High Temperatures of Russian Academy of Sciences
Email: sokotushenko@mail.ru

19 Universitetskaya St., Dubna, Moscow Region 141980, Russian Federation; 13 Bd.2 Izhorskaya St., Moscow 125412, Russian Federation

Candidate of Technical Sciences, Associate Professor of State University of Dubna, Leading Engineer of Joint Institute for High Temperatures of Russian Academy of Sciences

# References

- Z. P. Sklyarova, F. S. Sokolov, and V. S. Tkach, “Characteristics of the raw materials base of condensate [Kharakteristika syr’yevoy bazy kondensata],” Groups Gazprom. Conduct a gas science, vol. 18, no. 2, pp. 4-14, 2014, in Russian.
- R. I. Vyakhirev, A. I. Gritsenko, and R. M. Ter-Sarkisov, Development and operation of gas fields [Razrabotka i ekspluatatsiya gazovykh mestorozhdeniy]. Moscow: Nedra, 2002, 880 pp., in Russian.
- V. M. Zaichenko, I. L. Maikov, and V. M. Torchinskii, “Features of hydrocarbon mixtures filtration in a porous medium,” High Temperature, vol. 51, no. 6, pp. 776-784, 2013. doi: 10.1134/S0018151X13050222.
- V. M. Zaichenko, I. L. Maikov, V. M. Torchinskii, and E. E. Shpil’rain, “Simulation of processes of filtration of hydrocarbons in a gas-condensate stratum,” High Temperature, vol. 47, pp. 669-674, 5 2009. DOI: 10.1134/ S0018151X09050083.
- O. A. Lobanova and I. M. Indrupsky, “Modeling the mutual influence of hydroand thermodynamics processes in the filtration of hydrocarbon systems [Modelirovaniye vzaimovliyaniya gidroi termodinamicheskikh protsessov pri fil’tratsii uglevodorodnykh sistem],” Automation, telemechanization and communication in the oil industry, pp. 19-23, 10 2010, in Russian.
- O. A. Lobanova and I. M. Indrupsky, “Nonequilibrium of the phase behavior of hydrocarbon systems: modeling and scale effect [Neravnovesnost’ fazovogo povedeniya uglevodorodnykh sistem: modelirovaniye i masshtabnyy effekt],” in Abstracts of the IX All-Russian scientific and technical conference “Actual problems of development of the Russian oil and gas complex”, in Russian, 2012, pp. 96-97.
- A. V. Volokhova, E. V. Zemlyanaya, V. V. Kachalov, and V. N. Sokotushchenko, “Overview of the component enhancement methods in development of the gas condensate fields [Obzor metodov povysheniya komponentootdachi pri razrabotkakh gazokondensatnykh mestorozhdeniy],” The science. Innovations. Technologies, vol. 3, pp. 19-48, 2019, in Russian.
- K. Aziz and A. Settary, Petroleum reservoir simulation. London: Applied Science Publishers Ltd, 1979, 476 pp.
- I. N. Ponomareva and V. A. Mordvinov, Underground hydromechanics [Podzemnaya gidromekhanika]. Perm: Perm state technical university, 2009, 103 pp., in Russian.
- V. S. Mitlin, Underground hydromechanics of complex hydrocarbon mixtures [Podzemnaya gidromekhanika slozhnykh uglevodorodnykh smesey]. Moscow: VINITI, 1991, vol. 4, 154-222, in Russian.
- G. T. Bulgakova, T. A. Faizullin, and A. V. Zhiber, “Nonequilibrium two-phase filtration [Neravnovesnaya dvukhfaznaya fil’tratsiya],” Matematicheskoye modelirovaniye, vol. 18, no. 10, pp. 19-38, 2006, in Russian.
- D. O. Dil and A. M. Bubenchikov, “Diphasic filtration in a pipe filled with porous material [Dvukhfaznaya fil’tratsiya v trube, zapolnennoy poristym materialom],” Bulletin of Tomsk State University, Mathematics and Mechanics, vol. 5 (25), pp. 45-51, 2013, in Russian.
- A. L. Kovalev and Y. V. Sheberstov, “Numerical simulation of non-equilibrium local filtration in gas-condensate beds [Chislennoye modelirovaniye lokal’no-neravnovesnoy fil’tratsii v gazokondensatnykh plastakh],” News of gas science, no. 5 (37), pp. 164-171, 2018, in Russian.
- A. V. Volokhova, E. V. Zemlyanaya, V. V. Kachalov, V. N. Sokotushchenko, and V. S. Rikhvitskiy, “Numerical investigation of the gas-condensate mixture flow in a porous medium [Chislennoye issledovaniye fil’tratsii gazokondensatnoy smesi v poristoy srede],” Computer Research and Modeling, vol. 10, no. 2, pp. 209-219, 2018, in Russian.
- H. X. Vo, “Composition Variation During Flow of Gas-Condensate Wells,” Department of energy resources engineering of Stanford University, Tech. Rep., 2010.
- A. V. Volokhova, E. V. Zemlyanaya, V. V. Kachalov, V. S. Rikhvitsky, and V. N. Sokotustchenko, “Numerical modeling of dynamics of extraction of multicomponent gas-condensate hydrocarbon mixture in the mode of depletion of the filtration model of the formation [Chislennoye modelirovaniye dinamiki izvlecheniya mnogokomponentnoy gazokondensatnoy uglevodorodnoy smesi v rezhime istoshcheniya fil’tratsionnoy modeli plasta],” Geoinformatika, vol. 3, pp. 27-33, 2019, in Russian.
- A. N. Volkov, V. L. Lapshin, and A. V. Polyakov, “Simulation of gaz condensate system phase beghavior in porous media [Modelirovaniye fazovogo povedeniya gazokondensatnoy sistemy v poristoy srede],” Gaz Industry Magazine, no. 10, pp. 26-31, 2016, in Russian.
- V. G. Lysov and Y. G. Rykov, “On calculation of phase equilibrium in multicomponent filtration problems [O vychislenii fazovogo ravnovesiya v zadachakh mnogokomponentnoy fil’tratsii],” IPM Preprint im. M. V. Keldysh Russian Academy of Sciences, no. 94, 2014, in Russian.
- A. I. Brusilovsky, Phase transformations in the development of oil and gas fieldse [Fazovyye prevrashcheniya pri razrabotke mestorozhdeniy nefti i gaza]. Grail, 2002, 575 pp., in Russian.
- L. B. Director, V. V. Kachalov, I. L. Maikov, and S. N. Skovorodko, One-dimensional nonstationary model of two-phase filtration of a gascondensate mixture [Odnomernaya nestatsionarnaya model’ dvukhfaznoy fil’tratsii gazokondensatnoy smesi], 2-xs441. 2000, 45 pp., in Russian.