Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 10.6.1 - Characteristics of heat transfer in turbulent bubbly upflow in a vertical channel
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00265
 Material Information
Title: 10.6.1 - Characteristics of heat transfer in turbulent bubbly upflow in a vertical channel Multiphase Flows with Heat and Mass Transfer
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Tanaka, M.
Matsui, N.
Miyajima, Y.
Hagiwara, Y.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: DNS
bubbles
heat transfer
vortices
 Notes
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 PLIC-VOF 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 1-2mm in water were considered. It is found that the bubbles tend to collect on the near-wall regions of the channel and that vortices are activated by bubbles especially in the near-wall 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 wall-normal 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 heat-transfer enhancement. It is also found that the heat-transfer enhancement is more noticeable at higher Prandtl numbers.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00265
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1061-Tanaka-ICMF2010.pdf

Full Text



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 606-8585, 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 PLIC-VOF 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 1-2mm in water were considered. It is found that
the bubbles tend to collect on the near-wall regions of the channel and that vortices are activated by bubbles especially
in the near-wall 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 wall-normal
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 heat-transfer enhancement. It is also found that the heat-transfer
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 bubble-induced
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 heat-transfer 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 5-6mm, while it has
a peak in the core of the channel for the bubbles larger
than 5-6mm. 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 near-wall regions, while
strongly deformable bubbles tend to be expelled from
the near-wall 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
heat-transfer 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 heat-transfer 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, wall-normal 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 non-uniform 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 Navier-Stokes (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 phase-averaged 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 two-component 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 continuous-phase liquid. The cells of
0 < F < 1 include the interface. The time evolution
of F was estimated with a PLIC-VOF 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 mass-conserving scheme which also satisfies the con-
sistency condition (0 < F < 1) was applied for the
advection of the VOF function. The EI-LE (Eulerian
Implicit and Lagrangian Explicit) scheme (Scardovelli
& Zaleski (2003); Aulisa et al. (2003)), which is an ad-
vection scheme originally designed for two-dimensional


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 fractional-step method where
a pseudo-pressure 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 pseudo-pressure was solved
by using a multigrid method.
The 3rd-order 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 multi-dimensional CIP
method. The second-order 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 2nd-order improved Euler method was used for
the time-integration 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
so-called '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 single-phase turbulent
flow. In the present study, the channel Reynolds num-
ber was set at 3786, and the flow rate was kept constant.
The time-averaged friction Reynolds number was 127.2
in the present simulation in the case without bubbles.











We set the density inside the bubbles to be one-tenth
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 single-phase 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
single-phase flow if not otherwise specified.

Single-Phase 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,z-direction:Periodical
y-direction:Non-Slip


Table 3: Computational condition for the single-phase
flow
Number of grid points 96 x 192 x 48
Grid resolution Ax+ Az+ 4.16
Ay+ 0.47 ~ 2.08
Boundary condition x,z-direction:Periodical
y-direction:Non-Slip



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 near-wall 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 near-wall 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 single-phase 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 single-phase 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 upper-wall 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 mean-shear vorticity due to the bubbles.
Figs.12-14 show the profiles of the liquid velocity
fluctuations. The streamwise component decreases and
the wall-normal and spanwise components increase due
to the presence of bubbles in the near-wall 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 wall-normal and
spanwise components near the walls.
The wall-friction 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 right-hand side of Eq.(22) represent the vis-
cous shear stress, the Reynolds shear stress, the buoy-
ancy term, and the surface-tension 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 single-phase 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 time-averaged 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 near-wall regions as in the
case of single-phase flows. The surface-tension term has
large values in the regions of high void fractions (see
Fig.8). The bubbles are deformed by the mean shear in
the near-wall 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 single-phase 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
air-water 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
S-0.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 near-wall 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 time-averaged 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 single-phase 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
single-phase 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 single-phase 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

S-0.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 left-hand side of Eq.(23) represents the the total heat
flux in the wall-normal direction, and the first and the
second terms on the right-hand 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 left-hand side of
Eq.(23) for Pri = 2. Clear difference is not seen
between the single-phase and multiphase flows. The
change in the profile of heat-capacity 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 wall-normal 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 single-phase 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 air-water 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 near-wall 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 cross-sections 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 single-phase flow. e+ -
(K ,,l) is shown in the figures. Red and blue represent
the regions of + ( i) > 0 and 0+ (w+,) <
-30, respectively. The cross-sections of vortices are
represented by the black contour lines of the second in-
variant of velocity gradient tensor of Q+ -0.05.
In the single-phase turbulent flow, the vortices are lo-
cated only in the low-temperature regions away from the
walls. In the bubbly turbulent flow, on the other hand,
they are also located in the near-wall 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 heat-transfer en-
hancement depends on the liquid Prandtl number. As
is shown in Figs.25-26, the heat-transfer 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 near-wall 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 heat-transfer enhancement is more noticeable
at higher liquid Prandtl numbers.

References

Serizawa, A., Kataoka, I. and Michiyoshi, I., Turbu-
lence Structure of Air-Water Bubbly Flow- I. Measuring


C
.o

1.3
I)
z











Techniques, Int. J. Multiphase Flow, Vol. 2, pp.221-246,
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.99-113, 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.31-44, 1975

Tokuhiro A.T., Lykoudis PS, Natural Convection Heat
Transfer from A Vertical Plate-I. Enhancement with
Gas Injection, Int. J. Heat Mass Transfer, Vol.37, No.6,
pp.997-1003, 1994.

K.Kosuge,K Uchida,A Kitagawa,Y.Hagiwara, Enhance-
ment of Natural Convection Heat Transfer along a verti-
cal plate by micro-bubbles, 4th Japanese-European two-
phase flow group meeting, pp.1-8, 2006.

Rudman M., A volume-tracking method for incompress-
ible multifluid flows with large density variations, Int. J.
Numer. Meth. Fluids, Vol.28, pp.357-378, 1998.

Scardovelli R. and Zaleski S., Interface reconstruction
with least-square fit and split Eulerian-Lagrangian ad-
vection, Int. J. Numer. Meth. Fluids, Vol.41, pp.251-274,
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 area-preserving Volume-of-Fluid advection
method, J. Compt. Phys., Vol.192, pp.355-364, 2003. J.
Compt. Phys., Vol.213, pp.141-173, 2006.

Lorstad D. and Fuchs L., High-order surface tension
VOF-model for 3D bubble flows with high density ra-
tio, J. Compt. Phys., Vol.200, pp.153-176, 2004.

Francois M.M., et al., A balanced-force algorithm for
continuous and sharp interfacial surface tension models
within a volume tracking framework, J. Compt. Phys.
Vol. 213, pp.141-173, 2006.

Yabe T, Xiao F. and Utsumi T, The Constrained In-
terpolation Profile Method for Multiphase Analysis, J.
Compt. Phys., Vol.169, pp.556-593, 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.213-240,
1991.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs