Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Terminal velocity of a Taylor drop in a vertical pipe
Kosuke Hayashi, Ryo Kurimoto and Akio Tomiyama
Graduate School of Engineering, Kobe University
11, Rokkodai, Nada, Kobe, Japan
hayashi@mech.kobeu.ac.jp, kurimoto @cfrg.scitec.kobeu.ac.jp, tomiyama@mech.kobeu.ac.jp
Keywords: Taylor drop, Taylor bubble, Slug flow, Terminal velocity, Numerical simulation
Abstract
A correlation of Froude number, Fr, for Taylor drops rising through stagnant liquids in vertical pipes is proposed. The
fundamental functional form of the correlation expressed in terms of ReD, EoD, p* and L* is deduced by making use of a
scaling analysis based on the field equations of two phases and the jump conditions at the interface, where ReD is the drop
Reynolds number, EoD the Eotvos number, p* the density ratio, and u.* the viscosity ratio. Coefficients in the correlation are
determined using experimental and numerical data. Taking the limits, p* 0 and u 0, reduces the correlation to that for
Taylor bubbles in liquids. The coefficients relating with the effects of ReD and EoD are, therefore, determined by making use of
the knowledge on Fr of Taylor bubbles. Then, the effects of p* and uL* are investigated by using an interface tracking method.
The simulation showed that p* has no influence on Fr whereas the increase in u.* decreases Fr. By introducing the effect of u.*
into the correlation for Taylor bubbles, an extended correlation of Fr expressed in terms of ReD, EoD and uL* is deduced.
Comparisons between the extended Fr correlation, the measurements and the interface tracking simulations demonstrate that
the correlation is applicable to Taylor drops for a wide range of the dimensionless groups, i.e., 4.8 < EoD < 23, 1 < N< 14700,
12 < log < 4, 0.012 < u* < 20, and 0.0021 < ReD < 4959, where N and M are the inverse viscosity number and the Morton
number, respectively.
Introduction
Large gas bubbles in gasliquid twophase slug flows in
vertical pipes take bullet shapes. They are called Taylor
bubbles (Davies & Taylor, 1950). Since large drops in
vertical pipes also take bullet shapes, they can be also
referred to as Taylor drops.
Many studies have been carried out on the motion of a
Taylor bubble, and various models for the terminal velocity,
VT, have been proposed. A review and comparisons between
available VT models for a Taylor bubble can be found in
literature (Viana et al., 2003). Since the viscosity of the gas
phase has no influence on VT, these VT models do not take
into account the bubble viscosity. On the other hand, the
drop viscosity affects VT of the Taylor drop.
In contrast to the Taylor bubble, there are few studies on
the Taylor drop. Goldsmith & Mason (1962) deduced a VT
model for Taylor drops at low drop Reynolds numbers, ReD,
in which the effects of the drop viscosity on VT were
accounted for. The original expression of Goldsmith's VT
model can be rearranged as follows:
V [/f 1+( fi f2 f3 )(GD IC C) ApgD2
L 1+(f2+ f4) /wD ) I pC
(R h)2 R2 (Rh)2 h
8R 2 2(Rh)2 Rh
2 Rh 3 Rh
4
f4 =R h)
.2 R +
R h)
where Ap (=pc oD) is the density difference between the
two phases, p the density, g the magnitude of the
acceleration of gravity, D the pipe diameter, u. the viscosity,
andf, (i = 1 to 4) are functions of the pipe radius, R, and the
thickness, h, of the liquid film between the drop and the
pipe wall. The subscripts, C and D, denote the continuous
and dispersed phases, respectively. The drop Reynolds
number is defined by
pcVrD
ReD 
Equation (1) can be further rewritten in the following
nondimensional form:
Paper No
Fr _f 1+(flf2 + f3)*
+I +(f2 + f4)*
where Fr is the Froude number, N the inverse viscosity
number, and iL* the viscosity ratio:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Scaling analysis
Derivation of functional form of Fr. The local
instantaneous field equations and the jump conditions for
incompressible viscous fluids at a constant temperature are
given by
(4) V Vk = 0
Fr VT
VApgD/pc
N _pcApgD3
tLc
+vk +VkVVk
8t
(VD Vc) n =0
1 1
VPk +g+lv.t
Pk Pk
SLD
Lc
(PD PC KG)n + (T Tc) n = 0
They also carried out experiments of Taylor drops and
confirmed the validity of their model.
Mandal et al. (2007) measured terminal velocities of
Taylor drops at high ReD for various fluids. They confirmed
that the terminal velocity can be evaluated by using the
following VT model for sphericalcap bubbles in infinite
stagnant liquids (Joseph, 2003):
8 (1+8s) f 16s 32
Fr 8 (+8s) + 2 1s + 16 32 (1+8s)2 (7)
3 N 3 EOD N2
where s is the deviation of the interface from perfect
sphericity near the nose of a bubble, and EoD is the Eotvos
number defined by
SApgD2
EoD 
where the subscript k denotes the phase (k = D and C for the
dispersed and continuous phases, respectively), V is the
velocity, t the time, P the pressure, g the acceleration of
gravity, the viscous stress tensor, n the unit normal to the
interface, and K the curvature. Characteristic velocity, length
and time scales for the problem of concern are VT, D and
D/VT, respectively. As shown in Fig. 1, the orders of
magnitude of the viscous stress tensors can be estimated as
TD DV T /(R h) (13)
rc VcVT /h
Substituting Eqs. (13), (14) and the characteristic scales into
Eq. (10) yields
Pc CC1PCT2 +CC2PcgD+Cc3CLVTD/h2
where C is the surface tension. Their work implies that VT of
a Taylor drop at high ReD is independent of u.*. Mandal et al.
(2009) recently proposed another correlation of VT, which is
a combination of Brown's Fr model (Brown, 1965) and a
function of EoD used in Viana's correlation for a Taylor
bubble (Viana et al., 2003). Since the effect of u.* is not
taken into account, this model is also applicable only to
drops at high ReD.
These models, Eqs. (3) and (7), are applicable only to the
viscous force dominant (low ReD) regime and to the inertial
force dominant (high ReD) regime, respectively. In addition,
without the knowledge on the drop shape parameters, h and
s, they are of no use.
The purpose of this study is to develop a Fr correlation
for Taylor drops. A Fr correlation for Taylor bubbles has
been already developed by the authors (Hayashi et al., 2010)
by making use of a scaling analysis based on the local
instantaneous field equations and the jump conditions at the
interface (Tomiyama, 2004). By introducing the effect of u.*
into the correlation, the correlation is extended so as to be
applicable to Taylor drops. Coefficients in the extended
correlation are determined using a database of Froude
numbers obtained by using an interface tracking method
(Hayashi et al., 2006; Hayashi & Tomiyama, 2007, 2009).
Experiments of Taylor drops are also carried out to examine
the validity of the interface tracking simulation.
PD CCD1DT2 +CD2pDgD + CD3IDV D/(Rh)2 (16)
Pipe axis
Figure 1: Large drop in a pipe
Paper No
where Ck, (k = D or C and i = 1, 2, 3) are constants. The
curvature, K, in Eq. (12) is estimated as
K I1/(Rh) (17)
Substituting Eqs. (15)(17) into Eq. (12) gives the following
force balance:
D3 D3
CIpcV/2D2 +C2pDV2D2 +C VT 2+C DVT 2 h2
h (R h)
D2
+Cs +C6ApgD3 = 0
Rh
(18)
where C, (i = 1 to 6) are constants. This equation indicates
that the six different forces, i.e., inertial (F,c and FD),
viscous (Fac and FD), surface tension (F,) and buoyancy
(Fb) forces, govern the dynamics of a Taylor drop. The
relevant dimensionless groups, Fr, N, EoD, M, ReD, p* and
IL*, are represented by using these forces:
Fr= . = Fc (Froude number) (19)
Fb cApgD/pc
N = Fb cApgD3 (inverse viscosity number)(20)
V FP LC
EoD F ApgD (E6tv6s number) (21)
F, a
FF
F,3 2
4
gRC Ap
S(Morton number)
2 3
Pc0
1/) \
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
where c, (i = 1 to 5) are undetermined constants given by cl
= C6/C3, c2 = C1/C3, c3 = 2/C3, c4 = C4C3, and c5 = C5/C6.
The effects of It* vanish when the ReD is much larger than
the factor including It*. As described above, the viscosity
ratio does not affect Fr at high ReD (Mandal et al., 2007).
Hence the functional form of Eq. (26) is reasonable with
respect to the It* dependency.
The scaling analysis does not give the values of the
coefficients in Eq. (26). Hence the coefficients must be
determined using the knowledge on Fr.
Fr correlation for Taylor bubbles. It can be postulated
for a bubble in a liquid that p* z 0 and It* = 0. Then Eq.
(26) simplifies to
h )2
C D 1
Fr = C 2(h2 + 1+C5s
S(h2 R1 I R h) EO
\D) ReD
Let us consider two limiting cases, i.e., the inertial force
dominant (ReD co and EoD co) and the viscous force
dominant (ReD 0 and EoD co) cases to determine the
coefficients, cl and c2, in Eq. (27). In the former case, Eq.
(27) reduces to
Fr =
C2C
Since the Froude number in the inertial force dominant
regime is known to be 0.35 (Brown, 1965; Collins et al.,
1978; Dumitrescu, 1943; Griffith & Wallis, 1961; Harmathy,
1960; Wallis, 1969; White & Beardmore, 1962), we obtain
ReD = pcVD (drop Reynolds number) (23)
Fpc Itc
p*=D PD (density ratio) (24)
F, Pc
i* = =  (viscosity ratio) (25)
Fpc ftc
Only five independent dimensionless groups can be deduced
from the six forces. Furthermore the constraint given by Eq.
(18) reduces the number of independent dimensionless
groups from five to four. By selecting ReD, EoD, p* and It*
as the four independent variables, we can obtain the
following fundamental functional form of Fr from Eq. (18):
cC= 0.35
C2
Substituting a relation, ReD = FrN, into Eq. (27) yields
c2 Fr2 + =C l+c (30)
D) N D) 5 (R h) EOD (30)
According to Wallis (1969), Fr of a Taylor bubble in the
viscous force dominant regime is given by
Fr = 0.01N
Substituting the above equation into Eq. (30) and taking the
limit EoD oo yields
h 2 h 2
C2 ( Fr2 +0.01= cl
D D
Since Fr vanishes at ReD 0, Eq. (32) reduces to
hY2 =0.01
D)
h (c2 +(C3P*)+ 1+C4 h 2 2 *
D ReD (Rh)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
10 
S102
103 0 White & Beardmore (1962)
X Wallis(1969)
(0.01/(0.0816+ReD 1))05
1 0 4 I I I I I I I I
105104 103102 1010 101 102 103 104 105
ReD
Figure 2: Relation between ReD and Fr for EoD > 100
From Eqs. (29) and (33)
c2 j = 0.0816 (3
Substituting Eqs. (33) and (34) into Eq. (27) yields
Fr 0.01
Fr = G
S0.0816+
SReD
where
Figure 3: Function G
G =f 0.0816+
0.1 t ReD)
Figure 3 shows the values of G obtained by substituting
White & Beardmore's data (1962) into the above equation.
4) Though there is a small scatter in the data, G is well
correlated with the following function:
4 \ 463
+ 4 41
K EoD
By substituting Eq. (39) into Eq. (35), we obtain the
following correlation of Fr for Taylor bubbles:
D 1
G = 1+c5 E
Rh EoD
) Fr 463
0.01 41
(36) 0.0816+ReD EOD6
Then, taking the limit, EoD * o, in Eq. (35) gives
Fr= 0.01 (37)
0.0816+
ReD
Since the above equation is derived by making use of the
two limiting cases, ReD co and ReD 0, its applicability
to intermediate Reynolds numbers should be tested against
experimental data of high Eotvos numbers, EoD > 100. As
shown in Fig. 2, Eq. (37) agrees well with the experimental
data (White & Beardmore, 1962; Wallis, 1969) not only at
high and low Reynolds numbers, ReD > 103 and ReD < 1, but
also at intermediate Reynolds numbers, 1
Then the remaining term to be determined in Eq. (35) is
G. This term might be a function of two dimensionless
groups, e.g., ReD and EoD. As can be understood from Eq.
(36), the most relevant dimensionless number is, of course,
the E6tvos number. Hence, for simplicity, let us deduce an
empirical correlation of G in terms of EoD. Solving Eq. (35)
for G gives
Substituting the relations, ReD = FrN and N = (EoD3/l 14,
into Eq. (40) yields
] 1/2 1 463
0.0816Fr2 + 2Fr =0.01 1+ 41
Eo3/4 E96
EoD Eo ,
By taking the limit M 0, the above equation reduces to
( 232
Fr = 0.351+ 4196
y FoD J
This simplified equation is applicable to low viscosity
systems such as vapor bubbles in water. It should be noted
that the three limits, ReD co, N * o and M * 0, lead to
the same limiting case.
Equation (40) agrees with most of Fr data obtained by
White & Beardmore (1962) to within +10 % errors (Hayashi
et al., 2010). The accuracy is, however, not so good at
intermediate Morton numbers. This shortcoming can be
overcome by introducing a correction into Eq. (40) while
taking into account the nonlinear effects, e.g.,
Paper No
Paper No
0 .4 1 I I I 1
logM<7.7 4.9 4.3 2
o <7.
A4.9
0.3C 4.3
v 2.7
0 1.8
0.48
0.45
A 0.90
L0.20 2.0
v 2.9
S6.5
0.1
0
4 10 102
EOD
Figure 4: Comparison between Eq.
Beardmore's data (1962)
Fr \ 0.0089 K
SFr=31
0.0725 +1(10.11Re33)
ReD D
logM
0.45
0.90
2.0
2.9
(43) and White &
S463
41
1 96
FD
where the values of the coefficients, cl and c2, are also
slightly changed to improve the accuracy. This equation is
compared with White & Beardmore's experimental data
(1962) in Fig. 4. Good agreements between the data and Eq.
(43) are obtained for a wide range of the dimensionless
groups, EOD and M.
The coefficients, c3 and c4, in Eq. (26) must be
determined based on the knowledge on the effects of p* and
.L* on Fr of Taylor drops. There are, however, no databases
of Fr for various values of p* and ut*. The effects of p* and
uL* on Fr were, therefore, investigated by making use of an
interface tracking method.
Numerical method
The following mass and momentum equations based on
the onefield formulation for an incompressible twophase
flow are adopted (Hirt & Nicholds, 1981; Scardovelli &
Zaleski, 2000):
V = 0 (44)
,PM CV +V V
pM t (45)
= VP + V pM [VV + (VV)T ] + pg + CoKn6
where 6 is the delta function and the superscript T denotes
the transpose. The mixture density, pu, is given by
PM =(1 a)pD +apc (46)
where a is the cellaveraged volume fraction of the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
continuous phase. The harmonic mean of the viscosity
should be used rather than the arithmetic mean to correctly
deal with the continuity of the viscous stresses. Hence the
mixture viscosity, UtM, is given by
1 1a a
= +
RM RD PC
The density and viscosity of each phase are assumed to be
constant. The computational cell is filled with the
continuous phase when a = 1, and with the dispersed phase
when a = 0. A cell with 0 < a < 1 contains an interface. The
interface motion is calculated by solving the following
advection equation of a:
+V*.V = aV.V
Ot
Equations (44) and (45) are solved by using the modified
SOLA (solution algorithm) (Tomiyama & Hirano, 1994). An
interface tracking method based on a volumeoffluid
method, which we call the NSS (nonuniform subcell
scheme) (Hayashi et al., 2006), is used to solve Eq. (48).
The amount of fluid volumes transferred from an interface
cell to its neighbor cells is accurately calculated by making
use of nonuniform subcells which are temporarily
introduced only in the interface cells. A local level set
function is also computed in the interface cells, by which
the surface tension force in Eq. (45) is accurately evaluated.
The NSS can be implemented not only in Cartesian
coordinates but also in cylindrical and general curvilinear
coordinates. Since all the Taylor drops dealt with in this
study are axisymmetric, the twodimensional (r, z)
cylindrical coordinate system is used. The details of the
NSS can be found in Hayashi et al. (2006) and Hayashi &
Tomiyama (2007, 2009).
Validation
Experimental setup. Terminal velocities and shapes of
single Taylor drops were measured to obtain experimental
data for comparisons with simulations. The experiments
were carried out at atmospheric pressure and room
temperature (the temperature of the liquid and drops were
kept at 298 + 0.5 K). The experimental apparatus is shown
in Fig. 5. It consists of the vertical pipe and the square tanks.
Four pipes of D = 11.0, 20.1, 26.1 and 30.8 mm were used
to alter the values of EoD. The pipe length was 530 mm. The
pipes and the tanks were made of transparent acrylic resin.
Silicon oils of various viscosities (Sinetsu silicon,
KF9610, 30, 100, 300, 500, 1000 and KF965000) and
glycerolwater solutions were used for the dispersed and
continuous phases, respectively. The density, viscosity,
surface tension and temperature were measured using a
densimeter (Ando Keiki Co., Ltd., JIS B7525), a viscometer
(Rion Co., Ltd., Viscotester VT03E), capillary tubes (glass
tube, 1.02 mm i.d.) and a digital thermometer (Sato
Keiryoki MFG Co., Ltd., SK1250MC), respectively. Fluid
properties were measured at least five times. Uncertainties
estimated at 95 % confidence in measured p, ut and C were
0.021, 0.53 and 4.0 %, respectively. The measured values
of these fluid properties agreed well with the data given in
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
solution II
Silicon
drop,
Highspeed
video camera
Hemispherical
cup
Side view
Side view
 ipe
I Duct
E I
E LED
o
CO)
SPipe
SDuct
Highspeed
video camera
Crosssectional view
AA'
Figure 5: Experimental apparatus
nose
(a) Raw image (b) Binary image
Figure 6: Raw and binary images of a Taylor drop
literature (Ishikawa, 1968). The values in literature were,
therefore, used in calculations of the dimensionless numbers
and in numerical simulations. Various fluid systems shown
in Appendix were used.
A small amount of silicon oil was stored in the
hemispherical glass cup. A single drop was released by
rotating the cup. The visualization of the drop motion and
shape using a highspeed video camera (Integrated Design
Tools Inc., Motionscope M3, frame rate = 20 400 frame/s,
spatial resolution z 0.095 mm/pixel) confirmed that drops
were axisymmetric and rose rectilinearly along the pipe
axis.
Successive images of drops in the pipe were taken by
using the highspeed video camera at 350 mm above the
bottom of the pipe to obtain experimental data only at
terminal conditions. The pipe was enclosed with the acrylic
duct. The gap between the pipe and the duct was filled with
the same glycerolwater solution as that in the pipe to
reduce optical distortion of drop images. A LED light source
(Hayashi WatchWorks, LP2820) was used for backward
EOD
Figure 7:Froude numbers measured in the present study
and Fr data collected from literature (Goldsmith
& Mason, 1962; Mandal et al., 2007)
illumination.
The spherevolume equivalent diameters, d, of drops
were calculated from the drop images using an image
processing method (Tomiyama et al., 2002). An example of
drop images is shown in Fig. 6. The velocity of a drop was
computed from two binary images by measuring the
difference in elevation of the noses of the two drop images.
The measured VT showed that the drops reached their
terminal conditions at the measuring section.
Experimental results. All the Fr data obtained by 34
runs are shown in Fig. 7. The values of Fr and drop shapes
are given in Appendix. The experimental data obtained by
Goldsmith & Mason (1962) and Mandal et al. (2007) are
also plotted in the figure. Goldsmith & Mason dealt with
only low Reynolds number drops, in other words, low Fr
drops. The ranges of dimensionless groups examined by
them are 14 < EoD < 50, 1 < N < 12, logM = 0.76 and 3.4,
u* = 0.0012 and 0.075, 1.2 < p* < 1.5, 0.004 < ReD < 1.1,
and 0.004 < Fr < 0.091. On the other hand, Mandal's
experimental data were only for drops in low viscosity
systems. The ranges of dimensionless groups are 4.8 < EoD
< 228, 1430 < N< 14300, 11 < logM < 9.3, 0.69 < u <
1.2, 0.66 < p* < 0.88, 78 < ReD < 4960, and 0.05 < Fr <
0.349. Since the drop viscosity dissipates the kinetic energy,
Fr may decrease as uL* increases. However, their data
largely scatter and several data exceed the values of Fr
evaluated with Eq. (43) at logM= 11. Taylor drops in low
viscosity systems and those in large pipes tend to oscillate
due to the RayleighTaylor instability. This may be a reason
why Mandal's data show a large scatter. Figure 7 also shows
that the present experiments cover a wide range of
dimensionless groups, i.e., 7.2 < EoD < 92, 7.9 < N < 270,
4.7 < logAM < 0.4, 0.017 < u.* < 19, 0.74 < p* < 0.81,
0.029 < ReD < 74, and 0.0037 < Fr < 0.284, and therefore,
they fill the gap between Goldsmith's and Mandal's data.
Figures 810 are measured Fr at log = 0.45, 2.5 and
4.7, respectively. Equation (43) is also drawn in the figures
as the reference case of u.* = 0. The effects of u on Fr can
Paper No
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
EOD
Figure 8: Froude numbers at logM:
0.45
UA4
0 0.34
A 1.1
0.3 O 5.7
/ A Eq. (43)
0.2 /
0.1
0
4 10 102 1
EoD
Figure 9:Froude numbers at logM= 2.5
EOD
Figure 10: Froude numbers at logMA
0.2
logM=0.45
[]
[]
EOD
0 12
A 39
O 91
0.1 A uO
o o Q....O. o .O....C
102 101 1 10 102
1*
Figure 11: Froude numbers at various g*
be understood from these figures: (I) Fr monotonously
decreases with increasing *, (II) the influence of g*
vanishes when * is very small, and (III) * does not have
much influence on Fr at low M. The tendency (II) supports
the model of Goldsmith & Mason, and the term in the
square brackets in Eq. (3) must decrease with increasing g*
because of the tendency (I). Then the tendency (III)
qualitatively agrees with the observation of Mandal et al.
that Fr of a drop at low M is independent of g*.
The Froude numbers at EoD 12, 39 and 91 are plotted
against L* in Fig. 11 to make clear the effects of g*. Fr
decreases with increasing g*, and the magnitude of the
gradient, dFr/dg*, decreases as g* increases.
Since the density ratio is about 0.8 for all the experiments,
the effects of p* cannot be discussed. These are, therefore,
investigated by making use of the numerical simulation.
Simulation. Figure 12 shows a computational domain and
the initial condition. The drop was initially set at the center
of the computational domain. The initial drop shape was a
cylinder with two hemispheres attached to the top and
bottom of the cylinder. The drop length was L and the
thickness, h, of the liquid film was 0.2R. The ratio of the
spherevolume equivalent diameter, d, to the pipe diameter,
D, was 1.25. The domain size was (R, 4L). The top and
bottom boundaries were uniform inflow and continuous
outflow, respectively. The right boundary was a wall
moving downward at the drop terminal velocity. The left
was an axisymmetric boundary. Uniform computational
cells of Ar = Az = R / 32 were used. The total number of
cells was 18,816 (=32 x 588). The time step, At, ranged
from 1 x 103 to 1 x 102 ms which satisfied the condition of
numerical stability.
Numerical simulations with a higher spatial resolution, Ar
= Az = R / 64, were carried out at several conditions to
investigate the grid dependence. Terminal velocities
predicted with the higher resolution agreed well with those
predicted with the lower resolution. The relative error was
about 1 %. Hence the resolution of Ar = R / 32 was
sufficient for obtaining good predictions of VT.
Another domain size, (R, 7L), was also tested. The result
confirmed that the influence of the top and bottom
boundaries on predictions of VT was negligibly small even
with the domain size of (R, 4L).
Paper No
mn
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
V
SLiquid
i Pc Cc
1.5L
Drop
S PD RD
L
h=0.2R
Sh=0.________
S R
9 1.5L
Pipe axisi g
(Slip wall J
Moving wall
Outflow
Initial drop:
Length L
h=0.2R
d/D=1.25
Domain size:
Rx 4L
(32 x 588 cells)
Ar = Az
Figure 12: Computational domain
Measured Fr
Figure 14: Comparison between measured and predicted Fr
4 10 102
Figure 13: Comparison of Fr in the EoDFr diagram
Comparisons. Comparisons between measured and
predicted Fr are shown in Fig. 13. The simulations well
reproduce the measured Fr for the wide range of
dimensionless groups. The predicted Froude numbers are
directly compared with the measured Fr in Fig. 14, which
shows a very good agreement between them. Measured and
predicted drop shapes are compared in Fig. 15. The
predicted drop shapes also agree very well with the
experiments.
The good agreement assures that the NSS can be used to
investigate the effects of dimensionless groups on Fr.
Effects of p* and L*
Simulations were carried out at several values of p* to
investigate the effects of p* on Fr. The other dimensionless
groups were set at (logA, EoD, u*)=(1.0, 50, 1.0) and
(logM, EoD, u*)=(3.0, 20, 1.0). Predicted Froude numbers
are plotted against p* in Fig. 16, which clearly shows that
p* has no influence on Fr. Hence the coefficient c3 in Eq.
(26) must be zero, and the equation becomes a function of
ReD, EoD and Wt*:
Figure 15: Comparisons between measured and predicted
shapes (left: measured, right: predicted). (a)
(logM, EoD, *, d/D) = (0.45, 12, 0.11, 1.24);
(b) (0.45, 12, 1.1, 1.27); (c) (0.45, 39, 1.1,
1.07); (d) (2.5, 10, 1.1, 1.19); (e) (4.7, 26, 1.1,
1.01); (f) (4.7, 8, 3.8, 1.34); (g) (4.7, 8, 1.1,
1.24)
Paper No
(c)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.14
0.12
0.1 
n n00 8 Inn lo=_i n F=n = ,, *=1 n
S yv .I D u U
E logM=3.0, EoD=20, Ix*=1.0
0.06
0 0.2 0.4 0.6 0.N
P* (=PD/Pc)
Figure 16: Effect of p* on Fr
10
ig (=1:D/ctiH)
Figure 17: Function H
D)ReD
( h H ,
where
requirements to reproduce the experimental fact: (1)
H(iX*=0) = 1, (2) H increases with IL*, and (3) H becomes
(49) constant as * ( o. Referring to the numerical results and
to the functional form of Eq. (3), we adopted the following
functional form satisfying the above requirements:
D 1
1+ c5 h
5 (R h) EOD
H +c6C*
H +=7
1+c ci*
H = 1+c4 Rh *
(50) where c6 and c7 are constants and c6 > c7. These constants
were determined based on the data shown in Fig. 17, i.e. c6
= 1.75 and c6 = 0.27.
The Froude number defined by Eq. (19) can be regarded
as a combination of two dimensionless groups, VT/(gD)1 2
and p*, i.e. Fr = Vl/[(1p*)gD]1/2. Hence the dimensionless
terminal velocity, Vr/(gD)1/2 which is frequently used as the
Froude number for a Taylor bubble, is proportional to
(1p*)12, i.e. VT/(gD)12 = (1p*)12f(ReD, EOD, u'*).
Comparing Eq. (49) with Eq. (27) gives the following
expression
463
=Fr 0.0089 41 41
Fr = H+ EO196 (51)
0.0725 +R (10.11Re 33)
ReD
in which the effect of the drop viscosity is accounted for. As
can be understood from Eq. (50), H should be a function of
the viscosity ratio. Let us, therefore, deduce an empirical
correlation of H in terms of ut*.
Interface tracking simulations were carried out at various
values of M, EoD and ut*, i.e. 8 < logM< 4, 9.8 < EoD < 30
and 0 < ut* < 20, to obtain the values of H. Solving Eq. (50)
for H gives
/ 463
S Re 0.0089+ 41 
1 _O. 11Re33 Fr2 EO196)
101 eD D
0.0725 (52)
Substituting the numerical results into Eq. (52), we obtained
the values of H as shown in Fig. 17. Since computed H
scattered a lot due to the estimation error of Eq. (43), the
values of H obtained at the same value of u.* were averaged.
The functional form of H should satisfy the following
Froude number correlation for Taylor drops
Substituting Eq. (53) and the constants, c6 and c7, into Eq.
(51) gives the following correlation for Taylor drops:
463
9 41
0.008911+ 496 /
Fr = l EoD) (54)
0.0725+ 1+ (10.11Re )
S ReD 1+0.27It*J
This correlation is compared with the numerical predictions
in Fig. 18. The correlation well reproduces the tendencies of
Fr, i.e. Fr increases with N and EoD, while it decreases as
uL* increases. The agreement between Eq. (54) and the
predictions is good. Hence the function of H, Eq. (53), well
expresses the effects of u.* on Fr. Figure 19 shows Fr in
terms of EoD, M and u *. Since ReD is much larger than H
for EoD > 30 and logM = 8, the viscosity ratio does not
have much influence on Fr. In this case, Eq. (42) can be
used even if u.* / 0. Fr decreases with increasing u.* for
logM > 8 even at high Froude numbers. The effects of u.*
on Fr become significant as M increases because of the
decrease in ReD.
The correlation is compared with the simulations and the
experiments in Fig. 20, which includes 288, 34, 25 and 5
data of the simulations, the present experiments, Mandal's
data and Goldsmith's data, respectively. The ranges of
dimensionless groups are 4.8 < EoD < 23, 1 < N < 14700,
12 < logAM< 4, 0.012 < uL* < 20, and 0.0021 < ReD < 4959.
Most of the data obtained in this study lie within +10 %
Paper No
v.vW
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
1 10 102 103 104
(a) EoD = 10
(b) EoD = 20
N
(c) EoD = 30
Figure 18: Comparisons of Fr
(symbols) and Eq. (54)
between predictions
errors. The agreement between the correlation and
Goldsmith's data is good. Though several Mandal's data do
not agree with Eq. (54) due to the scatter of their data shown
in Fig. 7, Eq. (54) well reproduces most of their data.
Equation (54) reduces to Eq. (43) when ut* = 0. Hence Eq.
(54) is applicable not only to Taylor drops but also to Taylor
bubbles.
Figure 19: Equation (54) on the FrEoD plane
O 0 Predicted
+ Measured
Mandal et al. (2007)
Goldsmith & Mason (1962
0.1 0.2 0.3 0.4
Measured and predicted Fr
Figure 20: Comparisons of Fr between measurements,
predictions and Eq. (54)
Conclusions
A correlation of Froude number, Fr, for Taylor drops
rising through stagnant liquids in vertical pipes was
proposed. The fundamental functional form of the
correlation, Fr = f(ReD, EoD, p*, L*), was deduced by
making use of a scaling analysis based on the field
equations of two phases and the jump conditions at the
interface, where ReD is the drop Reynolds number, EoD the
Eotvos number, p* the density ratio and u.* the viscosity
ratio. Coefficients in the correlation were determined using
experimental and numerical data. Taking the limits, p* 0
and u* 0, reduces the correlation to that for Taylor
bubbles in liquids. The coefficients relating with the effects
of ReD and EoD on Fr were, therefore, determined by
making use of the knowledge on Fr of Taylor bubbles. Then,
the effects of p* and u.* were investigated by using an
interface tracking method. The simulation showed that p*
has no influence on Fr, whereas Fr monotonously decreases
with increasing ut*. The remaining terms in the correlation
Paper No
Paper No
were determined based on the numerical predictions. As a
result, the following correlation was obtained:
9 463
0.008 1+ 1 96
0.0725+1 (+1.75* (1 0.11ReD)
ReD 1+0.27t*) D
Comparisons between the correlation, the measurements and
the interface tracking simulations demonstrated that the
proposed correlation is applicable to Taylor drops for a wide
range of the dimensionless groups, i.e., 4.8 < EoD < 23, 1 <
N< 14700, 12 < logM< 4, 0.012 < L* <20, and 0.0021 <
ReD < 4959, where N is the inverse viscosity number and M
the Morton number.
Acknowledgements
This work has been supported by the Japan Society for
the Promotion of Science (grantinaid for scientific
research (B), No. 21360084). The authors would like to
express our thanks to Mr. Yasuki Yoshikawa for his
assistance in the experiments and simulations.
References
Brown, R. A. S. The Mechanics of Large Gas Bubbles in
Tubes. Canadian Journal of Chemical Engineering, 217223
(1965)
Collins, R., DeMoraes, F. F., Davidson, J. F. and Harrison, D.
The Motion of a Large Gas Bubble Rising through Liquid
Flowing in a Tube, Journal of Fluid Mechanics, 89, 497514
(1978)
Davies, R. M. & Taylor, G I. The Mechanics of Large
Bubbles Rising through Liquids in Tubes, Proc. R. Soc.
Lond. A, Vol. 200, 375390 (1950)
Dumitrescu, D. T. Stromung and Einer Luftbluse in
Senkrechten rohr. Z. Angew. Math. Mech., 23(3), 139149
(1943)
Griffith, P. & Wallis, G B. TwoPhase Slug Flow, Journal of
Heat Transfer, 83, 307320 (1961)
Harmathy, T. Z. Velocity of Large Drops and Bubbles in
Media of Infinite or Restricted Extent. A.I.Ch.E. Journal, 6,
281288 (1960)
Goldsmith, H. L. & Mason, S. G The Movement of Single
Large Bubbles in Closed Vertical Tubes, Journal Fluid
Mechanics, 10, 4258 (1962)
Hayashi, K., Sou, A. and Tomiyama, A. A Volume Tracking
Method based on NonUniform Subcells and Continuum
Surface Force Model using a Local Level Set Function.
Computational Fluid Dynamics Journal, 15, 225232 (2006)
Hayashi, K. & Tomiyama, A. Interface Tracking Simulation
of Bubbles and Drops in Complex Geometries. Multiphase
Science and Technology, 19, 121140 (2007)
Hayashi, K. & Tomiyama, A. Drag Correlation of Fluid
Particles Rising through Stagnant Liquids in Vertical Pipes
at Intermediate Reynolds Numbers, Chemical Engineering
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Science, 64, 30193028 (2009)
Hayashi, K., Kurimoto, R. and Tomiyama, A. Dimensional
Analysis of Terminal Velocity of Taylor Bubble in a Vertical
Pipe, Multiphase Science and Technology, (2010) in
I''"""
Hirt, C. W. & Nicholds, B. D. Volume of Fluid (VOF)
Method for the Dynamics of Free Boundaries. Journal of
Computational Physics, 39, 201225 (1981)
Ishikawa, T. Theory of viscosity of liquid mixture, Maruzen
(1968) in Japanese
Joseph, D. Rise Velocity of a Spherical Cap Bubble, Journal
of Fluid Mechanics, 488, 213223 (2003)
Mandal, T. K., Das, G. and Das, P K. Prediction of Rise
Velocity of a Liquid Taylor Bubbles in a Vertical Pipe,
Physics of Fluids, 19, 128109, 4 pages (2007)
Mandal, T K., Das, G. and Das, P. K. Liquid Taylor Bubbles
Rising in a Vertical Column of a Heavier Liquid: An
Approximate Analysis, Journal of Fluid Engineering, 131(1),
011303, 7 pages (2009)
Scardovelli, R. & Zaleski, S. Direct Numerical Simulation
of FreeSurface and Interfacial Flow. Annual Review of
Fluid Mechanics, 31, 567603 (2000)
Tomiyama, A. & Hirano, M. An Improvement of the
Computational Efficiency of the SOLA Method. JSME
International Journal, Series B, 37, 821826 (1994)
Tomiyama, A., Celata, G. P, Hosokawa, S. and Yoshida, S.
Terminal Velocity of Single Bubbles in Surface Tension
Force Dominant Regime. International Journal of
Multiphase Flow, 289, pp. 14971519 (2002)
Tomiyama, A. Drag, Lift and Virtual Mass Forces Acting on
a Single Bubble. on CDROM of 3rd International
Symposium on TwoPhase Flow Modelling and
Experimentation, Pisa (211 14)
Viana, F., Pardo, R., Yanez, R., Trallero, J. L. and Joseph D.
Universal Correlation for the Rise Velocity of Long Gas
Bubbles in Round Pipes, Journal of Fluid Mechanics, 494,
379398 (2003)
Wallis, G B. OneDimensional TwoPhase Flow.
McGrawHill (1969)
White, E. T & Beardmore, R. H. The Velocity of Rise of
Single Cylindrical Air Bubbles through Liquids Contained
in Vertical Tubes. Chemical Engineering Science, 17,
351361 (1962)
Appendix
Fluid systems used in the present experiments are
summarized in Table Al. Silicon oils of various viscosity
and pipes of D = 11, 20, 26 and 31 mm were used to deal
with various viscosity ratios and Eotvos numbers. The
ranges of dimensionless groups tested were 7.2 < EoD < 92,
7.9 < N < 270, 4.7 < logM < 0.4, 0.017 < t* < 19, and
0.74 < p < 0.81.
Table A2 and Figure Al show the measured Fr and drop
shapes obtained in the present experiments. These
experimental data might be useful for validation of interface
tracking methods.
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Table Al: Experimental conditions
Case 1 2 3 4 5 6 7 8 9 10 11 12 13 14
Pc [kg/m3] 1220 1190 1240
p, [kg/m3] 955 965 970 935 955 965 915 935 955 965 970 970 970 975
Lc [mPas] 85 25 262
toD [mPas] 29 97 485 9.4 29 97 4.6 9.4 29 97 291 485 970 4875
c [mN/m] 31 31 31 38 32 37 29 29 29 29 28 28 29 28
D [m] 11, 20,26 ,31 11, 20 ,26 11, 20 ,31
logM 2.5 4.7 0.45
p* 0.78 0.79 0.80 0.79 0.80 0.81 0.74 0.75 0.77 0.78 0.78 0.78 0.78 0.79
u* 0.34 1.1 5.7 0.37 1.1 3.8 0.02 0.04 0.11 0.37 1.1 1.9 3.7 19
EoD z10, 33,55 ,77 8, 26 ,44 z 12, 39 ,91
Table A2: Measured Froude numbers
Run EoD N logM u* ReD Fr
1 10 24 2.5 0.34 0.70 0.029
2 34 60 2.5 0.34 13.7 0.23
3 9.7 24 2.5 1.1 0.41 0.017
4 32 58 2.5 1.1 8.7 0.15
5 55 87 2.5 1.1 20.2 0.23
6 76 111 2.5 1.1 31.0 0.28
7 9.6 23 2.5 5.7 0.21 0.0088
8 32 58 2.5 5.7 5.3 0.092
9 54 86 2.5 5.7 13.8 0.16 Run 1 2 3 4
10 75 110 2.5 5.7 22.0 0.20
11 9.4 79 4.9 0.37 6.0 0.075
12 8.6 76 4.7 1.1 3.1 0.040
13 29 188 4.7 1.1 55.3 0.29
14 7.2 75 4.9 3.9 1.2 0.016
15 24 184 4.9 3.9 36.2 0.20
16 41 273 4.9 3.9 73.6 0.27
17 12 8.2 0.45 0.11 0.14 0.017
18 39 20 0.45 0.11 2.4 0.12
19 91 38 0.45 0.11 9.04 0.24 Run 7 8 9
20 13 8.8 0.40 0.017 0.23 0.027
21 12 8.5 0.44 0.036 0.19 0.022
22 11 8.1 0.45 0.37 0.096 0.012
23 38 20 0.45 0.37 1.7 0.086
24 11 8.0 0.42 1.1 0.064 0.0080
25 38 20 0.42 1.1 1.2 0.059
26 90 37 0.42 1.1 5.1 0.14
27 12 20 0.41 1.9 0.053 0.0067
28 39 37 0.41 1.9 0.98 0.050
29 91 24 0.41 1.9 4.3 0.12
30 11 8.0 0.46 3.7 0.041 0.0051 Run 12 13 14 15
31 37 20 0.46 3.7 0.83 0.042
32 88 37 0.46 3.7 4.0 0.11
33 11 7.9 0.46 19 0.029 0.0037
34 37 20 0.46 19 0.71 0.036
Run 17 18 19 20
Figure Al: Drop shapes
I 11I
I0
21
3 h
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Run 22 23 24 25 26 31
Run 32 33 34 27 28 29 30
Figure Al: Continued
