Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 17.4.1 - Three-dimensional Unit Slug in a Horizontal Pipeline
Full Citation
Permanent Link:
 Material Information
Title: 17.4.1 - Three-dimensional Unit Slug in a Horizontal Pipeline Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Febres, M.
Nieckele, A.O.
Fonseca Jr., R.
Azevedo, L.F.A.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: slug flow
horizontal pipeline
three dimensional
Abstract: Among all possible flow configurations that multiphase flow can take, the slug flow regime plays an important role because of its intermittent characteristics and its strong pressure gradients in pipelines, where this flow regime can be found. Due to its importance, several papers can be found in the literature. Although, the available works provide useful data, there are still missing relevant information of random gas-liquid characteristics, structures evolution and interaction along the pipeline. In the present work, a detailed numerical and experimental study of a slug flow in a horizontal pipeline is performed. The 3-D isothermal flow field was numerically determined with the finite volume method. The Volume of Fluid (VOF) model was employed to determine the interface of a single unit slug, obtained by employing a coordinate system coincident with the Taylor bubble. The turbulence was modeled with the RNG  model. The numerical predictions of the bubble shape, velocity profile and translation slug velocity were compared with the measured data, obtained with the PIV and Laser Induced Fluorescence techniques. Good agreement was obtained between the predicted and measured data.
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: VID00422
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1741-Febres-ICMF2010.pdf

Full Text

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

Three-dimensional Unit Slug in a Horizontal Pipeline

M. Febres*, A.O. Nieckele*, R. Fonseca Jr.t and L.F.A. Azevedo*

Department of Mechanical Engineering, PUC-Rio, Rio de Janeiro, RJ, 22453-900, BRAZIL
Petrobras, R&D Center, Rio de Janeiro, RJ, 21941-915, BRAZIL
mijail febres, nieckele, robertofonseca, lfaa@

Keywords: Slug flow; horizontal pipeline; three dimensional, VOF; PIV


Among all possible flow configurations that multiphase flow can take, the slug flow regime plays an important role because
of its intermittent characteristics and its strong pressure gradients in pipelines, where this flow regime can be found. Due to its
importance, several papers can be found in the literature. Although, the available works provide useful data, there are still
missing relevant information of random gas-liquid characteristics, structures evolution and interaction along the pipeline. In the
present work, a detailed numerical and experimental study of a slug flow in a horizontal pipeline is performed. The 3-D
isothermal flow field was numerically determined with the finite volume method. The Volume of Fluid (VOF) model was
employed to determine the interface of a single unit slug, obtained by employing a coordinate system coincident with the Taylor
bubble. The turbulence was modeled with the RNG K-s model. The numerical predictions of the bubble shape, velocity profile
and translation slug velocity were compared with the measured data, obtained with the PIV and Laser Induced Fluorescence
techniques. Good agreement was obtained between the predicted and measured data.


Multiphase flow is very common in the petroleum industry.
Not only different fluids, such as oil and water, can be found,
but these fluids can be present in different phases (liquid or
gas). Further, during its transport, depending on the phase's
superficial velocities, different flow patterns can be found,
such as sinooli stratified", "wavy stratified", "elongated
bubble", "slug", "annular", "chaotic", "bubble", etc.
Numerical simulation has became a very important tool to
predict multiphase flow inside pipelines due to the worldwide
increasing fuel demand for industry, as well as higher
standard requirements for equipment design for petroleum
processing and transport. For an efficient equipment or pipe
design, it is desired to determine the individual hydrodynamic
characteristics of each flow pattern, such as pressure drop and
flow intermittence parameters.
One of the most common flow patterns during operation
in Petroleum processing facilities is the "slug" pattern. There
have been many attempts in the literature to model a
two-phase slug flow. One of the first models, presented by
Wallis (1969), is based on the concept of the single unit slug
with a reference frame moving with the tip of the bubble.
Later, Taitel and Barnea (1989) using the same concept,
divided the single unit slug into two parts: the liquid slug and
the liquid film zones. Fagundes Netto et al. (1999) developed
a model based on mass and momentum conservation
equations to predict the shape of the tip and tail of the gas
bubble, as well as the liquid slug, and presented a
comparison with experimental data.
Taitel et al. (2000) showed that a negative pipe
inclination damps slug formation. Orell et al. (2'i 14) used a

sub-model of Taitel and Barnea (1990) model, and obtained
an increase of the pressure loss, obtaining good results in
agreement with experimental data. De Freitas et al. (2008)
applied the single unit slug concept through a vertical pipe
and studied the effect of gas expansion. Several codes have
been developed and are widely used, such as PLAC (Black et
al., 1990) and OLGA (Bendiksen et al., 1991).
A literature search reveals a large number of papers
presenting experimental investigations of the global
characteristics of slug flows. These experiments were
conducted to support the development and verification of
slug flow models in the horizontal configuration. In these
studies, the focus is on global measurements such as
pressure drop, overall void fraction and statistical parameters
such as slug length distribution, slug frequency and film
thickness (e.g., Bendiksen, 1984, Cook and Behnia, 2000,
Bertola, 2002). Probably due to the complex nature of the
slug flow pattern, relatively few papers were found reporting
results on detailed measurements of the flow field in the
liquid slug and film layer under the gas bubble. Detailed
information on the flow behavior is critical to the proper
understanding of the physical mechanisms governing the
Instantaneous full field measurements of the liquid
phase in slug flow were reported by Carpintero-Rogero et
al. (2006). The authors employed three simultaneous optical
techniques to produce time-averaged and turbulent data for
the liquid slug and film layer: PIV (Particle Image
Velocimetry), LIF (Laser Induced Fluorescence) and PST
(Pulsed Shadow Technique). These same techniques were
employed in the present paper.
Recently, Wang et al. (2007) experimentally showed

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

that the translational velocity varies with the distance from
the entrance. It was also shown that the slug frequency is a
function of the liquid superficial velocity.
Although, the available works provide useful data, there
are still missing relevant information of random gas-liquid
characteristics, structures evolution and interaction along the
Among the several techniques to numerically predict the
slug pattern, the Two Fluid Model, introduced by Ishii
(1975), is a very successful one, since it does not require the
introduction of transition criteria between flow patterns.
However, it does require the definition of interface transfers.
Issa e Kempf (2003) presented a one-dimensional
formulation based on the two fluid model to predict the slug
formation and Ortega and Nieckele (2005) investigated the
model performance with a conservative and non-conservative
Detailed slug flow information can be obtained with 2D
or 3D simulations, however those are very expensive. Thus,
to reduce the cost, most available works consider a single
slug unit. Cook et al. (2001) investigated an air bubble in
stagnant water using VOF, in symmetrical 3D domain with
several inclination angles (50 to 750), and obtained good
agreement with experimental data. Taha Taha et al. (2006)
and Ujang et al (2008) applied the VOF model to simulate
the slug flow with a reference frame moving with the bubble.
While the former considered a symmetrical 2-D and 3-D
vertical pipe, the later investigated a horizontal slightly
inclined slug flow. Ujang et al. (2008) showed a tendency of
the slug nose to move in direction of the pipe center line.
Vall6e et al (2008) using ANSYS CFX, simulated
air/water stratified and slug flow in a square pipe, obtaining
good qualitative results.
In the present work, the commercial CFD software
FLUENTTM was employed to solve numerically the
air/water two-phase slug flow with the single unit slug
approach. The slug translational velocity and film velocity
profiles were compared with measured data, presenting good


C empirical constant
D diameter (m)
F external force per unit volume (N/m3)
Fr Froude number
g gravitational constant (m/s)
G turbulent generation term (Pa/s)
L length (m)
n normal direction
p pressure (N im-)
S mean deformation rate (1/s)
u velocity vector (m/s)
V vertical velocity (m/s)
W axial velocity (m/s)

Greek letters
a volume fraction
e turbulent dissipation rate (m2/s3)
ic turbulent kinetic energy (n11- -
2 inverse of Prandtl number
pu viscosity (Pa s)
91 interface curvature (1/m)




specific mass (kg/m3)
surface tension (N/m)

superficial velocity

Mathematical Model

To determine the two-phase slug flow field inside a
horizontal pipeline, the VOF (volume of fluid) model (Hirt
and Nichols, 1981) was selected. The VOF model is based
on the solution of one single set of conservation equations of
mass and momentum. An auxiliary variable, named "volume
fraction" a,, is considered to identify the region occupied by
each phase. All variables and properties fields are shared by
both phases, and they represent average values. Properties
and variables on each cell are entirely representative of one
phase or a mixture of the phases depending on the value of
volume fraction, which are defined as follows: a, = 0 if the
cell does not contain fluid "i", a, = 1 if the cell is full of fluid
"i" and 0 < i < 1 if the cell contains the interface between the
phases. The sum of volume fractions is equal to 1, i.e.,
ag+ of =1.
Properties, like density p and viscosity p, which appear
in transport equations, are evaluated by taking in account the
amount of each phase in the cell, as follows: p =ag pg + cf pf;
S=ag p + a/ pa.
The interface tracking is determined by assuming a
material derivative of the interface equal to zero for a
referential on the interface, thus
a + iVag =0 (1)

where ii is the time average velocity vector.
The Reynolds average continuity and momentum
conservation equations can be written as

-(pu)+V (pu )=0 (2)

S(pi)+ V i u)= -Vp +p pg +
8t (3)
+V [/[ef 2S]+p +F

where p is the pressure, g is the gravity acceleration
vector, ef = /p + ut is the effective viscosity, p, is the
turbulent viscosity, S is the mean rate of strain tensor and

F is an external force which takes in account the effects of
surface tension (o-g) of two phases and it is based on the CSF
model (Continuum Surface Force) developed by Brackbill et
al. (1992)

S= (V +V7 ); F P= -ggp V- (4)
h 2 (pg + /2

where the interface curvature 91 = V i is determined as a

function of the volume fraction gradient in cells that contains
the interface (R = n / i Iwhere ii = Vag).
Turbulence was taken into account with the K-S RNG
model (Yakhot et al., 1992). For high Reynolds number, the
turbulent viscosity is

Pt = C/ p 2 / ; C, = 0.0845 (5)

where i is the turbulent kinetic energy and c its dissipation
rate, and are obtained from the following transport

t-(o -)+V(p A)= V[I PefV I] + G, pC

C-(e)+V.(* iE)= V[ae Pef VI]+

2 (7)
+ C, -G, C26P -
In these equations, G, represents the generation of
turbulence kinetic energy due to the mean velocity gradients

G, = S2 ; S

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

which should be sufficient for the formation of stable slugs.
Water from a reservoir was pumped in closed circuit
through the test pipe by a progressive cavity pump. A
centrifugal blower provided compressed air for the test
section. Calibrated rotameters were used to measure the
water and air flow rates. Air and water were mixed at a
Y-junction positioned at the entrance of the Plexiglas pipe.
After passing through the test pipe, the two-phase mixture
returned to the reservoir where a tangential inlet aided the
phase separation process. The green and blue arrows in
Fig. la indicate, respectively, the water and air flow paths
in the test section.

2 S, S,

.k and ), are the inverse Prandtl numbers for K and
c respectively, obtained from the following expressions:

(a) overview of test section

,8- -1.3929 06321
,o 1.3929

.,, + 2.3929 03679
o + 2.3929

The destruction of c depends on C, = C2, + R_ where

Cp r73(1-77/77) A:
R, = ; = S -
l+/ 77 C

The empirical constants are: Cg = 1.42, C, = 1.68, o = 1,
/= 0.012 and 7o = 4.38.

Boundary conditions: To predict a single unit cell, a
reference frame moving with the bubble was employed.
Therefore, at the pipe wall a slug translational velocity was
prescribed. At the entrance of the domain, the mixture
velocity was imposed as Wl=wsf+Wsg, where ws, and wsg are
the liquid and gas superficial velocities, relative to the wall.
The entrance turbulent quantities were

K- (W. )2 _=C 3/4 3/2/
2 P


with .=0.05 and I, =0.07 D. At the exit a constant pressure
was prescribed.

Experimental Facility

Figure 1 schematically presents the test section constructed
to conduct the experiments. A 24mm diameter, 10 meters
long Plexiglas pipe was mounted on a rigid steel frame that
could be rotated around a pivot to produce inclination
angles between 0 and +100 with the horizontal. In the
present paper only results for the horizontal cases will be
presented. The pipe length-to-diameter ratio was 400,

(b) Components of the measuring system
Figure 1: Experimental apparatus.

The measuring section was located at 350 diameters
from the pipe entrance. As it can be seen schematically in, the measuring section was specially prepared to
receive the components necessary for the implementation of
the three optical techniques employed in the experiments,
namely, PIV, LIF and PST. Light from a double cavity,
120-mJ, Nd-YAG laser manufactured by New Wave was
shaped into a plane sheet and directed to the meridional
vertical plane of the pipe. A set of cylindrical and spherical
lenses was used to form the light sheet that presented
approximate dimensions in the measuring region of 50 x 0.5
mm (width x thickness). A 450 mirror was also employed to
divert the laser beam that entered the test section from
below, as indicated in the figure. In order to minimize
reflections from the laser light at the curved pipe surface, a
rectangular glass box was mounted around the pipe at the
measuring area. The box was filled with water. The
maximum firing frequency of each laser cavity was 15 Hz.
The dual cavity laser system allows for very short pulse



intervals between laser firing. Typical pulse intervals of the
order of micro seconds can be easily achieved, which is
adequate for measuring liquid velocities significantly
higher than the levels encountered in the experiments
conducted. The laser was operated to emit green light at a
wavelength of 532 nm. The particle images were acquired
by a PIVCAM 10-30 digital camera manufactured by TSI.
The camera can be operated in the frame straddling mode,
which allows for the capture of pairs of consecutive images
with time intervals in the micro second range. Pairs of
consecutive images can be acquired by the camera at 15 Hz,
which is compatible with the laser firing frequency. The
camera offered a 1000 x 1000 pixel resolution. As indicated
in the figure, the camera was mounted orthogonally to the
light sheet plane. A 50-mm focal distance lens was used in
the experiments. The panel of LED's for the background
illumination required by the PST technique was mounted on
the opposite side of the camera. A sand-blasted diffuser
glass plate was positioned in front of the LED panel to
improve the spatial uniformity of the background
illumination. The high-pass optical filter required for the
LIF technique was fixed on the rectangular glass box in the
viewing path of the camera. In the present study seeding
was obtained with the use of 15-pm-diameter fluorescent
particles with peak excitation wavelength at 542 nm and
peak fluorescence emission at 612 nm. The particles density
was 1.05 g/cm3.
A photogate cell was installed around the pipe, upstream
of the glass box. This cell (not shown in the figure)
provided trigger signals for the velocity measuring system.
The cell output a high voltage signal when gas passes
through the cell, and a low voltage signal when water is
flowing through. The transition from a high to low voltage
signal is an indication of the passage of a gas-liquid
interface at the photogate position. A controlling circuit
triggered by the photogate signal outputs a trigger signal to
the velocity measuring system after a pre-set time delay. By
controlling the time delay, it was possible to always
measure the liquid velocity field at a pre-determined
position in relation to the gas bubble surface position. A
second photogate cell was installed at a known position
from the first one and was employed in the measurement of
statistical slug properties, such as length, frequency and
velocity distributions.

Numerical Scheme

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

Computational domain. The mesh of the circular pipe was
generated with GAMBITTM. Solution was obtained with
5100, 29000 and 232.000 control volumes in the cross
section, and the intermediate mesh illustrated in Fig. 2a was
selected, since the difference in the velocity profile was
inferior a 0.5%. Similar test was performed for the axial
position and 100 control volumes were selected with the
mesh concentrated in the region indicated in the Fig. 2c,
with length equal to 2D. At this region, where the tip of the
bubble was located, the mesh was approximated uniform in
the axial direction (Fig. 2b).

(a) cross section

-ID + ID
(b) detail at bubble tip.

(c) transversal section
Figure 2: Mesh distribution.

Initialization. The initial condition was defined according to
Ujang et al. (2008), in which water volume fraction was set
as "1" at first half of the pipe and equal to "0" at the second
half (Fig.3).

Wall Velocty W

Air Initial Interface
Figure 3: Boundary conditions and domain initialization.

The wall velocity was initialized according to
Bendiksen (1984) models to estimate the slug translational

WT =Co Wm +Wi,,,,

To solve the mathematical model, the commercial software
FLUENTTM was employed. It is based on the Finite Volume
technique which consists on dividing the computational
domain on small control volumes and integrating spatially
and temporarily the transport equations over them. First
order implicit technique was used for the time integration,
while the spatial integration was handled with the second
order "QUICK" scheme (Leonard, 1979).
The pressure-velocity coupling was solved with the
PISO algorithm (Issa, 1986) and PRESTO was selected as
the interpolation scheme for pressure.
To reconstruct the interface of the VOF methodology,
the "Geometric Reconstruction" technique was applied.
The resulting discretization algebraic equations were
solved using the AMG (Algebraic Aiih,, ,. Solver) option
(Hutchinson et al., 1986).

where Co stands for the distribution parameter of

Co =1.05 >Frc <3.5
C= 1.20F 5 Frc
Co = 1.20 Frc 2 3.5

WS/ ; C

0.54 (13)

g is the gravity acceleration, D the pipe diameter and Frc
stands for the critic Froude number.
The wall velocity obtained with the correlation given by
Eq. (12) and coefficients of Eq. (13) was a good estimated,
but it was corrected based on the flow field, so that the
bubble position was maintained fixed inside the
computational domain. Usually three iterations were
necessary to find a still bubble.
Figure 4 illustrates the bubble position inside the
domain at different time instants. Figure 4a corresponds to
the initial wall velocity prescribed from Eqs. (12) and (13),

W Ci -..Fg

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

and Fig. 4b illustrates the bubble position after the
converged wall velocity (translation bubble velocity) was

St= 0 20s

-t 0.40s
St= 0.60s


(a) Wall velocity. Eq. (12) (b) Converged wall velocity
Co=1.05, C,=0.54
Figure 4: Time evolution of Taylor bubble.

Results and Discussion

It was considered a pipeline with the same diameter as in the
experimental apparatus (D = 0.024m). The length was set as
L 25D, long enough to capture one single bubble and a
fully developed liquid slug). The working fluids were air and
water (air: pg=1.225 kg/m3, /g= 1.790x10-5 Pa s; water:
p,=998.2 kg/m3, p= 1.003 103 Pa s, 0g=0.073 N/m).
Six combinations of gas and liquid superficial velocities
were examined. Table 1 presents the superficial velocity of
each phase, the resulting mixture velocity, and the
superficial velocities ratio.
The initial slug translational velocity (Eqs. 12 and 13)
and the converged wall velocity were also included in Table
1. For all cases tested Frc was smaller than 3.5, and the
constants from Eq. (13) were: Co =1.05 and C1=0.54. The
slug translation velocity and the converged wall velocity
presented a variation from 5% to 11%. Therefore, a linear
regression of the numerical data for the converged slug
translation velocity was performed and the coefficients Co
and Ci of Eq. (12) were re-evaluated. The corrected
coefficients based on the numerical solution were:
Co =1.116 and C1=0.202 (Wdr,,=0.098 m/s). These values
are within the range observed experimentally by Bendiksen
(1984) for a pipe with 0.0242 m diameter (1.009 < Co <
1.188 and -11" l4111 s < Wdnft < 0.181m/s), indicating that the
present VOF solution was capable of predicting drift
velocity of a slug two phase flow.

Table 1: Test cases.
Cases w g ws Wm ws, Wr w
(m/s) (m/s) (m/s) wg/i (m/s)
1 0.475 0.295 0.770 1.610 1.071 0.956
2 0.475 0.393 0.868 1.209 1.174 1.067
3 0.475 0.516 0.991 0.921 1.303 1.205
4 0.788 0.295 1.083 2.671 1.399 1.308
5 0.788 0.393 1.181 2.005 1.502 1.416
6 0.788 0.516 1.304 1.527 1.631 1.552

The six cases selected to be tested were located inside a
pattern map (Fig. 5) where it can be seen that three of them
are in the elongated bubble region and the other three are at
the interface of the elongated bubble region with the slug

1 -s slug
elongated bubble l

0.1 anular

stratifed s
0.01 0.1 1 10
Figure 5: Map pattern

Bubble shape. To illustrate the bubble nose shape, Case 6
was selected. The volume fraction surface equal to zero is
plotted on top of the mesh (Fig. 6b) and it is compared with
the image obtained experimentally with the combination of
laser induced florescent and synchronized background
illumination (Fig. 6a). A similar configuration is obtained,
where it can be seen that the bubble nose is oriented
downward, agreeing with previous experimental
observation of Bendiksen (1984). However, the numerical
inclination is quite small. Further, the numerical simulation
predicted a more round nose than observed experimentally.

a) Experimental. b) Numerical
Figure 6: Bubble nose shape. Case 6

To visualize the bubble shape, it is shown in Fig. 7,
contours of the volume fraction at five cross section for
Case 1. The bubble presents an oval cross section, which
gets flatter along the bubble length. Note that the bubble is
squeezed and pushed to the upper part of the pipeline. It can
be also be seen that at z 0.8 D, the bubble touches the
wall and the liquid film above it disappears.

z=0.OD z=-0.2D z=-0.4D z=-0.6Dz=-0.8D
Figure: Volume fraction contours at different cross
section. Case 1

Velocity field. Figure 8 illustrates a comparison of the
measured velocity field with the numerical prediction, for
Case 1 and 6, where a reasonable agreement can be
observed. It can be seen the inclination of the bubble
towards the pipeline center, inducing a maximum velocity
dislocated to the center in front of the bubble.

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

a) Experimental. Case 1 b) Numerical. Case 1

a) Experimental. Case 6 b) Numerical. Case 6
Figure 8: Velocity field. Case 1 and Case 6

Axial velocity profiles were obtained along nine vertical
lines located at coordinates from -0.8D to 0.8D relative to
tip of the bubble (0.OD), as shown in Fig. 9. The
coordinates with negative signal will be called from now on
"upstream" and coordinates with positive signal will be
called "downstream".

(a)z= 0.0 D (b) z= 0.2 D
1 10

01 -100
00 03 06 09 12 00 03 06 09 12
W(m/s) W(m/s)

(c) z= 0.4 D

(d) z= 0.6 D

Figure 10: Axial velocity profile, upstream bubble nose.
Case 2

-0.8D -).6D -1.4D 4).2D (0.0D

0.2D o.4D 0.6D 0.8D

I I I -I I

I Ii

Figure 9: Data Acquisition positions


(a) z = 0.0 D

(b) z= 0.2 D

W(m/s) W(m/s)

(c) z= 0.4 D

Detailed results of selected cases are presented here.
Figures 10 through 12 show the velocity profiles upstream
of the bubble nose, corresponding to the selected Cases 2, 3
and 6, respectively, Figure 13 corresponds to the
downstream positions of the three cases. Experimental data
is presented with symbols while numerical results are
plotted with continuous lines.
Analyzing Figs. 10-12 it can be seen that there is no
experimental data in the gaseous phase, because tracing
particles were introduce only in the liquid phase. Due to the
smaller viscosity, the gas velocity presents a smaller
inclination at the wall, leading to much higher velocities
than the liquid phase. An excellent agreement between
numerical and experimental data was obtained in the liquid
region. However, it can be observed a local reduction on the
experimental liquid velocity profile near to the interface.
This behavior was not predicted numerically, where a
smooth transition from the liquid velocity to the gaseous
velocity profile is was obtained.
Analyzing Figs. 10-12, it can be seen that the maximum
gas velocity increases as one moves upstream, while the
maximum liquid velocity diminishes, agreeing with the
experimental data.

Figure 11: Axial velocity profile, upstream bubble nose.
Case 3


(a)z = 0.0 D

00 03 06 09 12 15 18
(c) z= 0.4 D


(b) z= 0.2 D

00 03 06 09 12 15 18
(d) z= 0.6 D

Figure 12: Axial velocity profile, upstream bubble nose.
Case 6

(d) z= 0.6 D


The agreement of the numerical results downstream of
the Taylor bubble is excellent, as it can be seen in Fig. 13 for
the three selected cases. According to several researches like
Taitel and Barnea (1990), Polonsky et al (1999) and Shemer
(2003), velocity profiles in the liquid slug ahead of the
bubble must be fully developed. For all cases tested, this
behavior was observed. This information is reinforced in Fig.
14, where velocity profiles for Case 1 corresponding to the
tip position and downstream positions from +0.2 D through
+0.8 D are plotted in Fig. 14a. It can be seen that at positions
equal to +0.4D, +0.6D and +0.8D the velocity profiles are
superimposed indicating a fully developed behavior for
velocity ahead of the bubble. Fig. 14b compares the numerical,
with experimental and the 1/7 law of fully developed velocity
profile, where a good agreement can also be seen.

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

again the numerical model predicted very well the velocity,
even though the vertical velocity components are quite
small, and more difficult to measure.

-010 000 010 -010 000 010

-010 000






(a) Case 4

-0 5

-10 C
00 03 06 09

(a) Case 2: z= + 0.2 D




00 03 06 09

(c) Case 3: z= + 0.2 D




-05 -

00 03 06 09 12

(e) Case 6: z= + 0.2 D

-0 08 000 0 08

-0 08 000 0 08

00 03 06 09

(b) Case 2: z= + 0.6 D


(d) Case 4: z= + 0.6 D


(b) Case 6

Figure 15: Vertical velocity profile, upstream bubble nose

-020 000 020

15 00 03 06 09 12 1

(f) Case 6: z= + 0.6 D


Figure 13: Axial velocity profile, downstream bubble nose.
101 --I -1 .--



00 03 06


04 06 08

(a) Case 1: several coordinates (b) Case 1: z= + 0.8 D

Figure 14: Axial velocity profile, downstream bubble nose.

A comparison between the numerical and experimental
velocity profile for the vertical velocity component can be
appreciated at Figs. 15 and 16, for Cases 4 and 6, for the
upstream and downstream positions, respectively. Once

-020 000 020

-020 000 020


1 00-
0,OD +0,2D +0,4D


(a) Case 4

-020 000 020 -020
-020 000 0 0
1 00-

1 00
0,OD +0,20

-020 000 020

-020 000 020


000 020
-0120 0

10 0 O

000 020


(b) Case 6

Figure 16: Vertical velocity profile, downstream bubble nose

An interesting result to analyze is the velocity field at
the pipe cross-section, in the liquid and gas region. Figure
17 illustrates the cross-section velocity vector distribution
for Case 1. It can be seen a strong recirculation in the
gaseous phase that begins on the interface, climbs to the top
of the pipe and returns on the sides. This recirculation is
probably responsible for the increase in the bubble diameter
as one move upstream of the bubble tip, as shown in Fig. 7.

(a) z= 0.2 D

(b) z= -0.4 D

....... ...;~i~
r%:,~..........~,,;s -
i ;---~----~---~ i'
1_-_---\~~~1 1111~--_~~
_-II I I I I I I-_
1 I ( I I I I
~ ( 1 I I I I ( ) )
1I 1 I

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

Case 1, 4 and 5 presented an error around 15%.
The vertical component is much smaller than the axial
component, thus, a more approximate definition of the error
was introduced, as

Z Vexp num
Error(%) = N V -, x 100 (15)
Vexpmnax exPminn

where Vexpx and Vexpm are the maximum and
minimum experimental vertical velocities. As expected the
relative error for these components were much higher. The
average error was around 20% for all cases.
One possible explanation of the better performance of
Cases 2, 3 and 6 can be explained by ratio of the superficial
velocities ratio, since these three cases presented the
smallest ratio Wsg w, This means that if the gas superficial
velocity is much higher than the liquid velocity, the
performance of the VOF model deteriorates.
Another reason of the discrepancies between the results
(although small) can also be attributed to an increase in the
dissipation at the interface, which was not modeled, that
reduced the liquid velocity at that region, as shown in Figs.


Numerical simulation of a single unit slug was performed
with the VOF model. The slug translation velocity was
obtained and the coefficients of the Bendiksen equation
were adjusted to fit all cases. The computed velocity was
compared with the experimental data obtained with PIV and
Laser Doppler anemometer. It can be said that a very good
agreement was obtained.


The authors thank CNPq and Petrobras for supporting the
development of this work.


(c) z= 0.6 D (d) z= 0.8 D
Figure 17: Velocity field at the cross-section, upstream of
bubble nose.

To quantify the difference between the velocities
measured and numerically determined, a relative error was
calculated as

S -exp W num
Error(%)= N x 100 (14)
N Wm

where N is the number of data, Wxp e Wnum, are
experimental and numerical velocities, respectively. The
velocity difference was normalized by the mixture velocity
Wm. Numerical results are available on each computational
nodes and do not necessarily coincide with position where
experimental data was measured, therefore, to estimate the
error, the numerical data was interpolated for the position
where the experimental data was available.
Cases 2, 3 and 6 presented the smallest error for the
axial component. At all nine axial positions indicated in Fig.
9, the average error was approximately equal to 5%, while

Barnea, D. & Taitel, Y. Transient Formulation Modes and
Stability of Steady-State Annular Flow. Chemical
Engineering Science, Vol. 44(2), 325-332 (1989).

Bendiksen, K. H. An Experimental Investigation o f the
Motion of long bubbles in inclined pipes. International
Journal of Multiphase Flow 10(4), 467-483 (1984).

Bendiksen, K.H.; Malnes, D.; Moe, R.; Nuland, S. The
dynamic two-fluid model OLGA: theory and application.
SPE Prod. Eng., Vol. 6, 171-180 (1991).

Bertola, V Slug velocity profiles in horizontal gas-liquid
flow, Experiments in Fluids, Vol. 32, pp. 722-727 (2002).

Black, P. S., Daniels, L. C., Hoyle, N. C., Jepson W. P
Studying Transient Multi-Phase Flow Using the Pipeline
Analysis Code (PLAC). Journal Energy Resource Technol.,
Vol. 112(1), 25 (1990)

Brackbill, J.U., Kothe, D.B., Zemach, C. A continuum
method for modeling surface tension. Journal of

Computational Physics, Vol. 100(2), 335-354 (1992).

Carpintero-Rogero, E., Kroess, B., Sattelmayer, T.
Simultaneous HS-PIV and shadowgraph measurements of
gas-liquid flows in a horizontal pipe. 13th Int. Symp. on
Applications of Laser Techniques to Fluid Mechanics

Cook, M. and Behnia, M. Slug length prediction in near
horizontal gas-liquid intermittent flow. Chemical
Engineering Science, Vol. 55, 2009-2018 (2000).

Cook, M. & Behnia, M. Bubble motion during inclined
intermittent flow. International Journal of Heat and Fluid
Flow, Vol. 22, 543-551 (2001).

Fagundes Netto, J.R., Fabre, J., P6resson, L. Shape of long
bubbles in horizontal slug flow. International Journal of
Multiphase Flow, Vol. 25, 1129-1160 (1999).

Hutchinson, B. R. and RAITHBY, G.D. A Multigrid
Method Based on the Additive Correction Strategy.
Numerical Heat Transfer 9, 511-537, (1986).

Issa, R.I. Solution of the implicitly discretized fluid flow
equations by operator-splitting. J. Comput. Phys. 62, 40-65,

Ishii, M. Thermo-Fluid Dynamic Theory of Two-Phase
Flow. Collection de la Direction des Etudes et Researches
d'Electricite de France. Eyrolles, Paris, France (1975).

Issa R.I. & Kempf M.H.W. Simulation of slug flow in
horizontal and nearly horizontal pipes with the two-fluid
model. International Journal of Multiphase Flow, Vol. 29,
69-95 (2003).

Leonard, B. P. A stable and accurate convective modeling
procedure based on quadratic upstream interpolation.
Comput. Method Appl. Mech. Eng. 19, 59-98 (1979).

Orell A. Experimental Validation of a simple model
gas-liquid slug flow in horizontal pipes. Chemical
Engineering Science, Vol. 60 (2205) 1371-1381 (2"'111

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

Ortega, A. J. & Nieckele, A. O.Simulation of horizontal
two-phase slug flows using the two-fluid model with a
conservative and non-conservative formulation,
Proceedings of COBEM Congresso Brasileiro de
Engenharia MecAnica. Ouro Preto, MG, Brasil (2005)

Polonsky, S., Shemer, L., Barnea, D. The relation between
the Taylor bubble motion and the velocity field ahead of it.
International Journal of Multiphase Flow 25, 957-975

Shemer, L. Hydrodynamic and statistical parameters of slug
flow. International Journal of Heat and Fluid Flow 24,
334-344 (2003).

Taitel, Y & Barnea, D. Two-phase slug flow. Adv. Heat
Transfer, Vol. 20, 83-132 (1990)

Taitel, Y, Sarica, C., Brill, J.P. Slug flow modeling for
downward inclined pipe flow: theoretical considerations.
International Journal of Multiphase Flow, Vol. 26, 833-844

Taha Taha, Z.F. Cui.CFD Modelling of slug flow in vertical
tubes. Chemical Engineering Science, Vol. 61, 676-687

Ujang, P.M., Pan, L., Manfield, P.D., Lawrence, C.J. and
Hewitt, G.F. Prediction of the translational velocity of
liquid slugs in gas-liquid slug flow using computational
fluid dynamics. Multiphase Science and Technology, Vol.
20(1), 25-79 (2008)

Vall6e, C., HOhne, T., Prasser, H., Siihnel, T. Experimental
investigation and CFD simulation of Horizontal stratified
two-phase flow Phenomena. Nuclear Engineering and
Design, Vol. 238, 637-646 (2008)

Wallis, G.B., 1969. One-dimensional Two-phase Flow,
McGraw-Hill, Inc, New York.

Wang X., Guo L. and Zhang X., An experimental study of
the statistical parameters of gas-liquid two-phase slug flow
in horizontal pipeline. International Journal of Heat and
Mass Transfer, Vol. 50, 2439-2443 (2007).

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