Proceedings of ENCIT 2016 16th Brazilian Congress of Thermal Sciences and Engineering

Copyright © 2016 by ABCM November 07-10th, 2016, Vitória, ES, Brazil

ASSESSMENT OF THE DYNAMIC PRESSURE EFFECT IN 1D TWO

FLUID MODEL FOR A VERTICAL ANNULAR GAS-LIQUID FLOW
Eric M. González, eric.gfont@gmail.com
Angela O. Nieckele, nieckele@puc-rio.com
Department of Mechanical Engineering, PUC-Rio, Rio de Janeiro, RJ, 22451-900 Brazil

João N. E. Carneiro, joao.carneiro@sintefbrasil.org.br

Instituto SINTEF do Brasil, Rio de Janeiro, RJ, 22290-906, Brazil

Abstract: Vertical annular flows are very common in many industrial applications. To predict this pattern of two-phase
flow, the Two Fluid Model can be applied. However, it is well known that the resulting system of equations is ill-posed.
Therefore, at the present work, a dynamic pressure term is included on the liquid phase to help increase the region at which
a grid-independent solution can be found. The conservation equations are discretized with the finite volume method, with
a first order time integration and second order TVD spatial discretization. The stabilizing effect of the dynamic pressure
is addressed. Flow parameters such as pressure drop, film thickness and wave characteristics numerically obtained are
compared against available experimental data, presenting good agreement.

Keywords: Dynamic pressure, Vertical Annular Flow, Gas-Liquid Flow, Two-Fluid Model 1D

1. NOMENCLATURE

Nomenclature usk phase superficial velocity, [m/s]

A pipe cross-sectional area, [m2 ] W viscosity number, [-]

C momentum flux parameter, [-] x spatial coordinate (distance from inlet), [m]

D internal diameter, [m] Re Reynolds number, [-]

Dhg gas hydraulic diameter, [m] Greek

F frictional force per unit volume, [P a/m] α volume fraction

g gravity acceleration, [m/s2 ] ρ density, [kg/m3 ]

L pipeline length, [m] µ dynamic viscosity, [Pa s]

Nσ surface tension factor, [-] σ surface tension, [N/m]

P average pressure, [Pa] τ shear stress, [Pa]

Rg universal constant gas, [m2 /s2 K] Subscripts

S wetted perimeter or interfacial contact length, [m] g gas phase

T temperature, [K] i interface

t time, [s] k concerning both gas and liquid phase

u phase velocity, [m/s] l liquid phase

ur wave velocity, [m/s] w wall

2. INTRODUCTION

Annular two-phase flow is encountered in many industrial applications: biological processes, nuclear industry, oil and
gas production systems. Annular pattern occurs typically at higher gas velocities. The liquid phase forms a film which
is distributed around the pipe circumference while the gas phase flows in the core. For vertical pipes the liquid film
distribution around the circumference can be considered as symmetric (Fowler and Lisseter, 1992; Inácio, 2012) unlike
horizontal pipes.
The high velocity gas disturbs the liquid film surface and forms waves, ripples and disturbance waves (Berna et al.,
2014; Zhao et al., 2013). When this difference is high enough, droplets are entrained from the liquid film surface into the
Proceedings of ENCIT 2016 16th Brazilian Congress of Thermal Sciences and Engineering
Copyright © 2016 by ABCM November 07-10th, 2016, Vitória, ES, Brazil

core of the gas stream. (Berna et al.,2014; Zhao et al.,2013).

To determine a multiphase flow, the Two-Fluid Model can be employed. Through an analysis of the characteristics of
the hyperbolic system of equations, it is possible to determine whether the characteristics are real and the system is well-
posed. It is well known that the classical 1D Two-Fluid Model is unconditionally ill-posed for vertical pipes because the
hydrostatic pressure term is negligible (Montini, 2011; Inácio, 2012; Carneiro, 2006). Therefore, in the present work, an
isothermal 1D Two-Fluid model is formulated including a dynamic pressure term, in order to render the equation system
well-posed for vertical pipes. Different dynamic pressure models are investigated, and the formulation was tested for an
air-water upwards (co-current) annular flow. The model equations are solved numerically with a finite volume technique.
Simulation results are compared with experimental data available in the literature.
In the next section, the mathematical modeling is presented, followed by the numerical method, simulation results and
finally conclusions.

3. MATHEMATICAL MODELING

Transport systems in oil/gas and nuclear industries are composed with kilometers of pipeline, in which the axial
variations of the flow field are more relevant than cross sectional variations. For this reason a one-dimensional form of
the governing equations was selected. The mathematical model employed to predict the vertical two-phase flow was the
Two-Fluid Model (Ishii and Hibiki, 2010).
The Two-Fluid Model consists of a set of conservation equations for each phase, obtained through a space average
over each phase. The volume fraction of each phase k, αk can be defined as
∀k
αk = (1)
∀
thus, the following restriction must apply
αg + αl = 1 (2)
The one-dimensional formulation was established through a space average over the cross-section. The mass and
momentum conservation equations for an upward vertical annular flow are:
∂(αg ρg ) ∂(αg ρg ug ) ∂(αl ρl ) ∂(αl ρl ul )
+ =0 ; + =0 (3)
∂t ∂x ∂t ∂x
∂(αg ρg ug ) ∂(Cg αg ρg u2g ) ∂Pgi ∂[αg (Pg − Pgi )]
+ = −αg − − αg ρg g − Fi (4)
∂t ∂x ∂x ∂x

∂(αl ρl ul ) ∂(Cl αl ρl u2l ) ∂Pli ∂[αl (Pl − Pli )]

+ = −αl − − αl ρl g − Fwl + Fi (5)
∂t ∂x ∂x ∂x
The liquid phase is taken as incompressible, while the density of gas phase is governed by the ideal gas law
Pg
ρg = (6)
Rg T
where T is the reference temperature of an isotherm flow.
To closure the system of equations, geometric parameters, momentum flux parameters, shear stress contributions, as
well as the pressure difference from the phase and interface must be defined, and are presented next.

3.1 Geometric parameters

In annular vertical flow the geometric parameters can be calculated from the configuration shown in Fig. 1.
 2
hl D √
αg = 1 − 2 ; hl = (1 − αg ) (7)
D 2
πD2
Sl = πD ; Si = πDg ; Dg = D − 2hl ; A= (8)
4
3.2 Momentum flux parameter

The momentum flux parameter is introduced to take into account that the velocity is not uniform along the cross
section. It is defined as
R 2
u dAk
Ck = R k 2 (9)
uk dAk
The recommended values by Fowler and Lisseter (1992) for the gas and liquid flux parameters are Cg = 1 and Cl = 1.33.
Proceedings of ENCIT 2016 16th Brazilian Congress of Thermal Sciences and Engineering
Copyright © 2016 by ABCM November 07-10th, 2016, Vitória, ES, Brazil

Figura 1: Geometric parameters

3.3 Shear stresses

In the eqs. (4) - (5) the frictional force per unit volume between liquid and wall, and at the interface are modeled as:
Sl τl Si τi
Fwl = ; Fi = (10)
A A
The shear stress for the liquid phase as well as at the interface were determined considering the flow as locally fully
developed, thus:
fl ρl |ul | ul fi ρg |ug − uf l |(ug − uf l )
τwl = ; τi = (11)
2 2
Here uf l is the liquid film velocity at the interface. Its definition is connected with each dynamic pressure model consi-
dered, and it will be presented in the dynamic pressure section.
The friction factor depends on the Reynolds numbers which can be calculated as:
ρl usl D ρg |ug − uf l | Dhg
Resl = ; Rei = (12)
µl µg
where the liquid superficial velocity is usl = αl ul and Dhg is gas hydraulic diameter defined as Dhg = 4Ag /Si .
There are many empirical correlations available in the literature to determine the friction factor. The correlations
selected to predict the friction factors for both laminar and turbulent regime have been employed by Alves et al. (2012);
Inácio (2012) and Berna et al. (2014), who reported reasonable friction factor prediction for both regimes.
The friction factors considered for laminar regime (Resl and Rei ≤ 2100) are:
"  0.33  #
24 16 ρl hl
fl = ; fi = 1 + 24 (13)
Resl Rei ρg D

and for turbulent regime (Resl and Rei > 2100) are:
"  0.33  #
0.0262 0.079 ρl hl
fl = ; fi = 1 + 24 (14)
(αl Resl )0.139 Re0.25
i ρ g D

3.4 Dynamic Pressure

The second pressure term on the right hand side of Eqs. (4) and (5) represent the difference between the bulk average
pressure and the interfacial average pressure.

∂[αk (Pk − Pki )] ∂(αk ∆Pki )

− =− (15)
∂x ∂x
The difference can be written based on the dynamic pressure closures proposed by different authors. They are described
below, with the interface liquid film velocity definition uf l :

3.4.1 Model 1

The first model considered was proposed by Fowler and Lisseter (1992). The gas interface pressure is considered
equal to the gas pressure (Pgi = Pg ). The difference bewteen the interface and average liquid velocity is the dynamic
pressure, which depends on the wave velocity of the liquid film interface ur
∆Pli = 0.02 ρl (ul − ur )2 ; ur = 2 ul (16)
In this model, the interface liquid film velocity definition, relevant for the friction closures, is equal to the wave
velocity, i.e., uf l = ur = 2 ul .
Proceedings of ENCIT 2016 16th Brazilian Congress of Thermal Sciences and Engineering
Copyright © 2016 by ABCM November 07-10th, 2016, Vitória, ES, Brazil

3.4.2 Model 2

The second model is based on the previous model, however, according to Berna et al. (2014), the wave velocity
depends on both phase flow conditions
√ √
2
ρg usg + ρl usl −0.38 0.16 −0.13
∆Pli = 0.02 ρl (ul − ur ) ; ur = 50 √ √ Resg Resl Nσ (17)
ρg + ρl

0.028 W −4/5

W ≤ 1/15 µl
Nσ = ; W = (18)
0.25 W > 1/15 q 1/2
ρl σ g(ρlσ−ρg )

The gas superficial velocity and Reynolds number are usg = αg ug and Resg = ρg usg Dhg /µg . The interface liquid
film velocity is equal to uf l = D/(4hl ) usl .

3.4.3 Model 3

The third model selected to be analysed was proposed by Bestion (1990). Here, there is a pressure difference for both
phases related to the dynamic pressure.
αl αg ρl ρg
∆Pki = 1.2 ρm (ul − ug )2 ; ρm = (19)
αg ρl + αl ρg
The interface liquid film velocity is equal to the liquid velocity uf l = ul .

4. NUMERICAL SCHEME

The finite-volume method was selected to discretize the system of equations. Scalar variables (volume fractions, den-
sities and pressure) were stored at the central nodal point, while a staggered mesh was employed for storing the velocities.
The equations were discretized using a high order TVD (Total Variation Diminishing) scheme for the convective term
using the Van Leer function as flux limiter (Versteeg and Malalasekera, 2007). For the time integration a fully implicit
Euler (first order) scheme was employed.
At each time step, the conservation equations are solved sequentially with a method based on the PRIME algorithm
(Ortega and Nieckele, 2005) to handle the velocity–pressure coupling. This algorithm consists of solving sequentially the
liquid momentum equation followed by the gas momentum equation in order to determine the respective phase velocities.
The global mass conservation equation was obtained by adding the normalized gas and liquid mass conservation equations.
Pressure is then determined by combining the global mass conservation equation with both phase momentum equations.
Each phase velocity was corrected with the new pressure field as well as their respectively explicit discretized momentum
equations (Ortega and Nieckele,2005; Simões et al.,2014). The system of algebraic equations for each variable is solved
with the TDMA algorithm (Patankar, 1980).

4.1 Time step and flow evolution

The time step must be small enough to capture the waves instabilities ( Simões et al., 2014; Inácio, 2012). The
Courant number establishes a relationship between the time step and the mesh spacing, thus:
∆t
Cnum = max(|u|) (20)
∆x
where max(|u|) represent the maximum velocity in the whole domain taking account both gas and liquid phases. The
Courant number value for all cases examined was maintained at 0.05.
During the flow evolution, fluctuations of velocities, pressure, and film thickness are observed. The time average
characteristics parameters of the annular flow where determined only after an initial transient, i.e., after the time mean
variables have reached an approximately constant value. Figure 2 shows the instantaneous and time average evolution of
the pressure gradient and spatial mean value of the film thickness and gas velocity.
During the initial transient strong variation of pressure gradient, and spatial mean film thickness and velocities are
observed. It can be seen in Figure 2(b) that the film thickness presented a quick drop toward its stabilized value. The time
average values take around 40 [s] to stabilize presenting variations less than 1%.

5. RESULTS AND DISCUSSION

To evaluate the performance of the three dynamic pressure models implemented, the experimental results of an air-
water flow of Zhao et al. (2013) was selected. It consists of an upward air-water flowing inside a vertical pipe with
D = 0.0345[m] and L = 2[m]. The fluids properties are presented in Tab. 1. The universal gas constant for air is
Proceedings of ENCIT 2016 16th Brazilian Congress of Thermal Sciences and Engineering
Copyright © 2016 by ABCM November 07-10th, 2016, Vitória, ES, Brazil

(a) Pressure gradient (b) Film thickness (c) Gas velocity

Figura 2: Evolution of instantaneous and time average pressure gradient, and spatial mean film thickness and gas velocity.

Tabela 1: Fluids properties.

Fluids ρ[kg/m3 ] µ[Pa s]
Air 1.18 1.79 × 10−5
Water 998.20 1.00 × 10−3

Rg = 287[N m/(kgK)]. The inlet gas and liquid superficial velocities are usg = 40.1[m/s] and usl = 1.75×10−2 [m/s],
respectively. The outflow pressure is the atmospheric pressure, and the operational temperature was equal to 298.15 [K].
To define the mesh, a grid test was performed, whereby an order of magnitude change in grid size was investigated,
with ∆x/D varying from 0.5 to 0.05. One consequence of an ill-posed system of equations is that a mesh independent
solution cannot be found. Thus, the grid test helps to identify which dynamic pressure models renders the system well-
posed.
The results for the mean film thickness, peak frequency from PSD (Power Spectrum Density) plot at x = 1.98[m]
and pressure gradient were plotted in Fig. 3, against dimensionless mesh spacing and are compared with the experimental
data. The uncertainty associated with the experimental data varied approximately from 5% to 7% (Zhao et al.,2013).

(a) Mean film thickness (b) PSD peak frequency (c) Pressure gradient
Figura 3: Grid test performed using the different models.

Analysing Fig. 3 corresponding to Model 1 of Fowler and Lisseter (1992), it can be clearly seen that the dynamic
pressure was not able to stabilized the model. Figs. 3(a) - 3(c) for the mean film thickness, power spectral density
peak frequency and pressure gradient show that these variables increase unboundedly with grid refinement, i.e., a mesh
independent solution was not found, indicating that the model is ill-posed. On the other hand, Model 2 (Berna et al.
(2014)) and Model 3 (Bestion (1990)) not only presented a mesh independent result for grid size inferior to ∆x/D < 0.1,
but also present good agreement with respect to the reference value, as can be observed in Figs. 3(a) - 3(c). For the case
analyzed, the system of equation became well-posed for Models 2 and 3, with the introduction of the dynamic pressure.

5.1 Axial distribution of flow parameters

Since Models 2 and 3 rendered well posed results, the distribution of flow parameters along the pipe are examined.
The mesh size employed was ∆x/D = 0.05. Figure 4(a) corresponds to the axial distribution of the time mean film
thickness, while Fig. 4(b) presents the power spectrum frequency obtained with Model 2 and 3.
Analyzing Fig. 4(a) it can be observed from the experimental data, an approximately constant film thickness for a
distance greater than 0.8 [m] from the inlet, where the flow can be considered as fully developed. Model 3 presents a
slightly better agreement with the experimental data than the Model 2, with a relative error around 19% with regard to the
Proceedings of ENCIT 2016 16th Brazilian Congress of Thermal Sciences and Engineering
Copyright © 2016 by ABCM November 07-10th, 2016, Vitória, ES, Brazil

(a) Mean film thickness (b) PSD peak frequency

Figura 4: Mean film thickness against axial distance from the inlet

reference in the fully developed region.

Figure 4(b) shows that both models presented a similar behavior for the variation of the dominant wave frequencies
with the axial distance. An approximate constant value was obtained for positions larger than 0.6 [m] for both models,
agreeing with the experimental data. In the fully developed region, reasonable agreement was obtained, with relative error
equal to 15% and 20% for Model 2 and Model 3, respectively.

6. CONCLUSIONS

It has been shown that the dynamic pressure closure for vertical annular flow may render the 1D Two Fluid Model well
posed, depending on the formulation chosen. The effect of the dynamic pressure model of Fowler and Lisseter (1992) was
not significant in this regard. Results were not grid independent, indicating ill-posedness. On the other hand, the model
by Bestion (1990) and the modified version of the model by Fowler and Lisseter (1992), including modifications inspired
by Berna et al. (2014), presented grid independent results for the mean film thickness, frequency and pressure gradient,
in good agreement with experiments.

7. ACKNOWLEDGEMENTS

The first two authors would like to acknowledge the supported granted during the development of this work by the
Brazilian Government agencies: CNPq and CAPES.

8. REFERENCES

Alves, M.V.C., Barbosa Jr, J. and Falcone, G., 2012. “Modeling the transient behavior of churn-annular flow in a vertical
pipe”. In 3rd Brazilian Conference on Boiling, Condensation and Multiphase flow. Vol. 400.
Berna, C., Escrivá, A., Muñoz-Cobo, J. and Herranz, L., 2014. “Review of droplet entrainment in annular flow: Interfacial
waves and onset of entrainment”. Progress in Nuclear Energy, Vol. 74, pp. 14–43.
Bestion, D., 1990. “The physical closure laws in the cathare code”. Nuclear Engineering and Design, Vol. 124, No. 3,
pp. 229–245.
Carneiro, J.N.E., 2006. Simulação Numérica de Escoamentos Bifásicos No Regime de Golfadas em Tubulações Horizon-
tais e Levemente Inclinadas. Master’s thesis, PUC-RJ.
Fowler, A. and Lisseter, P., 1992. “Flooding and flow reversal in annular two-phase flows”. SIAM Journal on Applied
Mathematics, Vol. 52, No. 1, pp. 15–33.
Inácio, J.d.R.G., 2012. Simulação do Regime Intermitente em Tubulações Verticais Utilizando o Modelo de Dois Fluidos
com Diferentes Relações de Fechamento. Master’s thesis, PUC-RJ.
Ishii, M. and Hibiki, T., 2010. Thermo-fluid dynamics of two-phase flow. Springer Science & Business Media.
Montini, M., 2011. Closure relations of the one-dimensional two-fluid model for the simulation of slug flows. Ph.D. thesis,
Imperial College London.
Ortega, A. and Nieckele, A., 2005. “Simulation of horizontal two-phase slug flows using the two-fluid model with a
conservative and non-conservative formulation”.
Patankar, S., 1980. Numerical heat transfer and fluid flow. CRC Press.
Simões, E.F., Carneiro, J.N. and Nieckele, A.O., 2014. “Numerical prediction of non-boiling heat transfer in horizontal
stratified and slug flow by the two-fluid model”. International Journal of Heat and Fluid Flow, Vol. 47, pp. 135–145.
Versteeg, H.K. and Malalasekera, W., 2007. An introduction to computational fluid dynamics: the finite volume method.
Pearson Education.
Proceedings of ENCIT 2016 16th Brazilian Congress of Thermal Sciences and Engineering
Copyright © 2016 by ABCM November 07-10th, 2016, Vitória, ES, Brazil

Zhao, Y., Markides, C.N., Matar, O.K. and Hewitt, G.F., 2013. “Disturbance wave development in two-phase gas–liquid
upwards vertical annular flow”. International Journal of Multiphase Flow, Vol. 55, pp. 111–129.

9. RESPONSIBILITY NOTICE

The following text, properly adapted to the number of authors, must be included in the last section of the paper: The
author(s) is (are) the only responsible for the printed material included in this paper.