7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Finite Volume Mlodeling of Fully Developed Vertical Annular Flow
F. M. Fernandez, A. Alvarez Toledo, E. E. Paladino* and J. A. Lima
Graduate Program in Mechanical Engineering, Universidade Federal de Rio Grande do Norte, Natal, RN
59072970, Brazil
emilioiiact.ufm.br
Abstract
This work presents a numerical model for predicting lwdrodynamic and heat transfer parameters for turbulent,
annular, vertical, gasliquid twophase flow. It solves mass and momentum transport differential equations for
both the gas core and the liquid film in their entire domains. Thus, the velocity and shear stress distributions
from the tube center to the wall are obtained, together with the film thickness and the pressure gradient, making
no use of empirical closure relations. The model was developed using the finite volume technique and an
algebraic turbulence model.
Keywords: Annular flow, Finite Volume Method, Turbulence Modelling
Nomenclature
e Entrained fraction of liquid
g gravitational constant (ms')
p pressure (Nm )
rradial coordinate
R Inner radius (m)
T Temperature (K)
U Velocity (m/s)
zAxial coordinate
Greek letters
8 Mean film thickness (m)
a Void fraction
SDynamic viscosity (Pa s)
SShear stress (Pa)
Subscripts
g Gas
gc Gas core
i Interface
1 Liquid
If Liquid film
w Wall
Introduction
Two phase vertical annular flow is a
common pattern encountered in the process, nuclear
and oil industries. This pattern is characterized by the
existence of a liquid film adjacent to the wall, a gas
core flowing in the center of the duct and a wavy
interface. The gas core can have liquid droplets
entrained and the liquid film can carry gas bublles.
As the gas mass rate increases, the morphology of the
interfacial waves changes, increasing the complexity
of the models.
After a review of the existing literature,
there were found different ways to model both phases
and their interaction. Regarding the liquid film, most
approaches make use of the triangular relationship
existing between pressure gradient, liquid film
thickness and liquid mass flow rate, in order to relate
the liquid film velocity profile with the interfacial
shear stress, which is obtained from empirical
correlations (Pu et al. 2006), (Peng 2008), (Okawa &
Kataoka 2005), (Fu & Klausner 1997). Other
researchers propose a more detailed description of the
interface considering the presence of disturbance
waves, and again calculating the interfacial shear
stress through empirical correlations (Jayanti &
Hewitt 1997). In some cases the velocity profile is
assumed to be known from the logarithmic law
(Kishore & Jayanti 2004) or 1/7 profile (Adecly &
Issa 2004). Then, the mass and momentum equations
in the gas core are solved incorporating the effects of
the liquid film through a friction factor with an
equivalent surface roughness, for the computation of
the interfacial stress. Some other approaches consider
a more sophisticated description of motion within the
liquid film by solving mass and momentum equations
for the film but, in all cases, these equations are
averaged along the radial direction (Adecly & Issa
(211114), (Antal et al. 1998). (Dobran 1983) considers
the film as to be divided into two sublayers: a
continuous liquid layer next to the wall and a wavy
one next to the liquidgas interface: the velocity
profile within both layers is obtained from lawof
thewall type relations similar to those of single
phase turbulent flow.
Concerning the gas core modeling, it is
usually (Dobran 1983), (Kishore & Jayanti (1***4);
treated as an homogeneous mixture of gas and
entrained liquid droplets, not considering slip
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(6) The total mass flow rates of both the gas core
and the liquid film are assumed to be known.
As usual, the proposed model has two
distinct phases (Fig. 1), liquid film (lf) and the gas
core (gc); the former is treated as a continuous, fully
turbulent layer having no entrained bubbles, with a
mean film thickness 6. The wavv nature of the
interface is taken into account through the turbulence
model.
The gas core is considered to have liquid
droplets entrained, and in view of assumption (4) and
the homogeneous model its properties are calculated
as:
between both phases. This means the liquid droplets
are assumed to travel at gas speed and allows for the
use of empirical correlations to predict the void
fraction E and the fraction of liquid entrained e which
are needed to compute the dropletladen gas core
density and viscosity.
The aim of this work is to aclueve an
accurate, simple and complete numerical
computation of all the hydrodynamic and heat
transfer parameters requiring just gas and liquid mass
flow rates as input data, and making no use of
empirical closure correlations to relate shear stresses,
pressure gradient and film thickness.
The work is in its early stage of
development and more complexity will be
progressively added as the models herein presented
prove effective. Regarding turbulence modeling, a
twoequation model is to be implemented. For the
moment, the algebraic turbulence model of
(Cioncolini et al. 2009) is successfully used to predict
eddy diffusivities in both phases. Yet simple, it
implicitly accounts for the effect of the interfacial
complex phenomena and its contribution to the
pressure gradient, as it requires calibration upon
experimental data.
The model is validated against experimental
results of (Wolf et al. 2001). Their measured data for
pressure gradient, film thickness and wall shear stress
is compared with those predicted by the model for
fullydeveloped conditions.
Model description
This model predicts pressure gradient, liquid
film thickness and velocity profiles for turbulent fully
developed, gas liquid, annular flow, using mass
fluxes as input values. As well, it solves de Nusselt
number for the thermal entry, calculating temperature
profiles for axial distances. It requires no closure
empirical correlations as it solves momentum and
mass conservation equations simultaneously for the
given mass flow rates.
Modeling is carried out under the following
assumptions:
(1) Fullydeveloped, steadystate, cocurrent,
vertical, annular flow.
(2) Flow is assumed to be axissymmetric and film
thickness is assumed to be constant in angular
coordinate.
(3) Entrainment and deposition processes have
reached fullydeveloped conditions, meaning that
there is no net mass transfer between both phases.
(4) There is no slip between the gas and the
entrained liquid droplets.
(5) Fluids thermodynamical and transport
properties are assumed not to vary with temperature
or pressure, even when the thermal entry problem is
solved.
(1 xg) xg
PI P,
9U, = (1 agc). pI + age AUg
Figure 1: Problem description
where a,, is the mass fraction
calculated as:
of the droplets,
apix,
go PI Pix g(1xg)
For predicting the entrained liquid fraction,
the entrainment correlation of (Ishii & Mishima
1989) was used, as it is for hydrodynamic
equilibrium conditions:
e = tanh(7.25 107 We 2 ReP 2)
where the Weber number We and the liquid Reynolds
number Rel are calculated as:
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Body forces are excluded from the linear
momentum equation, assuming that their effect on
the pressure gradient is negligible. Inertial forces due
to gas specific volume changes arising from pressure
variations are also neglected, that is to say, only
viscous stresses contribute to the pressure gradient. In
light of their experimental results, Wolf et al. (2000)
validate these assumptions. They present
experimental evidence which seems to justify the
commonly made assumptions in annular flow that the
wall shear stress is almost equal to the interfacial
shear stress and that the pressure gradient in annular
flow is dominated by interfacial friction.
Mathematical formulation
The timeaveraged linear momentum and
energy equations for hydrodynamically fully
developed flow are:
Id du,
pyr
r dr dr
1 d du,,
:
pyr
r dr dr
0
RG
8 (PI cp, u Ty) 1~( 8 d 8Tl
8z r Br B r
8z r Br Br
0
Figure 2: Grids for gas core and liquid film domains
Each of the equations above is integrated in
its corresponding domain, as shown in Fig. 2.
However, both fields are solved simultaneously, as
schematically shown in the equation below.
A :
The boundary conditions for these equations are
du_
dr
This is possible because a nonuniform
effective viscosity is used, which is a function of the
radial position. So dynamic viscosities can be written
as:
u = r = R
0
dr
0
RG
Pety~= Igc
Pet =P ~u
T (z, R) = Tw
T(0, r)= To
In the case o prescribed wall flux,
Special treatment is required at the interface,
in order to satisfy the momentum and energy balance
across it. Following (Patankar 1980), using an
aHRlOgy between heat transfer and momentum
diffusion, an equivalent viscosity at the interface is
calculated as explained below.
The governing equations are integrated in
the context of the finite volume method, at the
VOlumes at each side of the interface, shown
schematically in Figure 3, in which the east face of
the last volume in the gas side coincides with the
physical interface between both phases. So does the
'r Ir=R
The continuity of fluxes at interface gives,
We=P Jd pi pg
~'=" gs
Re, 
*
. U B
gc gc

Air Uzr Bzr
so OfR3 R3 1
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
s. p, .t~
west face of the first volume of the liquid side. Then,
for the east face of the last volume in the gas domain:
ae
INTERFACE
Arg arl
Figure 3: Finite volumes at the interface
This correlation was fitted by the authors
from a large database for gasliquid flow.
Regarding the gas core effective viscosity, it
varies linearly with wall distance and inversely
proportional to a constant a.
Pefgc Pgci
V=a
uE P
For the computation of the constant a, the
authors use the same experimental database used for
eq. 25 and 26, and propose an average value of 4.2 +
1.0 with an approximate standard deviation of 24%.
In the present study, its value is set to 4.3.
The thermal entry problem is described by
energy equations (9), (10), (13)(15) and (17). As
done with the velocity field, both temperature fields
are solved considering a variable conductivity. At the
interface a similar treatment as the one described for
viscosity is done for the conductivity.
Solution Algorithm
The flow chart that summarizes the scheme
of solution is shown in Fig. 2. The input data consists
of four major groups, namely, fluid properties of each
phase (9L, p, k, cp), numerical parameters such as
tolerance and mesh size, problem data such as mass
flows and tube dimensions and initial conditions and
guesses. The solution is achieved using a series of
subroutines as described below.
Subroutine MESH, creates and independent
meSh for each phase with its own number of finite
volumes. The physical interface is located in the east
face of the last volume of the gas core and the west
face of the first volume of the liquid film. MESH
takes the initial guess for the liquid film thickness
and remeshes each time the thickness is corrected.
Subroutine LIQENT uses an empirical
correlation to estimate fraction of liquid entrained
and calculates de properties of the gas core according
to the homogeneous model of multiphase flow. As
well, it determines the liquid film mass flow rate
which is used later.
Subroutine VISCOSITY takes the initial
guesses for pressure gradient and velocity field and,
together with the properties of the mixture recently
calculated, provides turbulent viscosity for both
phases.
Next, the first velocity field is estimated by
CALCVEL subroutine. As described before, a single
velocity field is obtained. As the pressure gradient is
where
ar = Ar + A
pe 1 J
,, t Eg
It is well known that, for flows with high
density ratios, as gas liquid flow, the interface is
wavy and the proposition above must be understood
in an average sense. The interfacial stress calculated
through equation (19), enforces the momentum
conservation at the interface. Nevertheless, the wavy
nature of the interface must be taken into account in
the calculation of gas and liquid effective viscosities
close to the interface using suitable turbulence
models as proposed by (Cioncolini et al. 2009),
which has been used in this work.
The mass flows of both phases are
calculated by the integration of the velocity profiles.
R6
Gc =2g 2xP u(r)rdr
Lf =: u (r)rdr
R6
These relations represent the global mass
conservation for the gas core and liquid film and are
used, together with momentum conservation, in the
solution algorithm for pressure gradient and film
thickness calculation.
As stated before, the effective viscosity is
calculated through the algebraic model proposed by
(Cioncolini et al. 2009), where the effective viscosity
only depends on the liquid film thickness.
#eyy = p 1+0.9x103 + 2 2
du
Ze = #Ue I
Figure 4: Numerical scheme flowchart
Results
Figures 5 and 6 below show the entire
velocity profile and shear stress distribution from the
center to the wall for a particular case.
The experiments by Wolf et al. (2000) were
used for model validation.
Table 1: Experimental mass fluxes kg/m2S
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
assumed to be uniform in each tube section, only the
gas mass flow is checked and the error is used to
correct the pressure gradient. The loop goes back to
the turbulent viscosity calculation feeding the
subroutine with the new values of velocity and
pressure gradient. Once a certain tolerance is
achieved, the liquid film mass flow is checked and
the error used to correct its thickness 8. With this new
predicted value, the mesh has to be redone and the
process is repeated until convergence.
LIQUID FILM
MASS FLUX
GAS ASS LUXTOTAL LIQUID
GAS MASS FLUX
71
71
97
97
124
124
154
154
154
29.8
54.2
10
14.3
10
21
10
12.7
19.7
Mass fluxes *1
gas 97 kg/mzs
liquid 20 kgimzs ,1
B 6
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016
r [mm]
Figure 5: Velocity profile
1E08 0.0000001 0.000001
Figure 6: Shear stress distribution
0.00001 0.0001
Log(R r)
0.001
0.01
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 8: Predicted vs. experimental wall shear
stress
As expected, there is a striking similarity
between Fig. 7 and 8, confirming that wall shear
stress is the main contributor to pressure gradient.
These experiments were conducted in a 10.8
m long tube of 31.8 mm internal diameter, with inlet
pressure kept at a constant value of 2.38 bar. Results
show that the different hydrodynamic parameters of
the annular flow compared in this work, namely,
pressure gradient, wall shear stress and film
thickness, do not take same distance from the inlet to
achieve fully developed conditions, so the set of data
to compare was carefully chosen.
Comparisons are presented below in Fig. 7 
9 for pressure gradient, wall shear stress and film
thickness, respectively, for air mass fluxes from 71 to
154 kg/nrs as shown in Table 1. There is excellent
agreement between predicted and experimental
values, with almost all points lying within an error of
+13%.
There seems to be some error in Fig. 9 in the
film thickness prediction for 71 kg/nrs. However,
according to Wolf et al. (2000), film thickness does
not reach the fully developed value within test tube
length for low mass flow rates, and this could explain
this behavior.
Gs asss
 71
O 97
+124
a154
'i
03
01 02 03
EXPERIMENTAL FILM THICKNESS [mm]
Figure 9: Predicted vs. experimental film thickness
Figures 10 and 11 show, respectively, the
variation of pressure gradient and mean film
thickness with liquid mass, flow rate for two gas
superficial velocities of 71 and 97 kg/m s. The
experimental data of Wolf et al. (2000) for fully
developed cases are also plotted showing very good
agreement, within an error of 17% at the most
1000 2000 3000 4000 5000
EXPERIMENTALPRESSURE GRADIENT[Pa/m]
Figure 7: Pressure gradient predicted vs.
experimental
10 20 30 40
PREDICTED WALL SHEAR STRESS [Pa]
771 kg/m2S
c97 kg/m2S
AWolf et al (2001) 71kg/m2S
eWolf et al (2001) 97 kg/m2S
5 15 25 35 45 55 65 75 85
9 19 29 39 49 59 69 79
LIQUID FILM MASS FLUX [kg/m's]
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1000
3500
3000
z 500
u 1500
v,1000
..500
LIQUID FILM MASS FLUX [kg/m's]
Figure 10: Pressure gradient vs. liquid mass flux
0.4
~71 kg/m2S
~97 kg/m2S
 Wolf et al (2001) 71 kg/m2S
eWolf et al (2001) 97 kg/m2S


0.35
0.3
**0.25
S0.15
0.1
0.05
Figure 11: Delta vs. liquid mass flux
In Fig. 12 the liquid film velocity profile
predicted by this model is compared to the one
provided by the universal velocity profile according
to the equations below:
Both profiles are predicted for a mass flow
rate of 97 and 20 kg/nrs for air and liquid
respectively, for which Wolf et al (2000)
experimental values of liquid film thickness and
pressure gradient show very little change towards the
test tube end (2% or less). These values were
predicted by the model within an error of about 1%
for the pressure gradient and less than 5 % for the
film thickness. With this in mind, it could be
concluded that the film thickness is accurately
predicted by the model. As for the universal velocity
profile, examining Fig. 6, the average velocity
predicted is clearly higher than the one for the
V' = v', O & 5
V' = 51n(y ) 3.05
I' = 2.5b2(y )+5.5
5
v>30
Table 2 shows some values of the pressure gradient
predicted compared with the ones given by the
expression above, using a firstorder approximation
for t .
It is recalled that pressure gradient is
calculated through the global mass conservation for
the gas core and the wall shear stress is calculated
through the liquid velocity gradient at the wall, and
equation (32) is not used in the solution algorithm.
Then, although these phenomena are physically
connected through the interfacial stress, and these
results are expected, they serve as a double check that
the velocity fields and pressure gradient predicted
through the proposed algorithm satisfy local and
global mass and momentum conservation equations.
Table 2: Comparison of pressure gradient values
GAS lASS FLUX TOTALLIQUID dp/dz dp _2cw
MASS FLUrX PRDTE dz R
71 40 1964.1 1964.3
97 20 2403.7 2403.8
124 40 3648.4 3648.5
154 20 4750.9 4751
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A model to predict hydrodynamic and heat
transfer parameters for fully developed vertical,
annular, twophase flow was implemented.
The hydrodynamic model is based on mass
and momentum conservation laws which are locally
and globally satisfied. The algorithm solves the well
known triangular relationship between pressure
gradient, liquid film thickness and liquid mass flow
rate with no empirical closure correlation or wall
laws for the liquid film velocity profile.
The liquid film and gas core effective
viscosities were calculated by algebraic turbulence
model proposed by Cioncolini et al. (2009), showing
excellent performance with little implementation
effort.
A twoequation turbulence model will be
implemented to provide a more general approach and
to be able to represent other situations such as non
developed flows with interfacial mass transfer. This
is object of current research.
Additionally, the validation of the model for
heat transfer problems is under run.
References
Adeclw, D. & Issa, R.I., 2004. Modelling of
annular flow through pipes and Tjunctions.
Computers & Fluids, 33, 289313.
Antal, S.P. et al., 1998. The Development of
Multidimensional Modeling Capabilities for Annular
Flows. In Lyon, France.
Cioncolini, A., Thome, J.R. & Lombardi, C.,
2009. Algebraic turbulence modeling in adiabatic gas
 liquid annular twophase flow. international
Journal of Addtiphase Flow, 35(6), 580596.
Available at:
http://dx.doi.org/10. 1016/j.ijmultiphaseflow.2009.02.
002.
Dobran, F., 1983. Hydrodynamic and Heat
Transfer Analysis of TwoPhase Annular Flow With
a New Liquid Film Model of Turbulence.
international Journal of Heat and Mass Transfer, 26,
11591171.
Fu, F. & Klausner, J.F., 1997. A separated
flow model for predicting twophase pressure drop
and evaporative heat transfer for vertical annular
flow. international Journal of Heat and Fluid Flow,
18(97), 541549.
Ishii, M. & Mishima, K., 1989. Droplet
entrainment correlation in annular twophase flow.
htt. 1 Heat Mass Transfer, 32, 18351846.
universal velocity profile, meaning that the latter
overpredicts film thickness.
O 2 4 6 8 10 12 14 16
D\IMENSIONLESS DISTANCE FROM\ WALL
Figure 12: Dimensionless velocity in the liquid film
Note that the eddy viscosity within the
liquid film in Cioncolini et al. (2009) turbulence
model is constant and so the velocity profile is linear.
A force balance in a slice of the tube
provides the widely used and simple relation between
pressure gradient and wall shear stress:
Conclusion
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Jayanti, S. & Hewitt, G.F., 1997.
Hydrodynamics and Heat Transfer in Wavy Annular
GasLiquid Flow: a Computational Fluid Dynamics
Study. Int. 1 Heat Mass Transfer, 40(10), 2445
2460.
Kishore, B.N. & Jayanti, S., 2004. A
multidimensional model for annular gasliquid flow.
Chem. Eng. Sci., 59, 35773589.
Okawa, T. & Kataoka, I., 2005. Correlations
for the mass transfer rate of droplets in vertical
upward annular flow. International Journal of Heat
and Mass Transfer, 48, 47664778.
Patankar, S.V., 1980. Numerical Heat
Transfer and Fluid Flow, Taylor & Francis.
Peng, S.W., 2008. Heat flux effect on the
droplet entrainment and deposition in annular flow
dryout. Commun. Nonlinear Sci. Numer. Simul., 13,
22232235.
Pu, F., Qiu, S. & Jia, D., 2006. An
investigation of flow characteristics and critical heat
flux in vertical upward round tube. Nuclear Science
and Techniques, 17(3), 170176.
Wolf, A., Jayanti, S. & Hewitt, G.F., 2001.
Flow development in vertical annular flow. Chem.
Eng. Sci., 56, 32213235.
