Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Twoway numerical study of bubbly upflows in a vertical channel using the
EulerLagrange model
M.J. Pang1, J.J. Wei1*, B. Yu2
1 Xi'an Jiaotong University, State Key Laboratory of Multiphase Flow in Power Engineering
Xi'an 710049, China
2China University of Petroleum, Beijing Key Laboratory of Urban Oil and Gas Distribution Technology
Beijinl02249, China
*Email: jjwei@mail.xjtu.edu.cn
Keywords: Bublly flow; EulerLagrange model; DNS; Turbulence modulation; Phase distribution
Abstract
Bubbly flows exist extensively in industrial processes, so it is very meaningful to study flow characteristics of them. This
paper introduces in detail a numerical method of the EulerLagrange twoway model for the airwater bubbly flows. The flow
field is simulated using direct numerical simulations (DNS) in Euler frame of reference, while the bubble dynamics are fully
analyzed by integration of Lagrangian equations of motion that takes into account interphase interaction forces including drag
force, shear lift force, wall lift force, pressure gradient force, virtual mass force, gravity force, buoyant force, and inertia force
in Lagrange fame of reference. The coupling between the phases is considered by regarding the interphase interaction forces
as momentum source terms of the continuous phase. Bubbles distribution and turbulent structures of the liquid phase under
the influence of bubbles are comprehensively analyzed. The results show that an overwhelming majority of bubbles cluster
near channel walls, and turbulent structures of the liquid phase are modified to some extent by addition of bubbles, namely,
the average streamwise velocity becomes enhanced, whereas the wallnormal and spanwise turbulent intensity and Reynolds
stress are dramatically reduced. It is also found that effects of bubbles on the liquid turbulent structure are similar to that of
the surfactant, but the mechanisms are much different in essence.
Introduction
Bubbly flows are widely used in a variety of industrial
processes such as smelting of metals, twophase heat
exchangers, aeration and stirring of reactors, fermentation
reactions, floatation of minerals, and bubble column
reactors, air conditioning engineering, wastewater treatment,
cavitating flows, transportation lines and so on. Bubbles
also play an important role in a number of natural
phenomena including the propagation of sound in the ocean,
the exchange of gases and heat between the oceans and the
atmosphere, and explosive volcanic eruption (Bunner and
Tryggvason, 2003). In particular, bubble columns are
frequently applied to the chemical industry due to their
excellent mass and heat transfer characteristics and simple
construction. Additionally, microbubbles as an additive, in
place of other additives such as polymer and surfactant, are
also employed to reduce the frictional drag force and save
energy due to their environmentally friendly characteristic.
One of numerous applications is a vertical upward bubbly
flow in a pipe (Guet and Ooms, 2006). This situation can,
for instance, be found in airlift reactors for enhancing
mixing, or for providing oxygen to microorganisms. It is
also encountered in mining technologies and wastewater
treatment. In any case, it is necessary for any application of
bubbly flows in practice that the fundamental properties of
the hydrodynamics should be fully understood. It is
reported that understanding the mechanisms that phases
distribute and the gasliquid phases interact with each is
crucial for a number of engineering processes, including
mixing, drag and heat transfer (TranCong et al., 2008). To
better understand the dynamics of bubbly flows, a number
of investigations have been done in the past. Although
many investigations on phase distribution and interaction in
bubbly flows were carried out in the past, perfect
mechanism models used to explain and predict physical
phenomena are still absent.
While experimental studies allow observations of large
systems under realistic conditions, it is generally difficult to
independently achieve full control of various effects and to
obtain very detailed data. Therefore, we investigate some
characteristics of the bubbly flow by using the numerical
method. The objective of the present work is to have a deep
Paper No
investigation on turbulence characteristics of bubbly flows
in an upward channel by using EulerLagrange twoway
model. Although bubbly flows can be simulated in
multiways (such as direct numerical simulation,
EulerLagrange trajectory model simulation and EulerEuler
twofluid model simulation), with different degrees of
approximation and assumption, different models have
different characteristics and application ranges. Comparing
with other models, the EulerLagrange model, in which the
motion of each bubble is individually tracked, provides the
most detailed insight into singlebubble dynamics as well as
into bubblebubble and bubblefluid interactions. It may be
used to evaluate the influence of basic physical and
geometric parameters like inertia, viscosity, surface tension,
bubble size and gas volume fraction on the evolution of
bubble swarms, the interaction of bubbles of different sizes,
and the induced flow structure of the liquid. It is well
known that interphase forces between bubbles and liquid
are numerous, but it is generally agreed that the drag force
is largely predominant over other contributions (Chen, 2004;
Cedric et al., 2009). Zhang et al. (2006) reported that the
virtual mass force has a small influence on the investigated
bubbly flow characteristic, however, the lift force strongly
influences the bubble plume dynamics and consequently
determines the shape of the vertical velocity profile. Tabib
et al. (2008) also considered that no significant contribution
of the added mass force on the simulation of local
hydrodynamics of bubble column is seen, but the drag and
lift forces have important influence on simulation.
Therefore, we take into account interphase forces including
drag force, shear lift force, wall lift force, pressure gradient
force, virtual mass force, gravity force, buoyant force and
inertia force, and other forces are neglected in the
simulation. In the present phase of the effort, attention is
focused exclusively on the void fraction and effects of
bubbles on the liquidphase turbulence. In order to reduce
errors caused by employing turbulence models, the
continuous phase is simulated by using direct numerical
simulations in Euler frame of reference. Trajectories of
bubbles' motion are obtained by solving the Newton
equations of motion in Lagrange frame of reference.
Because the void fraction is less than 10 interaction forces
between bubbles are neglected (Elghobashi, 1994).
Computational conditions
Computational domain and boundary conditions
We examine the bubbly flow in a vertical upward channel
between two parallel walls, where the streamwise,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
wallnormal, and spanwise directions are denoted by x, y,
and z, respectively. The flow geometric size and the
coordinate system are shown in Fig. 1. Gravity acts in the
negative x direction. Periodic boundary conditions are
imposed in both the streamwise and spanwise directions and
noslip boundary conditions are enforced at the right and
left walls. The computational domain size is 10hx2hx5h in
the streamwise, wallnormal, and spanwise directions,
respectively, where h is half the channel width. Our
calculations are made for a friction Reynolds number Re, =
150, which is based on the wall friction velocity u, and half
the channel width h. Uniform grids are used in the
streamwise and spanwise directions. Nonuniform grids are
generated in the normal direction by using the hyperbolic
tangent stretching function (It's defined by equation (1),
where N andj are the grid number and number index of grid
point in y direction) with denser mesh near the wall to
resolve smallscale eddies. The total number of grid points
of computational domain is 274625. Initial velocities of
bubbles are set to be zero. When steady turbulence is
statistically reached, a set of bubbles is injected into the
flow field.
h FIn39(
y h = tanh 1+2 (1)
0.95 L 2 Nj (
Figure 1: Flow geometry and coordinate system.
Physical and geometric parameters of liquid and
bubbles
The liquid phase is considered to be incompressible,
isothermal and has constant properties. Its thermophysical
properties data are chose as those of water at room
temperature. Liquids density is pf=1000 kg/m3, and liquids
kinematic viscosity is uf=106m2/s. Bubbles are considered
as spherical shape with the same diameter of 1.65 wall units,
and their density is b =1.3 kg/m3. The bubbles' deformation
is neglected, the amount of bubbles is M=19200, and the
average void fraction is ao=1.34x104. In the computational
process, Froude number is Fr = 0.0169. The related
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
parameters are listed in Table 1.
Table 1. Computational parameter of bubbles and liquid.
Computational parameter Specific values
Domain size (x, y, and z direction) 10hx2hx5h (=1500x300x750 wall units)
Resolution of mesh grid 64 x64 x64
Friction Reynolds number 150
Density (pb/pf) 1.3/1000.0
Diameter of bubbles (db) 0.01 h (1.65 wall units)
Number of bubbles (M) 19200
Kinematic visocity (vb/vf) 14.8 x 106/1.0 x 106
Average void fraction of bubbles (ao) 1.34 x104
Froude number (Fr) 0.0169
Governing equations Fn = F +FG +FD+Fm +Fs +Fw. (4)
Governing equations of the liquid phase The following specially presents the empirical correlations
According to the above assumptions, the nondimensional of the drag coefficient (CD) and the lift coefficient (CLF).
governing equations for the liquid phase can be written as The drag coefficient is calculated by the following
follows. empirical correlations (Lain et al. 2002).:
Dimensionless continuity equation: 161.5
RReb 1.5
au~
fi = 0.
Dimensionless momentum equation:
ap 1q _
P +11 2U f;. (3)
ax R Re,
14.9
S= Re078 1.5 Reb < 80) (5)
48 2.21 1 4756
(1 )+1.86x10 Re 80 Reb <1500
Reb Reb
2.61 1500 Reb
For the spherical bubble, the lift coefficient(Legendre and
Magnaudet 1998), CLF, is calculated by
Where all parameters related to length scale are normalized
by half the channel width (h), parameters having something
to do with velocity scale are normalized by the friction
velocity (u,), and variables related to time scale are
normalized by h/ut. The nondimensional values of velocity,
length, pressure, time and interfacial force are defined as
follows:
UZ
u h =h
w Z
U T'',
P p+ *
Pfu
t = t F =
h/u, u2 /h1
Note that the pressure field is decomposed into a fluctuating
pressure p+* and average value x* which drives the mean
flow. Here f, is the average interphase force, 8 is
Kronecker Delta.
Motion equations of bubbles
In view of assumptions for the flow, forces, such as FD the
drag force, Fn the inertial force, Fs the lift force, Fw the
wall lift force, Fm the virtual mass force, Fp the pressure
gradient force, FG the gravity force (including the buoyant
force), are applied to calculate trajectories of bubbles. A
force balance equation for the bubble, i.e., the constitutive
equation can be written as
CLF = K ,I"oRe (Reb,Srb)F + (gh Re (eb}2
where
6 Srb )0 5 2.255
CowRe(Reb,Srb) = (Reb Srb)05 2255
72 + 0.22 5
LhghRe R = 2 +16ReI eb exp(2.92d;221)
2 1+29/RebJ
Where Srb, Reb are nondimensional shear rate and bubble
Reynolds number. They can be respectively defined as
Srb = I4db/2luf b andReb = uIf bl db/ f
here (= Sb/Re.
Each force in equation (4) is replaced with specific
calculation expressions, and then it is normalized. Thus,
equation (4) is rewritten as
Paper No
tauf + 'u
7+u7 ,
Paper No
du+ Du, a1 1 3 )
P D 1+  t 1
Pf dr* Dt j 1 i 4 ab I f / I
du+ dut
C d + CLF Ukklm (U'
au
 U) f+
ax*
S+u + 2 mao 0.0064 0.0161
2 u u" max 0, + 
1 yfb )
where Re, is friction Reynolds number, Re, = u h/vf and Fr
is Froude number, Fr = uI/g .
Numerical methods
Direct numerical simulations of the liquid phase
(Kajishima et al., 2001; Yu, 2004)
The basic equations for fluid flow are discretized by using
the staggered grid in the Cartesian coordinate system.
Velocities of three directions are stored at the face center of
grid, and the pressure is stored at the center of grid. The
grid spacing in each direction is denoted by Ax, Ay, Az. The
grid spacing in y direction increases with the increase of
distance from walls. The continuity equation (2) is
discretized as
Sf = 0. (8)
Without considering the pressure term temporarily, the
equation of momentum is discretized as
T = gfdsc = H( fui, +  JJUjf, + (9)
t Re,
Where 6, represents the central finite difference operator
and the overline denotes the mean value between
neighboring terms. The second order central difference
operator is defined as
x Ifjk f+l,j,k 2f,'Jk + f1,j,k (10)
Ax
The timeadvancing scheme adopts the fractional step
method based on the AdamsBashforth method of the
second order accuracy in time. The intermediate velocity
components are expressed as
+\ + >(n) 3 + (n) 1 u+ (n 1)
uf, = f, +At* f dzsc 7Uf disc (11)
The index (n) symbolizes the current timestep, (nl) stands
for the previous timestep and () signifies the intermediate
step. Taking into account of the pressure term, the velocity
components for the next timestep (n+1) are expressed as
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
in a Poisson equation for the pressure
L L .." 1 Lu .... +'v ". a(,) ". (13)
a&x2 y*2 *2 At &* ay* z* I
After the pressure is solved from Eq. (13), then the velocity
field for the next timestep can be corrected by using Eq.
(12). When the effect of interphase force on the velocity
field of the liquid phase is considered, the velocity
components for the next timestep are rewritten as
+n+1 +, n+l
u l u At* Ap + At*f,*
Ax,
The Poisson equation is solved by using the multigrid
method. The MAC method is employed to couple velocity
and pressure fields.
Bubble motion equations
Liquid quantities at the bubble location
Because bubbles locations may be different from those used
to store velocity and pressure of the liquid phase, liquid
quantities at the bubbles locations should be firstly obtained
by the interpolation method before interphase forces related
to relative velocity between bubbles and liquids are
calculated. First of all, velocity values stored at the adjacent
interface center of cells are interpolated at vertices of cells
using a face interpolation. Liquid quantities at the bubbles
sites are calculated making use of liquid values of vertices
of cells by using a volume interpolation. The streamwise
velocity for a vertex of cell is given by
U i,,k = a1 U,,k + a2 'l,J+l,k +
a3 U,k,+l + a4* U,J+l,k+l
The steamwise velocity at a bubble location is given by
u",J ,k = l .U ,J,k +b2 1,Jl1,k + b3 .Ii.jk1 +
b4 .,J1,k + b5 ,u 1,,k +b6 .U1,jl,k + (16)
b7 lu 1,j,k1 +b8 nu1,1,k 
where al, a2, a3, a4, bl, b2, b3, b4, b5, b6, b7, b8 are
interpolation coefficients. Their values depend on the grid
system. Subscripts, i,j, k, are denoted as number indexes of
grid points. Figs. 2 and 3 provide illustration on two
interpolation methods. The interpolation method in the
other two directions is same to that of the streamwise
direction.
(n"+l) + ), ap n+l)
Uf, L Lxj
Substituting the divergence of Eq. (12) into Eq. (8) results
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
pfV 1"h
Figure 2: Face interpolation of liquid velocity.
Figure 3: Volume interpolation of liquid velocity
Bubble velocity and displacement
Bubble velocity equation is given by
U ( = + a At, (17)
where a is the nondimensional acceleration of bubbles,
a = du+ /dt
Bubble displacement is calculated through the integration
equation of velocity by using the second order
CrankNicholson method. Bubble displacement equation is
given by
*(n+1) *(n) At + (n+1) u+ (18)
Xb = X +b +U b
where X is the nondimensional displacement of bubbles.
Local and average interphase force
The resultant force acting on bubbles is translated into
liquids of the cell containing bubbles. This force is
calculated by the following expression.
/= yAm (a*).
f' pfAcAyAz j 'l
where a, is the acceleration of the ith bubble, m is the
amount of bubbles in a certain cell of which the number
index for the grid point is (i, j, k). The average interphase
force acting on liquids is calculated by using the following
equation.
where V is the volume of the computational domain.
Results and discussion
The lateral void fraction profile along the normalwall
distance is represented in Fig.4, and the local void fraction
(a) is normalized by the local maximum void fraction (amax).
Fig. 4 shows that the void fraction approximately displays
the symmetrical distribution along xz plane at y=0, the
peak values appear near the walls, and the place
corresponding to the maximum of the void fraction is about
one bubble diameter from the pipe wall under the influence
of the wall lift force as pointed out by many authors (Antal
et al., 1991; Lopez de Bertodano et al., 1994; Troshko and
Hassan, 2001). Presently, it is generally accepted that, for
upward bubbly flows, small bubbles tend to migrate toward
the pipe wall forming a wallpeak void fraction profile
under the action of the lift force, whereas large bubbles tend
to move toward the core of the pipe forming a corepeak
void fraction distribution because the direction of the lift
force is changed by the wake flow induced by large bubbles
(Kariyasaki, 1987; Ervin and Tryggvason, 1997; Lu and
Tryggvason, 2008). In the present investigation, we roughly
estimate that the bubble Reynolds number is in a range,
0
approximately less than 50, vortex shedding does not occur
(Koynov and Khinast, 2004) and only wake flow induced
by bubbles affects the liquid turbulence. Okawa et al. (2002)
considered that the influence of wake flow on the bubble
motion can also be neglected if the void fraction is
comparably low. Therefore, we can conclude that addition
of bubbles does not induce enough strong turbulence to
change the direction of the lift force in the present study.
And thus, a great number of bubbles migrate towards the
walls just under the influence of the shear lift force as
reported by many researches.
0os8
021
ne i i 
Figure 4: Void fraction profile along channel height.
Paper No
Paper No
Fig. 5 shows mean streamwise velocity profiles for the
bubble, singleliquid and liquid phases. Here, y+ is a
nondimensional wallnormal distance, Y+=yulvf. Fig. 5
shows that the result of the singleliquid phase is in good
agreement with the wellknown lawofthewall profile of
Newtonian fluids: the linear law u =y+ in viscous
sublayer (/<5) and logarithmic law u = 2.51ny+ + 5.0 in
the inner layer (y>30). One can clearly see from Fig. 5 that
bubbles flow faster than liquid in the whole channel height
under the influence of the buoyant force. In the present
Reynolds number range, the liquidphase velocity is slightly
high than that of the single liquid phase in the center region
of the channel, whereas the mean streamwise velocity
profiles for the liquid and single liquid phase almost overlap
near the wall. Similar results are obtained by Nagaya et al.
(2002) and Akoi et al. (2006). It should be noticed that this
phenomenon is not an essential characteristics of bubbly
flows through analysis for existing literatures. The influence
of bubbles on the liquidphase turbulence is associated with
many factors such as the global void fraction, the bubble
size and shape, the gas velocity, the Reynolds number, the
flow direction with respect to the gravity and so on.
Studying conditions of different literature are almost
different, so obtained results are also almost not same. And
thus, it also leads to different explaining mechanism to
different phenomena and even different mechanism to the
same phenomenon.
Y+
Figure 5: Mean streamwise velocity profiles of liquid and
bubbles.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
In order to ascertain the reason of the liquidphase velocity
change and fully understand the influence of bubbles on the
liquid turbulence, some turbulence statistics are further
investigated. Fig. 6 shows the turbulent intensity profiles of
liquid with and without bubbles. Injection of bubbles
reduces turbulence intensities in the wallnormal and
spanwise directions along the whole height of the channel
as shown in Fig. 6. It shows nonuniform suppression to the
turbulence intensity in the streamwise directions.
Suppression of the turbulence intensities in the wallnormal
and spanwise directions is very similar to the results of
Kitagawa et al. (2005) and Lu and Tryggvason (2008). Just
as the influence of bubbles on the mean steamwise velocity,
mechanisms on modulation of the liquid turbulence
intensities by bubbles, which are presented by different
researchers, are very different due to different study
conditions or other reasons. In the existing literature,
explanations on the effect of bubbles on the liquid
turbulence intensities are few and not unified. Kato et al.
(1999) concluded that the influence of bubbles on the liquid
turbulence intensities is just like the influence of solid
particles on the gas turbulence. Large bubbles (their
diameter is approximately larger than 1.5mm) enhance the
liquidphase turbulence, whereas small bubbles reduce the
liquid turbulence. This point is also accepted by Kim et al.
(2005). They pointed out that when there is no change in
volume of the dispersed flows which include bubbly flow,
droplet flow, and particleladen flow, the production and
dissipation of the turbulent continuous phase can be
understood by similar mechanisms. Kawamura and Kodama
(2002) pointed out that bubbles can influence the velocity
fluctuation in two ways. One is the socalled
"pseudoturbulence", that changes fluctuation of velocity
signals due to passage of bubbles, and the second is liquid
phase turbulence modulated by bubbles. However, the
increase or reduction of the turbulence intensities can not be
judged by this information. Mazzitelli et al. (2003)
considered that microbubbles accumulating in downflow
regions locally transfer momentum upwards. Consequently,
addition of microbubbles reduces the intensity of the
vertical velocity fluctuations and the turbulent kinetic
energy. Lu and Tryggvason (2008) believed that bubble
deformation has a great influence on the liquid turbulence
intensities. Nearly spherical bubbles reduce the liquid
turbulence intensities, whereas deformable bubbles enhance
the liquid turbulence intensities. In short, the mechanism of
modulation of the liquid turbulence by bubbles can not be
Figure 6: Turbulence intensity profiles.
Paper No
fully understood yet and the existing mechanisms can not
be widely accepted so far.
In the present study, the bubbles diameter is enough small,
and thus bubble deformation can be neglected due to small
Weber number (We =pU2Rb/a<
2004). Additionally, vortexes shedding can not also occur
due to the small bubble Reynolds number. Therefore, the
influence of bubble size and shape on the liquid turbulence
can be ignored. As is well known, the liquid fluctuations
correspond to a variety of different scaled turbulent vortices
in the turbulent flow. If the turbulent vortexes are
suppressed, the turbulence intensities will be reduced.
Compared with the singleliquid phase, a momentum source
of the interphase force appears in the momentum equation
of bubbly flows. It is easy to conclude that the interphase
forces exactly cause the decrease of the turbulence
intensities. In the present investigation, among all
interphase forces taken into account in computation, only
the lift force has the direct influence on the turbulence
vortexes. According to the definition of the lift force,
concurrence of gasliquid relative velocities and local
vorticities of fluids results in the lift force as visually shown
in Fig. 7. When the lift force, which is induced by the liquid
vortexes and the relative velocity between bubbles and
liquid, pushes bubbles toward the wall, an equal and
opposite force in the direction will exert on liquid vortexes.
This reaction force will decay the liquid vortexes, and even
move the whole liquid vortex away from the wall as
reported by Ferrante and Elghobashi (2' i4). And thus,
liquid vortexes suppressed by the lift force correspond to
the decrease of the turbulence intensities in the normal
direction. Experimental results of Kamp et al. (1995) and
numerical results of Chahed et al. (2002) also discovered
that bubbles have nearly no influence on the liquid
turbulence when the lift force was neglected in their
investigations. This fact further shows that the lift force has
a great effect on the turbulent intensities in the normal
direction. Besides the influence of the lift force on the
liquid vortexes, friction between phases and bubble
fluctuations are also able to consume part of liquidphase
energy, which indirectly reduces the liquidphase velocity
fluctuations.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Bubble b) b I
I\/ //L
s^ ^l I i \ \ ^\^ B,,,^t.
jF, y/ / trFj
Figure 7: The lift force acting on a bubble rising in a
spanwise vortex pushes it to the wall.
The total shear stress can be used to confirm if the
calculation has reached a statistically steady state. For the
fully developed turbulent flow, when a statistically steady
state is reached, the balance equation of shear stresses for
the singleliquid flow can be expressed as
 1 au+
Total = 1 = u'+ v''+ +
Re, Re, y* (21)
The righthand terms of equation (21) are respectively the
nondimensional Reynolds stress and viscous stress. Fig.8
shows the shear stress profiles for both the liquid phase and
the singleliquid phase. It can be seen that the addition of
bubbles reduces the liquidphase Reynolds stress, but it has
little impact on the viscous shear stress. For the dilute
bubbly flow, the total shear stress also maintains the linear
profile along the height of the channel. When the Reynolds
stress is reduced but the total shear stress and viscous shear
stress keep unchanged, there must be an additional force to
keep the stress balance. The force is no other than the
interfacial force, and therefore the shear stress balance
equation can be rewritten as
1 1 u y+
Total =1 R =" V + Re a
Re Re, dy
(22)
The last term on the right hand of equation (22) is the
interfacial stress term, which is equal to zero for the
singleliquid phase. The reason why the viscous shear stress
almost keeps unchanged is illustrated as follows. It is well
known that there are two factors influencing the viscous
shear stress: one is the liquid effective viscosity, and the
other is the velocity gradient. For a dilute suspension
system, the effective viscosity can be computed by the
following equation (Pan et al.,1999)
/t =/~lqud (1+2.5a +O(a0)). (23)
In the present investigation, the average void fraction is
equal to 1.34x104, so addition of bubbles has a very small
effect on the liquid viscosity. Additionally, compared with
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the singleliquid phase, the velocity gradient of the liquid
phase does not be greatly changed by bubbles. And thus, the
viscous shear stress almost keeps unchanged as shown in
Fig. 8.
Figure 8: Budgets of shear stress.
Comparing with the viscous shear stress, the Reynolds
shear stress is relatively complicated, and it is relatively
difficult to understand the reduction of it. In order to
provide the detailed information of the contribution to the
total turbulence production from various events occurring in
turbulent flow, and to provide information on the influence
of bubbles on such contributions, the quadrant analysis on
fluctuating velocity field (u, v) is performed. According to
Lu and Willmarth (1973) and Li et al. (2008), Fig. 9
schematically shows the categorization of fluid motions
according to the signs of u and v. The first quadrant,
corresponding to u>0 and v>0, denotes outward motion of
highmomentum fluid (Q1 motion for short); the second
quadrant, u<0 and v>0, shows ejections of lowmoment
fluid from the wall (Q2 motion for short). The third
quadrant, u<0 and v<0, depicts wallward motion of
lowmomentum fluid (Q3 motion for short); and the fourth
quadrant, u>0 and v<0, indicates sweeps of highmoment
fluid toward the wall (Q4 motion for short). The Q1Q4
depicted in Figs. 9 and 10 represent the Q1Q4 motions.
ui
ejection of lowmomentum
fluid from the wall
r>ylaO.QlImotion
outward interaction
WeII
u0f 0;:Q4motion
wallward interaction sweep of highmomentum
fluid toward the wall
Figure 9: Schematic diagram of the categorization of
turbulent fluid motions (quadrant motions)
Figure 10: Fractional contribution of Reynolds shear stress
from different quadrant motions.
Fig.10 shows the fractional contribution of the Reynolds
shear stress from different quadrant motions. It can be seen
that the addition of bubbles has a greater influence on the
second (Q2) and fourth quadrant (Q4) motions than on the
first (Q1) and third quadrant (Q3) motions. Bubbles
significantly suppress the ejection and sweep motions
occurring near the wall, which results in the reduced
Reynolds shear stress. The effect of bubbles on Reynolds
stress in Fig.10 is qualitatively similar to that of surfactant
additive as shown in the literature (Li et al., 2008).
Although the effects of the bubble on the liquidphase
turbulence are similar to surfactants, both of them are
different in essence. The influence of bubbles on the
liquidphase turbulence is an irreversible process with the
friction between the gasliquid phases translating a part of
energy into heat. However, the process happening to
surfactant is reversible, in which the surfactant solution
strongly resists deformation in shear flow but at the same
time also stores energy by deformation of the micellar
network, reducing the turbulence intensity. The deformed
elastic fluid elements spring back and releasing the energy
when external conditions are altered.
The budget of turbulent kinetic energy is very helpful to
further understand the liquid turbulence modulation induced
by bubbles. For a fully developed turbulent bubbly flow, the
budget terms of Reynolds stress (tu' u ) can be expressed
as
D =P,+T, +D +J E +E,, (24)
where the terms on the righthand side of the above
equation are identified as follows
p = k '+u' ' production rate;
k kx x
Paper No
Paper No
Ti = u, turbulent diffusion rate;
a 2
D, = +2u'+ u' vel.press.gradient term;
9x, ,.
nI = '+ a + u'+ ap viscous diffusion rate;
S x' xax
E = 2 dissipation rate;
dXk d k
E, = a(f, (u',) interfacial force term.
In equation (24), there is an extra term called the interfacial
force term El compared to the singleliquid phase. Here, a,
, are the local void fraction, the local interfacial
average force and the average velocity fluctuations,
respectively.
The budget terms of turbulent kinetic energy of three
velocity components for both the liquid phase and the
singleliquid phase are plotted in Fig. 11 as a function of the
dimensionless wall distance y It can be seen from Fig. 11
that there is a similar distribution trend for each term of
them, the addition of bubbles has little influence on each
term of budget of Reynolds stress in the streamwise
direction, whereas the addition of bubbles has a greater
influence on each term of Reynolds stress in the
wallnormal and spanwise directions. Compared with the
single liquid phase, both dissipation term and velocity
pressure gradient term in the wallnormal and spanwise
directions are reduced due to the addition of bubbles. It is
well known that no direct production term is available for
Reynolds stress in these two directions, as the velocity
pressure gradient correlations are usually considered to be
pseudoproductions. The decrease of velocity pressure
gradient term in these two directions means the reduction of
energy distributed to these two directions. This result
corresponds to the fact that the wallnormal and spanwise
turbulence intensities for the bubbly flow are reduced,
compared to the singlephase flow. In order to investigate
the effect of bubbles on the energy distribution from the
streamwise direction to the wallnormal and spanwise
directions in detail, Fig. 12 alone shows the profile of the
streamwise velocity pressure gradient terms for the
singleliquid phase and liquid phase. It can be clearly seen
that addition of bubbles reduces the velocity pressure
gradient terms in the streamwise directions. This implies
that the redistribution of turbulent energy from the
streamwise velocity components to the wallnormal and
spanwise velocity components is strongly suppressed due to
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
addition of bubbles. This part of spared energy is either
converted into useful work to enhance the mean streamwise
velocity or dissipated by friction between phases induced by
interphase forces. How to judge that the spared energy is
dissipated or converted into useful work is still scarce until
now, and it needs to be deeply explored in the next future.
 Singleliquid phase
0 4  Liqud phase
Production
0 3
0 2 / turbulent diffusion
S 0 1 v el pres gradient
0
'IF
(a) Budget of Reynolds stress u' u'
002
Singleliquid phase
 Liquidphase
Sel press gradient
0 01 
iscous diffusion
o
Turbulent diffusion
0 01 sspat
Dissipation
(b) Budget of Reynolds stress v' v'
(c) Budget of Reynolds stress w' w'
Figure 11: Budgets of Reynolds stress.
o 
t0 02
 Singlehliquid phase
 Llquid phase
50 100 150
Figure 12: Profiles ofvel. press.gradient of u' u' .
Paper No
Conclusions
A comprehensive numerical study of the bubbly flow in an
upward channel was carried out by using the
EulerLagrange twoway method. The following
conclusions can be obtained:
(1) A great number of bubbles tend to accumulate near the
walls and form a peakwall value distribution under the
action of the lift force, and the void faction shows a
symmetric distribution along y=0 plane.
(2)The addition of bubbles enhances the liquidphase
velocity at the core of the channel and reduces the
turbulence intensities in the wallnormal and spanwise
directions compared to those of the singleliquid phase.
Turbulence vortexes decayed by the lift force directly leads
to the decrease of the liquidphase turbulent intensities in
the wallnormal and spanwise directions along the channel
height.
(3)The addition of bubbles suppresses the strength of
bursting events, which results in reduction of Reynolds
shear stress. It also reduces velocity pressure gradient terms
in the streamwise directions, which means that the
redistribution of turbulent energy from the streamwise
velocity components to wallnormal and spanwise velocity
components is suppressed due to the addition of bubbles.
This further explains the reason why the wallnormal and
spanwise turbulence intensities of the liquid phase are
reduced compared to the single phase flow.
Acknowledgement
We gratefully acknowledge the financial support from the
NSFC Fund (No. 10602043, 50821064, 50876114), the
Specialized Research Fund for the Doctoral Program of
Higher Education of China (No. 20090201110002) and
Jiangsu Provincial Natural Science Foundation
(BK2009145).
References
Akoi K., Hishida K. & Kodama Y. Measurements of near
wall turbulent structure in a microbubble flow using a
highly magnifying telecentric PIV/PTV system, 13th Int
Symp on Applications of Laser Techniques to Fluid
Mechanics Lisbon, Portugal, 2629 June, 2006.
Antal S.P, Lahey Jr. R.T. & Flaherty J.E. Analysis of phase
distribution in fully developed laminar bubbly twophase
flow. Int. J. Multiphase Flow, Vol.7, 635652(1991)
Bunner B. & Tryggvason G Effect of bubble deformation
on the properties of bubbly flows. J. Fluid Mech., Vol.495,
77118(2003)
Chahed J., Colin C. & Masbernat L. Turbulence and Phase
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Distribution in Bubbly Pipe Flow Under Microgravity
Condition. J. Fluids Eng., Vol.124, 951956(2002)
Chen P, Sanyal J. & Dudukovic M.P CFD modeling of
bubble columns flows: implementation of population
balance. Chem. Eng. Sci., Vol.59, 52015207 l1 i4)
Cedric L.B., Faical L., Dromard N., et al. CFD simulation
of bubble column flows: Investigations on turbulence
models in RANS approach. Chem. Eng. Sci. Vol.64,
43994413(2009)
Elghobashi S. On Predicting ParticleLaden Turbulent
Flows. Appl. Sci. Res., Vol.52,309329(1994)
Ferrante A. & Elghobashi S. On the physical mechanisms of
drag reduction in a spatially developing turbulent boundary
layer laden with microbubbles. J. Fluid Mech.,
Vol.503,345355k 1" i4)
Ervin E.A. &Tryggvason G. The rise of bubbles in a vertical
shear flow. J. Fluids Eng., Vol.119, 443449 (1997)
Fujiwara A., Minato D. & Hishida K. Effect of bubble
diameter on modification of turbulence in an upward pipe
flow. Int. J. Heat Fluid Fl., Vol.25,4814XssX~,l I4
Guet S. & Ooms G. Fluid Mechanical Aspects of the
GasLift Technique. Annu. Rev. Fluid Mech., Vol.38,
22524'21"11,)
Kajishima T., Takiguchi S. & Hamasaki H. Turbulence
structure of particleladen flow in a vertical plane channel
due to vortex shedding. JSME International Journal Series
B, Vol.44,526535(2001)
Kamp A., Colin C. &Fabre J. The Local structure of a
turbulent bubbly pipe flow under different gravity condition.
Proceedings of the second international conference on
multiphase flow. Koto, A. Serizawa, T. Fukano, J. Bataille,
eds., 1995
Kariyasaki A. Behaviour of a single gas bubble in a liquid
flow with linear velocity profile, In: Proceedings of the
1987 ASMEJSME Thermal Engineering Joint Conference,
ASME, NewYork, NY, vol. 5, 261267,1987
Kato H., Iwashina T., Miyanaga M., et al. Effect of
microbubbles on the structure of turbulence on a turbulent
boundary layer. Journal of Marine Science and Technology,
Vol.4,155162(1999)
Kawamura T. & Kodama Y. Numerical simulation method
to resolve interactions between bubbles and turbulence. Int.
J. Heat Fluid Fl., Vol.23,6276' svil "'2)
Kim S., Bock Lee K. & Gu Lee C. Theoretical approach on
the turbulence intensity of the carrier fluid in dilute
twophase flows. Int. Commun. Heat Mass, Vol.32,
435444(2005)
Paper No
Kitagawa A., Hishida A.& Kodama Y. Flow structure of
microbubbleladen turbulent channel flow measured by PIV
combined with the shadow image technique. Exp. Fluids,
Vol.38,466475(2005)
Koynov A. & Khinast J.G., Effects of hydrodynamics and
Lagrangian transport on chemically reacting bubble flows.
Chem. Eng. Sci., Vol.59,3907'12 1'"4')
Lain S., Broder D., Sommerfeld M., et al. Modelling
hydrodynamics and turbulence in a bubble column using the
EulerLagrange procedure. Int. J. Multiphase Flow,
X6128,13x I14111i 1" 1i
Legendre D. & Magnaudet J. The lift force on a spherical
bubble in a viscous linear shear flow. J. Fluid Mech.,
Vol.369,81126(1998)
Li F.C., Kawaguchi Y. & Hishida K. Structural analysis of
turbulent transport in a heated dragreducing channel flow
with surfactant additives. Int. J. Heat Mass Tran.,
Vol.48,965973(2005)
Lopez de Bertodano M., Lahey Jr. R.T. &Jones O.C. Phase
distribution in bubbly twophase flow in vertical ducts, Int.
J. Multiphase Flow. Vol.20,805818(1994)
Lu J.C. & Tryggvason G., Effect of bubble deformability in
turbulent bubbly upflow in a vertical channel. Phys. Fluids,
Vol.20,040701(2008)
Lu S.S. & Willmarth W.W. Measurements of the structure
of the Reynolds stress in a turbulent boundary layer. J. Fluid
Mech., Vol.60,481(1973)
Mazzitelli I.M., Lohse D. & Toschi F. The effect of
microbubble on developed turbulence, Phys. Fluids,
Vol.15,685697(2003)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Nagaya S., Hishida K., Kakugawa A., et al. PIV/LIF
Measurement of Wall Turbulence Modification by
Microbubbles. Proceedings of the 3th International
Symposium on Smart Control of Turbulence, Tokyo,
Japan,6978, 2002
Okawa T., Kataoka I. & Mori M. Numerical simulation of
lateral phase distribution in turbulent upward bubbly
twophase flows. Nucl. Eng. Des., Vol.213,183197 (2002)
Pan Y, Dudukovic M.P & Chang M. Dynamic simulation
of bubbly flow in bubble columns. Chem. Eng. Sci.,
Vol.54,24812489(1999)
Tabib M.V., Roy S.A. & Joshi J.B. CFD simulation of
bubble columnAn analysis of interphase forces and
turbulence models. Chem. Eng. J. Vol.139,589614 (2008)
TranCong S., Marie J.L. & Perkins R.J. Bubble migration
in a turbulent boundary layer. Int. J. Multiphase Flow,
Vol.34,786807(2008)
Troshko A.A. & Hassan Y.A. A twoequation turbulence
model of turbulent bubbly flow. Int. J. Multiphase Flow,
Vol.27,19652000(2001)
Xu J., Maxey M. & Karniadakis G. Numerical simulation of
turbulent drag reduction using microbubbles. J. Fluid
Mech., Vol.468,271281(2002)
Yu B. & Kawaguchi Y. Direct numerical simulation of
viscoelastic dragreducing flow: a faithful finite difference
method. J. NonNewtonian Fluid, Vol.116, 4314,k1, 1r4)
Zhang D., Deen N.G. & Kuipers J.A.M. Numerical
simulation of the dynamic flowbehavior in a bubble column:
Astudy of closures for turbulence and interface forces.
Chem. Eng. Sci., Vol.61,75937608(2006)
