7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Modelling of Turbulence in Bubbly Flows
Eckhard Krepper, Martin Schmidtke, Dirk Lucas
Forschungszentrum DresdenRossendorf, POB 510119, D01314 Dresden,German
E.Krepperiafzd.de M. Schmidtkeiafzd.de D.Lucasiifzd.de
Keywords: CFD, turbulence, bubbly flow, model validation
Abstract
In the Euler/Eulerian approach simulating bubbly flow, the influence of the bubbles on the turbulence of the liquid has to be
modelled. Vice versa, the structures of the turbulent liquid flow influence the gas void fraction distribution, which is expressed as
a turbulent dispersion force. Reliable models for turbulence are an urgent precondition for the improvement of models describing
bubble coalescence and bubble breakup in any population balance model. In the present work the different approaches
simulating the influence of bubbles on the turbulence are revised and compared to measurements using ANSYSCFX. Models
for the turbulent dispersion force are validated.
Introduction
The concept proposed by Sato
In CFD modelling the concept of time averaging the Navier
Stokes equations proposed by Reynolds has proven
successful for the description of turbulent flows for many
industrial applications. Questions arise, when the influence
of bubbles on turbulence has to be considered. The correct
simulation of this influence is quite important since in
bubble size population balance models the liquid turbulence
influences the bubble fragmentation. Vice versa the liquid
turbulence flow structures influence the gas distribution.
The latter phenomenon is described as a turbulent dispersion
force.
In the present paper the common concept for modelling the
influence of bubbles on liquid turbulence quantities
proposed by Sato (1981) is described and analysed. Based
on this analysis a proposal is derived, describing this
influence by sources in the turbulence equations.
Simulations were done using ANSYSCFX and compared
with measured gas volume fraction distributions obtained
from the MTLoop test of FZD (Lucas et al., 2005) and
turbulence parameter obtained from measurements by
Shawkat et al. (2007, 2008).
In the ks turbulence model for a single phase flow the
turbulent viscosity prT is defined as a relation between
turbulent kinetic energy k and turbulent eddy dissipation e
bv:
k2
pT = C, pL
with 4ci as the liquid density and Cpu = 0.09.
Sato (1981) proposed for the case of two phase flow, to
consider the influence of bubbles by increasing the turbulent
viscosity by
P~s = CsPLn G UG UL
where G iS the gas volume fraction and UG TOSpective UL
the gas and liquid velocity. Cs = 0.6. The total viscosity is
the calculated from a laminar part p L, the turbulent viscosity
prT and the bubble induced contribution proposed by Sato pus.
Lu T Su u u
Experiments performed by Shawkat
This model approach has been used for long time as a CFD
standard and fields good agreement for the radial gas
fraction profiles between measurements and CFD
calculations considering the turbulent dispersion force,
which depends on the turbulent viscosity.
Simulations were done for Shawkats measurements (Figs.
14). A comparison for the turbulent kinetic energy shows a
reasonable agreement for the single phase case (Fig. 1, blue
asterisks and blue solid line) but a strong underestimation of
the turbulent kinetic energy in the pipe centre.
Unfortunately only few experiments on turbulence in two
phase flow in the interesting region of parameters were
published. Shawkat et al. (2007, 2008) have performed
measurements in a air/water flow of a vertical upward pipe
with diameter of 0.2 m and a height of 8 m. Besides radial
profiles of gas volume fractions velocities and velocity
fluctuations were measured. For the comparison with the
model here the axial and radial velocity fluctuations were
averaged and transformed in a profile of turbulent kinetic
energy.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The Sato model considers the increase of the turbulent
viscosity in bubbly flow, but the influence of the bubbles on
other turbulence parameters like the turbulent kinetic energy
and its dissipation rate are not modified in a physical based
way.
Consideration of turbulence modulation by
sources in the k and s equations
Bubble induced sources for turbulence kinetic energy
In the ks turbulence model for a liquid continuous phase the
presence of a dispersed gaseous phase can be considered in
the equation for the turbulence kinetic energy k:
~ (aILPLkL) V.OL L LkL)
(i' X L VkL L LS f (4)
and for the eddy dissipation E:
(aILPL L)V OLPL'L~L)
"(. L L L .SlePL 2eL L)+S (5)
SLk in equation (4) describes the bubble induced source in
the kequation. Assuming all frictional work of a rising
bubble is converted into turbulent kinetic energy tlus source
can be calculated by
o looo

,,.,,, ~ L=0.45 m/s
~ x X. X X JG=0.015 m/s
to.oloo G X0J.1 M /
S0.0001
0.00 0.02 0.04 0.06 0.08 0.10
[m]
Figure 1: Measured and calculated (SatoModel) profiles
for turbulent kinetic energy, measurements marked by
asterisks
"0.olooo
oE oo
oooo
JL=0.45 m/s
 JG=0
JG=0.015 m/s    
JG=0.1 m/s
0.00001
0.00 0.02 0.04 0.06 0.08 0.10
[m]
Figure 2: Calculated profiles (SatoModel) for the eddy
dissipation
JL=0.415 m/s
JG=P.015 m/s
S0.2 G~JG
0.0 '
0.00 0.02 0.04 0.06 0.08 0.10
[m]
Figure 3: Calculated profiles (SatoModel) for turbulent
eddy viscosity, the dotted lines mark the Sato contribution
CLs (see eq. 2)
0.60 
Where j is the drag force, 000 the gas volume fraction and
UG TOSpective UL the gas respective the liquid velocity.
For the Source in the Epsilon equation (5) is obtained:
SL" = Ce3 Sf
Bubble induced sources for turbulence eddy dissipation
Unfortunately there is no theoretical justification in the
literature for the corresponding source SLB (eq. 7). Most
approaches consider SLB proportional to St The relaxation
time z in most cases is calculated on a dimensional
background. Different concepts have been proposed.
Yao and Morel (2I r 14)assumed an additional dependency on
the bubble size da, because the wakes behind the bubbles are
assumed to be originally of the same diameter. The larger
these wakes, the larger will be the time for their energy to
cascade towards the smallest turbulence scales before to be
dissipated. So, if bubbles are quite large, the turbulence
generated is their wakes will have longer lifetime:
db 8)
Troshko and Hassan (2001) propose the consideration of the
 L=0.45 m/s
 JG=0.1 m/s
JG=0.015 m/s
 JG=0 m/s
0.40
S0.20
o.oo
0.000 0.020 0.040 0.060 0.080 0.100
[m]
Figure 4: Calculated profiles (SatoModel) for the liquid
velocity
S k =  aG G L~
2.5
2
1.5 *
1 *
0.5
0
0 0.05 0.1 0.15 0.2 0.25 0.3
co
Figure 7: Dependency of C,3 on the gas volume fraction
I~
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The FZD proposal for the Epsilon source
The concept
In the past lots of validation work was performed validating
the non drag forces (see e.g. Lucas et al 2007). Besides the
lift force and the wall force and the turbulent dispersion
force were considered. For the latter Burns et al. (21***4)
derived a Favre (mass weighted) averaged formulation:
FDISPlll TD4drcL LgradaG (10)
The turbulent dispersion force has a diffusive effect on the
gas distribution in the flow. Eq. 10 shows that the turbulent
dispersion force depends not only on the void fraction
gradient but also on the fluid viscosity Compp. eq. 3). For
these validations the bubble induced turbulence modulation
proposed by Sato (see eq. 2) was applied and good
agreement to measurements were found.
The idea aroused, to focus the consideration on the case of
homogeneous bubbly flow. For this case all radial
derivations and the turbulent viscosity (see eq. 1) are zero.
Then the sources in the turbulence equations Srk and SLB for
the case of homogeneous bubbly flows are adjusted in a way
that the turbulent viscosity obtained from eq. (1) equal to the
turbulent viscosity proposed by Sato (eq. 2) pt = Ps 
The assumption of a homogeneous bubbly flows reduces the
k equation (4) to
coefficients of virtual mass CIM and the drag force CD:
2Ca db
3CD UG UL
In the practical application of these approaches, the proposal
of Troshko and Hassan (eq. 9) causes some stability
problems because of the velocity difference in the
denominator. Both approaches need an additional
adjustment of the factor C,3 (see eq. 7).
JL=0~ ~.4 m/,J = 0.015 m/s
Morel, C,=0.1
 Morel, C,=0.45
 Morel, C,=1.0
 Sato
xxx exp
Ny 0.1000
E
S0.0100
0.0010
a
S0.0001
0.00 0.02
0.04 0.06 0.08 0.10
Figure 5: Calculated profiles for turbulent kinetic energy
(Morel model) and the Sato model Compp. Fig 1)
The application of the Morel approach (eq. 8) shows a much
better agreement to the Shawkat experiments than the Sato
approach. Fig. 5 shows the profiles of turbulent kinetic
energy and Fig. 6 of eddy viscosity for different values of
Cd. In Fig. 5 and 6 for comparison the corresponding values
according to Sato are added as red line. For C,3=1.0 a good
agreement of Morels approach for the turbulent kinetic
energy to Shawkat's measurements is found (see Fig. 5).
Furthermore this value shows an agreement to the turbulent
viscosity according to Sato's proposal (see Fig. 6).
0 = p, s, + Sfk
respective the E equation (5) to
6
S= , L 2PL L + SL
In combination with the Morel approach this results in an
adjustment for Cd as:
JL=0.45 m/s JG=0.015 m/s
Morel, C,=0.1
Morel, C,=0.45
Morel, C,=1.0
 Sato
IC C, 0.74
sE3 3 3
0.0 t
0.00
0.02 0.04 0.06 0.08 0.10
[m]
Figure 6: Calculated profiles for eddy viscosity for
different values of Cs, and the Sato model Compp. Fig 3)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Validation of models for turbulent dispersion force
Previous validation work concerning the non drag forces
compared measured and calculated gas volume fraction
profiles for fully developed flow situations. Therefore the
validation was directed only on the relative proportion of the
non drag forces. In the present work experimental data
gained at FZ DresdenRossendorf were used to compare
absolute values for the turbulent dispersion force by
investigating the gas distribution development near the gas
injection region (see Lucas et al. 2005).
Experiments
In the experiments air is inj ected in a vertical water flow (see
Fig. 11). Wire mesh sensor measurements of the cross
sectional distribution of the gas volume fraction were
performed in different tests at different distances from the
air injection. Assuming the same boundary conditions
during all tests of a series, the flow development can b
observed (see Lucas et al. 2005). The gas was injected by
nozzles distributed uniformly over the cross section of the
tube (see Fig. 12).
Fig. 7 shows the dependency of C, on the gas volume
fraction
Validation by Shawkat's experiments
Different tests were calculated applying this approach and
compared to Shawkat's measurements. The results are
shown in Figures 8 to 10.
JL=0.45 m/s
JG=0
0.1000 7;   =0.015 m/s
JG=0.1 m/s
0.0010
0.0010~1; x,
0.00 0.02 0.04 0.06 0.08 0.10
r
U)
N
E
x
a,
w
o
a,
Y
a,
o
a,
[m]
Figure 8: CFD simulation of the profiles for turbulent
kinetic energy, dotted lines denote the contribution of the
Sato model Compp. Fig 1)
0.3
E L J= 0.45 m/s
_IJG=0.015 m/s
0.0
0.00 0.02 0.04 0.06 0.08 0.10
[m]
Figure 9: Profiles for turbulent viscosity, dotted lines
denote the contribution of pus (see eq. 2)
~Water Flow
Figure 11: Arrangement of air injection and wire mesh
sensor during the MTLoop tests
JL=0.45 m/s
SJG=0
JG=0.015 m/s
in1.0000
~0.1000
~j0.0100
S0.0010
2 0.0001
,0.0000
0.0
)0 0.02 0.04 0.06 0.08 0.10
[m]
Figure 10: CFD simulation of the profiles for turbulent
eddy dissipation, dotted lines denote the contribution of the
Sato model Compp. Fig 2)
Figures 8 to 10 shows only a small illustrative part of the
performed validation calculations. All simulation show
similar tendencies.
Figure 12: Distribution of the 19 nozzles over the cross
section (compare Fig. 11 air injection)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
seen. This effect is less pronounced for a 50 % reduced
turbulent dispersion force at Cm=0.5 shown in Fig. 14 a)
than for Cm=1.0 shown in Figl4b. Also with reduced
turbulent dispersion force of Fig. 14 a) at z=0.23 m the trace
of the nozzle locations disappear completely.
oso
0.15
0.10
0.05
0.00
10.20
0.15
0.10
0.05
0.00
z = 0.08 m
z = 0.03 m
z = 0.08 m
z = 0.13 m
(0.20
0.1s
0.00
10.20
0.1s
0.10
o.os
U0.00
Figure 13: Measured gas void fraction distribution at
different distances z from the air injection for the test 074:
JL=1.067 m/s, k6=0.0368 m/s
In the measurements close to the injection plane a gas void
fraction maximum above each gas injection nozzle position
can be found (see Fig. 13, z=0.03 m). These maxima are
smoothened with larger distance from the injection plane.
Lucas et al. (2005) presents the corresponding gas volume
fraction distributions for further tests and further distances
from the gas injection. The smoothing process of the gas
fraction distribution with increasing distance from the
injection can be described by the diffusive effect of the
turbulent dispersion force (see eq. 10).
A turbulent dispersion coefficient less than 1 might be
reasonable since only a part of all eddies, namely eddies
larger than the bubble size contribute to the smoothing effect
of the turbulent dispersion force.
Validation by CFD calculations
CFD calculations were performed using the code
ANSYS/CFX. A 60 degree sector of the tube was simulated
with two side walls defined as symmetry planes. The liquid
was injected from the bottom via an inlet condition. The gas
was injected 0.5 m above the liquid inlet by point sources.
For the comparison tests with low gas fraction were selected
so that an Euler/Eulerian approach with a single
monodispersed bubble size could be applied. The drag was
simulated according to Grace. Besides the turbulent
dispersion force according to Burns (2l1~14 the lift force
according to Tomiyama (1998) and the wall force according
to Antal (1991) were considered. The bubble induced
influence on the parameters of liquid turbulence was
considered by sources in the k and a equations according to
the proposal described in the last section.
Fig. 14 shows the calculated gas volume fraction
distribution at different distances from the gas injection for
the test case 074 at a liquid superficial velocity of
JL=1.067 m/s. Gas was injected at a superficial velocity of
&6=0.0368 m/s. The smoothing effect on the gas volume
fraction distribution with increasing distance can be clearly
z = 0.13 m
z = 0.23 m
a) Cro=0.5
0o.20
D.15
0.os
o.oo
01.20
0.15
0.10
o.os
0.00
z = 0.03 m
z = 0.08 m
(0.20
0.1s
one
o.os
10.20
0.1s
0.00
z = 0.13 m
z = 0.23 m
b) CTD=1.0
Figure 14: Calculated gas volume fraction distributions at
different distances z from the gas injection.
Test 074: JL=1.067 m/s, k6=0.0368 m/s
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
072: z = 0.03 m
x x xExperiment
_ CTD =
 CTD = 0.5
074: z = 0.03 m
x x xExperiment
_, CTD
 CTD = 0.5
0.15
Z0.10
0.05
0.00
0.15
10.10
0.05
0.00 
0.000
0 000 0.005 0.010 0.015 0.020 0.025
0.005 0.010 0.015 0.020 0.025
[m]
074: z = 0.08 m
x x wExperiment
_ CTD1.0
CTD=0.5
072: z = 0.08 m
x x xExperiment
CTD=1.0 _
CTD=0.5
 
10.10
10.10
0 0 5 
0.00
0.000 0.005 0.010 0.015 0.020 0.025
[m]
0.20
074: z = 0.13 m
x x xExperiment
0.15 CTD=1.0
CTD=0.5
10.10
0.00 x
0.000 0.005 0.010 0.015 0.020 0.025
[m]
0.20
074: z = 0.23 m
x x xExperiment
0.15 __ __ CTD=0.5
CTD1.0
0.000 0.005 0.010 0.015 0.020 0.025
[m]
0 .2 0    
072: z = 0.13 m
x xx Experiment
0.15 _ CTD=1.0 
CTD=0.5
10.10 
0.00
0.000 0.005 0.010 0.015 0.020 0.025
[m]
0 .2 0     
072: z = 0.23 m
x xx Experiment
0.15 __ _ TD0. 5 
CTD=.0
0.05 x x x x x x
10.10
Txxx
x x x
r uXXXXXXXx
0.05
0.00
0.000 0.005 0.010 0.015 0.020 0.025 0.000 0.005 0.010 0.015 0.020 0.025
[m] [m]
Figure 15: Measured and calculated radial gas volume fraction profiles at different distances from the gas injection for the test
cases 072 (J, = 0.405 m/s, JG = 0.0368 m/s) and 074 (J = 1.067 m/s, JG = 0.0368 m/s)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Bumns, A.; Frank, T.; Hamill, I. and Shi, J.M. (21***4) The
Favre Averaged Drag Model for Turbulent Dispersion in
Eulerian MultiPhase Flows, 5th Intemnational Conference
on Multiphase Flow, ICMF'04, Paper No. 392
Grace, J.R., Wairegi, T. and Niguyen, T.H. (1976). Shapes
and velocities of simple drops and bubbles moving freely
through immiscible liquids, Trans. Inst. Chem. Eng., vol. 54,
pp. 167.
Lucas, D.; Krepper, E.; Prasser, H.M. (2005) Development
of cocurrent airwater flow in a vertical pipe, Int. J.
Multiphase Flow 31, 13041328
Lucas, D.; Krepper, E.; Prasser, H.M. (2007) Use of models
for lift, wall and turbulent dispersion forces acting on
bubbles for polydisperse flows, Chem. Eng. Sci., vol. 62,
pp. 4146 4157
Politano, M.; Carrica, P. & Converti, J. (2003) A model for
turbulent polydisperse twophase flow in vertical channels
International Joumnal of Multiphase Flow 29, 1153
Sato, Y. and Sekoguchi, K. (1975). Liquid velocity
distribution in two phase bubble flow, Int. J. of Multiphase
Flow, vol. 2, pp. 7995.
Sato, Y., Sadatomi, M., Sekoguchi, K. (1981) Momentum
and heat transfer in twophase bubble flowI, Int. J. of
Multiphase Flow, vol. 7, pp. 167177.
Shawkat, M.; Ching, C.Y. & Shoukri, M. (2007) On the
liquid turbulence energy spectra in twophase bubbly flow
in a large diameter vertical pipe, Int. J. of Multiphase Flow,
33, 300
Shawkat, M.; Ching, C. & Shoukri, M. (2008) Bubble and
liquid turbulence characteristics of bubbly flow in a large
diameter vertical pipe, Int. J. of Multiphase Flow, 34, 767
Tomiyama, A. (1998). Struggle with computational bubble
dynamics, in: Proceedings of Third International
Conference on Multiphase Flow, ICMF 98, Lyon, France,
June 812, 1998.
Troshko, A. A. & Hassan, Y. A. (2001) A twoequation
turbulence model of turbulent bubbly flows, Int. J. of
Multiphase Flow, 27, 1965
Yao, W., Morel, C. (2l1***4 Volumetric interfacial area
prediction in upward bubbly twophase flow, Int. J. of Heat
and Mass Transfer, 47, 307
In Fig. 15 the corresponding radial distributions are shown
for two test cases 072 and 074 to investigate the effect of
different liquid turbulence levels. Test 072 was performed at
a lower liquid superficial velocity JL=0.405 m/s with the
same gas superficial velocity.
Besides the two presented cases lots of other tests were
investigated. All show comparable tendencies. Commonly
at z=0.23 m the trace of the gas inlet nozzles disappears. The
lower distances show that for the two presented tests a
turbulent dispersion coefficient of Cm=0.5 seems to be too
small to reproduce the measured results.
Summary and Conclusions
Considering the effect of bubbles on liquid turbulence only
by increasing the turbulent viscosity like in the Sato
approach underestimates the turbulence production in a
vertical bubbly flow in the centre of the tube. The
introduction of corresponding source terms in the turbulence
equations results in a much better agreement to measured
radial profiles of the turbulent kinetic energy. Whereas the
source for turbulent kinetic energy can be derived
considering the work of a rising bubble, there exist no
satisfying confirmation for the corresponding source in the
equation of the turbulent dissipation. Both, the proposal of
(Yao and Morel, 2001), to consider the bubble size, (see
eq. 8) or the one described in the paper, to adjust the
changed viscosity for a homogeneous bubbly flow to the
Sato proposal (eq. 12) yield comparable good agreement to
measurements.
The liquid turbulence modified by gaseous bubbles has a
sensitive influence on the effect of the turbulent dispersion
force. Furthermore a sensitive influence on bubble
fragmentation and bubble coalescence has to be kept in
mind.
Acknowledgements
The work is funded by the German Federal Ministry of
Economics and Technology under the contract No. 150
1329.
References
Antal, S.P., Lahey, R.T., Flaherty, J.E. (1991). Analysis of
phase distribution in fully developed laminar bubbly
twophase flow, Int. J. Multiphase Flow, Vol. 17, pp.
635652.
