7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Characteristics of heat transfer in turbulent bubbly upflow in a vertical channel
M. Tanaka*, N. Matsui*, Y Miyajima*, and Y Hagiwara*
Department of Mechanical and System Engineering, Kyoto Inst. of Tech., Matsugasaki, Kyoto 6068585, JAPAN
Ililllli ki l ( kilic ip)
Keywords: Direct numerical simulations, Bubbles, Heat transfer, Vortices
Abstract
Direct numerical simulations have been conducted for turbulent bubbly upflow between two parallel heating walls in
order to clarify its heat transfer characteristics. A PLICVOF method was used for tracking the bubble interfaces, and
a balanced force algorithm was used to suppress spurious currents. The channel Reynolds number was set at 3768,
and the flow rate was kept constant. Air bubbles with a diameter of 12mm in water were considered. It is found that
the bubbles tend to collect on the nearwall regions of the channel and that vortices are activated by bubbles especially
in the nearwall regions. The wall shear stress and the Nusselt number are increased by a factor of about 1.6 and
1.3 respectively due to the effects of the bubbles. By examining the budget of the heat transfer in the wallnormal
direction, it is found that the enhancement of heat transfer is mainly caused by the increase in turbulent heat flux.
The vortices near the walls play a major role in the heattransfer enhancement. It is also found that the heattransfer
enhancement is more noticeable at higher Prandtl numbers.
Introduction
Turbulent bubbly flows have attracted a lot of attention
because of their importance for many practical applica
tions such as flows in chemical plants and nuclear power
plants. Enhancement of heat transfer by bubbleinduced
turbulence also attracts a lot of attention from the view
point of energy saving. Many studies have been con
ducted for the motion of bubbles and the characteristics
of heattransfer in turbulent bubbly flows.
The characteristics of bubbly flow strongly depends
on the motions of bubbles and resulting void fraction
distribution in the flow. Serizawa et al. (1975) found that
the local void fraction is high near the walls and is lower
in the core region of upflow in a pipe. Liu (1993) also
found in the experiments of turbulent bubbly upflow in a
vertical channel that the void fraction has peaks near the
walls for the bubbles smaller than 56mm, while it has
a peak in the core of the channel for the bubbles larger
than 56mm. Lu & Tryggvason (2008) also showed in
their direct numerical simulations of turbulent bubbly
upflow in a vertical channel that nearly spherical bub
bles tend to concentrate on the nearwall regions, while
strongly deformable bubbles tend to be expelled from
the nearwall regions. They also showed that the turbu
lence structures are changed by the motions of bubbles.
The detailed mechanism of turbulence modulation due
to the bubbles, however, has not been fully clarified yet.
Some experimental studies have been conducted for
heattransfer enhancement by the injection of bubbles.
Tamari & Niishikawa (1975) showed in their experi
ments of laminar natural convection heat transfer in wa
ter that the heat transfer is enhanced by the injection of
air bubbles. The enhancement of heat transfer by bub
ble injection was studied further in detail by Tokuhiro &
Lykoudis (1994) and Kosige et al. (2006). However, the
mechanism of the heattransfer enhancement has not yet
been fully clarified especially in turbulent flows.
In the present study, direct numerical simulations have
been conducted for turbulent bubbly upflow between
two parallel heating walls in order to clarify its heat
transfer characteristics.
Direct Numerical Simulations
Configuration. Fig.l shows the spactial configuration
considered in the present study. The x, y and z axes
were assigned to the streamwise, wallnormal and span
wise directions, respectively. Here, we consider the sys
tem where two parallel walls are heated and the mean
temperature increases linearly in the streamwise direc
tion. The wall heat flux, q,, is determined so that the
energy of the system is conserved. The wall heat flux is
kept nearly constant after the turbulence reaches a steady
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Y =bF + (1 F) /i. (6)
Here, p and p are the average values of the fluid den
sity and the viscosity in the grid cell, respectively. The
subscripts b and I denote the gas and liquid phases, re
spectively. p is the pressure, fs represents the volume
force associated with the interfacial surface tension of
the bubbles, and g denotes the gravitational acceleration.
Overlines represent the spatial average.
Since the mean temperature increases linearly down
stream, the temperature field is decomposed into the
mean and the variation parts as T 0 + Gxz, where
T is the temperature, O represents the temperature vari
ation, and G denotes the mean temperature gradient in
the streamwise direction.
The governing equation for the temperature variation,
O, is described by
Fig. 1: Spacial configuration.
state. It is assumed that the fluids are incompressible for
both the liquid and gas phases in the present study. The
buoyancy force originated from the nonuniform spatial
distribution of temperature is assumed to be negligibly
small. It is also assumed that the fluid properties do not
depend on the temperature.
Basic Equations. In this study, simulations of the
bubbly flow are conducted by using VOF (Volume of
Fluid) method. In the VOF method, the fraction of the
bubble gas, F, occupying each computational grid cell
is advected by using the equation
OF OF
at ax' (1)
where x = (x, x2, X3) = (x, y, z) denotes the position,
and u = (ul, u2, 3) (u, v, w) represents the fluid
velocity. In the present study, summation convention is
applied to repeated subscripts if not otherwise specified.
The motions of incompressible fluid are governed by
the NavierStokes (momentum) equations
ui + aui,
Kat xaj)
x+ + fsi + (P
Oxi Oxj
supplemented with the continuity equation
Oui
= 0,
where
Oui )
Ox
is the viscous stress, and
SpbF + (1 F) pi,
P) gi
(pCp) + U0 + Gul)
where
(pC,) =PbCpbF + (1 I
k =kbF+ (1 F)
a (K ae\
(7)
1 PCpi (8)
ki (9)
represent the phaseaveraged specific heat capacity and
heat conductivity, respectively. C, denotes specific heat.
Control Parameters. The parameters which con
trols the system considered here include the frictional
Reynolds number Re,, E6tvds number Eo, Morton
number M, Archimedes number Ar, Prandtl number
Pr, the density ratio Rp, the viscosity ratio R,, and the
ratio of specific heat Rcp. They are defined as
Re, = h,Eo ,M
Vi (PT p1(
A 2 7P
Ar = Pr
p '
Pb Cpb
RP i 0 Cp 1
VR
(k/pC,)'
where u,, h, v, do, and a denote the friction velocity,
(2) channel half width, kinematic viscosity, bubble diame
ter, and the surface tension coefficient, respectively. We
set h to be unity in the present study. Eotvos number
(3) Eo represents the ratio of the buoyancy and surface
tension forces, while Archimedes number Ar represents
the ratio of the buoyancy and viscous forces. Morton
(4) number M is not independent from Eotvos number Eo
and Archimedes number Ar and is expressed as M
Eo3/Ar2. The parameter B/(gAp), which is composed
of the mean pressure gradient B dp/dx+ pg, the gray
(5) itational acceleration g, the density difference between
0
C
FI
qwz const.
Y _1
1g
Tij ui
Oxy
two phases Ap = pi Pb, is also an important parame
ter. The buoyancy effect of bubbles may be large when
this parameter is small.
Hereafter, the velocity is normalized by the friction
velocity, which is expressed by
U , (11)
V Pi
where Tr denotes the wall shear stress. The time is nor
malized by the wall time unit defined by
The temperature is normalized by the friction tempera
ture defined by
0.7 q (13)
PlplU (13)
where qw denotes the wall heat flux.
The bulk mean temperature
(pCp)
is used as the representative temperature of the fluid in
the present study. The heat transfer coefficient can be
defined as
ht = (15)
Ow (b'
where O, is the averaged temperature variation at the
wall.
The enhancement of heat transfer can be estimated by
the Nusselt number which is defined by
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
incompressible velocity field, was applied for the three
dimensional incompressible velocity field. The exten
sion was conducted by decomposing the 3D velocity
field into three twocomponent velocity fields.
In the continuum surface force (CSF) method, the in
terfacial tension force is calculated as a body force,
f, = aln6s
where a is the coefficient of the interfacial tension, K is
the curvature, n is the normal to the interface, and 8s is
a delta function concentrated on the interface. The in
terfacial tension force in Eq.(2) was calculated by using
Eq.(17) and the relation n8s VF, which holds in
the continuum limit. The curvature, K, was calculated
with a high degree of accuracy by using height func
tions (Lorstad & Fuchs (21" i",1,1 which effectively elim
inates spurious currents for a static drop (Francois et al.
(2006)).
In a general VOF method, two interfaces inside the
same grid cell cannot be distinguished. Thus, coales
cence occurs when two bubbles are very close to each
other. To avoid this type of coalescence, separate VOF
functions are assigned to the bubbles. A repulsive force
is applied when two bubbles come very close to each
other to avoid the overlap of the bubbles.
Suppose Bubble 1 and Bubble 2 approach each other
and their VOF functions are given by Fi and F2. We
make use of their smoothed functions F1 and F2 which
are obtained by the convolution using the kernel
1 6(r)2 +6()3
K(r, y)= 2(3
0
'i. i.
( < 1)
( < < 1)
(1 < ).
The repulsive force
Numerical Methods. The position of the interface
was determined by the volume fraction F of the bub
ble. F 1 represents a cell filled with the gas of
the bubble, while F 0 indicates that the cell is
filled with the continuousphase liquid. The cells of
0 < F < 1 include the interface. The time evolution
of F was estimated with a PLICVOF algorithm (Rud
man (1998); Scardovelli & Zaleski (2003)) In PLIC
VOF algorithm, we take account of the slope of the in
terface. Young and Parkar's method was used to esti
mate the normal vectors of the interface from the val
ues of F in adjacent cells (Parker & Youngs (1992)).
A massconserving scheme which also satisfies the con
sistency condition (0 < F < 1) was applied for the
advection of the VOF function. The EILE (Eulerian
Implicit and Lagrangian Explicit) scheme (Scardovelli
& Zaleski (2003); Aulisa et al. (2003)), which is an ad
vection scheme originally designed for twodimensional
f,1 CFIVF2
is exerted on Bubble 1, where C is a positive constant.
Similarly, the force
f,2 CF2VFI
is exerted on Bubble 2. The integral of f, + f,, over
the whole computational domain is exactly zero. There
fore, this repulsive force does not affect the total mo
mentum in the system. In the present study, we set
7 = 3Ax.
The momentum and energy equations were solved
by using the finite difference schemes of the colloca
tion type. All of the velocity components and the pres
sure were defined at the center of the cell. The time
integration was based on a fractionalstep method where
a pseudopressure was used to correct the velocity field
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Computational conditions.
Edtv6s number(Eo ,. .., /a) 1
Archimedes number(N pf.. / ) 900
Density ratio(pb/pl) 0.02
Viscosity ratio(Pb/Pl) 0.02
Void fraction 0.06
so that the continuity equation was satisfied. A balanced
force algorithm developed by Francois et al. (2006) was
used to suppress unphysical flows (or spurious currents)
resulting from the unbalance between the interfacial ten
sion and the pressure difference across the interface.
Poisson's equation for the pseudopressure was solved
by using a multigrid method.
The 3rdorder CIP scheme (Yabe et al. (2001)) was
applied in the finite differencing of the convection terms
of momentum and energy equations. A directional split
ting method was used in the multidimensional CIP
method. The secondorder central difference scheme
was applied for the finite differencing of the viscous
terms of the momentum equations and the diffusion term
of the energy equation.
The 2ndorder improved Euler method was used for
the timeintegration of the convective and viscous (or
diffusion) terms (Rudman (1998)). The velocity field at
t" + At/2 was used to advect the VOF function.
Validation of Numerical Scheme. We have com
pared the rise velocity of a bubble in otherwise quies
cent fluid with that in Bunner & Tryggvason (2002).
The bubble rise velocity can be described by the bub
ble Reynolds number, Red ubdo/vi, where ub and
do denote the rise velocity and the bubble diameter, re
spectively. Numerical parameters employed here are the
same as that in Bunner & Tryggvason (2002). The com
putational conditions are summarized in Table 1. There
cases with different grid resolutions were examined with
32 x 32 x 32, 64 x 64 x 64, and 128 x 128 x 128 grid
points. For the three cases, the bubble diameter equals
to 15.5Ax, 31.1Ax, and 62.2Ax, respectively. Here,
Ax denotes the grid spacing. Periodical boundary con
ditions are imposed in all directions. Domain size was
set to be 27 x 27 x 27r.
Fig.2 shows the time evolution of the bubble rise ve
locity for the three cases. The horizontal straight line
of Red 31.74 represents the terminal velocity in the
case of d 38.9Ax in Bunner & Tryggvason (2002).
The tendency of convergence is clearly seen as the grid
resolution is increased. This indicates the validity of the
present numerical method.
It is found by an oscillation test that more than 20
grid points per bubble diameter are needed to simulate
oscillation motions of bubbles. A static drop test was
20
10
0
o
32
 64
1283
Red=31.74
0 5 10
t(g/d)o.5
15 20
Fig. 2: Time change of rise velocities for a single bubble.
also conducted under the same condition as in Francois
et al. (2006) to find that the magnitude of the spurious
current is as small as theirs.
Results and Discussion
Computational conditions. Since many grid points
are needed to resolve the fluid motions around and inside
the bubbles, simulations are limited to low Reynolds
numbers with small computational domains. We utilized
socalled 'minimal turbulent channel' in the present
study. Computational conditions are summarized in Ta
ble 2.
The simulations were performed on the domain of
7 x 2 x 7/2, as in Lu & Tryggvason (2008), with
192 x 192 x 96 rectangular grid cells. In the simula
tion of Lu & Tryggvason (2008), the constant pressure
gradient to drive the flow was set so that the friction
Reynolds number uh/v was 127.2. Therefore, the size
of the computational domain was 400 x 254.4 x 200
in wall units, and the domain was sufficiently large to
sustain turbulence, as was shown in Jimenez & Moin
(1991). The resulting channel Reynolds number was
3786 in their simulation of the singlephase turbulent
flow. In the present study, the channel Reynolds num
ber was set at 3786, and the flow rate was kept constant.
The timeaveraged friction Reynolds number was 127.2
in the present simulation in the case without bubbles.
We set the density inside the bubbles to be onetenth
of the liquid (Pb/pl 0.1), and we set the viscosi
ties to be equal (Pb/l = 1) to reduce the computa
tional cost as in Lu & Tryggvason (2008). Air bub
bles with a diameter of 1 2mm in water were con
sidered in the present study. EOtvds and Morton num
bers of the bubbles are 0.4 (Eo = ..i/ lo 0.4) and
3.27 x 10 10 (Mo gi/pr3 = 3.27 x 10 10), re
spectively Archimedes number is given by the relation
Ar Eo3/2/Mo1/2 = p ..I / and is 14000 in the
present study. These parameters correspond to a 1.7mm
bubble in the fluid whose viscosity is 1.9 times higher
than the water at room temperature.
We have imposed a constant temperature gradient in
the vertical direction (x direction) and uniform heat flux
from both walls. The instantaneous wall heat flux qw
was given by
2 Jo
so that the energy of the system was kept.
Prandtl numbers of gas and liquid phases were set at
Prb = 1.0 and Prl = 2.0, respectively. We set the
specific heat inside the bubbles to be half of the liquid
(Cpb/Cpi 0.5). A simulation of lower liquid Prandtl
number of Prl 1.0 was also conducted for compari
son. The ratio of heat conductivity kb/kl 1.0 for the
case of Prj = 2.0.
Fully developed singlephase turbulent fields were
used as initial conditions for the velocity and tempera
ture fields. Twelve bubbles with a diameter of 0.4 were
introduced randomly into the turbulent flow in the chan
nel. Hereafter, all quantities are normalized by the fric
tion velocity and the friction time in the corresponding
singlephase flow if not otherwise specified.
SinglePhase Turbulent Flow. For comparison,
a long simulation without bubbles has been also con
ducted at the same channel Reybolds number of 3786
and Prandtl numbers of Prl =1, 2. The computational
conditions are summarized in Table 3. Statistical quan
tities are obtained by taking averages over the period of
t+ 4200.
Fig. 3 shows the mean velocity profile. The numerical
result shows that the mean velocity agrees well with the
theoretical relation, (u+) y+, in the viscous sublayer,
while it has higher values than the empirical equation,
(u+) = 2.5 In y+ + 5.5, in the log layer. Here, () repre
sents the spatial average in the x and z directions. These
higher values may be attributed to the small computa
tional domain in the present simulation. Fig. 4 shows
the profiles of turbulence intensities. Solid lines repre
sent the DNS results of Lu & Tryggvason (2008). The
result shows a good agreement with theirs.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 2: Computational condition for the multiphase
flow
Channel Reynolds number
Domain size
Number of grid points
Grid resolution
Density ratio
Viscosity ratio
Diameter of bubbles
Number of bubbles
Edtv6s number
Morton number
Archimedes number
Prandtl number(liquid)
Prandtl number(bubble)
Ratio of specific heat
Void fraction
Boundary condition
3786
7 x 2 x 7/2
192 x 192 x 96
Ax+ Az+ 2.08
Ay+ 0.47 ~ 2.08
0.1
1.0
0.4( 24.4Ax )
(50.9 wall units)
12
0.4
3.27 x 10 10
14000
2
1
0.5
0.04
x,zdirection:Periodical
ydirection:NonSlip
Table 3: Computational condition for the singlephase
flow
Number of grid points 96 x 192 x 48
Grid resolution Ax+ Az+ 4.16
Ay+ 0.47 ~ 2.08
Boundary condition x,zdirection:Periodical
ydirection:NonSlip
As is shown in Fig.5, the Nusselt number exhibits
an oscillatory behavior due to intermittent generation of
vorticies in a small computational domain. The time
averaged values for Pr 1 and Pr = 2 are 15.6 and
20.2, respectively.
Turbulent Bubbly Flow. The turbulent bubbly flow
and the temperature field reached a fully developed state
about t 1000 after the introduction of bubbles. The
simulation was further conducted for about t+ 400
to obtain the statistical quantities of the turbulent bubbly
flow.
Fig.6 shows the typical example of the bubbles and
the vortical structures visualized by the second invari
ant of velocity gradient tensor Q = 0.05. It is
clearly seen in Fig.6(a) that the bubbles tend to collect
on the nearwall regions of the channel. The bubbles are
slightly deformed form the spherical shape. It is also
seen that vortices are activated by the bubbles especially
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Fig. 3: The mean velocity profile.
Fig. 4: The rms velocity fluctuations.
in the nearwall regions.
Fig.7 shows the time evolution of the locations of bub
bles. The bubbles rise along the walls almost all the
time. As is shown in Fig.8, the void fraction has sharp
peaks near the walls as a result of bubble accumulation
toward the walls.
Fig.9 shows the mean velocity profiles. The solid line
represents the liquid velocity in the singlephase flow,
and the broken line and circles represent the liquid and
gas velocities in the multiphase flow, respectively. Note
that the gas velocity is not the rise velocity of the bub
30 .
25 
20
S15 
10
 Pr 2
5  Pr 1
0 I I I I ,
4000 5000 6000 7000 8000
t+
Fig. 5: Time evolution of the Nusselt number.
(a) Multiphase
'ts
(b) Single phase
Fig. 6: The bubbles and vortices for the two simulations.
ble centroids, but just the velocity of the gas inside the
bubbles. The gas velocity takes higher values than the
liquid velocity near the walls, while it takes lower val
ues around the regions of y 0.25 and 1.75. The
bubbles considered in the present study are quite large
(db 0.4h), and they are exposed to high shear near the
walls. The balance between this shear stress and the in
terfacial surface tension leads to the higher gas velocity
near the walls and the lower gas velocity on the central
side of the bubbles.
The mean bubble velocity tu (z 12.0) is slightly
lower than the surrounding liquid velocity in spite of the
upward buoyancy forces acting on the bubbles. The rise
velocity of the bubbles considered here is about u+ = 5
in otherwise quiescent fluid, indicating the importance
'* r
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Fig. 9: The mean velocity profile.
0.06
I)
S0.05
C
.0.04
S0.03
>
0.02
5 0.01
R
Fig. 7: The path of bubbles.
0.15
0.1
O
0
> 0.05
0
Fig. 8: The mean void fraction profile.
of the buoyancy forces. In fact, the effects of interfacial
surface tension of bubbles dominate those of the buoy
ancy forces in this simulation as will be shown below.
The liquid velocity increases in the central region of the
channel due to the blockage effect of the bubbles near
the walls. The wall shear stress is increased by a factor
of 1.57 due to the effects of the bubbles.
Fig.10 shows the profile of streamwise vorticity
squared (w2). The amplitude of the streamwise com
ponent of vorticity has larger values compared with the
corresponding singlephase flow. There are two small
peaks at y/h 0.07(1.93) and 0.33(1.70) for each side
of the channel. These peaks indicate that the bubbles
rising along the walls enhance the generation of quasi
streamwise vortices. This is verified by the visualization
of the vortical structures around the bubbles. As shown
in Fig. 11, relatively small trailing vortices are seen on
0 0.5 1 1.5 2
Fig. 10: The streamwise vorticity squared profiles.
the wall side of the bubbles (e.g. a vortex marked by A).
Most of these vortices are generated by the tilting of the
vortex ring on the upperwall side of the bubble due to
the mean shear. Relatively large vortices are seen on the
side of the channel center (e.g. a vortex marked by B).
Most of these vortices are generated by the bending of
the background meanshear vorticity due to the bubbles.
Figs.1214 show the profiles of the liquid velocity
fluctuations. The streamwise component decreases and
the wallnormal and spanwise components increase due
to the presence of bubbles in the nearwall regions. The
fluid motions normal to the interfaces of bubbles are di
rectly suppressed by the presence of bubbles, leading to
the reduction of the velocity fluctuations. On the other
hand, the bubbles indirectly affect the fluid motions. The
turbulence is activated by the bubbles, and the redistribu
tion mechanism of Reynolds stresses is enhanced. This
may be one of the reasons for the decrease in the stream
wise component and the increase in the wallnormal and
spanwise components near the walls.
The wallfriction drag is increased due to the injec
tion of the bubbles (Note that the volume flow rate is
kept constant in the computation). Now, we examine the
mechanism for the increase of the drag by considering
i I
4 single phase I
S   multiphase 
\ I 
'\ I 
\ I
/
I 

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Fig. 11: The bubbles and vortices near the wall.
3
S single phase
ti  multiphase t
2'
0 0.5 1 1.5 2
y
Fig. 12: The rms velocity fluctuation profiles in stream
wise direciton.
the balance of forces in the channel. Taking the aver
age over time and the x and z directions of Eq.(2), and
integrating the averaged equation with respect to y, we
obtain the relation (Note that h 1.)
T(y) ',..,.' (y)
+ g f(p ) (y') dy'+ ( fs) (y') dy'
S T(1 y). (22)
Here, T, = T (0). The first, second, third, and fourth
terms on the righthand side of Eq.(22) represent the vis
cous shear stress, the Reynolds shear stress, the buoy
ancy term, and the surfacetension term, respectively.
The profiles of these four terms are drawn in Fig. 15.
They are normalized by the wall shear stress in the case
of the singlephase flow. The sum of the four terms and
the straight line of T") (1 y) are also plotted in the
figure, where T,) denotes the timeaveraged wall shear
stress in the multiphase flow and is 1.57 in this simula
tion. They agree well with each other, which indicates
that the overall balance of forces is satisfied. The viscous
0 0.5 1 1.5 2
Fig. 13: The rms velocity fluctuation profiles in wall
normal direciton.
single phase
  multiphase /
0 0.5 1 1.5 2
y
Fig. 14: The rms velocity fluctuation profiles in spna
wise direciton.
shear stress is dominant in the nearwall regions as in the
case of singlephase flows. The surfacetension term has
large values in the regions of high void fractions (see
Fig.8). The bubbles are deformed by the mean shear in
the nearwall regions and a restoring force due to the in
terfacial tension is acting to the liquid fluid. Since this
term has large values near the walls, it makes a major
contribution to the increase in the friction drag. The
buoyancy term has relatively large values in the core re
gion of the channel as well as the Reynolds shear stress.
Since this term is small near the walls, its contribution
to the increase of the friction drag is also small.
The profiles of viscous shear stress are plotted in
Fig.16. As was shown before, the wall shear stress
in the turbulent bubbly flow is higher than that in the
corresponding singlephase flow. The relative magni
tude of the viscous shear stress to the wall shear stress
(T (y) /Tw) is smaller in the bubbly flow, which means
that the effective viscosity of the fluid is increased lo
cally near the walls by the effect of the bubbles.
In the present study, we set the ratio of viscosities
(Pb/ml) as 1.0, which is higher than that in the actual
airwater situation. We examined the contribution from
1.5 ' Reynolds stress
S viscous stress
 surface tension
S0.5  buoyancy
0.5
1.5 
0 0.5 1 1.5 2
V
1.5
S1
0.5
0
S0.5
1
1.5
Fig. 15: Budget for shear stress.
single phase
  multiphase
^^ ~ 1
0 0.5 1 1.5 2
y
Fig. 16: The mean viscous shear stress profile.
the gas phase to the total viscous shear stress to find it
very small in the whole region of the channel (figures
not shown). This indicates that the effect of the larger
viscosity ratio is not large.
Fig. 17 shows the profiles of the Reynolds shear stress.
It is clearly seen that the magnitude of Reynolds shear
stress is increased in the nearwall regions due to the
presence of bubbles. This is because the momentum ex
change is enhanced by the vortices generated around the
bubbles. This increase in the Reynolds shear stress near
the walls contributes to the increase in wall friction.
Enhancement of Heat Transfer due to the Bubbles.
Fig. 18 shows the time evolution of the Nusselt numbers
in the cases of Pr = 1 and 2. The timeaveraged Nus
selt number for Pr 1 and Pr = 2 are 18.8 and 26.6,
which are 1.21 and 1.32 times higher than the corre
sponding values in the singlephase flow, respectively.
The profiles of the mean liquid temperature,
({ea, W+), are drawn in Fig. 19, where (OWii) is
the averaged temperature variation at walls. The tem
perature is normalized by the friction temperature of the
singlephase flow. The temperature difference between
the channel center and the walls is smaller in the multi
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
15 single phase
SI   multiphase
0.5
1.5
0 0.5 1 1.5 2
y
Fig. 17: The mean Reynolds shear stress profile.
30
20 .. 
10 
 Pr= 2
 Pr 1
0 100 200 300 400
t tt
Fig. 18: Time evolution of the Nusselt number.
phase flow than in the singlephase flow, which indicates
the enhancement of heat transfer due to the bubbles.
Now, we examine the mechanism for the enhance
ment of heat transfer by considering the energy balance
in the channel. The averaging of the energy equation
Eq.(7) over time and the x and z directions, and the inte
25
30 20
15 " / \
V10
 single phase
5  multiphase
0 0.5 1 1.5 2
y
Fig. 19: The mean temperature profiles for Pr 2.
. 1 1
K
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1 L.H.S. ofEq.(23)
1."  turbulent
S0.5  Molecular
= 0.5
0     4
1 
0 0.5 1 1.5 2
y
Fig. 20: Budget for heat transfer for Pr = 2.
1 single pahse
  multipahse
0.5
$0
S0.5
1 
0 0.5 1 1.5 2
y
Fig. 21: The profiles of L.H.S. of Eq.(23) for Pr = 2.
gration of the averaged equation with respect to y yield
q G G (ipCpu) (y') dy'
k o) (y) + (pCpv) (y) (23)
The lefthand side of Eq.(23) represents the the total heat
flux in the wallnormal direction, and the first and the
second terms on the righthand side represent the molec
ular heat flux and the turbulent heat flux, respectively. In
the figures below, each term is normalized by the wall
heat flux of each run.
The profiles of these three terms are drawn in Fig.20.
The sum of the molecular and turbulent fluxes is also
plotted in the figure. The sum agrees well with the left
hand side of Eq.(23), which indicates that the overall
balance of heat transfer is satisfied.
Fig.21 shows the profiles of the lefthand side of
Eq.(23) for Pri = 2. Clear difference is not seen
between the singlephase and multiphase flows. The
change in the profile of heatcapacity flow rate due to
the bubbles has an insignificant effect on the enhance
ment of heat transfer.
The profiles of the molecular heat flux for Prq = 2 are
  mumltpnase
~0.5
0.5
1
0 0.5 1 1.5 2
Fig. 22 The mean molecular heat flux profiles for Pr
Fig. 22: The mean molecular heat flux profiles for Pr
1
& 0.5
0.5
1
Fig. 23: The mean wallnormal turbulent heat flux pro
files for Pr = 2.
shown in Fig.22. The molecular heat flux in the bubbly
flow is small near the walls compared with the case of
the singlephase flow. This indicates that the effective
heat conductivity is increased due to the bubbles.
In the present study, we set the ratio of heat conduc
tivities (kb/ki) at 1.0, which is higher than that in the
actual airwater situation. We examined the contribution
from the gas phase to the total molecular heat flux to find
it very small in the whole region of the channel (figures
not shown). This indicates that the effect of the larger
heat conductivity ratio is small.
Fig.23 shows the profiles of turbulent heat flux for
Prl = 2. The turbulent heat flux is increased by the
effects of bubbles in the regions near the walls. As is
shown below, this is caused by the vortices whose gen
eration is activated by the bubbles. This increase in the
turbulent heat flux in the nearwall regions is responsible
for the increase in the effective heat conductivity of the
fluid near the wall, which enhances the heat transfer.
In summary, the enhancement of heat transfer in the
bubbly flow is caused by the increase of the turbulent
heat flux near the walls in the present simulation.
Fig.24(a) shows the instantaneous temperature field
and the crosssections of vortices in a horizontal (y z)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
4000 5000 6000 7000 80
Fig. 25: Nu(Pr =
phase flow.
(a) multiphase
2)/Nu(Pr = 1) for the single
1 I    i   i  I  
1.4
1.3
0 100 200 300 400
t+
(b) single phase
Fig. 26: Nu(Pr = 2)/Nu(Pr
flow.
1) for the multiphase
Fig. 24: Colormap of temperature in y
Pr 2.
z planes for
plane for Prl = 2. Those for the corresponding single
phase flow are also drawn in Fig.24(b) for comparison.
The temperature field is normalized by the friction tem
perature in the case of the singlephase flow. e+ 
(K ,,l) is shown in the figures. Red and blue represent
the regions of + ( i) > 0 and 0+ (w+,) <
30, respectively. The crosssections of vortices are
represented by the black contour lines of the second in
variant of velocity gradient tensor of Q+ 0.05.
In the singlephase turbulent flow, the vortices are lo
cated only in the lowtemperature regions away from the
walls. In the bubbly turbulent flow, on the other hand,
they are also located in the nearwall regions where the
temperature gradient is relatively high. It is obvious that
the vortices near the walls play a major role in the heat
transfer enhancement.
It is interesting to know how the heattransfer en
hancement depends on the liquid Prandtl number. As
is shown in Figs.2526, the heattransfer enhancement is
more noticeable at higher Prandtl numbers. Since the
thermal boundary layer is thinner at a higher Prandtl
number, the vortices near the walls (generated by the
bubbles) more effectively enhance the heat transfer at
a higher Pr.
Conclusions
Direct numerical simulations have been conducted for
turbulent bubbly upflow between two parallel heating
walls in order to clarify its heat transfer characteristics.
We have obtained the following results.
The bubbles accumulate toward the walls in the tur
bulent channel upflow.
The turbulence production is activated by the bub
bles in the nearwall regions.
The wall friction is increased by the injection of
bubbles. This is mainly caused by the interfacial
surface tension resulting from the deformation of
the bubbles due to the high shear near the walls.
The heat transfer is enhanced by the injection of
bubbles. This is because the turbulent heat flux is
activated by the generation of the vortices due to
the bubbles.
The heattransfer enhancement is more noticeable
at higher liquid Prandtl numbers.
References
Serizawa, A., Kataoka, I. and Michiyoshi, I., Turbu
lence Structure of AirWater Bubbly Flow I. Measuring
C
.o
1.3
I)
z
Techniques, Int. J. Multiphase Flow, Vol. 2, pp.221246,
1975.
Liu T.J., Bubble size and entrance length effects on void
development in a vertical channel, Int. J. Multiphase
Flow, Vol.19, No.l, pp.99113, 1993.
Lu J. and Tryggvason G., Effect of bubble deformability
in turbulent bubbly upflow in a vertical channel, Phys.
Fluids, Vol.20, 040701, 2008.
Tamari M., Niishikawa K., The Stirring Effect of Bub
bles Upon the Heat Transfer to Liquids,Trans, Japan
Soc. Mech. Engrs., pp.3144, 1975
Tokuhiro A.T., Lykoudis PS, Natural Convection Heat
Transfer from A Vertical PlateI. Enhancement with
Gas Injection, Int. J. Heat Mass Transfer, Vol.37, No.6,
pp.9971003, 1994.
K.Kosuge,K Uchida,A Kitagawa,Y.Hagiwara, Enhance
ment of Natural Convection Heat Transfer along a verti
cal plate by microbubbles, 4th JapaneseEuropean two
phase flow group meeting, pp.18, 2006.
Rudman M., A volumetracking method for incompress
ible multifluid flows with large density variations, Int. J.
Numer. Meth. Fluids, Vol.28, pp.357378, 1998.
Scardovelli R. and Zaleski S., Interface reconstruction
with leastsquare fit and split EulerianLagrangian ad
vection, Int. J. Numer. Meth. Fluids, Vol.41, pp.251274,
2003.
Parker B.J. and Youngs D.L., Two and three dimensional
Eulerian simulations of fluid flow with material inter
face, Technical Report,01/92, Atomic Weapons Estab
lishment, Aldermaston,Berkshire, February 1992.
Aulisa E., Manservisi S., Scardovelli R., Zaleski S., A
geometrical areapreserving VolumeofFluid advection
method, J. Compt. Phys., Vol.192, pp.355364, 2003. J.
Compt. Phys., Vol.213, pp.141173, 2006.
Lorstad D. and Fuchs L., Highorder surface tension
VOFmodel for 3D bubble flows with high density ra
tio, J. Compt. Phys., Vol.200, pp.153176, 2004.
Francois M.M., et al., A balancedforce algorithm for
continuous and sharp interfacial surface tension models
within a volume tracking framework, J. Compt. Phys.
Vol. 213, pp.141173, 2006.
Yabe T, Xiao F. and Utsumi T, The Constrained In
terpolation Profile Method for Multiphase Analysis, J.
Compt. Phys., Vol.169, pp.556593, 2001.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Bunner B. and Tryggvason G., Dynamics of homoge
neous bubbly flows Part. 1 Rise velocity and microstruc
ture of the bubbles, J. Comput. Mech., Vol.466, pp.17
52, 2002.
Jimenez J. and Moin P., The minimal flow unit in near
wall turbulence, J. Fluid Mech., Vol.225, pp.213240,
1991.
