Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Effects of Bubbles on Turbulent Flows in Vertical Channels
Shigeo Hosokawa and Akio Tomiyama
Kobe University, Graduate School of Engineering, Department of Mechanical Engineering
11 Rokkodai, Nada, Kobe, Hyogo, 6578501, Japan
hosokawa@mech.kobeu.ac.jp
Keywords: bubbly flow, turbulence, vertical channel, LDV, photobleaching molecular tagging velocimetry
Abstract
Turbulence structures in bubbly flows are very complex due to (1) bubbleinduced pseudo turbulence, (2) interaction between
bubbleinduced and shearinduced turbulences, (3) direct interaction between bubbles and turbulence eddies and (4)
modification of shearinduced turbulence due to the modulation of mean velocity distribution due to bubbles. Although many
studies have been carried out to understand the turbulence structure, our knowledge is still insufficient for developing accurate
models on the effects of bubbles on turbulence properties in bubbly flows. This report reviews our studies on turbulence
properties in airwater bubbly upflow in vertical channels. The authors proposed an eddy viscosity ratio, which is the ratio of
the bubbleinduced eddy viscosity (defined as the product of bubble diameter dB and relative velocity VR between bubbles and
liquid) to the shearinduced eddy viscosity (defined as the product of turbulence length scale It and turbulence velocity u'), to
correlate turbulence modification due to bubbles. Since the eddy viscosity ratio (dBVR/ltu') is regarded as the product of the two
ratios, dBl, and VR/u', turbulent bubbly flows could be classified by using these two ratios. Three dimensional turbulence
properties in bubbly flows with large dBl, in a vertical pipe were measured by using an LDV system. The results indicate that
(1) bubbleinduced pseudo turbulence augments liquid turbulence intensity when VR/u' is large and (2) turbulence attenuation
due to eddy breakup induced by bubbles takes place at small VR/u' where the bubbleinduced pseudo turbulence is not
prominent. Turbulence properties in bubbly flow in a vertical duct at small /, / and large VR/u' were measured by using a
photobleaching molecular tagging velocimetry which enables us to evaluate turbulence kinetic energy budget. The results
imply that (3) bubbleinduced pseudo turbulence is not prominent at small ,/ / (4) bubble passage through eddies inhibits the
growth and extension of shearinduced turbulence from the near wall region to the core region, and therefore, (5) it reduces the
production rate of turbulence kinetic energy in the core region. Enhancement of shearinduced turbulence due to modulation of
mean velocity distribution caused by momentum transfer between bubbles and liquid is also observed in the bubbly flow.
Introduction
Turbulence structure in bubbly flows is very complex
due to (1) bubbleinduced pseudo turbulence, (2) interaction
between bubbleinduced and shearinduced turbulences, (3)
direct interaction between bubbles and turbulence eddies
and (4) modification of shearinduced turbulence due to the
modulation of mean velocity distribution due to bubbles.
Many studies have been carried out to understand the
turbulence structure. Serizawa et al. (1983, 1990) and Wang
et al. (1987) measured turbulent intensities and Reynolds
shear stress in bubbly flows using hotfilm probes. Shawkat
et al. (2008) also measured void fractions and liquid
velocities using dual optical probes and hotfilm probes, and
reported turbulence characteristics and void distributions in
bubbly flows in a large diameter pipe. Turbulence energy
spectra in bubbly flows have also been measured by using
hotfilm probes (Lance & Bataille, 1991). Marie & Lance
(1983), So et al. (2002) and Hosokawa & Tomiyama
(i21" '4L applied LDV to twophase bubbly flows in vertical
pipes and measured turbulence intensities. Fujiwara et al.
_'""I4) measured turbulence properties in bubbly flows
using PIV, and they evaluated turbulence production and
dissipation from measured twodimensional velocity
distributions. Although these studies provide important
knowledge on the influence of bubbles on turbulence
property, our knowledge is still insufficient for developing
accurate models on the effects of bubbles on turbulence
properties in bubbly flows.
The authors proposed an eddy viscosity ratio, which is
the ratio of the bubbleinduced eddy viscosity to the
shearinduced eddy viscosity, to correlate turbulence
modification due to bubbles (Hosokawa & Tomiyama,
2004a, 2004b). The bubbleinduce eddy viscosity is defined
as the product of bubble diameter dB and relative velocity
VR between bubbles and liquid, and the shearinduced eddy
viscosity is defined as the product of turbulence length scale
It and turbulence velocity u'. Since the eddy viscosity ratio
(dBV / ,') is regarded as the product of the two ratios, dB/,t
and VR/u', turbulent bubbly flows could be classified by
using these two ratios. This report reviews our studies on
the effects of bubbles on turbulence properties in bubbly
flows from the point of view of the two ratios, ,/ / and
VR/u'.
Nomenclature
interfacial area concentration (1/m)
drag coefficient for a single bubble
drag coefficient for bubbly flows
dimensionless change in turbulence intensity
pipe diameter (m)
Paper No
Diffusion rate of k (m2/s3)
bubble diameter (m)
aspect ratio
Edtv6s number
focal length (m)
gravitational acceleration (in 1
turbulence kinetic energy (nr s)
turbulence length scale (m)
volumetric flux (m/s)
drag force (N/m3)
interfacial momentum transfer (N/m3)
number density (1/m3)
pressure (N in )
production rate of k (m2/s3)
distance from pipe center (m)
pipe radius (m)
Reynolds number
nondimensional shear rate
time (s)
irradiation time (s)
turbulence velocity (m/s)
frictional velocity (m/s)
mean velocity (m/s)
velocity (m/s)
relative velocity (m/s)
distance from wall (m)
dimensionless distance from wall (m)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
magnitude of the change. The critical parameter can be
regarded as the ratio of a characteristic length scale of
turbulence induced by the dispersed phase to that of
shearinduced turbulence. Since the eddy viscosity might be
one of the most fundamental quantities representing
turbulence characteristics, the ratio 4 of the eddy viscosity
Vd induced by the dispersed phase to the shearinduced eddy
viscosity v, would be a more appropriate indicator for the
index of turbulence modification.
The eddy viscosity v,, of an incompressible twophase
flow is introduced in the following averaged NavierStokes
equation based on the ensemble average and the socalled
eddyviscosity assumption (Delhaye, et al., 1981):
aV a VP+V a(v+v,,)(VV +V,)
Dt p (1)
M
+ +a g
P
where VL is the liquid velocity, t the time, p the liquid
density, P the pressure, aL the volume fraction of liquid
phase, v the kinematic viscosity, M, the interfacial
momentum transfer and g the gravitational acceleration. The
vT can be related with the eddy viscosity v, of
shearinduced turbulence, the eddy viscosity vB induced by
the bubbles, and v. Hence, it can be expressed as
Greek letters
a void fraction
aL liquid volume fraction
8t time interval (s)
Sk dissipation rate of k (m2/s3)
4) eddy viscosity ratio
X wave length (m)
v kinematic viscosity (m2/s)
p density (kg/m3)
T twophase multiplier for k
co magnitude of liquid velocity gradient (1/s)
Subsripts
B bubble
d dispersed phase
G gas
L liquid
t shearinduced turbulence
SP single phase flow
TP twophase flow
Indicator of Turbulence Modification
Gore & Crowe (1989, 1991) investigated turbulence
modification caused by the addition of particles in a gas
flow, and pointed out that the modification is well correlated
with the socalled critical parameter, o' / the ratio of a
particle diameter d to a turbulence length scale i,. They
applied this parameter to turbulence modification caused by
bubbles and drops and confirmed that ,/ is also applicable
to gasliquid dispersed flows. They also reported that the
critical parameter refers only to the question of increase or
decrease in turbulence intensity and does not relate to the
Since the magnitude of v is much smaller than those of vt
and Vd, one of the most significant dimensionless groups
involved in Eq. (2) would be
=
The magnitude of vB can be evaluated as the product of the
bubble diameter dB and the relative velocity VR between
bubbles and liquid (Sato & Sekoguchi, 1975):
VB ~ dBVR
When the void fraction a is low, the shearinduced eddy
viscosity in a twophase flow can be evaluated as the eddy
viscosity in a singlephase flow with the same volumetric
flux of the continuous phase as the twophase flow. Hence,
the magnitude of vt can be estimated as the product of the
turbulence length scale it and the turbulence velocity u,' for
the singlephase flow as follows:
v, u'l,
where u' is the turbulent velocity. Hence, the magnitude of
the eddy viscosity ratio 6 can be estimated as
4 dBVR
Srlt
u'/,
The eddy viscosity ratio 4 is regarded as the product of the
V TP = AVVB1V ***).
Paper No
length scale ratio, dB/l (Gore & Crowe's parameter) and the
velocity scale ratio VR/U'. Serizawa & Kataoka (1995)
reported that turbulence modification in gasliquid
twophase bubbly flows depends on the liquid volumetric
flux JL, in other words, on the intensity of shearinduced
turbulence. Due to the lack of u,' and VR, (, / cannot explain
these results, whereas 4 possesses a potential of correlating
the effects of u,' and VR on turbulence modification. It
should be also noted that 4 is regarded as the ratio of the
bubble Reynolds number ReB ( = VRdB/v ) to the turbulence
Reynolds number Re ( = u' l/v ).
Figure 1 shows the correlation between 4 and a
dimensionless change in turbulence intensity per unit
number density, Cn/Nd, for airwater bubbly upflows and
liquidsolid twophase upflows in vertical circular pipes.
Here, a dimensionless change in turbulence intensity CTI and
the number density of dispersed phase Nd are defined by
UT, Us
C, (7)
USP
S6ad
Nd
td'
where u denotes the RMS of axial fluctuating velocity,
U the axial mean velocity, Cd the volume fraction of
dispersed phase and d the diameter of dispersed phase. The
subscripts TP and SP indicate a twophase flow and a
singlephase flow with the same liquid volumetric flux as
the twophase flow, respectively. The pipe diameter D, the
diameter of dispersed phase d, the liquid volumetric flux JL
and the terminal velocity of dispersed phase were 20 30
mm, 1.0 4.9 mm, 0.5 1.0 m/s and 0.21 0.42 m/s
respectively. The liquid velocity was measured by using an
LDV system. The details on the experiments can be found in
Hosokawa & Tomiyama (2'1i"'4 It is clear that the
turbulence modification Cn/Nd is well correlated with the
eddy viscosity ratio ). This result indicates that (1) the eddy
viscosity ratio is an appropriate indicator to correlate
turbulence modification, (2) turbulence modification
depends not only on the length scale ratio ,/ / but also on
the velocity scale ratio VR/ut', and (3) turbulence
modification could be classified by using two ratios, ,, /
and VRI/U'.
Turbulence Modification at Large dBlit
Turbulent flows with millimetersize bubbles are
frequently encountered in industrial equipment. The authors
measured three dimensional turbulence properties in
turbulent flows with millimetersize bubbles in vertical
circular pipes by using the LDV system (Hosokawa et al.,
2000, 2004a, 2009a). Since the diameter of the pipe D was
20 or 25 mm, the turbulent length scale 1, ( 0.2 D ) was
about 4 5 mm. The mean bubble diameter was about 3.0 
4.5 mm, that is, the length scale ratio dB/lt was about unity.
Since the liquid volumetric flux JL ranged from 0.5 to 1.0
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
m/s, the turbulence velocity u' was roughly estimated as
0.025 0.15 m/s ( = 0.05 JL 0.15 JL) and the velocity scale
ratio VR/U' 1 10. To enable three dimensional velocity
measurements, the pipe was made of Fluorinated Ethylene
Propylene resin (FEP) which has almost the same refractive
index as water.
The measured mean velocities are plotted against
distance r from the pipe center normalized by the pipe
radius R in Fig. 2. The axial mean velocities in turbulent
singlephase flows (case 1 and 2) showed the power law
distributions. The axial mean velocities were flattened due
to the presence of bubbles. The radial and azimuthal mean
velocities (V and I) were almost zero under all flow
conditions.
3
2
z"
Co
1 '
Figure 1: Correlation between 4 and CnTNd
0.5
Case 1
SJL 0.5 m/s
So 0 JG = 0.018 m/s 0
0 0 0
0
 B
*uaaaa 
072 0 4 6 08
I I 'c se'2 I .
Case 2
JL=1.0 m/s
6 9 JG=0.036 m/s
8
Single phase flow
o U/JL
o W/JL
Bubbly flow
U/J L
A V/JL
W/JL
i 02 04 06 8
r/R
Figure 2: Distributions of mean velocity
The measured Reynolds stresses are shown in Figs. 3
and 4. Among the Reynolds stress components, the axial
component u'u' is the largest and about ten times larger
than the Reynolds shear stress u'v'. The magnitudes of the
radial component v'v' and azimuthal component w'w
are approximately half of that of u'u From the
axisymmetry of the flow field, u'w' and v'w' should
vanish over the cross section. We could confirm the
axisymmetry of void distribution from the measured u'w'
GasLiquid I
M JL=0 5m/s,4=0 017m/s(D=30 mm)
4 JL=0 5m/s,4=0 023m/s(D=30 mm) O
. JL= 7m/s,4=0 017m/s(D=30 mm)
V JL=0 7m/s,4=0 023m/s(D=30 mm) o
A JL=1 Om/s,4=0 017m/s(D=30 mm)
X JL=1 Om/s,4=0 023m/s(D=30 mm)
* JL=0 5m/s,6=0 015m/s(D=20 mm) O
" JL=0 9m/s,4=0 020m/s(D=20 mm)
LiquidSoli(JL=0 5m/s, D=30mm)
o d=4 Omm J=0 002m/s
 d=4 Omm :=0 004m/s ao
0 d=2 5mm J=O 002m/s
Sd=2 5mm =0O 004m/s Oo
Sd=1 Omm J=0002m/s %*
B d=1 Omm J=0 004m/so o y
a
Paper No
which was almost zero under all flow conditions. All the
Reynolds stress components were affected by the presence
of bubbles. The three components of fluctuation velocity
(u'u', v'v' and w'w') in case 1 are enhanced due to the
presence of bubbles over the crosssection. On the other
hand, u'u' and w'w' in case 2 are attenuated due to the
presence of bubbles except in the near wall region, whereas
v'v' does not change so much due to bubbles. Although the
Reynolds shear stress u'v' in case 1 is enhanced due to
bubbles, the change in u 'v' is much smaller than that in
u'u' in the core region. This indicates that bubbleinduced
pseudo turbulence augments liquid turbulence intensity
whereas it does not enhance Reynolds shear stress so much.
On the other hand, u'v' in case 2 is attenuated due to the
presence of bubbles except in the near wall region, and the
attenuation of u'v' is larger than that of u'u' in the core
region.
0.05
Case 1 Case 2
JL=0.5 m/s JL=1.0 m/s
"0.04 JG=0.018 m/s JG=0.036 m/&
S Single phase flow
u'u'/J 2
S0.03 A ''/JL2
o w'w'JL2
Bubbly flow
o 2
,0.02 A VV/JL2
'1 0 W'W'/JL2
0.01 2
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1
r/R r/R
Figure 3: Diagonal components of Reynolds stress
Case 1 Case 2
JL=0.5 m/s JL=1.0 m/s
S4 JG=0.018 m/s JG0.036 m/s
SX. "Single phase2flow
S3 o u'v'JL 2
SA U'W'/JL
o Bubbly flow o
SU'V'/JL2 0
c22 A U'W'/JL2
*0 c o
0 o 
1 0
0 02 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1
r/R r/R
Figure 4: Reynolds shear stress components
To understand the difference in turbulence structure between
turbulence enhancement and attenuation cases, examples of
the time series velocity data in the axial direction are shown
in Fig. 5. In the case of turbulence enhancement, the axial
velocity is clearly agitated by bubbles, therefore the
amplitude of fluctuating liquid velocity in the bubbly flow is
larger than that in the singlephase flow. The axial velocity
is accelerated around bubbles and the amplitude of the
fluctuation is in the order of the terminal rising velocity of a
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
bubble in stagnant water (about 0.2 m/s). On the other hand,
the agitation due to bubbles is not significant in the case of
turbulence attenuation in which the amplitude of fluctuating
velocity is larger than the bubble terminal velocity. This
clearly indicates that the velocity scale ratio VR/u' is an
important parameter of the turbulence modification in
bubbly flows.
(a) Turbulence enhancement
Velocity data
'0.5
Single phase flow
(JL=O 93m/s,r/R=0 7)
1.5Bubblyflo elcitydaa Phase distnbton function
(J =0 93m/s.J =0 018m/s.r/R=0 7) .(High gas, Low iquid) g
0 0.2 0.4
(b) Turbulence attenuation
Figure 5: Time series data of velocity in axial direction
105
104
103
o (Hz)
Figure 6: Turbulence power spectra
Paper No
The low frequency fluctuation observed in the singlephase
flow in Fig. 5 (b) becomes weak in the bubbly flow. This
property clearly appears in the turbulence energy spectra
shown in Fig. 6. The energy of large scale eddies, which
have a low frequency, decreases due to the presence of
bubbles in the case of turbulence attenuation, whereas the
turbulence energy in the case of turbulence enhancement
increases in the whole frequency range. From the time
series data and the spectra, it is supposed that the bubbles
tend to break up the large scale eddies containing higher
energy than the small eddies. When VR/u' is small, the
bubbleinduced pseudo turbulence is not significant and
the turbulence attenuation due to the eddy breakup induced
by bubbles occurs.
The above results imply that the turbulence properties
in bubbly flow strongly depend on VR/u', and therefore, an
accurate evaluation of VR must be important in numerical
simulations of turbulent bubbly flows. The authors
therefore carried out a multifluid simulation of the bubbly
flow in the pipe (Hosokawa & Tomiyama, 2009a). For an
accurate evaluation of VR, we focused on drag force
models in the simulation. The drag force MD in the
momentum equations is given by
1
MD =aC DSP V, V (9)
8
where CDS is the drag coefficient for bubbly flows and ajw
is the interfacial area concentration:
6a
aINT (10)
d,
The drag coefficient CDs (Tomiyama et al., 1995a) is give
by
Cs= CD2' (11)
DS= DOL
where CD is the drag coefficient for a single bubble and
aCL321 ( = 1.75 ) is the multiplier to account for the effects
of bubble swarm. Two types of drag correlations were
tested. One is an empirical correlation given by (Tomiyama
et al., 1995b)
CD =1max[min16 (1+0.15Re 6848 8 (12)
L Re ReJ 3Eo+4J
where Re is the bubble Reynolds number and Eo the
Edtvos number. The other is a theoretical correlation which
takes into account the effect of bubble aspect ratio E on CD
(Tomiyama et al., 2002a).
8 Eo
CF = EF 2 (13)
3 E2/3(1E2) 1Eo+16E43
sin 2 E E2
F sin 1 E (14)
IE2
Although the above equation can predict bubble terminal
velocities not only in an infinite liquid but also in a
stagnant liquid in a channel (Tomiyama et al., 2003), it
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
requires a correlation of E. The correlation of E can be
found in Hosokawa & Tomiyama (2009a).
Legendre & Magnaudet (1998) pointed out that the
presence of liquid velocity gradient increases the drag
force acting on a bubble, and proposed the following drag
multiplier ,,r:
sr 1+0.55Sr2
where Sr is the nondimensional shear rate defined by
Sr= dB
V
where co is the magnitude of the liquid velocity gradient.
Hence, the drag coefficient in Eq. (11) is given by
CDS = CDa 5 (1+0.55Sr2)
Numerical simulations were carried out using three sets of
drag force models; (a) Eq. (12), (b) Eq. (13), and (c) the
combination of Eqs. (13) and (17).
Figure 7 shows comparisons between measured and
predicted relative velocities. The predictions with Set (a)
overestimate VR. Although Set (b) well predicts VR in the
core region, it overestimates VR in the near wall region
especially in Case 2, in which the velocity gradient in the
near wall region is high. On the other hand, the predictions
with Set (c) agree well with the measured VR. Figure 8
shows comparisons between measured and predicted
turbulence kinetic energies with Sets (a) (c). The
accuracy of the prediction is improved especially in case 2
by using Set (c) instead of Set (a). This must be because
the turbulence kinetic energy depends on the relative
velocity, i.e. the drag force. This result implies the
importance of the effects of the bubble shape and shear
rate on the drag force for accurate predictions of turbulent
bubbly pipe flows.
Figure 9 shows measured and predicted distributions of
turbulence kinetic energy k for singlephase and twophase
flows. The Set (c) was used for CD in the predictions. The k
is augmented due to the presence of bubbles at low JL
(Cases 1), while it is attenuated in the core region at high
JL (Cases 2). The prediction gives the same trend for the
effect of JL on turbulence modification. These results
indicate that the RANS turbulence model (Lopez de
Bertodano et al., 1994) has a potential to predict not only
turbulence enhancement but also turbulence attenuation
resulting from the flattening of VL distribution caused by
the wall peaking void distribution at high JL. The broken
curves in Fig. 9 are the predictions using a simple
turbulence model based on the eddy viscosity ratio 4
(Hosokawa & Tomiyma, 2004a, 2004b), in which a
twophase multiplier T for the turbulent kinetic energy is
introduced:
k= y 2k,
where T is a function of the magnitude of the turbulence
modification Crz:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
=C, +1 (19)
As shown in Fig. 1, Cn/Nd is well correlated with the eddy
viscosity ratio 4. Hence, we proposed the following
empirical correlation of Cn/Nd (Hosokawa & Tomiyama,
2004b):
E
> 0.
r/R
r/R
Figure 7: Relative velocity between bubbles and liquid
r/R r/R
Figure 8: Turbulence kinetic energy in bubbly flows
Single phase flow
Measured
 Predicted
0.015 .
c"0.01
)(n
N
E
0.005
Twophase flow
o Measured
 Predicted (Lopez de Bertodano et al.)
 Predicted (Hosokawa & Tomiyama)
SCase 1
JL=0.5 m/s
JG=0.018 m/s

0 0 0 0
Case 2 a
JL=1.0 m/s
JG=0.036 m/s
0/
0 
, i , , , , , .
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8
r/R r/R
Figure 9: Turbulence modification due to bubbles
CT/Nd =2.54x107 ( 1.12)
4 Casel Case2
.JL=0.5m/s .JL=1.0m/s
JG=0.018m/s JG=0.036m/s
 .. .. .......
+34
2 0 0 00
o Measured + oy
1  Set (a)
S+Set (b)
 Set (c)
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8
Paper No
The predictions using this model agree well with the
experimental data and the predictions using the Lopez de
Bertodano's model. Although further studies on CTI are
required, the latter turbulence model has also a potential to
predict the turbulence modification in bubbly pipe flows.
From the point of view of accurate prediction, k in bubbly
flows is slightly underestimated in the core region at low JL,
and it is slightly overestimate in 0.4 < r/R < 0.8 at high JL.
These results indicate that improvement of bubbleinduced
pseudo turbulence and modeling of eddy breakup due to
bubbles are required for accurate simulations of bubbly
flows.
Turbulence Modification at Small ddlt and Large
VRIU'
The results at large ,U / in the former section indicate
that (1) bubbleinduced pseudo turbulence augments liquid
turbulence intensity when VR/u' is large and (2) turbulence
attenuation due to eddy breakup induced by bubbles takes
place at small VR/u' where the bubbleinduced pseudo
turbulence is not prominent. Turbulence modification at
small d/t1, and large VR/u' is discussed in this section. From
the point of view of eddy viscosity ratio 4, turbulence
attenuation is expected at small dB/t1 and large VR/u' since
small ) indicates turbulence attenuation as shown in Fig. 1.
On the other hand, large VR/u' might induce large velocity
fluctuation as shown in Fig. 5(a), and therefore, turbulence
enhancement is also expected at small o' / and large VR/u'.
Hence, we examined whether turbulence enhancement or
attenuation occurs at small dB/1 and large VR/u' (Hosokawa
et al., 2010). We also examined the turbulence kinetic
energy budget when the eddy breakup by bubbles occurs.
To understand turbulence kinetic energy budget, a
molecular tagging velocimetry based on photobleaching
reaction (PBMTV) was developed (Hosokawa &
Tomiyama, 2004c, 2009b) and applied to turbulent bubbly
flows with large VR/u' and small d/lt in a vertical square
duct (Hosokawa et al., 2010). The width W and length of the
duct were 50 mm and 1100 mm, respectively. The liquid
volumetric flux and the liquid Reynolds number Re were
0.06 m/s and 5000, respectively. The bubble diameters dB
was about 0.6 mm. Since the turbulence length scale I, can
be evaluated as 10 mm ( = 0.2 W ), the length scale ratio
dB/l, was about 0.06. The Kolmogorov scale of the turbulent
flow was about 200 itm, and therefore, the bubble size was
comparable to the Kolmogorov scale. The relative velocity
of the bubbles was about 0.06 m/s, and therefore, the
velocity scale ratio VR/u' was about 10 (u' is about 0.1 JL =
0.006 m/s.). Hence, this turbulent bubbly flow can be
classified as a small I / i lnd large VR/u' system.
The velocity components and velocity gradients were
measured at 1000 mm above the bottom of the vertical
duct using the photobleaching molecular tagging
velocimetry (PBMTV). Note that the flow in the
measurement position was fully developed. A brief
introduction of PBMTV is given below.
The photobleaching is an irreversible photodegradation
I
Paper No
of a fluorescent dye induced by intense illumination.
Figure 10 shows a schematic of PBMTV. Uranin
(fluoresceine sodium salt) was solved in water as a
fluorescent dye, which was photobleached by an intense
laser beam (wave length X = 488.0 nm). Its concentration
was about 106 mol/m3. The intense beam induces a
photochemical reaction which changes the fluorescent dye
to nonfluorescent dye. The intense beam passed through a
rectangular mask and focusing lens (focal length f = 76
mm) and then formed a tag on the plane of a laser sheet as
shown in Fig. 10. A low intensity laser sheet (L=514.5 nm)
was used to visualize the tags on the plane in a flow field,
and the image was recorded by the CCD camera
(DANTEC, HiSense MkII, 1280x1012 pixels) with a
macro lens (VS Technology Co., VSTC665). Tags
formed by the intense beam appeared as dark regions, i.e,
nonfluorescent regions, in the images as shown in the
right photographs in Fig. 10. Time sequences of the intense
laser beam and shutter of the camera are also shown in the
figure. First, the intense laser beam was irradiated for the
irradiation time T,, ( = 1 ms ) to form the tag. The
irradiation time was determined so as to obtain a high
contrast tag image and controlled by using the shutter.
After the irradiation of the intense laser beam, the CCD
camera recorded two consecutive images with the time
interval 6t ( 5 ms ). A delay generator (Stanford Research
Co. Ltd., DG535) was used to synchronize the shutter and
the CCD camera.
Mask
Laser sh e Laser mm 0
beam
Lens
Intense lose
laser
bea
Tag
Shutter CCD camera t
Camera Open
Close
ON
Intense laser beam
OFF IrT '
Time
Figure 10: Schematic of measurement system
Liquid velocities on the plane can be evaluated from
the displacement of the center of the tagged region
between the two consecutive images. The center positions
of tagged regions were calculated from the intensity
distributions in the tagged regions. Then, the displacement
vector of the center was divided by 6t to obtain the
twodimensional components (uo, vo) of the liquid velocity.
The deformation of the tagged region reflects the velocity
gradient tensor of the flow field. A simple searching
method based on crosscorrelation was adopted to evaluate
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the velocity gradients.
The measurements were carried out for xy and xz
planes. The size of the tag was 40 100 gLm and about 70
pixels in the recorded image. The tag size was smaller than
the Kolmogorov scale (about 200 gLm). The spatial
resolution of images was 1.1 gtm/pixel. The local mean
velocity and turbulence intensity were calculated from the
measured instantaneous velocities for about 200 tags at
each measurement point. The measurement uncertainties
for instantaneous velocity and velocity gradients were
0.312x104 m/s and 0.143 s', respectively.
The local axial mean velocity U is shown in Fig. 11
together with the void fraction a. Here, y and JG are the
distance from the wall and the gas volumetric flux,
respectively. The gradients of the axial velocity in the
bubbly flows are larger than that in the singlephase flow
in the near wall region. The axial velocity distribution in
the core region is flattened by the presence of bubbles in
the near wall region as many researches reported.
Figure 11: Axial velocity distribution
Figures 12 14 show the streamwise, wallnormal and
spanwise components of turbulence intensities u'2 v'2
w2 The turbulence intensities are enhanced in the near
wall region (y/W < 0.06) and they are attenuated in the
core region (y/W > 0.1). In the singlephase flow, the peak
positions and distributions are different among u'2, v'
w The difference among the horizontal distributions of
u'2 v2 w2 becomes small by addition of bubbles,
and the distribution is more or less flat in the core region.
This indicates that the wall peaking void fraction makes
the region of turbulence production move toward the wall.
The magnitudes of u'2 v'2 w2 in the core region of
the bubbly flows are almost the same, that is, the
turbulence is isotropic.
Figure 15 shows the Reynolds shear stress u'v'
Although the Reynolds stress is also enhanced due to the
presence of bubbles in the near wall region, the increase in
the peak value in the bubbly flow is not so large. To the
contrary, u'v' largely attenuates in the core region.
Since the above results indicate that strong turbulence
modification occurs in the near wall region, the turbulence
kinetic energy k is plotted against y+ ( = yu,/v ) in Fig. 16
together with the void fraction a. The frictional velocity u,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
n A
( = vOU ) was calculated by substituting the velocity
gradient OU/ v measured at the measurement point
nearest to the wall (the distance between the wall and the
measurement point was 250 rtm, y+ 2). The turbulence
kinetic energy k increases with JG in the near wall region
whereas it decreases due to the bubbles in the core region.
The y+ position of the maximum k shifts toward the wall by
addition of bubbles. This clearly indicates that the bubbles
induce velocity fluctuation in the near wall region, and that
the buffer region, where the active turbulence generation
occurs, extends toward the wall.
0 0.1 0.2 0.3
y/W
U ... .  u J
0.4 0.5
Figure 12: Streamwize component of turbulence intensity
U.0O.Ox i U
JG [mn/s 2 a
0.0x1 0
0.6 2.5x105 o
a 5.0x105 A 
7.5x105 
0.4 e
0 0
0.2 a
C' i
0 0.1 0.2 0.3
yNV
0.4 0.5
Figure 13: Wallnormal component of turbulence intensity
U)
E
S0.4
0.2
0 0.1 0.2 0.3
y/W
10.2
0.4 0.5
Figure 14: Spanwise component of turbulence intensity
1 0.3
E
o 0.2
x
I:
JG [m/st u'v' a
0.0x10 i
2.5x105 o
5.0x105 A
7.5x105 
o
00
0 0
o
a 0
"**
S.i i.
0.1 0.2 0.3
y/V
I
0.3
0.2
0.1
0.4 0.5
Figure 15: Reynolds shear stress
3 0.3
JG [m/s k a
0.0x10"
2.5x105 o 
5.0x10 A 
2 7.5x10 5  0.2
A A
%0
o *
1 < 0.1
.I*^.
1 0 1 .
Y
Figure 16: Turbulence kinetic energy
The turbulence production Pk, diffusion Dk and
dissipation Ek were evaluated by substituting the measured
local instantaneous velocities and velocity gradients into
the following simplified equations:
J x,
c9 1 _&+ k
Dk 8x 2 u pj'u +v 
Here, we assumed that the flow is fully developed and
symmetric with respect to y axis. Note that the convection
of turbulence is vanished in the fullydeveloped duct flow.
The evaluated production Pk, dissipation Sk and diffusion
Dk are shown in Figs. 17 19. Pk in the bubbly flow is
augmented fory+ < 20 due to the presence of bubbles, and
increases with JG. As JG increases, the position where Pk
takes a peak value moves from+ = 10 15 to y+ = 8 where
the void fraction takes a peak. Pk in the linear sublayer (y
< 2 3) takes low values and does not depend on JG so
much. The turbulence dissipation rate Ek is also enhanced
for y+ < 20 due to the presence of bubbles, and increases
with JG. The peak position of Ek also moves toward the
wall due to the presence of bubbles. The diffusion rate Dk
increases in the linear sublayer and decreases in the buffer
Paper No
JG [m/sl
0.0x10
 2.5x105
5.0x105
.Ap 7.5x105
o
a
'
9
lK^
Uf a
o 
A .. 0.2

JG [m/s] w ca
0 0.0x10
 2.5x105 o
5.0x105 A 
7.5x105 
A .
0,
O Og I
. . .. J
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
layer as JG increases. The position of the lowest peak of Dk
corresponds to the peak position of a. Since the positive
value of Dk indicates the influx of k, this result implies that
a part of k generated in the buffer layer is diffused by
bubbles and it dissipates in the linear sublayer.
1
0.8
(,)
' 0.6
x0.4
0.2
JG [m/s] P ac
0 0.0x10 5 *
2.5x10_5 o 
5.0x10s5 A 
7.5x105 a
D0 101 102
Figure 17: Production rate of k
JG [m/s] s ac
.0x105 *
2.5x105 o 
5.0x10'5 A 
7.5x10 5 
0
y
Figure 18: Dissipation rate of
o 101 102
Figure 18: Dissipation rate of k
0.4 JG [m/s D a
O.Oxl10
0 2.5x10 o 
m 0.2 5.0x105 
N 7.5x105 
E 0o
0
00.2
/ '\,
0.4 
100 101 102
y
H I I I
2 0.6
N`
E
0.4
0.2
X
<
0 0.02 0.04 0.06
a [%]
0.1
0.08 0.1
Figure 20: Relation between AP, As and a
n '
0.1
Figure 19: Diffusion rate ofk
0.6
Sn A
Figure 20 shows a relation between the modifications
of the local Pk and sk due to bubbles, APk and Ask, and the
local void fraction a in the near wall region (y+ < 20 ). APk
and Ask clearly increase with a, and there are strong
correlations between APk and a and between Ask and a.
The difference between the gradients, dAPJkda and
dAsk/da, is not large.
Figures 21 23 show the dimensionless production Pk,
dissipation Sk and diffusion Dk rates of k, which were
calculated by dividing Pk, Ek and Dk by u,/v. The
dimensionless Pk in the bubbly flow does not change due
to the presence of bubbles for y < 8. The dimensionless Pk
does not depend on JG under the present experimental
conditions. This result indicates that the bubbles increase
the frictional velocity due to the interfacial momentum
transfer, and this increase in the frictional velocity
enhances the shearinduced turbulence in the near wall
region. Hence, the main effect of the bubbles is the
increase in the liquid velocity gradient in the near wall
region, and the pseudo turbulence induced by the bubbles
is not significant. Since the bubble Reynolds number is
about 60 and the bubble size is comparable to the
Kolmogorov scale, the bubbleinduced turbulence
dissipates rapidly. The difference in the dimensionless e
between the singlephase and the bubbly flows in y < 8 is
not so large. The dimensionless Pk and sk decrease iny+ > 8.
This is not only due to the decrease in the gradient of axial
mean velocity but also due to the eddy breakup induced by
bubbles. That is, bubble passage through eddies inhibits
the growth and extension of shearinduced turbulence from
the near wall region to the core region, and therefore, and
it reduces the production rate of turbulence kinetic energy
in the core region. The dimensionless Dk also does not
change so much in y < 5, whereas it tends to decrease at
y  8 due to the presence of bubbles. Hence, the
bubbles might enhance the diffusion effluxx) of k due to
horizontal velocity fluctuation around the bubbles.
0.1
Figure 21: Dimensionless production rate of k
Paper No
JG [m/s] AP As
2.5x10 o *
5.0x105 A A
7.5x105 0
a
A
S A
A o
A0 *
~A 0
%
U I  .
JG [m/s] Pv/u a4
O.0x105
2.5x105 o 
5.0x105 A 
7.5x105 
0Do 101 + 102
y
Paper No
1
0.8
0.6
~ 0.4
0.2
0
1
..A. .. . ..
101
y
0.1
0
Figure 22: Dimensionless dissipation rate of k
4 JG [m/s] Dv/u,4 a
O.Oxl05 0
2 2.5x105 o
E2 5.0x105 A 
o 7.5x105 0 
0 I .N .,
o0
2
4
10 101 102
y
Figure 23: Dimensionless diffusion rate of k
In > 20, Dk is small, and therefore, Pk balances with sk in
the core region. These results indicate that the bubbles
modulate the mean velocity profile due to the interfacial
momentum transfer and enhance the frictional velocity and
shearinduced turbulence generated in the near wall region.
However, the bubbleinduced pseudo turbulence is not
prominent because the bubble diameter is comparable to
the Kolmogorov scale.
Conclusions
This report reviews our studies on turbulence
properties in airwater bubbly upflow in vertical channels.
An eddy viscosity ratio, which is the ratio of the
bubbleinduced eddy viscosity (defined as the product of
bubble diameter dB and relative velocity VR between bubbles
and liquid) to the shearinduced eddy viscosity (defined as
the product of turbulence length scale It and turbulence
velocity u') was proposed to correlate turbulence
modification due to bubbles. Since the eddy viscosity ratio
(dBVR/tu') is regarded as the product of the two ratios, dBlt/
and VR/u', turbulent bubbly flows was classified by using
these two ratios. Three dimensional turbulence properties in
bubbly flows with large dBlt in a vertical pipe were
measured by using an LDV system. Turbulence properties in
bubbly flow in a vertical duct at small dB/l and large VR/u'
were also measured by using a photobleaching molecular
tagging velocimetry which enables us to evaluate turbulence
kinetic energy budget. Several effects observed in the results
are summarized in Fig. 24, that is, (1) bubbleinduced
JG [m/s] sv/U4 a
O.Oxl105 0
2.5x105 o
5.0x105 A .
7.5x105
 /,'. ",, il
0 o
 I I I
00
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
pseudo turbulence augments liquid turbulence intensity
whereas it does not enhance Reynolds shear stress so much,
when VR/u' is large, (2) turbulence attenuation due to eddy
breakup induced by bubbles takes place at small VR/u' where
the bubbleinduced pseudo turbulence is not prominent, (3)
bubbleinduced pseudo turbulence is not prominent at small
dB/lt, (4) bubble passage through eddies inhibits the growth
and extension of shearinduced turbulence from the near
wall region to the core region, and therefore, (5) it reduces
the production rate of turbulence kinetic energy in the core
region, and (6) enhancement of shearinduced turbulence
due to modulation of mean velocity distribution caused by
momentum transfer between bubbles and liquid is also
observed in the bubbly flow.
S Bubbleinduced turbulence
Turbulence augmentation
Inhibition of growth of
Sshearinduced turbulence
S Reduction of production rate in core region
Breakup of turbulence eddies
due to bubbles
Modulation of mean velocity distribution
1 due to momentum transfer
between phases
Enhancement of
shearinduced turbulence
Figure 24: Several effects of bubbles on turbulent flow
These results indicate that the turbulence properties in
bubbly flows depend on ,' / and VR/u'. Figure 25 shows
rearrangement of the present experimental conditions and
the knowledge obtained so far on the ; I VR/u' plane.
Turbulence enhancement due to bubbleinduced turbulence
is remarkable when the eddy viscosity ratio 4 is large. To
the contrary, the effect of bubbleinduced turbulence is weak
in the case of small 4, and therefore, the other effects
become observable. When VR/u' is large, the effects of VR,
that is, the eddy breakup by bubbles and the modulation of
mean velocity distribution due to momentum transfer
between the phases, become strong. These effects have been
clearly observed in the experiments using the square duct at
large VR/u' and small dBlt. Although the effects of bubble
volume might be important at large dBlt, the experiments in
this region are still insufficient to understand the volume
effect. When both U' / and VR/u' are very small, the effects
of bubbles might be negligible. Although further
experiments are required for complete classification of
turbulent bubbly flows in the o; / VR/u' plane, the
importance of the classification to understand the structure
of turbulent bubbly flows has been demonstrated in this
report.
Paper No
Rich experimental data
dB/It
Figure 25: Classification of turbulent bubbly flows
Acknowledgements
The authors gratefully acknowledge the assistance of
experiments by Mr. Takashi Suzuki in Kobe University and
the financial support by Japan Society for the Promotion of
Science (grandinaid for scientific research (B) No.
21360084).
References
Delhaye, J.M., Giot, M. and Riethmuller, M.L.,
Thermohydraulics of TwoPhase Systems for Industrial
Design and Nuclear Engineering, A von Karman Institute
Book, Hemisphere (1981).
Fujiwara, A., Minato, D. and Hishida, K., Effect of bubble
diameter on modification of turbulence in an upward pipe
flow, International Journal of Heat and Fluid Flow, vol. 25,
pp. 481488 (2 '4).
Gore, R. A. and Crowe, C. T, Effect of particle size on
modulating turbulent intensity, Int. J. Multiphase Flow, Vol.
15 (2), pp. 279285 (1989).
Gore, R.A. and Crowe, C.T., Modulation of Turbulence by a
Dispersed Phase, Tans. ASME, Vol. 113, pp. 304307 (1991).
Hosokawa, S., Tomiyama, A., Takesaka, K. and Kondo, Y,
LDV Measurement of Reynolds Stresses in GasLiquid
TwoPhase Bubbly Flow in a Vertical Pipe, 2nd Japankorea
Symposium on Nuclear Thermal Hydraulics and Safety
(NTHAS2), pp.247252 (2000).
Hosokawa, S. and Tomiyama, A., Turbulence modification
in gasliquid and solidliquid dispersed twophase pipe
flows, International Journal of Heat and Fluid Flow, vol.
25, pp. 489498 (21i 114,
Hosokawa, S. and Tomiyama, A., Effects of Relative
Velocity between Phases on Turbulence Modification in
Dilute GasSolid TwoPhase Upflows in a Vertical Pipe,
Japanese Journal of Multiphase Flow, vol. 18, no. 3, pp.
255262 (2" '14b), in Japanese.
Hosokawa, S. and Tomiyama, A., Molecular tagging
velocimetry based on photobleaching reaction and its
application to flows around single fluid particles,
Multiphase Science and Technology, Vol. 16, Issue 4, pp.
335353 (2" 14c).
Hosokawa, S., and Tomiyama, A., MultiFluid Simulation of
Turbulent Bubbly Pipe Flows, Chemical Engineering
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Science, Vol. 64, pp. 5308 5318 (2009a).
Hosokawa, S. and Tomiyama, A, Application of
Photobleaching Molecular Tagging Velocimetry to
Turbulent Bubbly Flow in a Square Duct, Experiments in
Fluids, Vol. 47, pp. 745 754 (2009b).
Hosokawa, S., Suzuki, T and Tomiyama, A., Effects of
Bubbles on Turbulence Properties in a Duct Flow,
Multiphase Science and Technology, (2010) (In print).
Lance, M. and Bataille, J., Turbulence in the liquid phase of
a uniform bubbly airwater flow, Journal of Fluid
Mechanics, Vol. 222, pp. 95118 (1991).
Legendre, D. and Magnaudet, J., The lift force on a
spherical bubble in a viscous linear shear flow, Journal of
Fluid Mechanics, Vol. 368, pp. 81126 (1998).
Lopez de Bertodano, M., Lahey, R.T. and Jones, O.C., Phase
Distribution in bubbly TwoPhase Flow in Vertical Ducts,
International Journal of Multiphase Flow, Vol. 205, pp.
805818 (1994).
Marie, J.L. and Lance, M., Turbulence Measurements in
TwoPhase Bubbly Flow Using Laser Doppler Anemometry,
Measurement Techniques in Gas Liquid TwoPhase Flow
(Edited by Delhaye, J.M. and Cognet, G. ), pp. 141148,
Springer, Berlin (1983).
Sato, Y. and Sekoguchi, K., Liquid Velocity Distribution in
TwoPhase Bubbly Flow, Int. J. Multiphase Flow, Vol. 2, pp.
7995 (1975).
Serizawa, A., Tsuda, K. and Kataoka, I., Measurement
Techniques in GasLiquid TwoPhase Flows, IUTAN
Symposium, Nancy, France, Hemisphere (1983).
Serizawa, A. and Kataoka, I., Turbulence Suppression in
Bubbly TwoPhase Flow, Nuclear Engineering and Design,
vol. 122, pp. 116 (1990).
Serizawa, A. and Kataoka, I., Dispersed FlowI, Multiphase
Science and Technology, Vol. 7, pp. 125194 (1995).
Shawkat, W.E., Ching, C.Y. and Shoukri, M., Bubble and
liquid turbulence characteristics of bubbly flow in a large
diameter vertical pipe, International Journal of Multiphase
Flow, Vol. 34, pp. 767785 (2008).
So S., et al., Laser Doppler velocimetry measurement of
turbulent bubbly channel flow, Experiments in Fluids,
Vol.33, pp.135142 (2002).
Tomiyama, A., Kataoka, I., Fukuda, T and Sakaguchi, T,
Drag Coefficients of Bubbles: 2nd Report, Drag Coefficient
for a Swarm of Bubbles and its Applicability to Transient
Flow, Transactions of JSME, Vol. 61, Issue 588, pp.
28102817 (1995a), in Japanese.
Tomiyama, A., Kataoka, I. and Sakaguchi, T, Drag
Coefficients of Bubbles: 1st Report, Drag Coefficients of a
Single Bubble in a Stagnant Liquid, Transactions of JSME,
Vol. 61, Issue 587, pp. 23572364 (1995b), in Japanese.
Tomiyama, A., et al., Terminal velocity of single bubbles in
surface tension force dominant regime, International
Journal of Multiphase Flow, Vol. 28, pp. 14971519
(2002a).
Tomiyama, A., Nakahara, Y, Adachi, Y and Hosokawa, S.,
Shape and Rising Velocities of Single Bubbles rising
through an Inner Subchannel, Journal of Nuclear Science
and Technology, Vol. 40, Issue 3, pp. 136142 (2003).
Wang, S.K., Lee, S.J., Jones Jr., O.C. and Lahey Jr., R.T.,
3D turbulence structure and phase distribution
measurements in bubbly twophase flow, International
Journal of Multiphase Flow, Vol. 13, no. 3, pp. 327343
(1987).
