7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Threedimensional Unit Slug in a Horizontal Pipeline
M. Febres*, A.O. Nieckele*, R. Fonseca Jr.t and L.F.A. Azevedo*
Department of Mechanical Engineering, PUCRio, Rio de Janeiro, RJ, 22453900, BRAZIL
Petrobras, R&D Center, Rio de Janeiro, RJ, 21941915, BRAZIL
mijail febres @hotmail.com, nieckele @pucrio.br, robertofonseca @petrobras.com.br, lfaa@ pucrio.br
Keywords: Slug flow; horizontal pipeline; three dimensional, VOF; PIV
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 gasliquid 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 3D
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 Ks 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.
Introduction
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
twophase 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
submodel 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
flow.
Instantaneous full field measurements of the liquid
phase in slug flow were reported by CarpinteroRogero et
al. (2006). The authors employed three simultaneous optical
techniques to produce timeaveraged 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 30June 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 gasliquid
characteristics, structures evolution and interaction along the
pipeline.
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 onedimensional
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 nonconservative
formulation.
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 2D and 3D
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 twophase slug flow with the single unit slug
approach. The slug translational velocity and film velocity
profiles were compared with measured data, presenting good
agreement.
Nomenclature
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)
P
0
Subscripts
ef
g
m
s
t
specific mass (kg/m3)
surface tension (N/m)
effective
gas
liquid
mixture
superficial velocity
turbulent
Mathematical Model
To determine the twophase 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
8ag
a + iVag =0 (1)
a8t
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)
at
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 (og) 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 KS 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
equations:
(6)
t(o )+V(p A)= V[I PefV I] + G, pC
C(e)+V.(* iE)= V[ae Pef VI]+
2 (7)
+ C, G, C26P 
K K
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 30June 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
Yjunction positioned at the entrance of the Plexiglas pipe.
After passing through the test pipe, the twophase 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(177/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
(11)
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 lengthtodiameter 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
Fig.lb, 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,
120mJ, NdYAG 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
Uef
(10)
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 1030 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 50mm 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 sandblasted diffuser
glass plate was positioned in front of the LED panel to
improve the spatial uniformity of the background
illumination. The highpass 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 15pmdiameter 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 gasliquid
interface at the photogate position. A controlling circuit
triggered by the photogate signal outputs a trigger signal to
the velocity measuring system after a preset time delay. By
controlling the time delay, it was possible to always
measure the liquid velocity field at a predetermined
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 30June 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.
LZ
(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
Inlet
Pressure@
Water
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
velocity
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 pressurevelocity 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
distribution,
Co =1.05 >Frc <3.5
C= 1.20F 5 Frc
Co = 1.20 Frc 2 3.5
Wal
C1
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 30June 4, 2010
and Fig. 4b illustrates the bubble position after the
converged wall velocity (translation bubble velocity) was
obtained.
St= 0 20s
t 0.40s
St= 0.60s
1=0.80s
(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.790x105 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 reevaluated. 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
region.
1 s slug
elongated bubble l
0.1 anular
0.1
stratifed s
0.01
0.01 0.1 1 10
.rJm/ls
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 30June 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
I I
I Ii
I I I I
I I I I I I I I I I
I I I I I I I I I
I I I I I I I I I I
I I I I I I I I I I
I I II I I I I I I
I I I I I I I I I I I
I I I I I I I I I I I
Figure 9: Data Acquisition positions
W(m/s)
(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. 1012 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. 1012, 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
W(m/s)
(a)z = 0.0 D
00 03 06 09 12 15 18
W(m/s)
(c) z= 0.4 D
W(m/s)
(b) z= 0.2 D
00 03 06 09 12 15 18
W(m/s)
(d) z= 0.6 D
Figure 12: Axial velocity profile, upstream bubble nose.
Case 6
(d) z= 0.6 D
ii.
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 30June 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
1.00
S0.00
1.00
010
z(m)
(a) Case 4
0
0 5
10 C
00 03 06 09
W(m/s)
(a) Case 2: z= + 0.2 D
10
05
05
10
00 03 06 09
W(m/s)
(c) Case 3: z= + 0.2 D
10
05
00
05 
10
00 03 06 09 12
W(m/s)
(e) Case 6: z= + 0.2 D
0 08 000 0 08
0 08 000 0 08
00 03 06 09
W(m/s)
(b) Case 2: z= + 0.6 D
W(m/s)
(d) Case 4: z= + 0.6 D
z(m)
(b) Case 6
Figure 15: Vertical velocity profile, upstream bubble nose
020 000 020
15 00 03 06 09 12 1
W(m/s)
(f) Case 6: z= + 0.6 D
q
Figure 13: Axial velocity profile, downstream bubble nose.
101 I 1 .
05
y/R
1.1.
''I
05
+0.8D
.10
00 03 06
W(m/s)
05
04 06 08
W/Wmax
(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
100
000
1 00
0,OD +0,2D +0,4D
z(m)
(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
+0,6D
000 020
0120 0
020
L0
10 0 O
000 020
z(m)
(b) Case 6
Figure 16: Vertical velocity profile, downstream bubble nose
An interesting result to analyze is the velocity field at
the pipe crosssection, in the liquid and gas region. Figure
17 illustrates the crosssection 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
~t~s
~~"~
i~ir~T
u
....... ...;~i~
r%:,~..........~,,;s 
i ;~~~ i'
1__\~~~1 1111~_~~
_II I I I I I I_
\~~_~LI1111111~r~~
1 I ( I I I I
~ ( 1 I I I I ( ) )
1.)
1I 1 I
I I
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 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.
1012.
Conclusions
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.
Acknowledgements
The authors thank CNPq and Petrobras for supporting the
development of this work.
References
(c) z= 0.6 D (d) z= 0.8 D
Figure 17: Velocity field at the crosssection, 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 SteadyState Annular Flow. Chemical
Engineering Science, Vol. 44(2), 325332 (1989).
Bendiksen, K. H. An Experimental Investigation o f the
Motion of long bubbles in inclined pipes. International
Journal of Multiphase Flow 10(4), 467483 (1984).
Bendiksen, K.H.; Malnes, D.; Moe, R.; Nuland, S. The
dynamic twofluid model OLGA: theory and application.
SPE Prod. Eng., Vol. 6, 171180 (1991).
Bertola, V Slug velocity profiles in horizontal gasliquid
flow, Experiments in Fluids, Vol. 32, pp. 722727 (2002).
Black, P. S., Daniels, L. C., Hoyle, N. C., Jepson W. P
Studying Transient MultiPhase 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), 335354 (1992).
CarpinteroRogero, E., Kroess, B., Sattelmayer, T.
Simultaneous HSPIV and shadowgraph measurements of
gasliquid flows in a horizontal pipe. 13th Int. Symp. on
Applications of Laser Techniques to Fluid Mechanics
(2006).
Cook, M. and Behnia, M. Slug length prediction in near
horizontal gasliquid intermittent flow. Chemical
Engineering Science, Vol. 55, 20092018 (2000).
Cook, M. & Behnia, M. Bubble motion during inclined
intermittent flow. International Journal of Heat and Fluid
Flow, Vol. 22, 543551 (2001).
Fagundes Netto, J.R., Fabre, J., P6resson, L. Shape of long
bubbles in horizontal slug flow. International Journal of
Multiphase Flow, Vol. 25, 11291160 (1999).
Hutchinson, B. R. and RAITHBY, G.D. A Multigrid
Method Based on the Additive Correction Strategy.
Numerical Heat Transfer 9, 511537, (1986).
Issa, R.I. Solution of the implicitly discretized fluid flow
equations by operatorsplitting. J. Comput. Phys. 62, 4065,
(1986).
Ishii, M. ThermoFluid Dynamic Theory of TwoPhase
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 twofluid
model. International Journal of Multiphase Flow, Vol. 29,
6995 (2003).
Leonard, B. P. A stable and accurate convective modeling
procedure based on quadratic upstream interpolation.
Comput. Method Appl. Mech. Eng. 19, 5998 (1979).
Orell A. Experimental Validation of a simple model
gasliquid slug flow in horizontal pipes. Chemical
Engineering Science, Vol. 60 (2205) 13711381 (2"'111
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Ortega, A. J. & Nieckele, A. O.Simulation of horizontal
twophase slug flows using the twofluid model with a
conservative and nonconservative 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, 957975
(1999).
Shemer, L. Hydrodynamic and statistical parameters of slug
flow. International Journal of Heat and Fluid Flow 24,
334344 (2003).
Taitel, Y & Barnea, D. Twophase slug flow. Adv. Heat
Transfer, Vol. 20, 83132 (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, 833844
(2000)
Taha Taha, Z.F. Cui.CFD Modelling of slug flow in vertical
tubes. Chemical Engineering Science, Vol. 61, 676687
(2006).
Ujang, P.M., Pan, L., Manfield, P.D., Lawrence, C.J. and
Hewitt, G.F. Prediction of the translational velocity of
liquid slugs in gasliquid slug flow using computational
fluid dynamics. Multiphase Science and Technology, Vol.
20(1), 2579 (2008)
Vall6e, C., HOhne, T., Prasser, H., Siihnel, T. Experimental
investigation and CFD simulation of Horizontal stratified
twophase flow Phenomena. Nuclear Engineering and
Design, Vol. 238, 637646 (2008)
Wallis, G.B., 1969. Onedimensional Twophase Flow,
McGrawHill, Inc, New York.
Wang X., Guo L. and Zhang X., An experimental study of
the statistical parameters of gasliquid twophase slug flow
in horizontal pipeline. International Journal of Heat and
Mass Transfer, Vol. 50, 24392443 (2007).
