Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: P1.63 - Finite Volume Modeling of Fully Developed Vertical Annular Flow
Full Citation
Permanent Link:
 Material Information
Title: P1.63 - Finite Volume Modeling of Fully Developed Vertical Annular Flow Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Fernández, F.M.
Alvarez Toledo, A.
Paladino, E.E.
Lima, J.A.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: annular flow
finite volume method
turbulence modeling
Abstract: This work presents a numerical model for predicting hydro-dynamic and heat transfer parameters for turbulent, annular, vertical, gas-liquid two-phase 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.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00459
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: P163-Fernandez-ICMF2010.pdf

Full Text

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
59072-970, Brazil


This work presents a numerical model for predicting lwdro-dynamic and heat transfer parameters for turbulent,
annular, vertical, gas-liquid two-phase 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

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)

g Gas
gc Gas core
i Interface
1 Liquid
If Liquid film
w Wall


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 sub-layers: a
continuous liquid layer next to the wall and a wavy
one next to the liquid-gas interface: the velocity
profile within both layers is obtained from law-of-
the-wall 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
The gas core is considered to have liquid
droplets entrained, and in view of assumption (4) and
the homogeneous model its properties are calculated

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 droplet-laden 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
two-equation 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
fully-developed 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

(1) Fully-developed, steady-state, co-current,
vertical, annular flow.
(2) Flow is assumed to be axis-symmetric and film
thickness is assumed to be constant in angular
(3) Entrainment and deposition processes have
reached fully-developed 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 thermo-dynamical and transport
properties are assumed not to vary with temperature
or pressure, even when the thermal entry problem is

(1- xg) xg

9U,- = (1 agc). pI + age AUg

Figure 1: Problem description

where a,, is the mass fraction
calculated as:

of the droplets,

go PI Pix g(1-xg)

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 time-averaged linear momentum and
energy equations for hydro-dynamically fully
developed flow are:

Id du,
r dr dr

1 d du,,
r dr dr



8 (PI cp, u Ty) 1~( 8 d -8Tl
8z r Br B r

8z r Br Br

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


This is possible because a non-uniform
effective viscosity is used, which is a function of the
radial position. So dynamic viscosities can be written

u = r = R




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 OfR-3 R-3 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:



Arg arl
Figure 3: Finite volumes at the interface

This correlation was fitted by the authors
from a large database for gas-liquid flow.
Regarding the gas core effective viscosity, it
varies linearly with wall distance and inversely
proportional to a constant a.

Pef-gc Pgc-i

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 re-meshes each time the thickness is corrected.
Subroutine LIQ-ENT 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
Next, the first velocity field is estimated by
CALCVEL subroutine. As described before, a single
velocity field is obtained. As the pressure gradient is


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.

Gc =2g 2xP u(r)r-dr

Lf =: u (r)-rdr

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.

#ey-y = p 1+0.9x10-3 + 2 2

Ze = #Ue I

Figure 4: Numerical scheme flowchart


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.





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

1E-08 0.0000001 0.000001

Figure 6: Shear stress distribution

0.00001 0.0001
Log(R r)



7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

Figure 8: Predicted vs. experimental wall shear

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
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



01 02 03

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

Figure 7: Pressure gradient predicted vs.

10 20 30 40

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


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010




z 500

u 1500




Figure 10: Pressure gradient vs. liquid mass flux


~71 kg/m2S
~97 kg/m2S
- Wolf et al (2001) 71 kg/m2S
eWolf et al (2001) 97 kg/m2S









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


Table 2 shows some values of the pressure gradient
predicted compared with the ones given by the
expression above, using a first-order 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

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, two-phase 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
A two-equation 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.


Adeclw, D. & Issa, R.I., 2004. Modelling of
annular flow through pipes and T-junctions.
Computers & Fluids, 33, 289-313.

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 two-phase flow. international
Journal of Addtiphase Flow, 35(6), 580-596.
Available at: 1016/j.ijmultiphaseflow.2009.02.

Dobran, F., 1983. Hydrodynamic and Heat
Transfer Analysis of Two-Phase Annular Flow With
a New Liquid Film Model of Turbulence.
international Journal of Heat and Mass Transfer, 26,

Fu, F. & Klausner, J.F., 1997. A separated
flow model for predicting two-phase pressure drop
and evaporative heat transfer for vertical annular
flow. international Journal of Heat and Fluid Flow,
18(97), 541-549.

Ishii, M. & Mishima, K., 1989. Droplet
entrainment correlation in annular two-phase flow.
htt. 1 Heat Mass Transfer, 32, 1835-1846.

universal velocity profile, meaning that the latter
over-predicts film thickness.

O 2 4 6 8 10 12 14 16

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:


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
Gas-Liquid Flow: a Computational Fluid Dynamics
Study. Int. 1 Heat Mass Transfer, 40(10), 2445-

Kishore, B.N. & Jayanti, S., 2004. A
multidimensional model for annular gas-liquid flow.
Chem. Eng. Sci., 59, 3577-3589.

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, 4766-4778.

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,

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), 170-176.

Wolf, A., Jayanti, S. & Hewitt, G.F., 2001.
Flow development in vertical annular flow. Chem.
Eng. Sci., 56, 3221-3235.

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs