Journal of Petroleum Science and Engineering 205 (2021) 108923

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering

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

A fully implicit EOS based compositional two-phase transient flow

simulator in wellbores
Júlio César S. Nascimento a, Adriano dos Santos b, Carlos E. Pico Ortiz c, Adolfo Puime Pires c, *
a
Department of Materials Science and Technology, Universidade Federal da Bahia, Rua Professor Aristides Novis, 02, Federação, CEP:40210-630, Salvador, BA, Brazil
b
Petroleum Engineering Departament, Universidade Federal do Rio Grande do Norte, Campus Universitário, Lagoa Nova, CEP:59079-970, Natal, RN, Brazil
c
Laboratory of Petroleum Engineering and Exploration, Universidade Estadual do Norte Fluminense, CEP:27925-535, Macaé, RJ, Brazil

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

Keywords: Transient two-phase flow in pipes and wellbores appears in several processes in petroleum industry. Predicting
Two-phase transient flow pressure, temperature and composition behavior during offshore well tests can support investment decisions in
Two-fluid model exploration. The governing system of equations that model production (or injection) and build-up (or fall off) is
Two-phase segregation
composed by the mass conservation of each component, the transport equation for each phase and one overall
Fully implicit solution
Compositional flow
system energy equation. In this paper we develop a fully implicit equation of state based thermal compositional
simulator able to model two-phase multicomponent flow in wellbores and pipes. This simulator is fully
compositional and includes thermal effects. The thermodynamics calculations are based on the Peng-Robinson
equation of state, and the momentum equations are based on the Single-Pressure Two-Fluid model. The finite
volume method is used to discretize the non-linear system of equations, and time is discretized using Euler’s
implicit technique. The obtained numerical solutions for both co-current as well as countercurrent flow showed
excellent agreement when compared to benchmark solutions published in the literature. A well kick-off followed
by constant flow rate production and posterior shut-in was also simulated and the numerical results were dis­
cussed. It also models surge pressure and gravitational segregation effects in the wellbore. The presented
simulator may be applied in numerical pressure transient analysis and pressure behavior in flowlines connecting
wellheads and production platforms.

1. Introduction Several averaged models are available in literature to describe

multiphase flow in wellbores and pipelines, from a simplified drift-flux
Efforts for the development of numerical simulators for transient model (DFM) (Zuber and Findlay, 1965; Wallis, 1969) to a more
multiphase flow in pipes and wellbores accounting for heat and mass comprehensive two-fluid model (TFM) (Ishii, 1975; Ishii and Mishima,
transfer have grown. Modeling wellbore pressure and temperature 1984). The drift flux model considers only one momentum equation for
profiles is essential for well and topside facilities design and reservoir the mixture, and the coupling between phases is defined through an
performance forecast. Moreover, modeling transient pressure and tem­ expression based on two empirical parameters: the drift velocity and the
perature during production (injection) and build-up (fall-off) periods is distribution parameter. On the other hand, the two-fluid model uses a
important for numerical pressure transient analysis. Pressure waves momentum equation for each phase. The phase momentum transport
oscillation, shock waves propagation, phase appearance and disap­ equations are coupled via interfacial transport terms. Several models
pearance, and countercurrent flow with phases segregation are some of were proposed in the literature for the interfacial transport closure
the effects that occur during pressure transient tests (Almehaideb et al., equations of the two-fluid model, which are basically dependent on the
1989; Pourafshary et al., 2009; Han et al., 2013). Despite of the advances pipe geometry, fluid properties and the flow regime. Thus, two-fluid
in this research topic, the complex mathematical and numerical model presents a more detailed and realistic representation of the phe­
modeling of all the physical phenomena involved are still an open nomena. However, the two-fluid model is a more complex formulation
problem (Shirdel and Sepehrnoori, 2017; Raimondi, 2017). and computationally expensive (Ishii and Hibiki, 2006).

* Corresponding author.
E-mail addresses: jcsnascimento@ufba.br (J.C.S. Nascimento), adriano@eq.ufrn.br (A. Santos), capico@lenep.uenf.br (C.E. Pico Ortiz), puime@lenep.uenf.br,
adolfo.puime@gmail.com (A.P. Pires).

https://doi.org/10.1016/j.petrol.2021.108923
Received 22 November 2020; Received in revised form 4 May 2021; Accepted 5 May 2021
Available online 15 May 2021
0920-4105/© 2021 Elsevier B.V. All rights reserved.
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

There are two main approaches for the TFM formulation: the Two- state (Peng and Robinson, 1976). Then, transport equations for pressure,
Pressure Model (TPM) and the Single Pressure Model (SPM). If the temperature, velocities and volume fractions are solved. Flash calcula­
pressure of both phases is the same (SPM model), the problem becomes a tions are performed to determine the mole fractions in the oil and gas
mixed elliptic-hyperbolic type, due to the existence of regions where the phases from a given constant overall composition. Nevertheless, when
eigenvalues are complex. In this case, the problem is ill-posed, and small dealing with complex hydrocarbon mixtures in miscible floodings, or in
oscillations during the numerical iterations may generate instabilities the presence of high concentration contaminants (carbon dioxide,
and divergence issues (Ransom and Hicks, 1984; Dinh et al., 2003; hydrogen sulfide, etc.), the change in fluid composition must be taken
Prosperetti and Tryggvason, 2007). A Seven Equation Model was pro­ into account. Thus, it is necessary to solve a fully compositional equation
posed to avoid the aforementioned numerical instabilities and improve of state based system of equations, in which both the transport and
accuracy (Chao et al., 2018). It is based on the Two Pressure Two Fluid equilibrium equations must be satisfied simultaneously. In spite of the
Model and compares several reference cases published in the literature higher computational efforts necessary to solve the compositional
with good results. A similar approach using Godunov’s technique has formulation coupled with an equation of state, the results will be much
also been published (Ambroso et al., 2012). The main focus of these more reliable (Pourafshary et al., 2009; Forouzanfar et al., 2015).
works is the evaluation of the numerical aspects of the Seven Equation The flow pattern transitions during two phase flow may create dis­
Model for simple cases and do not present applications of the method to continuities in the momentum equations. To overcome convergence and
industrial problems. Due to the lack of realistic constitutive relations for stability issues due to these discontinuities, drift-flux model together
the interface momentum transfer, the SPM has been widely adopted for with the Shi et al. (2005) correlations have been used for the develop­
engineering applications (Städtke, 2006). Virtual mass and interface ment of fully implicit compositional wellbore simulators (Pourafshary
pressure regularization terms have been developed to eliminate the et al., 2009; Livescu et al., 2009; Forouzanfar et al., 2015; Xiong et al.,
elliptic regions that may appear in this model (Drew and Lahey, 1979; 2015; Shahamiri et al., 2015; Putcha and Ertekin, 2017). These models
Städtke, 2006; Bestion, 1990). have been applied successfully to solve co-current multi-phase flow
The two-fluid model has been analyzed from different aspects, both problems, but focused on the simulation of steady state and slow tran­
from the theoretical (Ishii, 1975; Delhaye and Achard, 1976; Drew and sients of gas-oil production.
Lahey Jr, 1979; Ishii and Mishima, 1984; Städtke, 2006), as well as from Although the drift-flux model allows slippage between phases, it
the numerical point of view (Suarel and Abgrall, 1999; Evje and Flåtten, leads to incorrect interphase momentum transfer in separated (stratified
2003; Garcia-Cascales and Paillère, 2006; Munkejord et al., 2009; and annular), fast transient and counter-current flows, in which the
Shekari and Hajidavalloo, 2013; Nascimento et al., 2018). This model phases are weakly coupled (Ishii and Hibiki, 2006). According to Ishii
has been extensively used in the development of different applications in and Hibiki (2006) the two-fluid model is more suitable than the
nuclear (Ransom, 1982; Micaelli, 1987; Borkowski and Wade, 1992) and drift-flux for the above-mentioned applications. Recently, Raimondi
petroleum enginering (Almehaideb et al., 1989; Stone et al., 1989; (2017) developed a semi-implicit compositional two-phase flow using
Bendiksen et al., 1991; Bonizzi and Issa, 2003; Issa and Kempf, 2003; the two-fluid model (SPM), which is able to simulate fast transients
Shirdel and Sepehrnoori, 2012, 2017, 2017; Krasnopolsky and Lukya­ flows. However, this model was only applied to analyze oil and gas
nov, 2018; Sondermann et al., 2019). depressurization in pipelines.
Several wellbore models using the two-fluid (SPM) approach coupled In this paper we present a fully implicit one-dimensional thermal
with a reservoir simulator have been developed. These simulators solve compositional two-phase flow simulator accounting for counter-current
the transport equations in the reservoir and wellbore, simultaneously, flow and gravitational segregation in wellbores using the single pressure
fully-implicitly (Stone et al., 1989; Almehaideb et al., 1989) or two-fluid model (SPM). The non-linear system of equations is discretized
semi-implicitly (Shirdel and Sepehrnoori, 2012, 2017). Stone et al. through the finite volume method, considering staggered grids and the
(1989) developed a three-phase (oil, water and gas) flow model in order first order upwind scheme. Time is discretized using Euler’s implicit
to simulate steam-assisted gravity drainage (SAGD) recovery process. In technique. The thermodynamic equilibrium is solved simultaneously
their wellbore model, five conservation equations (mass and momentum with the transport equations through the Newton-Raphson iterative
for liquid and gas, and energy for mixture) were considered. Further­ method.
more, the liquid phase consists of a mixture of oil and water with no In order to avoid discontinuity issues in the momentum equations, a
slippage. Almehaideb et al. (1989) proposed an isothermal three phase simplified flow pattern map is applied (RELAP5). Robust constitutive
flow model to investigate the effect of phase segregation on pressure relations are adopted, such as the used by RELAP5 - Reactor Excursion
responses during a build-up test. Assuming simplified constitutive and Leak Analysis Program (RELAP5) and TRACE - TRAC/RELAP
equations for the wall shear stress and interfacial friction, the afore­ Advanced Computational Engine (US-NRC, 2007), to calculate the wall
mentioned authors showed that if phase segregation occurs, pressure shear stress and interfacial friction, respectively. The obtained numeri­
increases quickly due to the increase of the hydrostatic column and gas cal solutions were compared to reference solutions for benchmark phase
compression at the bottom and at the top of the wellbore, respectively. segregation problems. Finally, a well kick-off followed by constant flow
Shirdel and Sepehrnoori (2012, 2017) improved the previous models by rate production and posterior shut-in was simulated. It is important to
introducing more comprehensive constitutive equations for interfacial highlight that to the best of our knowledge, the proposed fully-implicit
momentum transfer and phase-behavior calculation. Moreover, a slip­ solution and the obtained results for phase segregation considering the
page between oil and water modeled through a drift-flux approach was variation of temperature and fluid composition along the wellbore are
suggested (Shirdel and Sepehrnoori, 2017). not available in the literature.
Fluids properties also play an important role in the calculation of
pressure, temperature, saturation and concentration profiles during 2. Governing equations
transient two-phase flow in pipes. Many of the developed transient
wellbore models (TFM and DFM) are based on the black-oil model In this section, the governing equations for two-phase (liquid (l) and
(Stone et al., 1989; Almehaideb et al., 1989; Livescu et al., 2010; Pan vapour (g)) flow of a mixture consisting of Nc components are discussed.
and Oldenburg, 2014) or the pseudo-compositional model (Shirdel and The thermal compositional 2-phase flow model in pipes is composed by
Sepehrnoori, 2012, 2017). Both approaches are similar, but differences the mass conservation equation of each component, the momentum
appear in the phase behavior calculation. In the pseudo-compositional conservation equation of each phase, and the overall energy equation for
approach (Shirdel and Sepehrnoori, 2012, 2017) thermodynamic the 2-phase system. Assuming instantaneous equilibrium between pha­
properties (e.g., density, enthalpy, viscosity, gas solubility and forma­ ses, the mass balance equation for component c is expressed as:
tion volume factor) are calculated using the Peng-Robinson equation of

2
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

( Np
) ( Np
)
∂ ∑ ∂ ∑
(1) ∑
Nc
αp xc,p ξp + αp xc,p ξp vp = qnc,p ,
∂t p=1
∂z p=1
xc,p = 1, p = l and g (8)
c

for c = 1, …, Nc and p = 1, …, Np , where αp , ξp , vp are the volume

fraction, molar density and velocity of phase p, respectively; xc,p and qnc,p
Np
∑
αp = 1. (9)
are the mole fraction and internal mass source term of component c in p
phase p, respectively; Nc is the number of components, and Np denotes
In summary, the resulting Np × (Nc +2) + 2 independent equations
the total number of phases. The internal mass source term of component
c in phase p given by (differential and algebraic constraints, Eqs. (1)–(9)) are solved in order
to determine the unknowns: pressure (P), temperature (T), Np × Nc mole
Np
∑ ( ) fractions (xc,p ), Np velocities (vp ) and Np volume fractions (αp ).
qnc,p = − W λp xc,p ξp PRp − P , (2)
The constitutive equations necessary for obtaining a closed system of
p
equations are presented in the next subsections.
where W is the well index, λp is the phase mobility, and P and Pp R are the
area averaged pressure in the wellbore and the phase pressure in the grid
2.1. Interfacial pressure
cell of a coupled reservoir simulator in which the well segment is
located, respectively.
In this work, the following correlation for interfacial pressure ΔPi is
The momentum equations for each phase using the single pressure
assumed (Bestion, 1990):
two-fluid model approach are given by
α g α l ρg ρl ( )2
∂( ) ∂( ) ∂P ∂α ΔPi = (P − Pi ) = σ ( ) vg − vl . (10)
αρv + α ρ v v = − αp − αp ρp g sin θ − τwp − ΔPi p + MpD , αg ρl + αl ρg
∂t p p p ∂z p p p p ∂z ∂z
(3) The aforementioned expression, which has been used to model non-
for liquid (l) and vapour (g) phases, where P is the area averaged stratified flow pattern in CATHARE code (Barre and Bernard, 1990), is
pressure; τwp is the wall shear stress; ΔPi (the pressure correction term) used to guarantee the hyperbolicity of the system (1)–(9). For this
purpose, Paillère et al. (2003), Evje and Flåtten (2003, 2005, 2006),
represents the average interfacial pressure, and MDp is the drag force at
Shekari and Hajidavalloo (2013) and Yeom and Chang (2011) proposed
gas-liquid interface per unit of volume. In the right hand side of Eq. (3),
1 < σ ≤ 2.0, and we adopted σ = 1.2.
the last three terms are flow pattern dependent and are evaluated by
means of semi-empirical correlations.
The hypothesis of identical temperature T(x, t) for all phases at each 2.2. Interfacial momentum transfer
position results in the following total energy equation:
[N ( ) ] [N ( ) ] The interfacial momentum transfer is the most important term in the
∂ ∑ p
1 2 ∂ ∑ p
ρp v2p two-fluid formulation because it controls the degree of mechanical
αp ξp hp + ρp vp − P + α p ξp h p + vp =
∂t p 2 ∂z p 2 coupling between phases (Ishii and Hibiki, 2006). Therefore, an accu­
(4) rate correlation is required in order to guarantee the effectiveness of the
Np
∑ numerical model. Generally, only interfacial shear forces are taken into
− αp ρp vp g sin θ − Qw + Qh , account in separated flows (stratified and annular flow), and only drag
p
forces are taken into account in dispersed (bubbly, slug and mist) flows.
The models available in the literature describing the interactions among
where hp is the molar enthalpy of phase p and Qh is the internal energy
phases were chosen aiming flexibility and general use. In this paper, the
source term:
approach proposed in RELAP5 was chosen to model the interfacial drag
Np
∑ ( ) force:
Qh = − W hp λp ξp PRp − P . (5)
p MgD = − Ci vr |vr |, (11)

In addition, the heat transfer between the wellbore and the sur­
where Ci and vr are the interfacial drag coefficient and relative velocity,
rounding formations (Qw ) is expressed in terms of the overall heat
respectively. For liquid phase MDl = − MDg .
transfer coefficient Ut as follows:
( ) For dispersed flows, Ci and vr are computed by using the drift-flux
Qw = 2πrto ΔzUt Tf − Te . (6) approach (RELAP5), as follows:
( )
The overall heat transfer coefficient depends on the wellbore features αg α3l g ρl − ρg
(casing and tubing thermal conductivities), formation (thermal con­ Ci = , (12)
v2D
ductivity) and flowing fluids properties (convection coefficients). In
( )
addition, Tf is the flowing fluid temperature, Te is the surrounding for­ 1 − C0 αg
vr = vg − C0 vl , (13)
mation temperature, rto and Δz are the external tubing radius and cell 1 − αg
length, respectively. Furthermore, we assume that Ut is a given constant
or determined by using Hasan and Kabir’s model (Hasan and Kabir, where vD and C0 are the drift-velocity and distribution parameters,
1994). which are calculated by using Kataoka and Ishii’s correlations (Kataoka
Additional closure equations include phase equilibrium conditions, and Ishii, 1987).
given by: For mist flow, Ci is given by

fc,l = fc,g , (7) 1

C i = C D Ai ρf , (14)
8
where fc,l and fc,g are the fugacities of component c in the liquid and gas
phase, respectively. Finally, the following molar and volume fractions where f indicates the continuous phase (f = l for bubbly and f = g for
constraints are required: mist flow), CD is the drag coefficient; and Ai is the interfacial area, which
can be written as

3
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

3.6αf wall for the bubbly, slug, transition slug/annular and annular flows (i.e,
Ai = , (15)
dd the gas shear stress is set to zero). Finally, for the slug/annular transition
regime, the shear stress is calculated by linear interpolation, as sug­
where dd is the diameter of dispersed phase. In addition, the following gested by TRACE (US-NRC, 2007).
drag coefficient (CD ) and diameter of dispersed phase (dd ) are consid­
ered: CD = 0.44 and dd = 1.0 mm. 2.4. Two-phase flow map
For slug flow the interfacial drag coefficient Ci is obtained by the
contribution of the drag force in two regions: liquid slug region and As shown above, the constitutive equations for interfacial mo­
Taylor bubble region. The complete description of the calculation of Ci mentum transfer and wall shear stress depend on the definition of the
for the slug flow can be found in RELAP5. flow pattern, which indicates how the fluid distribution affect the
In annular flow the gas phase flows at the center of the pipe (gas coupling between phases. Usually, flow patterns are identified from a
core), which may contain liquid droplets, while the liquid is pushed to flow pattern map obtained from experimental observations or theoret­
the pipe wall and flows as a thin film (Shoham, 2006). The interfacial ical models. For example, the maps described in Shoham (2006) are
friction term, MDg , includes the shear stresses between the liquid film and based on the physical modeling of the specific (assumed) underlying
the gas and the drag force between trapped liquid droplets and the gas transition mechanisms. Despite the extensive experimental verification
core. and agreement with different mechanistic approaches and empirical
Neglecting liquid droplets dispersed in the gas phase, MDg can be correlations, those types of flow pattern maps are too complex to be
written as included in a fully implicit numerical solution method. From the point of
view of a numerical simulator for the two fluid model, it is more
MgD = − τi Ai , (16) convenient and consistent to use transition criteria based on a geometric
parameter than those based on volumetric fluxes (Ishii and Mishima,
where τi is the interfacial shear stress, defined as 1984). The chosen formulation has three major transitions: bubbles to
1 ( )⃒( )⃒ churn/slug (BS), churn/slug to annular/mist (SA) and annular/mist to
τi = fi ρg vg − vl ⃒ vg − vl ⃒. (17) mist (AM). These transitions were determined using the volumetric
2
fraction as a criterion and an adaptation of the models presented in
and fi is the interfacial friction factor obtained from Wallis (1969) RELAP5 (RELAP5). This map can be applied to any pipe inclination
model, given by angle between 60 and 90◦ upward, downward and for counter-current
[ ( )] flows.
fi = 0.005 1 + 75 1 − αg . (18)

In addition, if the volume fraction of the liquid droplets in the gas 3. Numerical method
core is much smaller than gas volume fraction (αl,d ≪αg ), the interfacial
area Ai for annular flow can be written as: The system of governing equations (Eqs. (1)–(9)) is discretized using
√̅̅̅̅̅ the Finite Volume Method (Patankar, 1980; Versteeg and Malalasekera,
αg
Ai = 4 . (19) 2007), as shown in Fig. 1. The pipe length is divided into N computa­
D
tional cells of uniform size Δzk ≡ zk− 12 − zk+12 , where z is the spatial co­
Finally, the interfacial friction, MDg , for the slug/annular transition ordinate, subscript k denotes cell index, and Ak and θk are the cross
flow pattern is obtained by linear interpolation, as follows: sectional area and pipe inclination with horizontal direction, respec­
( ) ( ) tively. In order to avoid the commonly observed oscillations on pressure
αg,A − αg αg − αg,SA
D
Mg,SA = D
Mg,S + D
Mg,SA , (20) when co-located grids are used (Prosperetti and Tryggvason, 2007;
αg,A − αg,SA αg,A − αg,SA
Versteeg and Malalasekera, 2007), a forward staggered-grid scheme is
where the subscripts S, SA and A represent slug, transition slug/annular adopted. In this grid, the scalar variables such as pressure, temperature
and annular flow regimes, respectively. and volume fractions are stored and solved at cells centers zk , while
velocities are stored and computed at grid cells faces zk+12 . The mass and
The liquid interfacial momentum transfer is computed as MDl = −
MDg . energy equations are spatially integrated across the main grid volumes,
while the momentum equations are integrated across volumes of the
staggered grid (see Fig. 1).
2.3. Wall shear stress To prevent stability problems due to the strong nonlinearity of the
system (1)–(9), time integration was performed using a fully implicit
Due to its simplicity, to the weak dependence on the flow regime and first-order backward Euler scheme (Chen, 2007). Consequently, the
to the similarity with the single-phase expression, an empirical model of proposed numerical solution is unconditionally stable allowing large
the wall friction force, in the form of a two-phase multiplier, was chosen time-steps (i.e., the time-step is not restricted by CFL condition in the
to compute the wall friction pressure gradient. So, the wall shear stress is form of Δx Δt
max(λp ) ≤ 1, where λp are the eigenvalues associated with
modeled according to the TRACE (US-NRC, 2007) approach: wave speeds) (LeVeque, 2004; Prosperetti and Tryggvason, 2007).
( )2 ⃒ ⃒ It is important to highlight that transport of mass and energy (scalar
( )2 αp fp ρp vp ⃒vp ⃒
τwp = Φp , (21) properties) at the faces of the main control volume (zk±12 ), and mo­
2D
mentum at faces of the staggered grid (zk,k+1 ) need to be interpolated
where D is the hydraulic diameter and fp is the Darcy-Weisbach friction using their values at centers and edges of the main grid cell, respectively
factor computed using the correlation developed by Zigrang and Syl­ (see Fig. 1). In order to avoid oscillations and obtain smooth solutions
vester (1982). In Eq. (21), Φp is the two-phase multiplier (Lockhart and near discontinuities, a first-order upwind scheme was considered.
Martinelli, 1949), defined according to the flow pattern as: The nonlinear algebraic system resulting from the discretized equa­
⎧ tions written in residual form (see Appendix A) for all grid cells is solved
⎨ α− 1.8 , if ​ bubbly, ​ churn ​ or ​ slug simultaneously by applying the Newton-Raphson method, as follows:
(22)
p
Φp =
⎩ α−p 2.0 , if ​ annular ​ or ​ mist.
δX(ν+1) = − [J(Xv )]− 1 R(Xv ), X(ν+1) = X(ν) + δX(ν+1) , ν = 0, 1, …. (23)
It is considered that only the liquid phase is in contact with the pipe

4
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Fig. 2. Algorithm to solve one time step using Newton-Raphson linearization.

Fig. 1. Schematic of one-dimensional control volume and forward staggered

grid scheme.

where ν is the Newton iteration, R is the residual vector, δX is the vector

of increments of the unknown variables X, and J is the Jacobian matrix,
defined as
[ ](ν)
∂R
J(ν) = , (24)
∂X nR ×nX

where nR is the number of residual equations. The unknown vector is

obtained by iterating system (23) until the maximum norm of residual
vector R(ν+1) and increment vector δX(ν+1) , defined as:
⎧
⎪ max|Rν | < ε, and
⎨
⃒ ⃒
⃒ δX ν+1 ⃒ (25)
⎩ max⃒⃒
⎪ ( )⃒ < ε,
max 1, Xν+1 ⃒

reach a specified tolerance ε. The linear system (23) is solved using

PARDISO (MKL, Intel Fortran). Finally, in order to improve accuracy
and performance, the Jacobian matrix (J) is calculated analytically.
A schematic representation of the algorithm for one time step inte­
gration in the proposed model is shown in Fig. 2. For the first time step,
the initial condition is used as a guess for the unknown vector. Other­
Fig. 3. Algorithm to perform a phase appearance and disappearance test.
wise, the solution of the previous time step is used for the unknown
vector at first Newton iteration. Then, the stability/flash calculations
(Fig. 3) are performed for every grid cell at each time step in order to achieved within the specified maximum number of Newton iterations,
compute the number of phases. the time step size is divided by two and the process is repeated until
According to the number of phases present in a grid cell, the set of convergence. In this paper, the tolerance ε = 10− 5 and a maximum
unknowns is chosen for each cell and the residual vector and Jacobian number of Newton iterations equal to 50 are adopted.
matrix are formed. Then, the linear system (Eq. (23)) is solved and the Each grid cell in the domain may contain up to 2Nc + 6 primitive
unknown X(ν+1) vector is updated. If the specified tolerance is not variables if there are two phases present (P; T; vp , αp and xc,p , for liquid

5
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

and vapour) or Nc + 3 primitive variables if the system is single-phase,

liquid or gas (P; T; vp and xc,p , for liquid or gas). The number of pha­ 4. Numerical tests
ses is set as 1; then a stability test is performed followed by a two-phase
Different scenarios of two-phase flow in vertical wells are investi­
flash calculation when the original single phase is unstable. If a phase
gated in this section. Both co-current and countercurrent flow, consid­
disappears in an originally two-phase cell, a recombination procedure is
ering phase appearance and disappearance, were studied in the
performed. In this work we follow the approach presented in Cao (2002)
aforementioned scenarios. The obtained numerical solutions, which
to determine if a phase disappears during the iterative solution process
were compared to solutions available in the literature, are presented in
(Fig. 3).
the next subsections.
After every iteration the appearance/disappearance of a phase is
checked. In a two-phase grid cell, if one of the volumetric fractions
become negative the phase condition of this grid cell is updated to single 4.1. Case 1: co-current two-phase flow
phase. In such cases, the overall composition of the remaining phase is
recalculated recombining the composition and volumetric fractions of Three scenarios proposed by Pourafshary et al. (2009), here named
the existing phases in the previous iteration. In a single phase grid cell, a as Case 1A, 1B and 1C, were studied. The tubing data is presented in
stability test is performed using the updated pressure, temperature and Table 1, and the fluid composition used in each simulation is shown in
composition. If this phase is thermodynamically unstable, a new phase Table 2. Different fluids were considered: black-oil, volatile oil and
appears, and a flash calculation is performed to determine the compo­ retrograde gas for cases 1A, 1B and 1C, respectively. The boundary
sition of the new system. conditions for all three cases are: well head constant pressure (see
The residual vector R is given by Table 1); total bottom hole mass flow rate equal to 680.39 kg-moles/day
[ ]T (1500 lb-moles/day); and constant bottom hole global fluid composition
→ → → →
R = R 1 , R 2 , …, R k , …, R N , (26) (see Table 1).
We considered that the production pipe is initially filled up with a
mixture (liquid and gas) at rest (vg = vl = 0.0 m/s). Besides that, it is also
→
where R k denotes the vector of residual equations for the grid cell k and considered that the fluid is in thermal equilibrium with the surrounding
→
N is the total number of grid cells. For the k − th grid cell, R k , is defined rock. The phase (liquid or gas) and composition at each depth was
as determined through a flash calculation using the global composition and
the hydrostatic equilibrium.

⎧⎡ ⎤T
⎪
⎪ [ ]N [ ] [ ]Np.k
⎪
⎪⎢ c 1
αp ⎥
⎪
/
Nc h P p,k+ 2 fc fNc xc,p
⎪
⎪ ⎣Rk , ..., Rk , Rk , Rp,k , Rk , ..., Rk , Rp.k , Rk ⎦ , if ​ two − phase
⎪
⎪ p=1 p=1
⎪
⎪
⎨
→
Rk = ⎡ ⎤T (27)
⎪
⎪
⎪
⎪
⎪
⎪ ⎢ c [ ]Np,k+ xc,p ⎥
⎪ 1/
⎪
⎪ ⎢R , ..., RNk c , Rh , RP p=1 , Rp.k
2
⎥ , if ​ single − phase
⎪
⎪ ⎣ k k k ⎦
⎩

Fig. 4 compares the solution for the three cases using the proposed
where Np,k and Np,k+12 denote the number of phases at cells centers and solution and the steady-state solution of Pourafshary et al. (2009) and
edges, respectively. In order to reduce the computational cost of the Putcha and Ertekin (2017). The results shown in Fig. 4a, b and c were
stability/flash calculation at grid cell edges, Np,k+12 is obtained as follows determined using a grid containing 20, 40 and 100 cells, respectively;
⎧ the initial time step adopted for all cases was 2.5 s; the time step was
⎨ Np,k , if ​ Np,k = 2 and Np,k+1 = 1 doubled at every iteration up to a maximum value of 50 s. Note that our
Np,k+12 = Np,k+1 , if Np,k = 1 and Np,k+1 = 2 , (28) transient results converge to the published steady-state solutions.
⎩
Np,k , otherwise

Similarly, the unknown vector X is given by

[ ]T Table 1
→ → → Tubing data.
X = X 1 , X 2 , …, X nX , (29)
Parameter Case 1A Case 1B Case 1C

→ Depth (m) 1524.00 2438.40 3657.60

where X k denotes the unknown vector for the grid cell k. For the k− th Wellhead pressure (MPa) 2.068 6.895 10.342
→ Bottom-hole temperature (K) 338.7 363.7 349.8
grid cell, X k is defined as
Geothermal gradient (K/m) 0.026 0.010 0.012
⎧⎡ ⎤T
⎪ Inner tubing diameter (m) 0.076200
⎪
⎪⎢ [ ] [ ] Outer tubing diameter (m) 0.076384
⎪
⎪ [ ]Np,k+ ( 1 ) c p,k ( ) p,k ⎥
N ,N N
⎪
/

⎪⎢
⎪ ⎥
⎣Pk , Tk , vp,k p=1 , xc,p k c=1,p=1 , αp k p=1 ⎦ , if ​ two − phase Roughness (m) 1.82 × 10− 4
2

⎪
⎪
⎪
⎪ Heat conductivity (W/m.K) 60.55
⎨
→ Inner diameter (m) 0.08665
Xk = ⎡ ⎤T .
⎪
⎪ Outer diameter (m) 0.19226
⎪
⎪
⎪ Heat conductivity (W/m.K) 60.55
⎪
⎪ ⎢ [ ]Np,k+ [( ) ] Nc ⎥
⎪ Wellbore diameter (m) 0.17
1/

⎪
⎪ ⎢Pk , Tk , vp,k p=1 , xc,p k
2
⎥ , if ​ single − phase
⎪ ⎣ c=1 ⎦
⎪
⎩ Heat conductivity (W/m.K) 6.96
Heat conductivity (W/m.K) 2.42
(30) Heat diffusivity (m2/s) 1.3845 × 10− 6

6
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Table 2 along the pipe.

Initial and inlet fluid composition. A sketch of phase distributions along the tube during the initial,
Component Case 1A. Blackoil Case 1B. Volatile oil Case 1C. Retrograde gas transient and stationary states are shown in Fig. 5. When t > 0, liquid is
accelerated and flows downward while gas accelerates and flows in the
C1 0.3 0.55 0.8
opposite direction, leading to countercurrent flow. During the transient
C3 0.12 0.10 0.04
period, two interfaces (named interfaces A and B, see Fig. 5) separate
iC4 0.12 0.10 0.04
two single phase regions (liquid at the bottom and gas at the top) from
iC5 0.12 0.10 0.04
one intermediate two-phase region. When the steady state flow is
C7 0.17 0.075 0.04
reached, liquid and gas are completely segregated and their velocities
C8 0.17 0.075 0.04
are both zero.
At the beginning of the simulation a 2 m length vertical pipe is filled
Table 3 presents the average deviation between our proposed technique with a fluid at rest, and the pressure and temperature are assumed to be
and literature results. For all cases, they are smaller than 10%. constant along the tube (1.0 MPa and 325.15 K). Molar fraction of each
component and volumetric fraction of each phase are determined after
performing a series of flash calculations. The grid was composed by 400
4.2. Case 2: phase segregation cells and time step Δt was set equal to 5 × 10− 4 s. Fig. 6a, b and c show
that the numerical solution developed in this work and the reference
The benchmark proposed by Coquel et al. (1997) analyzes the fric­ solution present close agreement. It is important to point out that,
tionless immiscible and isothermal two-phase flow in a closed vertical regardless of the numerical diffusion introduced by the upwind scheme,
pipe. In this subsection we compare the results of an isothermal, small the discontinuities positions are correctly predicted by the proposed
miscibility two-phase flow, composed by 30% methane (C1 ) and 70% solution scheme.
decane (C10 ), with the analytical results published in Nascimento et al. The numerical spikes observed in the superficial velocities (see
(2018). It is important to point out that in this problem, the counter­ Fig. 6b and c), which occur near the region of phase disappearance, are
current flow will be considered, in which phases appear and disappear due to discontinuities in the volumetric fraction inside the staggered

Fig. 4. Comparison between the transient pressure profiles and the steady state solutions proposed by Pourafshary et al. (2009) and Putcha and Ertekin (2017): (a)
Case 1A black-oil; (b) Case 1B volatile oil; (c) Case 1C retrograde gas.

7
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Table 3
Relative deviation (%) between proposed solution and literature data.
Case 1A Case 1B Case 1C

Depth Pourafshary et al. Putcha and Ertekin Depth Pourafshary et al. Putcha and Ertekin Depth Pourafshary et al. Putcha and Ertekin
(m) (2009) (2017) (m) (2009) (2017) (m) (2009) (2017)

152.40 0.792 5.284 243.84 0.800 1.896 365.76 1.296 1.117

304.80 0.366 6.498 487.68 0.085 2.568 731.52 1.717 1.934
457.20 0.154 7.747 731.52 0.531 3.199 1097.28 4.313 2.657
609.60 1.669 8.543 975.36 0.792 3.746 1463.04 6.671 3.419
762.00 0.659 8.787 1219.20 0.862 4.044 1828.80 8.851 4.271
914.40 2.203 8.461 1463.04 1.438 4.221 2194.56 10.379 5.242
1066.80 1.078 7.607 1706.88 0.150 4.083 2560.32 11.858 6.274
1219.20 0.748 6.386 1950.72 0.969 3.774 2926.08 13.033 7.419
1371.60 1.423 4.546 2194.56 1.446 3.245 3291.84 11.839 8.645
1524.00 0.612 2.955 2438.40 2.637 2.453 3657.60 10.002 10.115
Average 0.882 6.074 0.883 3.021 7.269 4.645

P(t) Z(t) ng (t) Vg0

= 0 , (31)
P0 Z n0g Vg (t)

where the superscript “0” indicates initial condition, Zg and Vg are the
gas compressibility factor and volume, respectively. Because ng and Vg
are constant for immiscible solutions, the gas region pressure is constant
and equal to the initial pressure (P(t) = P0 ). During miscible flow,
( )〈
ng (t) Vg0
because Z is constant and n0g Vg (t)
1, pressure decreases (Fig. 9).

4.3. Case 3: well kick-off, production and shut-in

In this section, kick-off, production, and shut-in periods for a vertical

well are simulated. The obtained results show that the proposed simu­
lator is able to predict pressure oscillatory effects during the kick-off
period, countercurrent flow, phase segregation and flow pattern tran­
sition. It is considered that initially there is only liquid hydrocarbons in
Fig. 5. Sketch of the phase segregation process. the wellbore at hydrostatic and thermal equilibrium. The boundary
conditions consist of prescribed: (a) wellhead pressure and (b) mass and
volume (see Fig. 1). Consider that the interface “A” (see Fig. 5), whose energy sources (as a function of pressure, temperature and produced
position is not predicted by the model (A.1)-(A.6), is located somewhere fluid composition, see Eq. (2)). Main simulation data and fluid proper­
below the face k + 12. In this case, the linear interpolation between (αpk ) ties are presented in Tables 4 and 5.
⎛ ⎞ When t = 0, the well is suddenly opened and the wellhead pressure is
⎜ ⎟
and (αpk+1 ) will produce an unreal ⎝αpk+1 ⎠ used in the momentum instantaneously changed from 7 MPa to the prescribed pressure (5 MPa).
2 When t = 6 hours, the wellbore is shut-in, and phase segregation begins.
Pressure and liquid flow rate at the center of the first cell at the
equation, resulting in excessive acceleration in the staggered sub- beginning of the production period are presented in Fig. 10. Note that
volume between the faces k and k + 12. For this reason, the magnitude immediately after the beginning of the production period pressure drops
of the aforementioned spikes is amplified if friction is neglected. abruptly. As a consequence, the fluid suddenly expands and both pres­
Fig. 7 shows the methane global composition along the pipe for sure and flow rate oscillate close to the wellhead. As the pressure os­
different times. As expected, when steady state develops, the gas at the cillations disappear due to friction, the liquid flow rate also stabilizes.
top of the pipe is composed by methane only; while at the bottom the Liquid velocity profile evolution indicates fluid depressurization at early
liquid composition is 88% decane and 12% methane. times (see Fig. 11).
Finally, Fig. 8 compares the steady state pressure solution developed Fig. 12 presents the pressure, temperature, liquid hold-up, and phase
in this work and the reference solution. Besides that, in order to stress velocity profiles for different times during the production period.
the importance of the miscibility in the solution, Fig. 8 also presents the Fig. 12a shows that the Well Head Pressure (WHP) equals the prescribed
immiscible case, where there is no mass transfer. Note that the pressure value (boundary condition) when t = 0.05s. At this point, the heavier
behavior is similar for both cases; nevertheless, the pipe pressure is fluid (initially in the wellbore) starts to be replaced by the lighter
smaller when mass transfer takes place (miscible solution). It is also reservoir fluid (see Table 5). As a consequence, bottomhole pressure
important to point out that the immiscible case solution matches the decreases (see Fig. 12a) leading to the release of gas at t = 100s
reference solution (also immiscible). (Fig. 12b). The gas released creates a discontinuity in the liquid hold-up
For miscible solutions, the number of mols and the gas phase volume profile that moves upward and reduces the mixture density and pressure
(ng (t) and Vg (t), respectively) decrease with respect to their initial values gradient until it reaches the top, establishing two-phase flow along the
n0g and Vg0 (Fig. 9). Considering constant temperature (T(t) /T 0 = 1), the wellbore (see Fig. 12b when t = 100, 400 and 600s). When t = 1000s
ratio P(t)/P0 can be found from the real gas law: steady state flow is fully developed. These results show that the pro­
posed model is able to calculate phase appearance during countercur­
rent flow in wellbores.
Liquid and gas superficial velocities are presented in Fig. 12c and d.

8
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Fig. 6. Proposed (dashed lines with markers) and reference (continuous lines) solutions for: (a) liquid volume fraction; (b) superficial liquid velocity and (c) su­
perficial gas velocity.

Fig. 8. Comparison of proposed numerical solution (miscible and immiscible)

Fig. 7. Overall methane composition calculated in this work. and reference solution for steady state pressure along the pipe.

9
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

compressibility, phase appearance at the bottom and the fixed reservoir

pressure boundary condition, both bottomhole pressure (BHP) and
liquid flow rate oscillate. These effects are no longer present when
steady state develops (Fig. 13b).
Fig. 14 presents WHP and BHP after the shut-in. Due to the well
known water hammer effect, pressure increases at the upper part of the
pipe. Depending on the pipe length and the mixture compressibility,
BHP may drop sharply (Nascimento et al., 2018). In this case, we
consider reservoir production (mass and energy source terms at the

V0 n
Fig. 9. PP0 , ZZ0 , Vgg and n0g relations at the top (x = 0, t) of the pipe, where super­
g
script denotes initial value.

Table 4
Wellbore configuration, discretization model and boundary conditions.
Wellbore configuration Value Numerical parameters

Depth (m) 1000.0 Number of grid 250

cells
Inner diameter (m) 0.10 Δt initial (s) 0.01
Roughness (m) 2.4 × 10− 4 Δt minimum (s) 10− 4

Inclination (◦ ) 90.0 Δt maximum (s) 10.0 Fig. 10. Pressure and liquid flow rate behavior at the beginning of produc­
Geothermal gradient (K/m) 0.02 Max NR iteration 25 tion period.
Bottom-hole temperature 333.15 Δt multiplier 2.0
Overall heat transfer coeficient (Ut ) 20.0 Simulation time 1
(W/m2.K) (days)
Boundary conditions
Wellhead pressure (MPa) 5.0
Reservoir pressure (MPa) 17.0
Reservoir temperature 333.15
Well index W (m3) 9.868 ×
10− 14
Oil mobility λo (1/Pa.s) 20.268 ×
103
Gas mobility λg (1/Pa.s) 10.134 ×
105

As expected, liquid velocity increases in the upward direction due to

both pressure decrease and gas interfacial drag. Due to the higher
temperature of the produced fluid, temperature increases at the bot­
tomhole (see Fig. 12e).
As shown in Fig. 13a, liquid depressurization in the wellbore is
responsible for most of the liquid production up to 20s. Due to fluid Fig. 11. Liquid velocity profile along tube depth.

Table 5
Fluid composition and component critical properties.
Components Initial wellbore fluid Produced fluid Critical pressure Critical temperature Molecular weight (g/ Critical volume (l/ Acentric
composition composition (atm) (K) mole) mole) factor

N2 1.67 × 10− 2
8.17 × 10− 2 33.500 126.200 28.013 0.090 0.040
C1 10.0 × 10− 2
32.5 × 10− 2 45.400 190.600 16.043 0.099 0.008
CO2 0.20 × 10− 2
0.20 × 10− 2 72.800 304.200 44.010 0.094 0.225
C2 7.53 × 10− 2
7.52 × 10− 2 48.200 305.400 30.070 0.148 0.098
C3 4.74 × 10− 2
4.74 × 10− 2 41.900 369.800 44.097 0.203 0.152
iC4 0.05 × 10− 4
0.05 × 10− 4 36.000 408.100 58.124 0.263 0.176
nC4 4.12 × 10− 2
4.12 × 10− 2 37.500 425.200 58.124 0.255 0.193
iC5 0.05 × 10− 4
0.05 × 10− 4 33.400 460.400 72.151 0.306 0.227
nC5 2.97 × 10− 2
2.97 × 10− 2 33.300 469.600 72.151 0.304 0.251
nC6 1.38 × 10− 2
1.38 × 10− 2 32.460 507.500 86.000 0.344 0.275
C7+ 67.3 × 10− 2
36.3 × 10− 2 27.563 583.849 115.685 0.328 0.450

10
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Fig. 12. Solution for the production period: (a) pressure; (b) liquid fraction; (c) liquid velocity; (d) gas velocity and (e) temperature.

bottomhole, see Eqs. (2) and (5)) even after the well is shut-in, until the afterflow period may last long. Otherwise, high liquid flow rates lead to
bottomhole pressure is equal to the original reservoir pressure. This short afterflow periods, but may generate pressure shocks (Han et al.,
phenomenon is known as “afterflow”, and is associated with wellbore 2013).
storage. If the fluids are compressible and the flow rate is small, the For a long shut-in period, gravitational segregation becomes very

11
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Fig. 13. Pressure and liquid flow rate versus time: (a) during startup process and (b) up to steady-state.

and two-phase countercurrent flow in the intermediate region. This

configuration lasts until steady state is reached. Fig. 15b shows the
different pressure gradient for each region separated by the interface
where liquid hold up solution is discontinuous.
Gravitational segregation causes the increase of bottomhole pressure
due to the gas accumulation at the top and to the hydrostatic column of
the liquid. This effect lasts until the steady state develops (Fig. 16a). In
addition, due to the afterflow effect liquid continues to accumulate in
the well up to 1000s (Fig. 16b). Comparison between the proposed and
reference solutions (Fair, 1981) is presented in Fig. 16a.
Liquid and gas superficial velocities are presented in Fig. 17a and b,
respectively. Note that, although the proposed model considers wall and
interfacial friction, the velocities at the interfaces that separate gas and
two-phase regions presented numerical spikes similar to those observed
in Fig. 6b and c. The reason for the appearence of these spikes is dis­
cussed in section 4.2.
Fig. 14. WHP and BHP profiles after shut-in. Temperature solution is presented in Fig. 18. Note that for t < 1000 s,
fluid temperature increases near the bottomhole due to fluid compres­
clear (Fig. 15). Note that, before complete segregation (time around sion and the warmer source fluid flowing in during afterflow. As ex­
1000s), liquid moves downward and accumulates in the lower part of pected, when afterflow effects are over, the temperature profile tends to
the wellbore, while gas moves upward (countercurrent flow). Besides the geothermal gradient. Spurious temperature oscillations appear at the
that, the solution presents three different regions: single-phase flow at gas-liquid interface region after the development of the steady state
the upper and lower regions of the pipe (gas and liquid, respectively) flow, probably caused by the velocity oscillation (Fig. 17a and b).

Fig. 15. Solution for the shut-in period: (a) liquid volume fraction and (b) pressure.

12
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Fig. 16. This work solution and reference solution (Fair, 1981) for BHP for shut-in period.

Fig. 17. Solution for the shut-in period: (a) liquid velocity and (b) gas velocity.

Finally, the global composition of the most representative compo­

nents at the top, at 750 m and at the bottomhole at initial conditions and
at t = 52000 s are presented in Fig. 19. Note that the fluid composition is
almost constant in the wellbore at the shut-in time. At 750 m part of the
nitrogen (N2 ) and of the methane (C1 ) condense, while at the top, due to
the low pressure and temperature, the gas phase is mostly composed by
light reservoir fluid fractions; such as methane, ethane and nitrogen.

5. Conclusions

A fully-implicit Equation of State based numerical model for thermal

compositional 1D two-phase flow in pipes using the single pressure two-
fluid model was developed. The governing equations were discretized by
applying the finite volume method, considering forward staggered grids
and upwind scheme for spatial discretization. The resulting non-linear
system of equations was solved using the Newton Raphson method
with the Jacobian calculated analytically. The simulator was validated
Fig. 18. Temperature solution after shut-in. by comparing the obtained solutions to different literature and reference
solutions. The proposed model was able to handle the phase appearance
and disappearance during phase transitions along the pipeline and to

13
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Fig. 19. Overall composition (zc ) for the shut-in period: (a) shut-in time (t = 0.0 s) and (b) t = 52000 s after shut-in.

take into account the effect of composition and temperature on the Investigation, Writing – original draft, Writing – review & editing,
pressure gradient. Those phenomena are directly involved in artificial Adriano dos Santos: Methodology, Formal analysis, Investigation,
lift processes and well testing procedures in reservoirs of complex hy­ Writing – review & editing, Carlos Pico Ortiz: Methodology, Software,
drocarbon mixtures with contaminants. Formal analysis, Investigation, Writing – review & editing, Adolfo Pires:
Frictionless low miscibility problem simulations showed excellent Methodology, Formal analysis, Investigation, Writing – review & editing
agreement with reference benchmark solutions. Nevertheless, the
miscible solution did not capture perfectly the pressure shock that ap­
pears at the interface between the liquid single phase region and the Declaration of competing interest
two-phase region.
A test consisting of a well kick-off followed by production with The authors declare that they have no known competing financial
subsequent shut-in was proposed and the numerical results were dis­ interests or personal relationships that could have appeared to influence
cussed. In summary, it was shown that the proposed numerical simulator the work reported in this paper.
is able to predict surge pressure effects, and co-current/countercurrent
flow at the beginning of the production and shut-in periods. Acknowledgments
To the best of our knowledge, a compositional fully-implicit solution
using the two-fluid model, has not been published for the cases pre­ The authors are grateful for the financial support provided by PET­
sented here. ROBRAS (grant No. 2013/00029–9). J. C. Nascimento acknowledges the
PhD scholarship provided by CAPES/CNPq (grant no. 235263/2014–1).
Credit author statement We also thank Prof. Albert C. Reynolds (University of Tulsa) and Prof.
Santos Alberto E. Remigio (Universidade Federal de Uberlândia) for
Júlio Nascimento: Methodology, Software, Formal analysis, useful discussions and suggestions.

Nomenclature

Latin
A cross sectional area [m2 ]
CD drag coefficient
Ci interfacial drag
D inner diameter [m]
dd diameter of dispersed phase [m]
g acceleration gravity [m/s2 ]
f fugacity [Pa] ; friction factor
h molar enthalpy [J/mol]
J Jacobian matrix
K vapor-liquid equilibrium ratio
MD drag force per volume [N/m3 ]
qn mass source term [mol/m3 .s]
n number of moles
Np number of phases

14
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

Nc number of components
P pressure [Pa]
R residual vector
Qh energy source term [J/m3 .s]
Qw radial heat transfer [J/m3 .s]
rto outer tubing radius [m]
Re Reynolds number
t time [s]
T temperature [K]
Ut overall heat transfer coefficient [W/m2 .K]
v velocity [m/s]
V volume [m3 ]
x mole fraction
X unknown vector
W Well index [m3]
z spatial coordinate; global composition
Z compressibility factor

Greek
α volume fraction
ξ molar density [mol/m3 ]
ρ mass density [kg/m3 ]
θ pipe inclination
σ Bestion damping coeficient
λ phase monility [1/Pa.s]
τ shear stress [N/m3 ]
Φp friction two-phase multiplier
Δt time-step [s]
Δz gridblock length [m]
ΔPi pressure correction term [Pa]
δX increment of unknown vector

Subscripts
c, p component c in phase p
d dispersed phase
f continuous phase
i interfacial
g gas phase
k gridblock index
l liquid phase
p phase
w wall
A annular
S slug
s superficial
SA slug - annular transition

Superscripts
n+ 1 new time
n old time
ν iteration index
R reservoir

Appendix A. Discretized governing equations

The residual equation corresponding to the mass balance equation for component c (given in Eq. (1)) is given by

⎡ ⎤ 1
Np,k+ n+1
(AΔz)k ⎣
Np,kn+1
∑ ( )n+1 Np,kn
∑ ( )n ∑2 ( )n+1
Rck = n+1
Np,k αp xc,p ξp k
− n
Np,k αp xc,p ξp k ⎦ + Ak+12 n+1
Np,k+ 1 αp xc,p ξp vp
k+1
Δt p p p
2 2

1 (A.1)
Np,k− n+1
∑2 ( )n+1 ( )n+1
n+1
− Ak− 1 Np,k− 1 αp xc,p ξp vp k− 1 − qnc,p ,
2 2 2 k
p

Similarly, the residual equation corresponding to the energy mixture equation is given by

15
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

⎡ n+1 ( )2 )n+1 Np,k ( )2 )n ⎤

Np,k ( ∑ ( )n (
n
(AΔz)k ⎣ ∑ ( )n+1 vp vp ⎦+
Rhk = αp k ξp hp + ρp − αp k ξp hp + ρp
Δt p
2 k p
2 k

N n+1 1 ( )2 ) ]n+1 N n+1

( )2 )]n+1
p,k+ [ ( p,k− 1 [ (
∑ 2
vp ∑ 2
vp
Ak+12 αp ξp hp + ρp vp − Ak− 12 αp ξp hp + ρp − (A.2)
p
2 k+ 1
p
2 k− 1
2 2

⎛ ⎞
N n+1
(AΔz)k ⎜ n+1 ⎟ ∑
p,k
[ ]n+1
⎝Pk − Pnk ⎠ − (AΔz)k αp ρp vp k g sin (θ)n+1 n+1
k+12 − (Qw )k − (Qh )n+1
k .
Δt p

Based on the staggered cell control volume (Fig. 1), the residual momentum equation for phase p momentum is written as
⎡ ⎤
P
(AΔz)k+1 ⎢( )n+1 ( )n ⎥ ( )n+1
Rp,k = 2
⎣ αp ρp vp k+1 − αp ρp vp k+1 ⎦ + Ak+1 αp ρp vp vp k+1
Δt 2 2

( )n+1 ( )n+1 [ ] ( )n+1 (A.3)

− Ak αp ρp vp vp k + Ak+12 αp k+1 Pn+1 k+1 − Pk
n+1
+ (AΔz)k+1 τwp k+1
2
( )n+1
2 2
( )n+1
n+1 [( )n+1 ( )n+1 ]
+(AΔz)k+1 αp ρp k+1 g sin (θ)n+1
k+1 + Ak+12 (ΔPi )k+1 αp k+1 − αp k − (AΔz)k+1 MpD 1 ,
2 2 2 2 2 k+2

Finally, the equilibrium and constraints residuals equations are expressed as follows:
( )n+1 ( )n+1
Rfkc = fc,l k − fc,g k , (A.4)

∑
Nc
( )n+1
(A.5)
x
c,p
Rp,k = xc,p k
− 1,
c

n+1
Np,k
∑ ( )n+1
(A.6)
α
Rk p = αp k
− 1.
p

Appendix B. Thermodynamic model

In this paper, the Peng-Robinson (PR) equation of state (Peng and Robinson, 1976) was used to determine the properties of the hydrocarbon phases
and equilibrium calculation:
RT a(T)
P= − (B.1)
v − b v(v + b) + b(v − b)

where a and b for a pure component are given by:

R2 Tc2
a = 0.45724α(T) , (B.2)
Pc

RTc
b = 0.07780 , (B.3)
Pc
{ [ ( )12 ]}2
T
α(T) = 1 + m 1 − (B.4)
Tc
{
0.3764 + 1.54226ωc − 0.26992ω2c , se ωc ≤ 0.491
m= (B.5)
0.379642 + 1.48503ωc − 0.1644ω2c + 0.01667ω3c , se ωc > 0.491

where Pc and Tc are critical pressure and temperature, respectively, and ωc is the acentric factor.
For a multicomponent mixture a and b are obtained by van der Waals’s mixing rule as follows:
∑
Nc ∑
Nc
√̅̅̅̅̅̅̅̅( )
a= xi,p xj,p ai aj 1 − δij , p = l, g (B.6)
i=1 j=1

∑
Nc
b= xi,p bi , p = l, g. (B.7)
i=1

where δij is the binary coefficient for components i and j.

Fluid mass and molar densities are given by:

16
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

PMWp
ρp = , (B.8)
Zp RT

ρp
ξp = , (B.9)
MWp

where MWp is the molecular weight, R is the universal gas constant and Zp is the compressibility factor.
In addition, it follows from (B.1) and (B.8) that
( ) ( )
Z 3 + (B − 1)Z 2 + A − 3B2 − 2B Z + B2 + B3 − AB = 0 (B.10)

where A and B coefficients are given by:

aP bP
A= , B= . (B.11)
(RT)2 RT

Eq. (B.10) has three roots. If only one root is real, this root defines the compressibility factor of the single phase. In the case where there are three
real positive roots, the largest and smallest define the compressibility factor of the vapor (gas) and liquid phases, respectively.
Fluid enthalpy in terms of PR EOS is given by:
∑
hp = xc,p hc,p , (B.12)
c

[ ( √̅̅̅ ) ]
( ) T ∂∂Ta − a Zp + 1 + 2 B
hc,p = hIG
c,p + RT Zp − 1 + √ ̅̅
̅ ln ( √̅̅̅ ) , (B.13)
2 2b Zp + 1 − 2 B

hIG 2 3 4 5
c,p = a0 + a1 T + a2 T + a3 T + a4 T + a5 T . (B.14)

where hc,p is the enthalpy of a pure component and hIGc,p is the enthalpy at the ideal gas state. The values of a0 -a5 coefficients are obtained in Passut and
Danner (1972).
Furthermore, the Lohrenz-Bray-Clark (LBC) (Lohrenz et al., 1964) correlation is used to determine the fluid viscosities.
Finally, the flowchart showing the two-phase flash calculation steps in the proposed simulator is shown in Fig. B.20. The procedure consists in
determining the phase compositions (xc,l and xc,l , c = 1, …, Nc ) and the number of mols of liquid (L) and gas (V). In this paper, the successive sub­
stitution method is employed to solve the thermodynamic equilibrium relation Eq. (7).

Fig. B.20. Algorithm to perform a two phase flash calculation.

17
J.C.S. Nascimento et al. Journal of Petroleum Science and Engineering 205 (2021) 108923

References Lockhart, R.W., Martinelli, R.C., 1949. Proposed correlation of data for isothermal two-
phase, two-component flow in pipes. Chem. Eng. Prog. 45 (1), 39–48.
Lohrenz, J., Bray, B.G., Clark, C.R., 1964. Calculating viscosities of reservoir fluids from
Almehaideb, R.A., Aziz, K., Pedrosa, O.A., 1989. A reservoir/wellbore model fo
their compositions. J. Petrol. Technol. 16 (10), 1171–1176.
multiphase injection and pressure transient analisys. In: SPE Middle East Oil
Micaelli, J.C., 1987. CATHARE an Advanced Best-Estimate Code for PWR Safety
Technical Conference and Exhibition. Society of Petroleum Engineers.
Analysis. SETH/LEML-EM.
Ambroso, A., Chalons, C., Raviart, P.-A., 2012. A Godunov-type method for the seven-
Munkejord, S.T., Evje, S., Flatten, T., 2009. A musta scheme for a nonconservative two-
equation model of compressible two-phase flow. Comput. Fluids 54, 67–91.
fluid model. SIAM J. Sci. Comput. 31, 2587–2622.
Barre, F., Bernard, M., 1990. The CATHARE code strategy and assessment. Nucl. Eng.
Nascimento, J.C.S., dos Santos, A., Pires, A.P., 2018. A fully-implicit solution for the
Des. 124, 257–284.
single-pressure two-fluid model with sharp discontinuities. Comput. Fluids 175,
Bendiksen, K.H., Malnes, D., Moe, R., Nuland, S., 1991. The dynamics two-fluid model
214–229.
OLGA: theory and application. SPE Prod. Eng. 6, 171–180.
Paillère, H., Corre, C., Cascales, G.J.R., 2003. On the extension of the AUSM+ schemes to
Bestion, D., 1990. The physical closure law in the CATHARE code. Nucl. Eng. Des. 124,
compressible two-fluid models. Comput. Fluids 32, 891–916.
229–245.
Pan, L., Oldenburg, C.M., 2014. T2well an integrated wellbore–reservoir simulator.
Bonizzi, M., Issa, R., 2003. A model for simulating gas bubble entrainment in two-phase
Comput. Geosci. 65, 46–55.
horizontal slug flow. Int. J. Multiphas. Flow 29 (11), 1685–1717.
Passut, C.A., Danner, R.P., 1972. Correlation of ideal gas enthalpy, heat capacity and
Borkowski, J.A., Wade, N.L., 1992. TRAC-BF1/MOD1 Models and Correlations. NREG/
entropy. Ind. Eng. Chem. Process Des. Dev. 11 (4), 543–546.
CR-1391.
Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. McGrAW-HILL Books.
Cao, H., 2002. Development of Techniques for General Purpose Simulators. Ph.D. thesis.
Peng, D.Y., Robinson, D.B., 1976. A new two-constant equation of state. Ind. Eng. Chem.
Stanford University. Department of Petroleum Engineering.
Fundam. 15, 59–64.
Chao, F., Liu, D., Shan, J., Gou, J., Wu, P., 2018. Development of temporal and spatial
Pourafshary, P., Varavei, A., Sepehrnoori, K., Podio, A., 2009. A compositional wellbore/
high-order schemes for two-fluid seven-equation two-pressure model and its
reservoir simulator to model multiphase flow and temperature distribution. J. Petrol.
applications in two-phase flow benchmark problems. Int. J. Numer. Methods Fluid.
Sci. Eng. 69 (1), 40–52.
88 (4), 169–192.
Prosperetti, A., Tryggvason, 2007. Computational Methods for Multiphase Flow.
Chen, Z., 2007. Reservoir simulation: mathematical techniques in oil recovery. In: CBMS-
Cambridge University Press.
NSF Regional Conference Series in Applied Mathematics. Society for Industrial and
Putcha, V.B., Ertekin, T., 2017. A fast and robust compositional, multi-phase, non-
Applied Mathematics.
isothermal wellbore hydraulics model for vertical wells. In: SPE Annual Technical
Coquel, F., Amine, K.E., Godlewski, E., Perthame, B., Rascle, P., 1997. A numerical
Conference and Exhibition. Society of Petroleum Engineers.
method using upwind schemes for the resolution of two-phase flows. J. Comput.
Raimondi, L., 2017. Compositional simulation of two-phase flows for pipeline
Phys. 136, 272–288.
depressurization. SPE J. 22, 1–242.
Delhaye, J.M., Achard, J.L., 1976. On the averaging operators introduced in two-phase
Ransom, V.H., 1982. RELAP5/MOD1 Code Manual Volum 1: Code Structure, System
flow modelling. In: CSNI Meeting in Transient Two-phase Flow. CSNI.
Models and Numerical Method. NUREG/CR-1826.
Dinh, T.N., Nourgaliev, R.R., Theofanous, T.G., 2003. Understanding the ill-posed two-
Ransom, V.H., Hicks, D.L., 1984. Hyperbolic two-pressure models for two-phase flow.
fluid model. In: The 10th International Topical Meeting on Nuclear Reactor Thermal
J. Comput. Phys. 53, 124–151.
Hydraulics (NURETH-10).
RELAP5, February 2001. RELAP5-3D Code Manual Volume 1: Code Structure, System
Drew, D., Lahey Jr., R.T., 1979. Application of general constitutive principles to
Models and Solution Method. Idaho National Engineering and Environmental
derivation of mutidimensional two phase flow equations. Int. J. Multiphas. Flow 5,
Laboratory.
243–264.
Shahamiri, A., Heidari, M., Buchanan, W., Nghiem, L., et al., 2015. A CFD compositional
Evje, S., Flåtten, T., 2003. Hybrid flux-splitting schemes for a commom two-fluid model.
wellbore model for thermal reservoir simulators. SPE Annual Technical Conference
J. Comput. Phys. 192, 175–210.
and Exhibition. Society of Petroleum Engineers.
Evje, S., Flåtten, T., 2005. Weakly implicit numerical schemes for two-fluid model. SIAM
Shekari, Y., Hajidavalloo, E., 2013. Apllication of osher and price-c schemes to solve
J. Sci. Comput. 26, 1449–1484.
compressible isothermal two-fluid models of two-phase flow. Comput. Fluids 86,
Evje, S., Flåtten, T., 2006. CFL-violating numerical schemes for a two-fluid model. J. Sci.
363–379.
Comput. 29 (1), 83–114.
Shi, H., Holmes, J.A., Durlofsky, L.J., Aziz, K., Diaz, L.R., Alkaya, B., Oddie, G., 2005.
Fair, W.B., 1981. Pressure buildup analysis with wellbore phase redistribution. Soc.
Drift-flux modeling of two-phase flow in wellbores. SPE J. 10 (1), 24–33.
Petrol. Eng. J. 21, 259–270.
Shirdel, M., Sepehrnoori, K., 2012. Development of a Transient Mechanistic Two-phase
Forouzanfar, F., Pires, A.P., Reynolds, A.C., 2015. Formulation of a transient multi-phase
Flow Model. Society of Petroleum Engineers, pp. 942–955.
thermal compositional wellbore model and its coupling wih a thermal compositional
Shirdel, M., Sepehrnoori, K., 2017. Development of transient mechanistic three-phase
reservoir simulator. In: SPE Annual Technical Conference and Exhibition. Society of
flow model for wellbores. Society of Petroleum Engineers 374–378.
Petroleum Engineers.
Shoham, O., 2006. Mechanistic Modeling of Gas-liquid Two-phase Flow in Pipes. SPE.
Garcia-Cascales, J.R., Paillère, H., 2006. Application of AUSM schemes to multi-
Sondermann, C.N., Baptista, R.M., de Freitas Rachid, F.B., Bodstein, G.C.R., 2019.
dimensional compressible two-phase flow problems. Nucl. Eng. Des. 236,
Numerical simulation of non-isothermal two-phase flow in pipelines using a two-
1225–1239.
fluid model. J. Petrol. Sci. Eng. 173, 298–314.
Han, G., Ling, K., Khor, S.H., Zhang, H., Thakur, R.K., et al., 2013. Simulation of
Städtke, H., 2006. Gasdynamics Aspects of Two-phase Flow. Wiley-VCH.
multiphase fluid-hammer effects during well startup and shut-in. Oil and Gas
Stone, T.W., Edmunds, N.R., Kristoff, B.J., 1989. A comprehensive wellbore/reservoir
Facilities 2, 68–77.
simulator. In: SPE Symposium on Reservoir Simulation. Society of Petroleum
Hasan, A.R., Kabir, C.S., 1994. Aspects of wellbore heat transfer during two-phase flow.
Engineers.
SPE Prod. Facil. 9, 211–216.
Suarel, R., Abgrall, R., 1999. A multiphase Godunov method for compressible multifluid
Ishii, M., 1975. Thermo-fluid Dynamics of Two-phase Flow. Eyrolles.
and multiphase flows. J. Comput. Phys. 150, 425–467.
Ishii, M., Hibiki, T., 2006. Thermo-Fluid Dynamics of Two-phase Flow. Springer.
US-NRC, 2007. TRACE V5.0 Theory Manual: Field Equations, Solution Methods and
Ishii, M., Mishima, K., 1984. Two-fluid model and hydrodynamic constitutive relations.
Physical Models. United States Nuclear Regulatory Commission.
Nucl. Eng. Des. 82, 107–126.
Versteeg, H.K., Malalasekera, W., 2007. An Introduction to Computational Fluid
Issa, R.I., Kempf, M.H.W., 2003. Simulation of slug flow in horizontal and nearly
Dynamics. Pearson.
horizontal pipes with the two-fluid model. Int. J. Multiphas. Flow 29, 69–65.
Wallis, G.B., 1969. One-Dimensional Two-phase Flow. McGraw-Hill Book Company.
Kataoka, I., Ishii, M., 1987. Drift-flux model for large diameter pipe and new correlation
Xiong, W., Bahonar, M., Chen, Z., et al., 2015. Development of a thermal wellbore
for pool void fraction. Int. J. Heat Mass Tran. 30 (9), 1927–1939.
simulator with focus on improving heat loss calculations for SAGD steam injection.
Krasnopolsky, B.I., Lukyanov, A.A., 2018. A conservative fully implicit algorithm for
In: SPE Canada Heavy Oil Technical Conference. Society of Petroleum Engineers.
predicting slug flows. J. Comput. Phys. 355, 597–619.
Yeom, G.-S., Chang, C.H., 2011. Flux-based wave decomposition scheme for an isentropic
LeVeque, R.J., 2004. Finite-Volume Methods for Hyperbolic Problems. Cambridge
hyperbolic two-fluid model. Numer. Heat Tran., Part B: Fundamentals 59 (4),
University Press.
288–318.
Livescu, S., Durlofsky, L.J., Aziz, K., Ginestra, J.C., 2009. Development and application of
Zigrang, D.J., Sylvester, N., 1982. Explicit approximations to the solution of Colebrook’s
a fully coupled thermal compositional wellbore flow model. In: SPE Western
friction factor equation. AIChE J. 28 (3), 514–515.
Regional Meeting. Society of Petroleum Engineers.
Zuber, N., Findlay, J., 1965. Average volumetric concentration in two-phase flow
Livescu, S., Durlofsky, L.J., Aziz, K., Ginestra, J.C., 2010. A fully-coupled thermal
systems. J. Heat Tran. 87 (4), 453–468.
multiphase wellbore flow model for use in reservoir simulation. J. Petrol. Sci. Eng.
71, 138–146.

18