Journal of Petroleum Science and Engineering 185 (2020) 106609

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering

journal homepage: http://www.elsevier.com/locate/petrol

Transient simulation of gas-condensate two-phase flow in pipes

Di Fan, Jing Gong *, Shengnan Zhang, Guoyun Shi, Qi Kang, Changchun Wu
National Engineering Laboratory for Pipeline Safety/ MOE Key Laboratory of Petroleum Engineering/ Beijing Key Laboratory of Urban Oil and Gas Distribution
Technology, China University of Petroleum-Beijing, Fuxue Road No. 18, Changping District, Beijing, 102249, PR China

A R T I C L E I N F O A B S T R A C T

Keywords: An algorithm for transient gas-condensate two-phase flow in pipes is developed in this paper. It adopts a two-
Gas-condensate fluid model for hydraulic simulation and a single energy model for thermal simulation. For ease of calcula­
Transient simulation tion, a temperature-based energy equation is derived from the law of conservation of energy. In order to consider
Two-fluid model
the influence of the phase transition, a new mass transfer model is developed in terms of total mass flow rate.
Mass transfer
Based on these models, a SIMPLE(Semi-Implicit Methods for Pressure-Linked Equations)-like algorithm is pro­
posed to simulate the transient gas-condensate two-phase flow. It consists of four steps: pressure and velocity
correction, volume fraction solution, temperature calculation and property parameters update. Case study shows
that the proposed algorithm could well simulate gas-condensate stratified flow with low liquid mass fraction and
high liquid mass fraction. In addition, the algorithm could provide a more consistent result for transient gas-
condensate two-phase flow.

1. Introduction equation and an empirical slip relationship between phase velocities.

Unfortunately, the slip relationship is usually related to the flow pattern.
Gas-liquid two-phase flow is often encountered in chemical engi­ At present, most closure relationships are constructed under mixed flow
neering (Yue et al., 2007; Chen et al., 2011; Li et al., 2018), petroleum regimes (Hibiki and Ishii, 2003; Evje and Flåtten, 2007; Liapidevskii and
industry (Bendiksen et al., 1991; Goldszal et al., 2007; Krasnopolsky Tikhonov, 2019). For stratified flow, the closure relationships are so
et al., 2016) and nuclear technology (Vernier and Delhaye, 1968; Ishii, limited that the two-fluid model is more prevalent. Although Wallis
1975; Bestion, 1990). In the exploitation of offshore condensate gas proposed a two-fluid model earlier (1969), the widely accepted form
fields, the gas-liquid mixture is usually transported to the onshore pro­ was established by Drew et al. (1979). Compared to the drift-flux model,
cessing terminal through a single pipeline. It leads to the unique the biggest difference is that gas and liquid have separate momentum
gas-condensate two-phase transient flow in the petroleum industry. equations in the two-fluid model. Despite this, they are not completely
Since the transport distance usually reaches tens or even hundreds of independent but interact with each other. To describe their interaction,
kilometers, both pressure and temperature will change greatly so that the interfacial shear stress has been proposed in a lot of studies (Cohen
neither of them could be ignored. In addition, due to the change of and Hanratty, 1968; Taitel and Dukler, 1976; Andritsos and Hanratty,
pressure and temperature, phase transition will occur between the gas 1987; Andreussi and Persen, 1987; Xiao et al., 1990). However, the
phase and the liquid phase. To better simulate the transient flow, the two-fluid model only contains two continuous phases: the gas phase and
influence of mass transfer caused by phase transition has to be consid­ the liquid phase. If another dispersed phase such as a bubble or a droplet
ered in the algorithm. appears, its performance will be affected a lot. Therefore, many efforts
In terms of the hydraulic model, there are two main types: the drift- had been made to introduce dispersed phase into the two-fluid model
flux model and the two-fluid model. The concept of drift-flux was pro­ (Bendiksen et al., 1991; Danielson et al., 2005; Goldszal et al., 2007;
posed by Zuber and Findlay (1965). In this model, the phase fraction is Bonizzi et al., 2009). In these models, the dispersed phase also had a
affected by the relative velocity between the gas and the liquid. Mean­ mass equation and even a momentum equation. To a certain extent, they
while, the distribution of each phase was uneven in the cross section. were no longer the traditional two-fluid model which only considered
Later, Ishii (1975) established a general form of drift-flux model. It the continuous phase. It was more appropriate to call them multi-field
consists of two mass conservation equations, a mixing momentum model.

* Corresponding author.
E-mail address: ydgj@cup.edu.cn (J. Gong).

https://doi.org/10.1016/j.petrol.2019.106609
Received 29 July 2019; Received in revised form 16 October 2019; Accepted 19 October 2019
Available online 21 October 2019
0920-4105/© 2019 Elsevier B.V. All rights reserved.
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

The study of the thermal model is far less attractive than that of the 2. Models
hydraulic model. In most studies, investigators adopted a single energy
equation to describe the thermal process (Bendiksen et al., 1991; Pau­ 2.1. Governing equations
chon et al., 1993; Henriot and Pauchon, 1997; Henriot et al., 2002; Abili
and Kara, 2014). Despite this, their specific forms were not the same Although multi-field model introduces a description of the dispersed
owing to the different hydraulic models. For drift-flux model, mixing phase, it needs more closure relationships, e.g., droplet entrainment
velocity had been obtained through the mixture momentum equation. rate, drag on droplet and droplet deposition rate. These relationships are
Therefore, it was easy to give a mixing energy equation. Unfortunately, not unique and their choice may affect the final simulation results. In
in two-fluid model, the momentum equations were separate for gas and order to minimize the influence of closure relationships, a two-fluid
liquid, which makes it difficult to give such a mixing energy equation. model is adopted in this paper since it requires only a few closure re­
Bendiksen et al. (1991) proposed a solution: they first gave the energy lationships. In this model, the mass balance equation for each phase is
equation for each phase. Afterwards, assuming that they have the same written as:
temperature, the two energy equations could be combined into a single
∂ � ∂ �
energy equation. Since the two-fluid model was suited for separate en­ α g ρg þ αg ρg vg ¼ ϕg (1)
∂t ∂x
ergy equations, some researchers had abandoned the assumption of
equal temperature (Bestion, 1990; Tiselj and Petelin, 1997; Zhou and ∂ ∂
Podowski, 2001; Raimondi, 2017). In their thermal models, each phase ðα ρ Þ þ ðα ρ v Þ ¼ ϕl (2)
∂t l l ∂x l l l
had an energy equation. Whether it is a single energy equation or
The momentum conservation equation for each phase is written as:
multiple energy equations, almost all of them are given in the form of
internal energy and enthalpy. As all property parameters are tempera­ ∂ � ∂� � ∂p τwg Sg τi Si
αρv þ α ρ v2 ¼ αg αg ρg g sin θ þ ϕg vi (3)
ture dependent, it is more convenient to take the temperature as the only ∂t g g g ∂x g g g ∂x A A
independent variable in energy equations.
There are two main mass transfer models: one based on the total heat ∂
ðαl ρl vl Þ þ
∂ �
αl ρl v2l ¼ αl
∂p
αl ρl g sin θ
τwl Sl
þ
τ i Si
þ ϕl vi (4)
transfer and the other based on the phase fraction. In nuclear technol­ ∂t ∂x ∂x A A
ogy, since only water and steam exist in the cooling pipe, a mass transfer In Eqs. (1)–(4), α, ρ, p and v are the volume fraction, density, pressure
model is established using the total heat transfer and phase transition and velocity. τwg and τwl are the wall shear stress while τi is the interfacial
enthalpy (Bestion, 1990; Tiselj and Petelin, 1997). However, because of shear stress. Sg, Sl and Si are the wetted perimeters of the gas, liquid, and
the complex composition of gas-condensate, it is difficult to determine the interface. ϕ is the mass transfer rate between the gas phase and the
the phase transition enthalpy. Therefore, the phase fraction is intro­ liquid phase. A is the pipe cross-sectional area. vi is the velocity before
duced to calculate mass transfer in the petroleum industry (Bendiksen the phase transition. Subscripts g and l indicate the gas phase and the
et al., 1991; Krasnopolsky et al., 2016). Several parameters could liquid phase, respectively.
describe the phase fraction, e.g. Gas-Oil Ratio (GOR), Gas-Liquid Ratio For the volume fraction α and the mass transfer rate ϕ, they obey the
(GLR), Equilibrium Gas mass fraction (EGF) etc. GOR and GLR are relation shown in Eqs. (5) and (6). Eqs. (1)–(6) form a compete two-fluid
usually obtained by measurement in the field while EGF is usually ob­ hydraulic model. There are two significant differences compared to
tained by the equation of state (Soave, 1972; Patel and Teja, 1982; models that ignore phase transition. For the mass conservation equation,
Stryjek and Vera, 1986). For gas-condensate, EGF is usually sued to mass transfer occurs as a mass source so that the right side of Eqs. (1) and
describe the gas mass fraction because of its relatively clear composi­ (2) is no longer zero. For the momentum conservation equation, mass
tion. That’s why OLGA proposed such a mass transfer model based on transfer also brings momentum exchange. Therefore, there is one more
EGF and total fluid (Bendiksen et al., 1991). However, this model ne­ term on the right side of Equations (3) and (4).
glects the velocity difference between gas and liquid so that the mass
transfer is usually overestimated. αg þ αl ¼ 1 (5)
Several numerical methods have been used to simulate transient gas-
ϕg þ ϕl ¼ 0 (6)
liquid two-phase flow, e.g. explicit Godunov-type schemes (Saurel and
Abgrall, 1999), Semi-Implicit Methods for Pressure-Linked Equations In this paper, a single energy model is adopted to simulate the
(SIMPLE) (Morales-Ruiz et al., 2012; Krasnopolsky et al., 2016) and fully transient thermal process. It could be written as:
implicit Newton-Raphson-based iterative methods. However, due to the
X∂ � � �� X �
∂
� ��
characteristics of gas-condensate two-phase flow in pipes, a complete 1
ρk αk Ek þ v2k þ gz þ
1
ρk αk vk Hk þ v2k þ gz ¼ Q (7)
algorithm should include mass transfer estimation, temperature calcu­ k¼g;l
∂t 2 k¼g;l
∂x 2
lation and property parameters update. Therefore, none of the existing
algorithms could be directly used for the transient simulation of Here, Q is the total heat transfer between the fluid and the environment.
gas-condensate two-phase flow in pipes. E is the internal energy while H is the enthalpy. They have the following
To overcome above problems, an algorithm for the transient simu­ relationship:
lation of gas-condensate two-phase flow is developed in this paper. For p
thermal state, it adopts a temperature-based energy model. In order to Ek ¼ Hk (8)
ρk
better consider the influence of mass transfer, a new model is derived in
terms of total mass flow rate. Based on them, a complete algorithm is Eq. (7) is not suitable for simulation because the temperature is
proposed to simulate the transient gas-condensate two-phase flow in implicitly included. Therefore, a transformation is required so that the
pipes, which includes pressure and velocity correction, volume fraction temperature appears explicitly in the energy equation. For simplicity,
solution, temperature calculation and property parameters update. we assume that the mass transfer has little effect on energy so that it
This paper is organized as follows: the models are introduced in could be ignored. Substituting Eqs. (1)–(4) and (8) into Eq. (7) yields an
section 2, including the governing equations, shear stress model, mass energy equation containing only enthalpy H.
transfer model and heat transfer model. In section 3, we introduce the
numerical scheme, which is composed of the pressure equation, dis­
cretization method and correction method. In section 4, two cases, one
low liquid mass fraction pipe and one high liquid mass fraction pipe, are
presented separately. Conclusions are drawn in section 5.

2
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

8 34σ
>
> N N � 0:005
> ρg v2l we μ
<
ε¼ (15)
>
> 0:3
: 170σðNwe Nμ Þ Nwe Nμ > 0:005
>
ρg v2l

ρg v2l ε
Nwe ¼ (16)
σ

Fig. 1. A control volume with I-I as inlet section and II-II as outlet section. μ2l
Nμ ¼ (17)
ρg σε
X� ∂ ∂
� X ∂ �
∂p τwk Sk vk
�
ðρk αk Hk Þ þ ðρ α v H Þ ¼ Q þ ðpαk Þ þ αk vk þ
k¼g;l
∂t ∂x k k k k k¼g;l
∂t ∂x A 2.3. Mass transfer
τi Si vg τi Si vl
þ
A A Since the transient process in gas-condensate pipe is usually very
(9) slow, an assumption is made that the flow is in steady state. Considering
a control volume of length dx shown in Fig. 1, the inlet gas mass through
The enthalpy H could be calculated using Eq. 10 I-I during time dt can be calculated by Eq. (18) while the outlet gas mass
Hk ¼ Cpk T (10) through II-II can be calculated by Eq. (19). In these two formulas, Xe is
the equilibrium gas mass fraction (EGF) which can be determined using
Substituting Eq. (10) into Eq. (9) yields the final form of energy the equation of state.
equation in terms of temperature, as follows: � � �
min ¼ ρg vg αg A þ ρl vl αl A Xe dt (18)
X� ∂ � ∂ �
�
ρk αk Cpk Tk þ ρk αk vk Cpk Tk ¼ Q � � � � �
k¼g;l
∂t ∂x � ∂ ρg vg αg A þ ρl vl αl A Xe
X� ∂ � mout ¼ ρg vg αg A þ ρl vl αl A Xe þ dx dt (19)
∂p τwk Sk vk τi Si vg τi Si vl ∂x
þ ðpαk Þ þ αk vk þ þ (11)
∂t ∂x A A A
k¼g;l Due to the influence of mass transfer, the inlet gas mass flow does not
Eqs. (1)–(6) and (11) constitute the system of differential equations equal to that of the outlet even in the steady state. Thus the mass transfer
used in this paper. During simulation, it still needs several closure re­ per unit time per unit volume could be modeled as follows:
lationships which are introduced in section 2.2-2.4. � � � � � �
mout min ∂ ρg vg αg A þ ρl vl αl A Xe dxdt ∂ ρg vg αg þ ρl vl αl Xe
ϕg ¼ ¼ ¼
Adxdt ∂xAdxdt ∂x
2.2. Shear stress (20)

In the governing equations, the shear stress, τ, can be calculated by 2.4. Heat transfer
Eqs. (12) and (13).
1 Still using the control volume shown in Fig. 1, the total heat transfer
τwk ¼ fk ρk v2k (12)
2 between the fluid and the ambient environment during time dt is given
by:
1
τi ¼ fi ρg v2g (13) ΔQ ¼ KðT Te ÞπDdxdt (21)
2
Where K is the total heat transfer coefficient from the fluid to the
where the friction coefficient, f, can be computed by Eq. (14) (Xiao et al.,
ambient environment. T is the fluid temperature while Te is the envi­
1990).
ronmental temperature. Therefore, in terms of unit time, unit cross-
� �
1 2ε 9:35 sectional area and unit pipe length, the heat transfer from the fluid to
pffiffi ¼ 3:48 4lg þ pffiffi (14)
f D Re f the environment can be calculated as follows:

In Eq. (14), ε is the absolute roughness. Since there are two kinds of KðT Te ÞπDdxdt KðT Te ÞπD
Q¼ ¼ (22)
friction coefficient, ε is determined through different ways. For the Adxdt A
friction coefficient between the fluid and the pipe, ε is the roughness of
the pipe wall. However, for the interfacial friction coefficient, an initial 3. Numerical scheme
value will be calculated by Eq. (15). Meanwhile, ε should be bounded
between the pipe’s absolute roughness and 0.25(h/D). Where the Weber 3.1. Pressure equation
number Nwe and the Liquid Viscosity number Nμ are defined by Eqs. (16)
and (17), respectively. We adopt a SIMPLE-like algorithm to simulate the transient gas-

Fig. 2. Schematic diagram of calculation unit.

3
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

condensate two-phase flow in pipes. As we all know, the SIMPLE algo­ Table 1
rithm is initially proposed for single-phase flow, where the mass balance Length and elevation data of the field pipe.
equation is used to correct pressure (Patankar and Spalding, 1972). For a Length(km) 0 27 47 54 64 69 74
gas-liquid two-phase flow, a similar equation can be obtained by adding Elevation(m) 974 973 1009 1000 1024 1026 1040
the gas mass equation (Eq. (1)) and the liquid mass equation (Eq. (2)) Length(km) 84 94 101 111 126 131 144
Elevation(m) 1026 1055 1052 1068 1014 1015 965
together (Moukalled and Darwish, 2004a,b). However, in order to avoid
the effect of different orders of magnitude, the two equations are divided
by their densities respectively before the addition. Therefore, the pres­ on a staggered grid as shown in Fig. 2. The scalar variables (pressure,
sure corrected equation is as follows: temperature, volume fraction, density) are arranged at the cell centers
while the vector variables (velocity) are arranged at the cell edges. Other
Xαk ∂ρ X1 ∂ Xϕ
k
þ ðαk ρk vk Þ ¼ k
(23) property parameters e.g. viscosity, heat capacity, liquid surface tension
ρ
k¼g;l k
∂t ρ
k¼g;l k
∂x ρ
k¼g;l k etc. are also located at the cell centers. If a variable is required but not
defined at that location, it will be calculated by the first order upwind
It should be noted that due to mass transfer, the right side of Eq. (23)
scheme, e.g.
is no longer zero. This is different from single-phase flow and gas-liquid �
two-phase flow without considering phase transition. It shows the ρ 1 if vg;j 1=2 > 0
ρg;j 1=2 ¼ ρg;j otherwise (24)
impact of mass transfer on pressure calculations. g;j

In the discretization, two momentum equations (Eqs. (3) and (4)) are
performed on the momentum control volume and other equations (Eqs.
3.2. Discretization method
(1), (2), (11) and (23)) are performed on the mass control volume, as
follows:
The governing equations are discretized by a finite volume method

1h t *
i 1 h t # # * i
α #
ρ
1=2 g;j 1=2 vg;j 1=2 αtg;j t
ρ t
1=2 g;j 1=2 vg;j 1=2 þ α ρ v v αtg;j 1 ρ#g;j 1 v#g;j 1 v*g;j ¼
Δt g;j Δx g;j g;j g;j g;j 1

(25)
p#j p#j 1 1� # �
α t
g;j 1=2 α t #
ρ
g;j 1=2 g;j 1=2 g sin θ τ t
1=2 Sg þτ # t
i;j 1=2 Si þ ϕtg;j 1=2 vtg;j 1=2
Δx A wg;j

1h t *
i 1 h t # # * i
α #
ρ
1=2 l;j 1=2 vl;j 1=2 αtl;j ρt t
1=2 l;j 1=2 vl;j 1=2 þ α ρ v v αtl;j 1 ρ#l;j 1 v#l;j 1 v*l;j ¼
Δt l;j Δx l;j l;j l;j l;j 1

(26)
p#j p#j 1 1� # t t
�
t
αtl;j 1=2 αtl;j #
1=2 ρl;j 1=2 g sin θ τ 1=2 Sl τ#i;j 1=2 Si þ ϕl;j
t
1=2 vl;j
Δx A wl;j 1=2

� � � �
1 tþΔt # 1
αg;j ρg;j αtg;j ρtg;j þ αtþΔt ρ# v# αtþΔt ρ# #
1=2 vg;j ¼ϕtg;j
Δt Δx g;jþ1=2 g;jþ1=2 g;jþ1=2 g;j 1=2 g;j 1=2

(27)
� � � �
1 tþΔt # 1
α ρ αtl;j ρtl;j þ αtþΔt ρ# v# αtþΔt #
l;j 1=2 ρl;j
#
1=2 vl;j 1=2 ¼ ϕtl;j
Δt l;j l;j Δx l;jþ1=2 l;jþ1=2 l;jþ1=2
(28)

Table 2
Composition of the condensate gas.
number component mol% mole weight(g/mol)

1 C1 86.37 16.043
2 C2 6.59 30.07
3 C3 1.34 44.097
4 i-C4 0.32 58.124
5 n-C4 0.42 58.124
6 i-C5 0.18 72.151
7 n-C5 0.18 72.151
8 C6 0.24 86.178
9 C7þ 1.9 171.21
10 CO2 0.16 44.01
11 N2 2.3 28.014

Fig. 3. The algorithm flow chart.

4
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

1150 55
OLGA
1125 50
TNSP
1100 TNSP-N 45

Geometry (m)
Geometry
1075 40

1050 35
30
1025
25
1000
20
975
15
950
0 20 40 60 80 100 120 140
Length (km)

Fig. 4. Outcomes of the field pipe: (a) Pressure; (b)Temperature; (c) Gas velocity; (d) Liquid velocity; (e) Holdup; (f) Mass flow rate.

Xαtk;j ρtþΔt
k;j ρtk;j
t
ρ
k¼g;l k;j
Δt
where the superscripts t and t þ Δt represent the value of the current
X 1 αtk;jþ1=2 ρtþΔt tþΔt
k;jþ1=2 vk;jþ1=2 αtk;j tþΔt
ρ tþΔt
1=2 k;j 1=2 vk;j 1=2
Xϕtk;j time step and the value of the next time step, respectively. The super­
þ ¼ (29)
script * and # represent the value to be solved during the current iter­
t t
ρ
k¼g;l k;j
Δx ρ
k¼g;l k;j
ation and the value of the last iteration, respectively. It seems that two

X� 1 � �#
ρ#k;j αtþΔt t tþΔt
k;j Cpk;j T k;j ρtk;j αtk;j Ctpk;j T tk;j þ
k¼g;l
Δt

X 1 �� � � ��
ρ#k;jþ1=2 αtþΔt (30)
# t tþΔt
k;jþ1=2 vk;jþ1=2 Cpk;jþ1=2 T k;jþ1=2 ρ#k;j tþΔt
α # t tþΔt
1=2 k;j 1=2 vk;j 1=2 C pk;j 1=2 T k;j 1=2 ¼ QtþΔt
k;j þ
k¼g;l
Δx

X� 1 � �
1
� �
τ#wk;j Stk;j v#k;j
�
τ#i;j Sti;j v#g;j τ#i;j Sti;j v#l;j
p# tþΔt
k;j αk;j ptk;j αtk;j þ αtþΔt #
k;j vk;j p# p# þ þ
k¼g;l
Δt Δx k;jþ1=2 k;j 1=2
A A A

5
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

Fig. 5. Mass flow rate at different times (t ¼ 0 represents the time when the mass flow change at the inlet has just been completed).

unknowns (α and T) appear simultaneously in Eq. (30). In fact, since the � �

* *
solution is sequential, only the temperature is unknown while the vol­ a#j 3=2 vg;j 3=2 þ a#j 1=2 vg;j 1=2 þ a#jþ1=2 v*g;jþ1=2 ¼ dtj 1=2 p#
j p# #
j 1 þ bj 1=2
ume fraction α is known. It will be described in detail in section 3.4. (31)

3.3. Correction method For gas velocity, pressure and gas density at the next time step, they
can all be written as the sum of the current value and the correction
Taking the gas momentum equation (Eq. (25)) as an example, using value, as follows:
the first order upwind scheme, it could be rewritten as: vtþΔt ¼ v*g þ v’g (32)
g

6
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

Fig. 5. (continued).

ptþΔt ¼ p# þ p; (33) second-order correction, the pressure corrected equation yields:

a# ’
j 1 pj þ a# ’ # ’ #
j pj þ ajþ1 pjþ1 ¼ bj (40)
∂ρ#g ’ 1
ρ tþΔt #
¼ρ þ ρ ¼ ρ þ ’
p #
(34)
Using this equation, the pressure correction value could be
g g g g
∂p
tþΔt tþΔt tþΔt computed.
Since the v , p and ρ also satisfy the momentum equation, it
yields
� � � � � � 3.4. Simulation process
a#j 3=2 v*g;j 3=2 þ v’g;j 3=2 þ a#j 1=2 v*g;j 1=2 þ v’g;j 1=2 þ a#jþ1=2 v*g;jþ1=2 þ v’g;jþ1=2
h� � � �i Since the composition of gas-condensate is very complex, the prop­
¼ dtj 1=2 p# p# þ b# erty table method is adopted to update related parameters. It is assumed
’ ’
j þ pj j 1 þ pj 1 j 1=2

(35) that the composition of fluid is constant along the pipe. Consequently,
fluid properties e.g. viscosity, heat capacity, liquid surface tension et al.,
Because v* is solved from Eq. (31), it could be eliminated from Eq. only vary with temperature and pressure. With this assumption, a table
(35). We have: containing all property parameters could be calculated under a wider
� � range of pressure and temperature in advance. In simulation, all pa­
a#j 3=2 v’g;j 3=2 þ a#j 1=2 v’g;j 1=2 þ a#jþ1=2 v’g;jþ1=2 ¼ dtj 1=2 p’j p’j 1 (36)
rameters could be updated by interpolation from this table. Under the
For simplicity, the influences of node j-3/2 and j þ 1/2 are ignored, assumption of constant composition, another important parameter, Xe,
the velocity correction value could be computed as follows: is also only related to temperature and pressure. Therefore, it is regarded
as a special kind of property parameter and updated in the same way as
dtj � � other property parameters.
(37)
1=2
v’g;j ¼ p’j p’j
1=2
a#
j 1=2
1
The whole algorithm is shown in Fig. 3. When the simulation begins,
the initialization is executed first. It usually involves reading the basic
Eqs. (34) and (37) show that both the gas density correction value parameters of the pipe, reading the property table etc. Then, the
and the gas velocity correction value could be represented by the pres­ correction process is carried out until the pressure correction value is
sure correction value. For liquid phase, there are also two similar small enough. Currently, an error criterion of 10 3 is adopted in our
relationships. program. Afterwards, the volume fraction of each phase is computed by
mass conservation equation. With the new pressure, velocity and vol­
∂ρ#l ’
ρtþΔt
l ¼ ρ# ’ #
l þ ρl ¼ ρl þ p (38) ume fraction, the algorithm will calculate the temperature at the next
∂p
time step. Finally, using new pressure and temperature, all the property
d’j � � parameters are updated from the property table for the next calculation.
(39) Based on this algorithm, we developed a Transient Natural gas-
1=2
v’l;j 1=2 ¼ p’j p’j 1
a#
j 1=2 condensate Simulation Program (TNSP) using Java.
Substituting Eqs. (34) and (37)–(39) into Eq. (29) and neglecting the

7
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

Fig. 6. Pressure at different times: (a) Calculated by TNSP; (b) Calculated by OLGA; (t ¼ 0 represents the time when the mass flow change at the inlet has just
been completed).

1150 56 1150 56
0 0.25h 0.5h 0 0.25h 0.5h
1125 1h 2h 4h
48 1125 1h 2h 4h
48
8h 16h 32h 8h 16h 32h
Geometry (m)

Geometry (m)
1100 64h 128h Geometry 1100 64h 128h Geometry
40 40
1075 1075
1050 32 1050 32

1025 24 1025 24
1000 1000
16 16
975 975
8 8
950 950
0 20 40 60 80 100 120 140 0 20 40 60 80 100 120 140
Length (km) Length (km)

(a) (b)
Fig. 7. Temperature at different times: (a) Calculated by TNSP; (b) Calculated by OLGA; (t ¼ 0 represents the time when the mass flow change at the inlet has just
been completed).

4. Case study pressure is 105.3 bar. For thermal state, the inlet fluid temperature is
50 � C and the ambient temperature is 5 � C. It is simulated by both TNSP
In the simulation of gas-condensate two phase flow, two boundary and OLGA and the outcomes are shown in Fig. 4. For comparison, the
conditions are often encountered: constant boundary condition (CBC) results of neglecting mass transfer (TNSP-N) are also given by our
and variable boundary condition (VBC). For the CBC case, both the inlet program.
mass flow and the outlet pressure are independent of time. In contrast, if As can be seen from Fig. 4 (a), the pressure results given by the two
either of them is time dependent, it is the VBC case. To test the perfor­ programs are very close. For our program, the pressure of TNSP is a little
mance of TNSP under both boundary conditions, two gas-condensate lower than that of TNSP-N. As shown in Fig. 4 (c) and (d), the TNSP-N’s
pipes are simulated in this section. One is a low liquid mass fraction gas velocity is greater than that of TNSP while their liquid velocities are
pipe and another is a high liquid mass fraction pipe. During simulation, almost the same. Since the pressure drop is related to the phase velocity,
the grid size is 200 m and the max time step is 1 s. For comparison, these a larger gas velocity generates a larger friction pressure drop. Therefore,
cases are also simulated by OLGA using the same grid size and time step the TNSP-N yields the greater pressure drop. As for the pressure differ­
with Lookup table method. All of these cases are performed on a com­ ences between our program and OLGA, it is likely due to their different
puter equipped with an Intel(R) Core(TM) i3-8100 CPU and 8.00 GB of friction coefficients. For temperature shown in Fig. 4 (b), both outcomes
RAM. Meanwhile, to avoid the influence of property parameters, both of TNSP-N and TNSP are similar with that of OLGA. Since the mass
the property table for OLGA and the property table for TNSP used in transfer is neglected in the energy model of our program, it shows that
these cases are obtained by the PR equation of state in the commercial the approximation is acceptable in terms of CBC cases. Another very
software PVTsim. important parameter presented in Fig. 4 (e), the liquid holdup, indicates

4.1. Application in a low liquid mass fraction pipe

Table 3
Composition of the condensate gas.
The first example is a real field pipe. The length of the pipe is 144 km
and the inner diameter is 575 mm. More detailed length and elevation number component mol% mole weight(g/mol)

data are shown in Table 1 (Wang, 2015). For internal pipe roughness, a 1 C1 59.04 16.043
fixed value of 0.02 mm is assumed. The total heat transfer coefficient is 2 C2 10.42 30.07
3 C3 15.12 44.097
1.0 W/m2⋅K 1. The composition of the condensate gas is listed in
4 i-C4 2.39 58.124
Table 2. At 105.3 bar and 5 � C (Outlet pressure and ambient tempera­ 5 n-C4 7.33 58.124
ture), the liquid mass fraction is 20.13%. 6 i-C5 2.00 72.151
7 n-C5 1.72 72.151
4.1.1. CBC case 8 C6 1.18 86.178
9 C7þ 0.80 171.21
In this case, the inlet total mass flow is 82.25 kg/s while the outlet

8
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

Fig. 8. Outcomes of the virtual pipe: (a) Pressure; (b)Temperature; (c) Gas velocity; (d) Liquid velocity; (e) Holdup; (f) Mass flow rate.

that the results of TNSP are closer to that of OLGA when mass transfer is 4.1.2. VBC case
considered. In the processing of gas-condensate, liquid flow has an During the production of gas-condensate, some gas wells may be
important influence on the efficiency of the separator. Fig. 4 (f) shows temporarily shut down for equipment maintenance, resulting in a sharp
that if mass transfer is ignored, the predicted liquid flow at the outlet drop in mass flow at the inlet. In this section, we assume that due to the
will have a large deviation. Instead, if mass transfer is considered, the adjustment of gas field production, the inlet mass flow rate changes from
liquid flow provided by TNSP is almost the same as that of OLGA. All of 82.25 kg/s to 40 kg/s in 1000 s. Other parameters are the same as those
the outcomes indicate that TNSP could deal with CBC cases in section 4.1.1. Through this case, the effect of the new mass transfer
appropriately. model under variable boundary conditions is tested. Since the phase

9
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

Fig. 9. Outcomes at different times (t ¼ 0 represents the time when the mass flow change at the inlet has just been completed).

velocity and the phase fraction interact with each other, for ease of that of TNSP in the first hour. However, their liquid mass flow rates are
comparison, the mass flow is used instead of presenting velocity and almost the same. It is likely due to the fact that the two programs use
holdup separately. The outcomes are shown in Fig. 5, Fig. 6 and Fig. 7. different mass transfer models. In commercial software OLGA, the mass
OLGA takes 58 min 27.28 s to simulate this case while TNSP takes transfer rate is computed as follows (Bendiksen et al., 1991):
49 min 33.61 s. "� � � � � � � � #
The zero point of time in the figure is the moment when the mass ϕg ¼
∂Xe ∂p
þ
∂Xe ∂p ∂x
þ
∂Xe ∂T
þ
∂Xe ∂T ∂x
mg þ ml
�
flow change at the inlet has just been completed. From Fig. 5 (a)–(d), it ∂p T ∂t ∂p T ∂x ∂t ∂T p ∂t ∂T p ∂x ∂t
can be seen that the gas mass flow rate given by OLGA is greater than (41)

10
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

50 0 0.25h 50 0 0.25h
0.5h 1h 0.5h 1h
2h 4h 2h 4h
40 40
8h 16h 8h 16h

30 30

20 20

10 10

0 0
0 10 20 30 40 50 60 0 10 20 30 40 50 60
Length (km) Length (km)

Fig. 9. (continued).

where mg and ml are the gas phase mass and the liquid phase mass in the 4.2. Application in a high liquid mass fraction pipe
pipe section, respectively. Since Xe is calculated in a stationary state, this
model neglects the velocity difference between the gas phase and the In this section, a virtual pipe is presented to test TNSP under high
liquid phase. Actually, because the gas phase always flows faster than liquid mass fraction. The composition of fluid is shown in Table 3. At
the liquid, the gas mass fraction in a pipe section is always smaller than 50 bar and 5 � C (Outlet pressure and ambient temperature), the liquid
that of equilibrium gas mass fraction Xe. Therefore, the OLGA model mass fraction is 60.67%. The pipe is 54 km long with an inner diameter
tends to overestimate the amount of mass transfer. That’s why the gas of 409.4 mm. The internal pipe roughness is also 0.02 mm. The total heat
mass flow during the first hour is much larger than that of TNSP. transfer coefficient is 4.0 W/m2⋅K 1. Moreover, in order to avoid the
Afterwards, the change of OLGA’s gas mass flow rate begins to influence of topographic relief, the pipe is set as horizontal.
accelerate. Fig. 5 (e) shows the gas mass flow rates given by two pro­
grams are almost the same from the second hour. For liquid mass flow, 4.2.1. CBC case
there is always a pipe section that OLGA gives a lower result than that of For hydraulic state, the inlet total mass flow rate is 50 kg/s while the
TNSP (shown in Fig. 5 (e)–(j)). It appears that a large amount of liquid is outlet pressure is 50 bar. For thermal state, the inlet fluid temperature is
used to replenish the vaporized fluid in these sections. Therefore, the 50 � C and the ambient temperature is also 5 � C. It is simulated by OLGA,
liquid mass flow rates are much smaller than those in TNSP. When this TNSP and TNSP-N. The outcomes are shown in Fig. 8.
process is finished, TNSP and OLGA reach the same state as shown in As shown in Fig. 8 (a) and (b), both pressure and temperature of
Fig. 5 (k). TNSP and TNSP-N are similar to those of OLGA. For gas velocity, liquid
Fig. 6 shows that during this process, the pressure trends given by the velocity and holdup (Fig. 8 (c)–(e)), when mass transfer is taken into
two programs are almost the same. However, for temperature shown in consideration, the trends given by TNSP are consistent with those of
Fig. 7, although they finally reach the same state at 128 h, their change OLGA. However, if it is neglected, there will be a great deviation.
processes are a little different. The result given by TNSP has more Because of the interaction among these three variables, it makes no
fluctuations than that of OLGA. It is probably because the mass transfer sense to compare their specific values separately. Therefore, the mass
is neglected in our energy model. Since the differences of most param­ flow rate is compared here. As shown in Fig. 8 (f), by applying the mass
eters are small, TNSP performs well in the transient process. transfer model proposed, both gas mass flow rate and liquid mass flow
rate are almost the same as those of OLGA. It is worth noting that the
flow in the front of the pipe is dominated by gas while the rear becomes
liquid dominated. That means even in this complex flow, TNSP could

11
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

achieve good results. h Liquid film height, m

m mass of fluid, kg
4.2.2. VBC case Nwe Weber number
In this case, the inlet total mass flow rate evenly decreases from Nμ Liquid viscosity number
50 kg/s to 25 kg/s within 1000s while other parameters are the same as p Pressure, Pa
those in section 4.2.1. The outcomes are shown in Fig. 9. OLGA takes Q Heat transfer, w/m3
2 min 50.65 s to simulate this case while TNSP takes 3 min 17.72 s. Re Reynolds number
In general, TNSP gives a more consistent result than that of OLGA. As S Wetted perimeter or interfacial contact length, m
can be seen from Fig. 9 (a), (c) and (e), when the change of inlet total t Time, s
mass flow is finished, the gas mass flow, the liquid mass flow and the v Phase velocity, m/s
holdup along the pipe begin to fluctuate synchronously. But the results vi Velocity before the phase transition, m/s
given by OLGA are a little complex. Only the liquid mass flow and the x Mileage coordinate, m
holdup change synchronously (shown in Fig. 9 (d) and (f)). As for the gas Xe Equilibrium gas mass fraction
mass flow, the fluctuation has completely dissipated. It is also caused by z Height coordinate, m
the mass transfer model in OLGA. Since the mass transfer is over­
estimated, a large amount of fluid evaporates into the gas. Consequently, Greek Symbols
the gas mass flow rate changes slowly. Just from these three parameters, α Volume fraction
TNSP reaches the new steady state much earlier than OLGA. However, ε Absolute roughness, m
for pressure presented in Fig. 9 (g) and (h), it shows that OLGA changes θ Pipe inclination, rad
more quickly than TNSP. This phenomenon may also be caused by the μ Viscosity, Pa⋅s
different mass transfer models and a detailed research will be conducted ρ Density, kg/m3
in the future. For temperature shown in Fig. 9 (i) and (j), their trends are σ Surface tension, N/m
almost the same. Although the change process is slightly different, both τ Shear stress, Pa
of them reach the new steady state at the 16th hour. Therefore, TNSP ϕ Mass transfer rate, Kg/(m3⋅s)
could simulate liquid-dominated VBC cases very well.
Subscripts
5. Conclusions g Gas
i Interface
A transient algorithm for the simulation of gas-condensate two-phase in inlet
flow in pipes is presented in this paper. It adopts a two-fluid hydraulic l Liquid
model and a single thermal model, which include two mass equations, out outlet
two momentum equations and one energy equation. Due to the char­ w Wall
acteristics of gas-condensate two-phase flow in pipes, both thermal state
and phase transition need to be considered during the simulation. Appendix A. Supplementary data
Therefore, a temperature-based energy equation is derived in this paper.
In order to account for the effects of the phase transition, a new mass Supplementary data to this article can be found online at https://doi.
transfer model is developed in terms of total mass flow. Afterwards, a org/10.1016/j.petrol.2019.106609.
complete algorithm is proposed, which consists of pressure and velocity
correction, volume fraction solution, temperature calculation and References
property parameters update. Conclusions are drawn as follows:
Abili, N., Kara, F., 2014. Compositional fluid tracking: an optimised approach to subsea
fluid sampling. Int. J. Oil Gas Coal Technol. 8 https://doi.org/10.1504/
(i) The temperature-based energy equation proposed in this paper is
IJOGCT.2014.064421.
more convenient in simulation. Compared with OLGA, it achieves Andreussi, P., Persen, L.N., 1987. Stratified gas–liquid flow downwardly inclined pipes.
almost the same outcomes. Int. J. Multiph. Flow 13, 565–575. https://doi.org/10.1016/0301-9322(87)90022-
(ii) The mass transfer model derived in this paper could well simulate X.
Andritsos, N., Hanratty, T.J., 1987. Influence of interfacial waves in stratified gas–liquid
CBC cases. For VBC cases, it could provide a more consistent flows. J. AIChE 33, 444–454. https://doi.org/10.1002/aic.690330310.
result. Bendiksen, K.H., Maines, D., Moe, R., Nuland, S., 1991. The dynamic two-fluid model
(iii) The algorithm developed in this paper is effective for both low OLGA: theory and application. Soc. Pet. Eng. J. https://doi.org/10.2118/19451-PA.
Bestion, D., 1990. The physical closure laws in the CATHARE code. Nucl. Eng. Des. 124,
liquid fraction pipe and high liquid fraction pipe. 229–245. https://doi.org/10.1016/0029-5493(90)90294-8.
Bonizzi, M., Andreussi, P., Banerjee, S., 2009. Flow regime independent, high resolution
Acknowledgments multi-field modelling of near-horizontal gas–liquid flows in pipelines. Int. J. Multiph.
Flow 35, 34–46. https://doi.org/10.1016/j.ijmultiphaseflow.2008.09.001.
Chen, J.F., Chen, G.Z., Wang, J.X., Shao, L., Li, P.F., 2011. High-throughput microporous
This work was supported by the National Natural Science Foundation tube-in-tube microreactor as novel gas–liquid contactor: mass transfer study. AIChE
of China (51874323 and 51534007), and the National Major Science J. 57 (1), 239–249. https://doi.org/10.1016/j.ces.2006.12.057.
Cohen, L.S., Hanratty, T.J., 1968. Effect of waves at a gas–liquid interface on a turbulent
and Technology Projects of China (2016ZX05028004-001, air flow. J. Fluid Mech. 31, 467. https://doi.org/10.1017/S0022112068000285.
2016ZX05066005-001 and 2016ZX05037005), all of which are grate­ Danielson, T.J., Bansal, K.M., Hansen, R., Leporcher, E., 2005. LEDA: the Next
fully acknowledged. Multiphase Flow Performance Simulator, 12th International Conference on
Multiphase Production Technology. BHR Group, Barcelona, Spain, p. 16.
Drew, D., Cheng, L., Lahey, R.T., 1979. The analysis of virtual mass effects in two-phase
Nomenclature flow. Int. J. Multiph. Flow 5, 233–242. https://doi.org/10.1016/0301-9322(79)
90023-5.
Evje, S., Flåtten, T., 2007. On the wave structure of two-phase flow models. SIAM J.
A Pipe cross-sectional area, m2
Appl. Math. 67, 487–511. https://doi.org/10.1137/050633482.
Cp Specific heat capacity at constant pressure, J/(kg⋅K) Goldszal, A., et al., 2007. Ledaflow 1D: Simulation Results with Multiphase Gas/
D Diameter, m Condensate and Oil/Gas Field Data, 13th International Conference on Multiphase
E Internal energy, J/kg Production Technology. BHR Group, Edinburgh, UK, p. 15.
Henriot, V., Duret, E., Heintz�e, Eric, Courbot, A., 2002. Multiphase production control:
f Friction coefficient application to slug flow. Oil Gas Sci. Technol. 57 https://doi.org/10.2516/ogst:
H Enthalpy, J/kg 2002006.

12
D. Fan et al. Journal of Petroleum Science and Engineering 185 (2020) 106609

Henriot, V., Pauchon, C., 1997. Tacite: Contribution of Fluid Composition Tracking on Raimondi, L., 2017. Compositional simulation of two-phase flows for pipeline
Transient Multiphase Flow Simulation, Offshore Technology Conference. Offshore depressurization. Soc. Pet. Eng. J. https://doi.org/10.2118/185169-PA.
Technology Conference, Houston, Texas, p. 10. https://doi.org/10.4043/8563-MS. Saurel, R., Abgrall, R., 1999. A multiphase godunov method for compressible multifluid
Hibiki, T., Ishii, M., 2003. One-dimensional drift-flux model and constitutive equations and multiphase flows. J. Comput. Phys. 150, 425–467. https://doi.org/10.1006/
for relative motion between phases in various two-phase flow regimes. Int. J. Heat jcph.1999.6187.
Mass Transf. 46, 4935–4948. https://doi.org/10.1016/S0017-9310(03)00322-3. Soave, G., 1972. Equilibrium constants from a modified Redlich-Kwong equation of state.
Ishii, M., 1975. Thermo-fluid Dynamic Theory of Two-phase Flow. Eyrolles, Paris. Chem. Eng. Sci. 27, 1197–1203. https://doi.org/10.1016/0009-2509(72)80096-4.
Krasnopolsky, B., Starostin, A., Osiptsov, A.A., 2016. Unified graph-based multi-fluid Stryjek, R., Vera, J.H., 1986. PRSV: an improved Peng—Robinson equation of state for
model for gas–liquid pipeline flows. Comput. Math. Appl. 72 (5), 1244–1262. pure compounds and mixtures. Can. J. Chem. Eng. 64 (2), 323–333. https://doi.org/
https://doi.org/10.1016/j.camwa.2016.06.020. 10.1002/cjce.5450640224.
Li, W.L., et al., 2018. CFD analysis of gas–liquid flow characteristics in a microporous Taitel, Y., Dukler, A.E., 1976. A model for predicting flow regime transitions in
tube-in-tube microchannel reactor. Comput. Fluids 170, 13–23. https://doi.org/ horizontal and nearly horizontal gas–liquid flow. J. AIChE 22, 47–55. https://doi.
10.1016/j.compfluid.2018.04.022. org/10.1002/aic.690220105.
Liapidevskii, V.Y., Tikhonov, V.S., 2019. Enhancement of the drift-flux model for gas- Tiselj, I., Petelin, S., 1997. Modelling of two-phase flow with second-order accurate
liquid slug flow in a long vertical pipe. Int. J. Multiph. Flow 50–58. https://doi.org/ scheme. J. Comput. Phys. 136, 503–521. https://doi.org/10.1006/jcph.1997.5778.
10.1016/j.ijmultiphaseflow.2018.08.009. Vernier, P., Delhaye, J.M., 1968. General two-phase flow equations applied to the
Morales-Ruiz, S., Rigola, J., Rodriguez, I., Oliva, A., 2012. Numerical resolution of the thermodynamics of boiling nuclear reactors. Rev. Energy Primaire 4, 1–43.
liquid–vapour two-phase flow by means of the two-fluid model and a pressure based Wallis, G., 1969. One Dimensional Two Phase Flow. McGraw-Hill, New York, NY, USA.
method. Int. J. Multiph. Flow 43, 118–130. https://doi.org/10.1016/j. Wang, Z., 2015. Online Simulation Study for Subsea Oil-Gas Gathering and
ijmultiphaseflow.2012.03.004. Transportation System (In Chinese): Beijing, China.
Moukalled, F., Darwish, M., 2004. Pressure-based algorithms for multifluid flow at all Xiao, J.J., Shonham, O., Brill, J.P., 1990. A comprehensive mechanistic model for two-
speeds—part i: mass conservation formulation. Numer. Heat Transf., Part B: phase flow in pipelines. In: SPE Annual Technical Conference and Exhibition. Society
Fundamentals 45 (6), 495–522. https://doi.org/10.1080/10407790490430651. of Petroleum Engineers, New Orleans, Louisiana, p. 14. https://doi.org/10.2118/
Moukalled, F., Darwish, M., 2004. Pressure-based algorithms for multifluid flow at all 20631-MS.
speeds—part ii: geometric conservation formulation. Numer. Heat Transf., Part B: Yue, J., Chen, G., Yuan, Q., Luo, L., Gonthier, Y., 2007. Hydrodynamics and mass transfer
Fundamentals 45 (6), 523–540. https://doi.org/10.1080/10407790490437997. characteristics in gas–liquid flow through a rectangular microchannel. Chem. Eng.
Patankar, S.V., Spalding, D.B., 1972. A calculation procedure for heat, mass and Sci. 62 (7), 2096–2108. https://doi.org/10.1016/j.ces.2006.12.057.
momentum transfer in three-dimensional parabolic flows. Int. J. Heat Mass Transf. Zhou, J., Podowski, M.Z., 2001. Modeling and analysis of hydrodynamic instabilities in
15, 1787–1806. two-phase flow using two-fluid model. Nucl. Eng. Des. 204, 129–142. https://doi.
Patel, N.C., Teja, A.S., 1982. A new cubic equation of state for fluids and fluid mixtures. org/10.1016/S0029-5493(00)00364-2.
Chem. Eng. Sci. 37, 463–473. https://doi.org/10.1016/0009-2509(82)80099-7. Zuber, N., Findlay, J.A., 1965. Average volumetric concentration in two-phase flow
Pauchon, C., Dhulesia, H., Lopez, D., Fabre, J., 1993. TACITE: a comprehensive systems. J. Heat Transf. 87, 453–468. https://doi.org/10.1115/1.3689137.
mechanistic model for two-phase flow. In: Proceedings of the BHRG Conference on
Multiphase Production, Cannes, France.

13