Numerical evaluation of phase behavior properties for gas condensate under non-equilibrium conditions

Fuel 226 (2018) 675–685

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

Numerical evaluation of phase behavior properties for gas condensate under non-equilibrium conditions

T

⁎

Hanmin Tua, Ping Guoa, , Na Jiab, Zhouhua Wanga a b

State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China Program of Petroleum Systems Engineering, Faculty of Engineering and Applied Science, University of Regina, Regina, Saskatchewan S4S 0A2, Canada

A R T I C LE I N FO

A B S T R A C T

Keywords: Non-equilibrium effect Interphase mass transfer Retrograde condensation saturation Gas condensate

It is generally known that for a gas reservoir if its pressure changes or flow rate is smaller than the phase transition or mass transfer rate, non-equilibrium effect may arise. Commonly, this situation is ignored when the approaches based on the local equilibrium assumption were applied. Previously there are only a few works focused on this subject and now it attracts more and more attentions now. In this work, a mathematical model was formulated to simulate the gas-related properties of Constant Composition Expansion (CCE) tests under nonequilibrium conditions. In this method, the condensed liquid or condensate is assumed to form a continuous film when brought into contact with gas condensate that was depressurized incrementally at a constant temperature. The whole process is controlled by interphase mass transfer and is completed when abandonment pressure or thermodynamic equilibrium is reached. The elapsed time to reach the final state is determined by the pressure dropping rate and the pressure difference between the initial and final states. The results obtained using the current method match satisfactory with the experimental data. The results demonstrate that the larger the pressure drop rate, the lower the retrograde condensate saturation, and the stronger degree of non-equilibrium effect. It proved that a higher pressure dropping rate can not only alleviate the retrograde condensate phenomenon but also improve the recovery of gas condensate.

1. Introduction Gas condensate reservoirs, many of them with industrial exploitation value play a significant role in natural gas industries. For retrograde condensate, the stratum pressure of the reservoir usually decreases below the dew point pressure during the production process. The condensed oil will not flow and therefore a region of high saturation will build up near the wellbore unless its saturation exceeds a threshold value. It makes the phase behavior of such system (gas and condensate) quite complex and difficult to be simulated, and the more crucial issue is that it will cause the condensate blockage problem and reduce the recovery of both gas and oil phases [1]. In contrast to the conventional natural gas reservoirs, numerous studies were focused on condensate oil which may cause lots of formation problems. In these studies, some of them derived a simulation model that considers the radius of condensate blockage as a function of time [2]; some findings demonstrated that condensate blockage would reduce well deliverability [3,4]; a few studies introduced the first gas rate equation to describe the condensate blockage [5]; some other studies focused on Equation of State (EoS) compositional simulation

⁎

[6–8]. All of the above studies concentrated on the macroscopic well deliverability and they in fact alleviated the effect of condensate blockage on gas production to some extent. However, how to attenuate the condensate blockage and explore and develop gas condensate reservoirs economically and effectively are still critical issues. In this work, our attention is focused on the better understanding of the mechanisms of condensate of gas condensate reservoir. A.X [9] pointed out that when the pressure of the condensate gas reservoir decreased below the dew point pressure with a constant rate, it was difficult for both gas phase and condensate liquid phase to reach thermodynamic equilibrium at high speed gas flow rate. This phenomenon was henceforth denoted as a non-equilibrium effect, namely that when the pressure or temperature change rate is large enough to supersede the phase transition rate (mass or heat transfer), non-equilibrium effect will arise. Based on this rationale, it was observed that the amount of condensate liquid gradually decreased with increasing pressure dropping rate. And a correlation was proposed to describe the relationship between the condensate liquid and pressure dropping rate under the non-equilibrium conditions. At present, it is generally believed that the non-equilibrium effect

Corresponding author. E-mail address: [email protected] (P. Guo).

https://doi.org/10.1016/j.fuel.2018.04.048 Received 5 December 2017; Received in revised form 23 February 2018; Accepted 10 April 2018 Available online 24 April 2018 0016-2361/ © 2018 Published by Elsevier Ltd.

Fuel 226 (2018) 675–685

H. Tu et al.

Nomenclature

L V xi yi n x iI

yiI SV0 SL0 Ni Nt (F ) PiV

PiL JiV JiL CtV CtL D B Γ Pd T

Vl Vt D̵ δij γi k· k l Ξ Φ Θ QiI Ki Psi P I (χ ) Δt A Ni j j Npositive (i )

liquid phase gas phase mole fraction of component i in bulk liquid phase mole fraction of component i in bulk gas phase number of components interfacial mole fraction of component i on the liquid phase side Interfacial mole fraction of component i on the gas phase side initial mole amount of gas phase,% initial mole amount of liquid phase,% mole flux of component i, mol/(m2·s) total molar flux, mol/(m2·s) vector of Independent equations mass transfer rate equations of component i in the gas phase mass transfer rate equations of component i in the liquid phase molar diffusion flux of component i in gas phase, m/s molar diffusion flux of component i in liquid phase, m/s total mole density of gas phase, mol/m3 total mole density of liquid phase, mol/m3 matrix of effective diffusion coefficient, m2/s A square matrix matrix of thermodynamic factor dew point pressure, MPa simulated temperature, K

measured volume of condensed liquid, m3 measured total volume of the system, m3 matrix of Max-Well diffusion coefficient, m2/s Kronecker delta function activity coefficients of component i matrix of finite flux mass transfer coefficient, m/s mass transfer coefficient, m/s thickness of the film, m matrix of correction factor matrix of Mass transfer rate factor matrix of thermodynamic correction factor interfacial equilibrium equation of component i equilibrium factor of component i saturation vapor pressure of component i, MPa system pressure, MPa unit matrix vector solution variables Time step, hour Effective area of continuous film, m2 mole amount of component i at time step j, mol positive mole amount of component i at time step j, mol

j Nnegative (i) negative mole amount of component i at time step j, mol

NLj (i) NVj (i)

(J ) (Ef )

mole amount of component i in the liquid phase at time step j, mol mole amount of component i in the gas phase at time step j, mol matrix of Jacobi vector of deviation parameter

dynamic flow parameters, and then a new mass transfer rate named normalized component mass transfer rate was obtained. All of these correlations are used in the compositional flow equations. Actually, as gas condensate flows from formation to the wellbore, due to the drastic pressure declining rates, the gas velocity increases and reaches a maximum near the wellbore. When the gas velocity is extremely low, the system has enough time to reach equilibrium; however as velocity increases, the system has no sufficient time to reach equilibrium thus it shifts to a non-equilibrium state. The situation described above has been studied experimentally in most of the aforementioned works, and the process was unpredictable by using a commercial compositional simulator since they are all based on the local equilibrium assumption. Consequently, the non-equilibrium mechanism for gas condensate reservoirs is still obscure and the required gas related properties are also insufficient at different pressure stages since the whole experimental process is at a non-equilibrium state. It is

has a significant effect on gas condensate within the reservoir and it needs to be taken into account to accurately model this process. At this moment, only a few works were proposed in the literature to model non-equilibrium effect on gas related properties of gas condensate, especially with the assistance of theoretical research. On one hand, experimental studies indicated that the produced gas and oil are in non-equilibrium state for both ‘condensing’ and ‘vaporing’ gas drives [10] and there is a precipitation lag when a constant pressure dropping rate is maintained below the dew point pressure [11–14]. In the process of depressurization, portions of the formed condensate liquid moved along with the gas phase until they are deposited. This nonequilibrium phenomenon leads to a higher content of liquid condensate in the gas phase compared with those at equilibrium state situation [15,16]. Other researches also confirmed that the ultimate condensed liquid recovery had approximately increased from 36.15 wt% to 42.55 wt% due to the inadequate contact between the gas and liquid phases [17]. Results show that the higher the pressure dropping rate, the more serious of the non-equilibrium effect, the more condensate fluid produced with gas, the lower amount of the retrograde condensate remained in the formation, and the greater deviate degree from the equilibrium state [18]. On the other hand, theoretical research [19] stated that for the nonequilibrium state of a two-phase mixture, mass transfer took place at the interface. The above process could be expressed by using an efficiency factor which was an exponential function of the contact time and the mass transfer coefficient. Later on, a new K-factor for a system at non-equilibrium has been proposed [20]. It stated that these equilibrium and non-equilibrium K-factors could be related by a variable so called the degree of non-equilibrium of the system (E). And a semiempirical equation was proposed to obtain the value of E based on the fact that the mass transfer rate not only depends on the equilibrium Kfactor and the contact time but also is inversely proportional to the effective mobility ratio of the two phases. Following that, several studies were focused on the mass transfer rate [21–23]. They introduced a coefficient C which is a function of reservoir characteristics and

l y iV

y

Condensate gas

y iI

Condensate liquid

Interphase

x iI

x

x iL

Fig. 1. The schematic of interphase mass transfer. 676

Fuel 226 (2018) 675–685

H. Tu et al.

Table 1 The required initial iterative value for calculation. Initial parameters

Pressure

Temperature

Compositions of bulk gas/liquid phase

Saturation of liquid/gas phase

Mole amount for component i

The thickness of film

Initial value

Pd

T

yi / xi

SL0 /SV0

N L0 (i)

l

Fig. 2. Algorithm flow chart.

[24]. Herein, the interphase mass transfer phenomena occurred. Currently, mass transfer theory is commonly used in the chemical industry. Numerous researchers calculated the mass and energy transfer rate based on film models [25–37]. Among those models, some of them are concentrated on the condenser, which is a widely used instrument in chemical engineering industry. In a typical condenser, condensation occurs when a gas mixture is brought into contact with a cold surface.

reasonable to establish a mathematical model with the capability to predict the related properties at the non-equilibrium state. From a thermodynamic point of view, when pressure of gas condensate reduced to below its dew point pressure, the system will display complicated physical and chemical responses in a gas-extraction process. The condensed liquid is more likely to form a continuous film. Drop-wise condensate may also be formed under certain circumstances 677

Fuel 226 (2018) 675–685

H. Tu et al.

Viscometer

Separator

Hydrometer

Gasometer Separator of Ground

Intermediate Container

Chromatograph

PVT cell Automatic Pump Fig. 3. Schematic diagram of the experimental apparatus [41]. 0.016

Component

Composition of gas phase (mol %)

Composition of liquid phase (mol %)

C1 C2 C3 C4 C5 C6 C7+

80.95 11.90 1.94 1.08 1.11 1.12 1.87

60.45 13.22 2.97 2.02 2.78 3.99 14.58

Mass transfer flux of C1 (mol/(m2h))

Table 2 The compositions of gas and liquid phase slightly below the dew point pressure.

Table 3 The calculated average gas diffusion coefficients (× 10−8 m2/s).

5MPa/h

0.012 0.008

0.004 0

0

7

C2

C3

C4

C5

C6

C7+

C1 C2 C3 C4 C5 C6 C7+

0.000 1.994 1.970 1.830 1.688 1.532 1.482

1.994 0.000 1.878 1.745 1.609 1.461 1.412

1.970 1.878 0.000 1.409 1.300 1.180 1.141

1.830 1.745 1.409 0.000 1.218 1.105 1.069

1.688 1.609 1.300 1.218 0.000 1.037 1.003

1.532 1.461 1.180 1.105 1.037 0.000 0.947

1.482 1.412 1.141 1.069 1.003 0.947 0.000

C1

C2

C3

C4

C5

C6

C7+

C1 C2 C3 C4 C5 C6 C7+

0.000 2.444 2.389 2.376 2.374 2.372 1.740

2.444 0.000 4.786 4.760 4.752 4.755 3.486

2.389 4.786 0.000 3.697 3.692 3.694 2.708

2.376 4.760 3.697 0.000 3.224 3.227 2.365

2.374 4.752 3.692 3.224 0.000 2.928 2.147

2.372 4.755 3.694 3.227 2.928 0.000 1.920

1.740 3.486 2.708 2.365 2.147 1.920 0.000

28

35

28

35

Pressure (MPa)

0.002 5MPa/h

Table 4 The calculated average liquid diffusion coefficients (× 10−9 m2/s). Component

21

Fig. 4. The mass transfer flux of C1.

Mass transfer flux of C4 (mol/(m2h))

C1

14

-0.004 -0.008

Component

2MPa/h

2MPa/h

0.0015

0.001

0.0005

0 0

-0.0005

7

14

21

Pressure (MPa) Fig. 5. The mass transfer flux of C4.

problem and facilitate theoretical analysis, only film-wise condensate is considered. In this method, mass transfer mechanism on the base of two-film theory was utilized. The model focuses on the saturation of retrograde liquid, the mass transfer flux, the interface and the bulk phase compositions. A serials of material balance equations are solved for the calculation of condensate liquid growth. With a view to facilitate the entire calculation, the calculating method only considered the existed condensed liquid with the assumption of continuous film. The interphase mass transfer is assumed to occur in the close proximity of the interface (as shown in Fig. 1).Let’s consider the two gas and liquid phases L and V as shown in Fig. 1 where typical composition profiles are presented. The bulk phase compositions are yiV (i = 1,2,3⋯n) and x iL (i = 1,2,3⋯n) for gas and liquid

This isobaric process is somewhat similar to the isothermal process of the retrograde condensation phenomenon of gas condensate system during CCE processes due to the fact that they all formed condensate liquid when the single variable P or T changed. Therefore, some new idea appears to establish a mathematical model based on this mass transfer theory to simulate the gas related properties under non-equilibrium effect.

2. Mathematical model A new model is developed using a kinetic mass transfer approach to model the non-equilibrium state of gas condensate. To simplify the 678

Fuel 226 (2018) 675–685

H. Tu et al.

0.1

0.004

Liquid phase

2MPa/h

0.08

0.003

Mole fraction of C4

0.002

0.001

0.06

0.04

2MPa/h 0.02

0 0

7

14

21

28

0

1

Components

Liquid phase

AAD =

1 M

C1

C4

C7+

−0.0006 21.4874

0.0001 20.4372

0.0003 23.4739

Mole fraction of C7+

0.8

AAD RAADb (%) M

∑i = 1 (y5i −y2i ) , where AAD is the absolute average deviation, M is

the number of pressure point, y5i is the mole fraction of component i in gas phase at the pressure dropping rate of 5 MPa/h and y2i is the mole fraction of gas phase at 2 MPa/h. b

RAAD (%) =

1 M

M

∑i = 1

abs (y5i − y2i ) y2i

14

21

28

35

Fig. 8. The mole fraction of C4 in the two interface phases.

Table 5 The deviation between the pressure dropping rate of 5 MPa/h and 2 MPa/h.

a

7

Pressure (MPa)

Fig. 6. The mass transfer flux of C7+.

Parameters

Gas phase

0

35

Pressure (MPa)

a

5MPa/h

Mole fraction (%)

Mass transfer flux of C7+ (mol/(m2h))

5MPa/h

0.6

0.0012

5MPa/h 2MPa/h

0.0008 0.0004

Gas phase

0 0

7

14

21

28

35

Pressure (MPa)

0.4 5MPa/h 0.2

2MPa/h Gas phase

× 100 , where RAAD is the relative average

0 0

deviation.

5

10

15

20

25

30

35

Pressure (MPa) 1

Fig. 9. The mole fraction of C7+ in the two interface phases.

Mole fraction of C1

0.8

Gas phase

6) No mass transfer resistance is on the interface.

0.6

0.4

The number of independent equations is 3n , these equations are: (n−1) mass transfer rate equations for the gas phase; (n−1) mass transfer rate equations for the condensate liquid phase; n equilibrium equations for the interface; two equations for gas and liquid mole fractions. The 3n equations can be used to solve 3n variables, which are n mass transfer molar fluxes Ni , 2n mole fractions on each side of the interface yiI and x iI .The following is the set of independent equations being represented by a vector (F ) ,

5MPa/h 2MPa/h

0.2 Liquid phase 0

0

7

14

21

28

35

(F ) = (P1V ,P2V ,⋯PNV − 1,P1L,P2L,⋯PNL− 1,Q1I ,Q2I ,⋯QNI ,RV ,RL)

Pressure (MPa) Fig. 7. The mole fraction of C1 in the two interface phases.

Where the PiV and PiL represent the mass transfer rate equations of component i in the gas and liquid phase respectively as given below,

phases, respectively. n is the number of components. The interfacial compositions for component i are denoted as x iI on the liquid phase side, and yiI on the gas phase side of the interface.The initial calculation is starting from the dew point pressure. Under this condition, the gradient of molar flux takes place in a direction with respect to the interface. Assuming that the initial saturation of gas condensate is SV0 and the condensed liquid is SL0 (which is deemed to be quite small). The following assumptions are applied, 1) 2) 3) 4) 5)

(1)

V V V ⎧ F1 = P1 = J1 + Nt y1 −N1 = 0 ⎪ = V = V + J2 Nt y2V −N2 = 0 ⎪ F2 P2 ⎪⋮ ⎪ ⎪ Fn − 1 = PnV− 1 = JnV− 1 + Nt ynV− 1−Nn − 1 = 0 ⎨ Fn = P L = J L + Nt x L−N1 = 0 1 1 1 ⎪ ⎪ Fn + 1 = P2L = J2L + Nt x 2L−N2 = 0 ⎪⋮ ⎪ ⎪ F2n − 2 = PnL− 1 = JnL− 1 + Nt x nL− 1−Nn − 1 = 0 ⎩

The calculation is carried out under isothermal condition; The gas-liquid interface is in equilibrium state; Mutual diffusion occurs between gas and liquid phases; Diffusion coefficient varies with pressure; The molar fluxes on each side of the gas-liquid interface must be equivalent;

(2)

Here, the transfer direction from the gas phase to the liquid phase was considered as leading to a positive flux, the opposite situation leading to a negative flux.QiI represents the interfacial equilibrium equations, 679

Fuel 226 (2018) 675–685

H. Tu et al.

0.67

Mole fraction of bulk liquid phase

Mole fraction of bulk gas phase

0.85

0.83

0.81

0.79

5MPa/h

2MPa/h

0.77

0

5

10

15

20

25

30

5MPa/h

2MPa/h

0.65

0.63

0.61

0.59

35

0

5

10

Pressure (MPa)

15

20

25

30

35

25

30

35

25

30

35

Pressure (MPa) Fig. 10. The mole fraction of C1 in the bulk phases. 0.036

5MPa/h

0.018

Mole fraction of bulk liquid phase

Mole fraction of bulk gas phase

0.02 2MPa/h

0.016

0.014

0.012

0.033

0.03 5MPa/h

2MPa/h

0.027

0.01 0

5

10

15

20

25

30

0

35

5

10

15

20

Pressure (MPa)

Pressure (MPa) Fig. 11. The mole fraction of C4 in the bulk phases. 0.149

5MPa/h

Mole fraction of bulk liqquid phase

Mole fraction of bulk gas phase

0.039

2MPa/h

0.03

0.021

0.144

0.139

0.134 5MPa/h

2MPa/h

0.129

0.012

0

5

10

15

20

25

30

35

0

5

10

15

20

Pressure (MPa)

Pressure (MPa) Fig. 12. The mole fraction of C7+ in the bulk phases. I I I ⎧ F2n − 1 = Q1 = K1 x1 −y1 = 0 ⎪ I I I F2n = Q2 = K2 x 2 −y2 = 0 ⎨⋮ ⎪ F = QnI = Kn x nI −ynI = 0 ⎩ 3n − 2

interface, n

F3n − 1 = RV =

∑

yiI −1 = 0

i=1

(4)

(3) n

F3n = RL =

Ki factor can be defined as Ki = γi Psi/ P , where γ is the activity coefficient calculated by Van Laar model with van der Waals mixing rules [38].RV and RL denote the mole fraction summation equations at the

∑

x iL−1 = 0

i=1

The vector of the solution variables defined as follows, 680

(5)

Fuel 226 (2018) 675–685

H. Tu et al.

8

Parameters

Retrograde condensate saturation %

Table 6 The deviation of components’ distribution in bulk gas and liquid phases at different pressure dropping rates of 5 MPa/h and 2 MPa/h. Component C1

C4

C7+

AAD

Gas phase Liquid phase

−0.0039 0.0028

0.0006 −0.0005

0.0013 −0.0006

RAAD (%)

Gas phase Liquid phase

0.4745 0.4498

4.5541 1.4167

6.9290 0.4115

Note: The definition of AAD and RAAD is the same as before.

4

2 Experiment Calculated

0

10

Retrograde condensate saturation %

6

0

7

14

21

28

35

Pressure MPa 8

Fig. 15. Comparison between experimental and calculated retrograde condensate saturation at 5 MPa/h.

6

Table 7 The comparison of components’ molar fraction between non-equilibrium with equilibrium states at two pressure dropping rates of 2 MPa/h and 5 MPa/h.

4 Equilibrium

2

Non-equilibrium(2MPa/h)

16

20

24

28

31.03

28.97

26.00

22.07

18.00

14.00

RAD2d

18.32 31.41 21.96 27.06

66.46 60.35

28.80 34.55

10.38 16.49

4.52 10.57

3.31 8.98

(%) RAD5 (%) RAAD2 (%) RAAD5 (%)

Non-equilibrrium(5MPa/h)

0 12

P (MPa)

32

Pressure MPa

d

Fig. 13. Experimental data of the SL at different pressure dropping rate.

RAD2 (%) =

abs (y2i − yeqi )

× 100 , RAD5 (%) =

yeqi

abs (y5i − yeqi ) yeqi

× 100 where RAD is

the relative absolute deviation.

Retrograde condensate saturation (%)

8

Table 8 Summary of the difference between experimental and simulated retrograde condensate saturation.

6

P (MPa) RAD (%)

4

RAAD (%) SDe

2

28.97

26.00

22.07

18.00

14.00

9.55 16.59 5.34 5.49 0.12 0.07

4.67 5.63

12.03 6.39

1.55 0.08

0.31 0.66

3.93 3.61

Experiment Calculated

0

0

7

e

14

21

28

35

Fig. 14. Comparison between experimental and calculated retrograde condensate saturation at 2 MPa/h.

(χ )T = (N1,N2,⋯Nn,x1I ,x 2I ,⋯x nI ,y1I ,y2I ,⋯ynI )

(6)

(7)

j j j j j Here, if Ni j > 0,Npositive (i ) = Ni , else Ni < 0,Nnegative (i) = Ni . Based on the material balance equation, for the liquid phase, j j NLj +(i)1 = NLj (i) + Npositive (i ) + Nnegatuve (i)

(8)

And for the gas phase,

=

j j NVj (i) −Npositive (i ) −Nnegatuve (i)

2 ∑iM = 1 (xi − x ) M−1

,where SD is the standard deviation. {x1,x2,…,xM } are the

Table 9 The compositions of gas condensate examples.

Taking Δt as the time step in the process of reducing pressure, A as the effective area of continuous film, Ni j as the mole amount of component i at time step j.

Ni j = Ni ·Δt·A

SD =

observed non-equilibrium experimental values, x is the mean value of experimental and calculated values.

Pressure (MPa)

NVj +(i1)

2 MPa/h 5 MPa/h 2 MPa/h 5 MPa/h 2 MPa/h 5 MPa/h

31.03

(9)

Hence, the saturation of liquid phase is defined as, 681

Samples

Sample 1

Sample 2

Component

Gas phase (mol %)

Liquid phase (mol %)

Gas phase (mol %)

Liquid phase (mol %)

Gas phase (mol %)

Liquid phase (mol %)

C1 C2 C3 C4 C5 C6 C7+ Content of condensate oil T and Pd

81.76 8.17 4.26 2.72 1.12 0.70 1.27 80 g/m3

50.36 7.56 6.37 5.52 3.09 2.69 24.41

85.64 7.65 2.52 1.07 0.48 0.96 2.18 150 g/m3

61.48 7.86 3.47 2.18 1.09 1.75 22.16

85.03 8.25 2.63 1.21 0.51 0.40 1.97 235 g/m3

63.34 8.77 3.49 2.00 1.06 1.01 20.33

88.3 °C/23.05 MPa

Sample 3

78 °C/30.29 MPa

80.6 °C/35.30 MPa

Fuel 226 (2018) 675–685

H. Tu et al.

Retrograde condensate saturation(%)

2.4

SLj + 1 (%) = 1.8

∑ NLj +(i)1 ∑ NLj +(i)1

+ ∑ NVj +(i1)

× 100 (10)

Newton-Raphson method is recommended by Taylor and Krishna [24] to solve all of the model independent equations represented by a vector (F). The detailed model can be seen in the supplementary Appendix A.

1.2

Experiment(Equilibrium) Calculated in this work(1MPa/h) Calculated in this work(2MPa/h) Calculated in this work(5MPa/h)

0.6

2.1. Calculation procedure At the specific pressure and temperature, the independent equations can be calculated, and an independent equation deviation parameter (Ef ) can be defined as,

0

0

5

10

15

20

25

Pressure (MPa/h)

(Ef ) = |(χ j + 1 )−(χ j )|

Fig. 16. The retrograde condensation curve within the content of 80 g/cm3.

(11)

(χ )T

= The solution vector of variables can be obtained if the discrepancy function gives a minimum value (⩽ 10−5) of Ef for the above equation at each pressure.All of the parameters displayed in the above equations must be specified and calculated. The required initial input data for the calculations are listed in Table 1. The value of the initial molar amount for component i is given by N L0 (i) = x i0 SL0 . For the purpose of accelerating convergence, at each round of calculation, the initial iterative values for mole fractions of gas and liquid at the interface are made to equal to the mole fractions of the bulk gas and liquid phases. The initial iterative molar fluxes are estimated by definition. Here, we define that N > 0 represents a positive flux that components are transferred from the gas phase to liquid phase, while N < 0 denotes the opposite case. This set of initial iterative values is simple and easy to reach convergence at any pressure.To calculate the diffusion coefficient, various correction factors such as [B], [Γ], [Ξ], etc. are estimated. Ki factors in Eq. (3) can be calculated through the relationship between the activity coefficient and the saturated vapor pressure. As the pressure is reduced at a constant rate, and the time step can be determined by pressure dropping rate. The value of the effective area was decided to be the bottom area of the PVT cell. The whole calculation must be completed when the abandonment pressure (2 MPa) has been reached or the difference of calculated two adjacent mass transfer fluxes is quite small. The final calculated interface compositions and the amount of liquid phase should be consistent with those estimated through a single-stage flash calculation by using an EoS. An algorithm for calculating the gas-related properties under non-equilibrium effect is given in Fig. 2.

Retrograde condensate saturation(%)

6 5 4

3 2

Experiment(Equilibrium) Calculated in this work(1MPa/h)

1

Calculated in this work(2MPa/h) Calculated in this work(5MPa/h)

0

0

7

14 21 Pressure (MPa)

28

35

Fig. 17. The retrograde condensation curve within the content of 150 g/cm3. 10

Retrograde condensate saturation(%)

(N1,N2,⋯Nn,x1I ,x 2I ,⋯x nI ,y1I ,y2I ,⋯ynI )

8

6

4 Experiment(Equilibrium)

Calculated in this work(1MPa/h)

2

3. Experimental method

Calculated in this work(2MPa/h) Calculated in this work(5MPa/h)

0

0

10

20

30

The experimental procedure is similar to a conventional CCE test [39,40], and the experimental setup is presented in Fig. 3 in which a mercury-free DBR PVT cell is one of the most important experimental apparatus. At the beginning, the required volume of the gas condensate at a pressure higher than its dew point pressure was charged into the PVT cell at the preset temperature of 351.95 K. The pressure was then decreased gradually in quite small steps till the first condensate liquid dropped out at the dew point pressure. This was the initial state of the non-equilibrium PVT measurement and the liquid composition at this point was treated as the initial composition for subsequent simulation. Pressure decreases further at a constant pressure dropping rate and two phases were formed. The only difference between the non-equilibrium PVT measurement and the traditional CCE test is that the system pressure drops at a constant rate. In this work, two different pressure dropping rates of 2 and 5 MPa/h were applied. At each pressure, the condensate liquid volume Vl was measured and the corresponding retrograde condensate saturation SL was calculated for the CCE as,

40

Pressure (MPa)

Fig. 18. The retrograde condensation curve within the content of 235 g/cm3. Table 10 The deviation between the equilibrium experimental data and the calculated retrograde condensate saturation. Depletion pressure rate

1 MPa/h 2 MPa/h 5 MPa/h

Samples Sample 1 RAD (%)

Sample 2

Sample 3

6.05 9.34 18.10

9.69 13.11 21.90

8.65 10.67 19.09

SL (%) = 682

Vl × 100 Vt

(12)

Fuel 226 (2018) 675–685

H. Tu et al.

an equilibrium state, for the same type of fluid sample, the mole fraction curves in the interface gas and liquid phases (shown in Figs. 7–9) coincide no matter what the pressure dropping rate is. As can be seen, the shape of the curves for each component changes with the component characteristics. For light components, the mole fraction curves gradually decrease with decreasing pressure. For heavy components, the mole fraction curves show an opposite situation and they constantly rise as the pressure decrease. For intermediate components, the change of the curves is between the above two situations. This particular transition is evident for the liquid components shown in Fig. 8, the curve first increases to maximum and then decreases as the pressure decreases. The situation of liquid phase is mainly predominated by heavy hydrocarbons while the majority gas phase is governed by light hydrocarbons. In order to evaluate the mole fraction of both bulk phases, material balance equations were used (as seen in Appendix A). From the previous calculation, it was known that light components tend to be condensed and then re-evaporated while heavy components continue to precipitate with decreasing pressure. Figs. 10–12 confirm that the bulk oil phase is becoming increasingly heavy and the bulk gas phase is becoming continuously light. And these growth and reduction rates arrive a maximum at the maximum retrograde condensate pressure. At the lower pressure, the curves of lightest components tend to level off, while the curves of heavy components do not. At different pressure dropping rate, the component mole fraction within bulk phases has a completely opposite characteristics: light components in the bulk liquid phase is larger at a higher pressure dropping rate compared to a quite opposite situation in the bulk gas phase (see Table 6, the data of AAD). While the situation of heavy components is exactly reversed. This illustrates the different characteristic of components since light components are easier to evaporate while heavy components are more readily to condensate. Fig. 13 presents the comparison of the experimental data of retrograde condensate saturation at two pressure dropping rates with the data at the equilibrium state. It is obvious that all curves show the same trend and the difference between two pressure dropping rates and equilibrium state is gradually decreased with decreasing pressure (it is also can be seen in Figs. 14 and 15). Based on those phenomena, we can speculate that at different developing stage of gas condensate reservoir, the recovery of condensed liquid at non-equilibrium state is much larger than that at equilibrium state. And the greater the degree of nonequilibrium, the higher the recovery rate is. At the same time, it also can be seen that the growth rate of condensate liquid increases rapidly at high pressure stage, while the growth rate slows down or almost equals to that at equilibrium state when the system pressure is less than the maximum retrograde condensate pressure. This phenomenon can be owned to the driving force of the concentration gradient. At high pressure stage, there was a large amount of condensate oil in gas condensate. As pressure decreased, the system mainly focused on condensation and more condensed liquid was carried away by gas condensate stream. This phenomenon is particularly significant in nonequilibrium state. Whereas when the pressure is below the maximum retrograde condensate pressure, the content of condensate oil is less than that at high pressure stage (which means the components of condensate gas are constantly changing), the system was dominated by evaporation. Consequently, at later development period of gas condensate reservoir, the non-equilibrium effect is too weak to be ignored (Table 7). The retrograde condensation saturation was obtained through the Eq. (10). Figs. 14 and 15 show a comparison between the experimental values and the simulated results. The calculated results match satisfactory with the measured retrograde condensation saturation. The average errors are 5.34% and 5.49% for 2 MPa/h and 5 MPa/h, respectively (see Table 8). It confirms that the proposed method can accurately estimate retrograde condensation saturation at different pressure dropping rates. The maximum amount of retrograde condensation

4. Results and discussion To study the effect of non-equilibrium on gas-related properties, one gas condensate sample with 231 g/cm3 condensate content was chosen. The simulated gas condensate sample was divided into 7 pseudo components. The initial compositions of gas and liquid phases were calculated through a flash calculation at a condition that pressure slightly is below its dew point pressure. This is to ensure the estimated retrograde condensate saturation is small enough and coincides with the experiment. The calculated initial compositions at 351.95 K are listed in Table 2.As the interphase mass transfer continues to carry on, the difference between phases gradually reduces, and as long as the time is long enough the system would reach a dynamic equilibrium state finally. The following diffusion coefficient corrections were used to determinate the effective diffusion coefficients: Wilke-Lee correction [42] for gases and Wilke-Chang correction [43] for liquids. As it is known that the diffusion coefficients vary little with pressure, it is reasonable to consider constant diffusion coefficients. As the variables of molar fluxes needed to be solved are all quite small, in order to accurately simulate the results, the diffusion coefficients in this work were calculated at each pressure. Here, Dij = Dji . The calculated average diffusion coefficients of gas and liquid phases are presented in Tables 3 and 4. It is observed from these two tables that the diffusion coefficients decrease with the increasing carbon number, and the gas phase diffusion coefficients are almost one order of magnitude larger than those of the liquid phase. To analyze the process of interphase mass transfer, the mole fluxes of C1, C4 and C7+ are plotted versus pressure (Figs. 4–6). Overall, the curves are smooth, and they descend with the reduction of pressure and intersect with X axis. The intersection location moves toward the pressure shrinking side of X axis with increasing carbon number till there is no intersection point with the axis for heavy components. The nature change of the mass transfer before and after the intersection point indicates the directional change of mass transfer. That is to say, for a gas condensate fluid sample, the heavy components are not substantially evaporated and extracted, large portions of light and intermediate components cause phase transition from the initial gas state to liquid state and then back to gas state, these can also be observed from the gas condensate phase envelop. This explains why a lighter gas phase and a heavier liquid phase remained at the end of the experiment. From the observation of the shape of the curves, it can be seen that the tails do not level off. The primary reason is the system has not achieved a balanced state when the abandonment pressure was reached. It indicated that a non-equilibrium state has been maintained during the whole process of pressure depletion. Furthermore, the variations of the curves are not the same for two specified pressure dropping rates. As depicted in these figures, the calculated mole fluxes curves of 5 MPa/h are above the curves of 2 MPa/h on the right side of the intersection point, and an opposite trend is observed on the left side of the intersection point. It inferred that the system has a strong ability to extract and evaporate fluid components at a larger pressure dropping rate. The main reason is that there is no sufficient time for the gas and liquid phases to attain equilibrium during the non-equilibrium process. Thus, the greater the pressure dropping rate, the stronger the instability of the system, and the larger the difference between gas and liquid phases. A summary of the deviation results between the 5 MPa/h and 2 MPa/h is shown in Table 5. It is demonstrated that there is an obviously difference between these two pressure dropping rates, and it’s also evident that a higher pressure dropping rate can leads to a higher fluxes. Study [44] has shown that the concentration profile in the region close to the gas-liquid interface undergoes a rapid change. There will be a considerable concentration difference between the interface and the bulk phases, therefore the calculation of mass transfer and diffusion process cannot be ignored. Due to the assumption that the interface is in 683

Fuel 226 (2018) 675–685

H. Tu et al.

problem at high pressure stage while at a lower pressure stage this effect is too weak to be ignored. Under the circumstances where there are no other special problems for the gas condensate reservoir, a larger mining speed is beneficial to improve oil and gas recovery.

saturation of 2 MPa/h is 7.36% and 6.85% for 5 MPa/h. It demonstrates that a larger pressure dropping rate will strength the non-equilibrium effect and a relatively lower values for the retrograde condensation saturation is obtained. The standard deviations of non-equilibrium experimental data and the calculated values are also shown in Table 8. A lower standard deviation 0.07 of 5 MPa/h indicates the data points tend to be closer to the calculated value than that of 2 MPa/h. But overall, both the calculated data are in an acceptable range compared with experimental results.

Acknowledgment This work is supported by the National Natural Science Foundation of China (No. 51374179), “Theory and molecular dynamics study on CO2-crude oil non-equilibrium diffusion considering capillary pressure and adsorption” and the China Scholarship Council (No. 201708510114)..

5. Examples In the following study, it estimates the retrograde condensation saturation for different gas condensate samples through a series of pressure depletion rates. These examples have been selected from three different gas condensate reservoirs. The initial compositions of samples for both gas and liquid phases are presented in Table 9. The sample 1 is a condensate gas with low condensate oil content, 80 g/m3. The resulting retrograde condensation saturation has been plotted in Fig. 16. The curves start at the dew point pressure of 23.05 MPa, the condensate saturation then slowly increases to reach a maximum value at 11 MPa, and subsequently the condensate amount slightly declines. The trends of the curves are the same as the equilibrium experimental data. The maximum retrograde condensation saturations of 1 MPa/h, 2 MPa/h and 5 MPa/h are 1.925%, 1.857% and 1.708%, respectively. The sample 2 is a gas condensate with condensate oil content of 150 g/m3.The curve variation is basically the same as the sample 1. The condensate saturation first keeps ascending and then declines slowly. The condensate oil content of sample 3 is 235 g/m3. The trajectory of the curve is different from those of sample 1 and sample 2 where an obvious re-vaporization is observed leading to a significant downward trend in the curves (see Figs. 17 and 18). Table 10 is the relative absolute deviation between the equilibrium experiment data and the different depletion pressure rates. It is obviously that the calculated results are in a good agreement with the equilibrium experimental data. It is further indicated that the greater the depletion rate, the stronger the degree of deviations from the equilibrium state.

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.fuel.2018.04.048. References [1] Ivind F, Whitson CH. Modeling gas-condensate well deliverability. SPE Reservoir Eng 1996;11(4):221–30. [2] Hubbert MK. Physical principles of oil production. J Geol 1950;58(6):655–60. [3] Eilerts CK, Sumner EF, Potts NL. Integration of partial differential equation for transient radial flow of gas-condensate fluids in porous structures. SPEJ 1965:141–52. [4] Kniazeff VY, Navill SA. Two-phase flow of volatile hydrocarbons. SPEJ 1965:37–44. [5] O'dell HG, Miller RN. Successfully cycling a low-permeability high-yield gas condensate reservoir. J Pet Technol 1967;19(1):41–7. [6] Jones JR, Raghavan R. Interpretation of flowing well response in gas-condensate wells. SPE Form Eval 1988;3(3):578–94. [7] Jones JR, Vo DT, Raghavan R. Interpretation of pressure-buildup responses in gascondensate wells. SPE Form Eval 1989;4(1):93–104. [8] Fussell DD. Single-well performance predictions for gas condensate reservoirs. J Pet Technol 1973;25(7):860–70. [9] Mиpэзжaнзaцe AX. Natural gas extraction process. Translated by Zhu EL, Wei Z, Wang ZW; 1993. [10] Al-Wahaibi YM, Muggeridge AH, Grattoni CA. Gas-oil non-equilibrium in multicontact miscible displacements within homogeneous porous media. J Pet Sci Eng 2009;68(1–2):71–80. [11] Downar-Zapolski P, Bilicki Z, Bolle L, Franco J. The non-equilibrium relaxation model for one-dimensional flashing liquid flow. Int J Multiph Flow 1996;22(3):473–83. [12] Civan F, Rasmussen ML. Analysis and Interpretation of Gas diffusion in quiescent reservoir, drilling, and completion fluids: equilibrium vs. non-equilibrium models. In: Presented at the SPE annual technical conference and exhibition, Denver, Colorado; 2003. [13] Civan F. Including non-equilibrium relaxation in models for rapid multiphase flow in wells. In: Presented at the SPE Annual Technical Conference and Exhibition, Houston, USA; 2004. [14] Civan F. Modeling well performance under non-equilibrium deposition conditions. In: Presented at the SPE Production and Operations Symposium, Oklahoma; 2001. [15] Qi MM, Lei ZD, Kang XD, Xue XM, Lv SX. Study on inflow performance relationship of condensate reservoir considering high velocity flow effect. Oil Drill Prod Technol 2006;28(3):74–7. [16] Liu YJ, Li XF, Kang XD. Determination of reasonable pressure difference for condensate gas reservoir. Acta Petrol Sin 2006;27(2):85–8. [17] Wang CQ, Tang Y, Du ZM, Chen L, Sun Y, Pan Y, et al. Phase behaviors of condensate gas with vaporous water and liquid production characteristics in a nonequilibrium pressure drop process. Acta Petrol Sin 2006;34(4):740–6. [18] Wei DC, Li YJ, Yao L. Study on non-equilibrium effect of condensate gas system with water vapor under HT and HP. J Chongqing Univ Sci Technol: Nat Sci Ed 2013;15(3). [19] Murphree EV. Rectifying column calculations-with particular reference to N component mixtures. Ind Eng Chem 1925;17(7):747–50. [20] Saha S, Peng D. A non-equilibrium phase behaviour model for compositional reservoir simulation. In: Presented at the Annual Technical Meeting, Calgary, Alberta; 1992. [21] Indrupskiy IM, Lobanova OA, Zubov VR. Non-equilibrium phase behavior of hydrocarbons in compositional simulations and upscaling. Comput Geosci 2017;21(5–6):1173–88. [22] Long N, Sammon P. A non-equilibrium equation-of-state compositional simulator. In: Presented at the SPE reservior simulation symposium, Dallas Texas; 1997. [23] Lobanova OA, Indrupskiy IM. Non-equilibrium and scale effects in modeling phase behavior of hydrocarbon mixtures. Neft Khoz (Oil Industry) 2012(6):18–23. [24] Taylor R, Krishna R. Multicomponent mass transfer. New York: Wiley; 1993. [25] Taitel Y, Tamir A. Film condensation of multicomponent mixtures. Int J Multiph Flow 1974;1(5):697–714. [26] Estrin J, Hayes TW, Drew TB. The condensation of mixed vapors. AIChE J 1965;11(5):80–2.

6. Conclusions A simple method based on the film theory was presented to estimate the amount of retrograde condensation saturation at different pressure dropping rates. The model introduces the essential interphase mass transfer to better know about the characteristics for each component, and to predict mole fluxes, compositions of each component at the interface and in the bulk phase with regard to the variation of the pressure. The present model requires no equations of state calculation but only uses mass transfer and material balance equations to depict the non-equilibrium effect. This assumption is reasonable since most of the equations of state are aimed for equilibrium calculation. A physical method has been implemented to simulate a gas condensate with seven pseudo-component components which shows satisfactory calculation of the retrograde condensate saturation with 2 MPa/h and 5 MPa/h experimental data. The results demonstrate that the positive and negative quantities of the mole fluxes can well explain the retrograde condensation and re-vaporization process of gas condensate fluid. As the differences of each component characteristics, the gas phase tends to be lighter while the condensed liquid phase tends to be heavier. In general, the greater the pressure drop rate, the more serious the non-equilibrium effect. The so called “condensate hysteresis” occurred due to the insufficient time for both gas and condensed liquid phase to contact and the condensate liquid will be carried away by gas phase. The proposed method was confirmed from the calculation results of three samples with different condensate content. To summarize, for a gas condensate reservoir, non-equilibrium effect can relieve condensate blockage 684

Fuel 226 (2018) 675–685

H. Tu et al.

condensation by a film model. Chem Eng Sci 1982;37(1):117–9. [37] Furno JS, Taylor R, Krishna R. Condensation of vapor mixtures. 2. Comparison with experiment. Ind Eng Chem Process Des Dev 1986;25(1):98–101. [38] Kontogeorgis GM, Economou IG. Equations of state: from the ideas of van der Waals to association theories. J Supercrit Fluids 2010;55(2):421–37. [39] Drohm JK, Trengrove RD, Goldthorpe WH. On the quality of data from standard gas-condensate PVT experiments. In: Presented at the SPE Gas Technology Symposium, Dallas, Texas; 1988. [40] Aily ME, Khalil MHM, Desouky SM, Batanoni MH, Mahmoud MRM. Experimental studies on constant mass-volume depletion of gas-condensate systems. Egypt J Pet 2013;22(1):129–36. [41] Wang ZH, Tu HM, Guo P, Zeng FH, Sang TY, Wang ZD. Fluid behavior of gas condensate system with water vapor. Fluid Phase Equilib 2017;438:67–75. [42] Wilke CR, Lee CY. Estimation of diffusion coefficients for gases and vapors. Ind Eng Chem 1955;47(6):1253–7. [43] Wilke CR, Chang P. Correlation of diffusion coefficients in dilute solutions. AIChE J 1955;1(2):264–70. [44] Riazi MR. A new method for experimental measurement of diffusion coefficients in reservoir fluids. J Pet Sci Eng 1996;14(3–4):235–50.

[27] Sage FE, Estrin J. Film condensation from a ternary mixture of vapors upon a vertical surface. Int J Heat Mass Transf 1976;19(3):323–33. [28] Krishna R, Panchal CB. Condensation of a binary vapor mixture in the presence of an inert gas. Chem Eng Sci 1977;32(7):741–5. [29] Krishna R, Standart GL. Mass and energy transfer in multicomponent system. Chem Eng Commun 1979;3(4–5):201–75. [30] Wallace JL, Daison AW. Condensation of mixed vapors. Ind Eng Chem 2002;30(5768):948–53. [31] Sardesai RG. Studies in condensation [PhD. Thesis]. Manchester, England: University of Manchester Institute of Science and Technology; 1979. [32] Tamir A, Merchuk JC. Verification of a theoretical model for multicomponent condensation. Chem Eng J 1979;17(2):125–39. [33] Bandrowski J, Kubaczka A. On the condensation of multicomponent vaporous in the presence of inert gases. Int J Heat Mass Transf 1981;24(1):147–53. [34] Taylor R, Noah MK. Simulation of binary vapor condensation in the presence of an inert gas. Lett Heat Mass Transfer 1982;9(6):463–72. [35] Webb DR. Heat and mass transfer in condensation of multicomponent vapors. In: Presented at the seventh international heat transfer conference, Munich, Germany; 1982. [36] Webb DR, Taylor R. The efficient estimation of rates of multicomponent

685