7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Bubble Breakup in EulerLagrange Simulations of Bubbly Flow
Y.M. Lau, N.G. Deen* and J.A.M. Kuipers
Fundamentals of Chemical Reaction Engineering, IMPACT, University of Twente, Enschede, The Netherlands
Corresponding author: N.G.Deen@tnw.utwente.nl
Keywords: EulerLagrange simulations, bubbly flow, breakup
Abstract
A basic model for bubble breakup was developed and implemented in an EulerLagrange model. In this approach,
bubbles are tracked by solving Newton's second law of motion. Bubblebubble collisions are accounted for by means
of a hardsphere model developed by Hoomans et al. (1996), and coalescence is modeled with the model proposed by
Darmana et al. (2001). Bubble breakup is based on comparing the bubble Weber number with a critical value, which
depends on the properties of the surrounding liquid and the bubble shape. Simulations were performed for a square
bubble column, which was studied earlier experimentally by Deen (2001), with superficial gas velocities of 5 mmn/s.
The results of the lower velocity showed reasonable good agreement compared with measurements for liquid and
bubble velocities.
Nomenclature
Roman symbols
C coefficient (m)
D diameter (m)
E bubble aspect ratio ()
E6 Eitvos number ()
F force vector (N)
IPI interferometric particle imaging
M Morton number
PIV particle image velocimetry
Q rate of change of the distribution function ()
R radius (m)
V volume (m3)
We Weber number ()
m number of daughter bubbles ()
n normal vector ()
n number of bubbles ()
g gravitational constant (nms1)
h film thickness (m)
p pressure (Nm 2)
t time (s)
u liquid velocity vector (nms 1)
v bubble velocity vector (nms 1)
x location vector (m)
Greek symbols
Q breakup frequency (s 1)
ID interphase transfer source term (Nm 3)
3 daughter size distribution ()
e fraction ()
p density (kgmn3)
a surface tension (Nm 1)
T stress tensor (Nm 2)
C critical Weber value correction ()
Subscripts
D drag
G gravity
L lift
P pressure
S source term for number density
VM virtual mass
W wall
b bubble
break breakup
coal coalescence
coll collision
crit critical
1 liquid
ph phase change
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Introduction
equations of motion are written as:
Due to the computational costs industrial bubbly flows
are usually modeled with an EulerEuler model, also
known as the TwoFluid Model (TFM). Here the liq
uid and dispersed phase are both treated as continuous
phases with interfacial closures. However, typical two
fluid models are less suited to handle complex phenom
ena such as coalescence and breakup, because only one
characteristic size of the dispersed phase is used. For
that reason a population balance equation (PBE) is con
sidered for handling the variations in the size and shape
distributions of the dispersed phase statistically.
On the contrary, with an EulerLagrange model, the
socalled Discrete Bubble Model (DBM), the bubble
size distribution can be determined in a straightforward
manner, provided that coalescence and breakup mod
els are incorporated. Darmana et al. (2001) already
considered coalescence, while Van den Hengel et al.
(2005) applied stochastic models for both coalescence
and breakup mechanisms. To improve the modeling of
the breakup of bubbles a fundamental breakup model is
developed and described in this paper. The model is sub
sequently tested for the case of a square bubble column.
Discrete Bubble model
In the DBM each individual bubble is treated in a La
grangian manner, while the liquid phase motion is de
scribed on an Eulerian grid, taking into account the cou
pling or interaction between the gas and the liquid phase.
Liquid phase hydrodynamics. The liquid phase is rep
resented by the volumeaveraged NavierStokes equa
tions, which consist of the continuity and momentum
equations:
(Cfii) +V cipu 0
at
09
(e6P1U)+V.e p quu = etVP Vte1T+ep1g+<
(2)
The presence of the bubbles is reflected by the liquid
phase volume fraction et and the total interphase mo
mentum transfer i due to action of the interface forces
between the liquid and the bubbles. The liquid phase
flow is assumed to be Newtonian and a subgridscale
model by Vreman (2004) is employed for the turbulence.
Bubble dynamics. The bubbles are tracked by solving
the 2nd law of Newton for each individual bubble. In
terfacial forces are taken into account by the net force
E F, experienced by each individual bubble. Then the
dxb dv
dt v P9dt
with Xb the location, v the velocity and Vb the volume
of the bubble.
The net force acting on each individual bubble is cal
culated considering all the relevant forces (see table 1),
which are separate and uncoupled distributions originat
ing from far field gravity, pressure, drag1, lift2, virtual
mass3 and wallinteraction4:
SF F F + F FD+FL FVM+ FW (4)
Collision/Coalescence. Collision and coalescence of
bubbles are treated with a model described by Darmana
et al. (2001), who combined the hard sphere collision
model by Hoomans et al. (1996) and the coalescence
model by Sommerfeld et al. (2003).
For a given bubble collision pair a and b, the film
drainage time for coalescence to occur is calculated as:
tab pi 11
16o
h}
hf
where the ho is the initial film thickness and hf the final
film thickness just before film breakage. The equivalent
bubble radius for a system of two different sized bubbles
is obtained from:
a1 1b
R2 R
The contact time between two bubbles is calculated by
assuming that it is proportional to a deformation distance
divided by the normal component of the collision veloc
ity:
t CcRab
tab v, v (7)
where Cc represents the deformation distance normal
ized by the effective bubble radius.
The criteria for the occurrence of coalescence or col
lision (and bounce off) is determined by the ratio of the
contact time and the film drainage time. When the con
tact time is less than the film drainage time (tab < tab)
coalescence will not occur and the bubbles will bounce
off each other. In the other case (tab > tab) coales
cence will commence and a new coalesced bubble is cal
culated.
'Ishii & Zuber (1979)
2Auton (1983)
3Auton (1983)
4Tomiyama et al. (1995)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Overview of forces acting on a bubble.
Force Closure
FG pggVb
Fp VbVP
FD Cv,. IV UI(v U) CD = uEO65
FL CLPVb(V u) x (V x u) CL 0.5
FVM  (cvMflVbDt"u+ CVMPlVb(V u ) CVM = 0.5
S d[ 1 1 i. f exp(0.933Eo+ 0.179) 1< E6 < 5
C 2 (Ly j u) n O2y 0.007E6i+0.04 5 < E6 < 33
Breakup models in literature
Before deriving a breakup model for the DBM, breakup
models in literature are evaluated on the basis of suitable
implementation in the DBM.
Breakup has several causes: thermodynamic mech
anisms, such as cavitation collapse, and fluid dy
namic mechanisms, such as high relative veloci
ties/acceleration, resonance and turbulence. Most
breakup models in literature, as described in reviews by
Lasheras et al. (2002) and Liao (2009), are statistical
models and consider breakup due to turbulence. These
statistical breakup models are only applicable in TFM,
where the continuous phase is described in the Eulerian
framework with the k E turbulence model and where
a PBE is solved for the dispersed phase. This PBE de
scribes the temporal and spatial rate of change of a dis
tribution function n(D, x, v, t), which is defined as the
probable number of particles (bubbles) with diameter D,
located at x, with a velocity v, at time t:
otn + Vx (vn) + Vv (Fn) = (8)
Sph + Scoll + Scoal + Sbreak
where the rates of change of n with time due to breakup,
coalescence and collision are denoted by Sbreak, Scoal
and Sco11. The force per unit mass acting on a particle
(bubble) is denoted by F, and the rate of change with
time of its diameter due to evaporation, condensation or
dissolution is given by Sph. The breakage kernel Sbreak
is usually modeled as:
Sbreak (D)
j m(Do)",(D, Don)(Do)n(D,t)dD
(D)n(D,t) (9)
where m(Do) represents the mean number of bubbles
resulting after the breakage of a mother bubble hav
ing the diameter Do, 3(D, Do) is the distribution of
the daughter bubbles formed after the fragmentation of
a mother bubble with diameter Do and (2(D) is the
breakup rate of bubbles with diameter D.
The breakup rate in the breakup kernel of a PBE is
expressed as the breakup fraction multiplied by the fre
quency. The breakup fraction is usually modeled in
terms of a socalled activation complex, e.g. in the
Prince & Blanch (1990) model the activation complex
is a function of the ratio between surface energy and the
mean turbulent kinetic energy. The frequency is deter
mined by the collision rate between the bubbles and the
turbulent eddies.
For a complete breakup model, one must include
an expression for the size distribution 3(D, Do) of the
daughter bubbles resulting from breakup of a mother
bubble of size Do in addition to a model for (2(D). This
term can be formulated with: statistical models, phe
nomenological models based on the change in surface
energy of a breaking bubble and hybrid models which
are based on a combination of both. In all cases, both
the shape of the daughter probability density function
(pdf) and the number of daughter bubbles formed by a
breakup event, m(Do), must be determined. The latter
is done by either assuming a given number of daughter
bubbles a priori, or by deriving an empirical relation for
the number of daughter bubbles from available experi
mental data.
Breakup model for the DBM
From the breakup models in literature none is really suit
able for implementation in the DBM, because all models
are derived for the PBE and breakup is determined statis
tically rather than deterministically. Here a fundamental
breakup model is derived for the DBM with the same
structure as the ones from literature: a breakup criteria
and a daughter bubble distribution.
Breakup criteria. Hinze (1955) is one of the first to
establish a breakup theory by introducing a dimension
less ratio between the force which causes deformation
and the surface tension which tends to restore the bubble
Table 2: A few experimentally derived critical Weber
values in literature.
Flow conditions Critical Weber value
Turbulent pipeflow 1.1 (t)
Microgravity conditions 5
Balancing drag and buoyancy 0.25
Turbulent upward jet 2.5
sphericity. In turbulent flows the deformation is induced
by inertia and the dimensionless number is known as the
Weber number:
We = p 2d (10)
where 6u2 is the mean square velocity difference over
a characteristic distance. Assuming that only velocity
fluctuations over a distance close to the bubble diameter
are capable of causing large deformations, the character
istic distance is taken equal to the bubble diameter.
The theory implies that there exists a critical Weber
value above which the bubble breaks up. Various exper
imental investigations (see table 2) have been conducted
in search for this critical value, e.g. pipeflow5, micro
gravity6, dragbuoyancy balance7, upward jet8, etc.
For the DBM a critical value is derived from the force
balance between external stresses from the continuous
phase, which attempt to destroy the bubble, and the sur
face stress of the bubble plus the viscous stress of the
fluid inside it, which restore its form. Since the surface
stress (surface energy) is dependent on the surface area
of the bubble, this critical value incorporates a correction
factor, which is a function of the bubble aspect ratio. The
critical Weber value is as follows:
We = 12 ((E)
(1+ 2E1 6075) 02213
3E10717 )
and from Clift et al. (1978):
1
1 + 0.1631. "
(E6 < 40, M < 10 6)
Daughter bubble distribution. Binary breakup of bub
bles is assumed and the breakup fraction, that deter
mines the volumes of the daughter bubbles is derived
5Hesketh et al. (1991)
6Risso (2000)
7Senhaji (1993)
8Sevik & Park (1973)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
from a random number distribution, for which an U
shaped beta probability density function (see Figure 1)
is chosen. This function assumes that the probability of
equal volume breakup is smaller than breakup in a small
and a large daughter bubble. As for the locations, the
larger daughter is placed at the same location as the par
ent bubble, with the smaller daughter bubble randomly
placed around it (see figure 2).
10
9
8
7
6
LL
o 5
0
4
3
2
1
0
0 0.2 0.4 0.6 0.8
Relative daughter bubble volume
Figure 1: Beta probability density function of the
daughter bubble size distribution.
/
Figure 2: Illustration of the binary breakup process of
a parent bubble in two daughter distribution, with the
larger bubble at the original location and the smaller
bubble randomly placed around it.
Results and discussion
Liquid and bubble velocities. Simulations were per
formed for the experimental setup of Deen (2001),
who performed PIV measurements in a square bubble
column with a perforated bottomplate and a superficial
gas velocity of 5 mm/s. The inlet bubble diameter was
set to 1.5 mm and a snapshot of the DBM simulation is
illustrated in figure 3, where various data (gasfraction,
bubble radii, locations and velocities) are plotted. To
evaluate the DBM results, liquid and bubble velocities of
experiment and simulation are compared at three differ
ent heights (0.15m, 0.25m and 0.35m) within the bub
ble column. The averaged vertical liquid velocities are
plotted in figure 4, along with its standard deviation in
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Sasfroction
0.44
0.37
0.29
0.22
0.15
0.07
0.00
radius (m)
5.2e03 ..;
4.4e03
3.6e03
2.8e03
2.0e03
1.2e03
3.7e04
rI ..r, (m/s)
:0 39
': J49
0.359
0.269
0.180
0.090
0.000
Figure 3: Simulation results for a superficial gas velocity of 5 mm/s: (a) gas fraction; (b) bubble locations and radii;
and (c) liquid velocity flow field.
(a) height = 0.15m (b) height = 0.25m (c) height = 0.35m
Figure 4: Timeaveraged vertical liquid velocities at different heights.
+ + + + + + +
018+ +
016
014 +
012
S 01
.8 oo
006
004
002 operment +
simulaton 
+ + + ++++
percent +
simulation 
02 04 06 08 1 0 02 04 06 08 1 0 02 04 06
XD XD XD
(a) height = 0.15m (b) height = 0.25m (c) height = 0.35m
Figure 5: Standard deviation of the vertical liquid velocities at different heights.
+ + +
+ + + +
005
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
mula t35 n
(a) height = 0.15m (b) height = 0.25m (c) height = 0.35m
Figure 6: Timeaveraged vertical bubble velocities at different heights.
18000
16000
14000
12000
10000
8000
6000
4000
2000
0
experiment
simulation
0.002 0.004 0.006 0.008 0.01 0.012
diameter[m]
Figure 7: Comparison between bubble size distributions from experiment and simulation at a superficial gas velocity
of 5 mm/s.
figure 5 and the averaged vertical bubble velocities in
figure 6. The averaged liquid velocities at the bottom of
the column show good agreement with the experiment,
while in the higher region of the column the values are
underpredicted. This also holds for the standard devia
tion, but differs from the averaged vertical bubble veloc
ities. Here there is overprediction in the lower region,
while in the higher region the velocities at the centerline
show good agreement. A reason for the differences be
tween simulation and experiments could be due to the
applied closure relations, which are strictly only valid
for single bubbles. Using closure relations derived for
heterogeneous bubbly flows could yield a better predic
tion of the liquid and bubble velocities.
Bubble size distribution. Hansen (2009) used the in
terferometric particle imaging technique (IPI) to mea
sure in the same setup with similar settings the bubble
size distribution. The experimental and simulation re
sults are plotted in figure 7. The experimental data con
tains more smaller bubbles than the simulations, imply
ing that the simulation might require a better definition
of the boundary condition (inlet bubble diameter). Be
sides this difference, the bubble size distribution of the
simulation does show the same qualitative trend with the
experiment.
Conclusions
A fundamental breakup model has been successfully for
mulated and implemented. The obtained simulation re
sults shows that the DBM extended with a coalescence
and a breakup model can be used to study the bubble size
distribution. For the case of a low superficial gas veloc
ity the trend of the bubble size distribution agrees well
expsr,%; __
34 06
with the experiments. Velocities of liquid and bubble
show some agreement. However, simulations with im
proved closure relations are still needed for better eval
uation. For cases of high superficial gas velocity more
investigations are still required.
Acknowledgements
This project is part of the Industrial Partnership Program
"Fundamentals of Heterogeneous Bubbly Flow", which
is funded by FOM, AkzoNobel, Corus, DSM and Shell.
References
Auton, T.R., The dynamics of bubbles, drops and par
ticles in motion in liquids., PhD Thesis, University of
Cambridge, Cambridge, U.K., 1983.
Clift, R., Grace, J.R. and Weber, M.E., Bubbles, drops
and particles., Academic Press, New York, 1978
Darmana D., Deen N.G. and Kuipers J.A.M., Paral
lelization of an EulerLagrange model using mixed do
main decomposition and a mirror domain technique:
Application to dispersed gasliquid twophase flow,
Journal of Computational Physics, Vol. 220(1), pp. 216
248, 2006.
Deen N.G., Solberg T. and Hjertager B.H., Large Eddy
Simulation of the Gasliquid Flow in a Square Cross
sectioned Bubble Column, Chemical Engineering Sci
ence, Vol. 56, pp. 63416349, 2001.
Hansen, R. Computational and experimental study of
bubble size in bubble columns., PhD Thesis, Aalborg
University Esbjerg, Denmark, 2009.
Hengel E.I.V. van den., Deen N.G. and Kuipers J.A.M.,
Application of coalescence and breakup models in a dis
crete bubble model for bubble columns, Industrial and
Engineering Chemistry Research 44: 52335245, 2005.
Hesketh, R.P, Etchells, A.W. and Russel, T.W.F., Ex
perimental observations of bubble breakage in turbulent
flow., Ind. Eng. Chem. Res., 30, 845, 1991.
Hinze, J.O., Fundamentals of the hydrodynamic mech
anism of splitting in dispersion processes., AIChE J. 1,
289a295, 1955.
Hoomans B.P.B., Kuipers J.A.M., Briels W.J. and van
Swaaij W.P.M., Discrete Particle Simulation of Bubbles
and Slug Formation in a Twodimensional Gasfluidized
Bed: a Hardsphere Approach, Chemical Engineering
Science, Vol. 51(1), pp. 99118, 1996.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Ishii, M. and Zuber, N., Drag coefficient and relative ve
locity in bubbly, droplet or particulate flows., AIChE J.
25, pp. 843855, 1979.
Lasheras, J.C., Eastwood, C. MartfnezBazAn, C and
Montai6s, J.L., A review of statistical models for the
breakup of an immiscible fluid immersed into a fully
developed turbulent flow., International Journal of Mul
tiphase Flow, Vol. 28, pp. 247278, 2002.
Liao, Yixiang and Lucas, Dirk., A literature review of
theoretical models for drop and bubble breakup in tur
bulent dispersions., Chemical Engineering Science, Vol
ume 64, Issue 15, pp. 33893406, 2009.
Prince, M.J. and Blanch, H.W., Bubble coalescence and
breakup in airsparged bubble columns., AIChE J., vol.
36, pp. 14851499, 1990.
Risso, F., The mechanisms of deformation and breakup
of drops and bubbles., Multiphase Science and Technol
ogy, Vol.12, pp.150, 2000.
Senhaji, R., Qualification global du fractionnement
d'une phase dispersed de faible viscosit6 en function des
propri6t6s turbulentes de lh6coulement exteme, Thesis,
Ecole Centrale de Nantes, 1993.
Sevik, M. and Park, S. H. The splitting of drops and bub
bles by turbulent fluid flow., J. Fluids Engng 95, pp. 53
60, 1973
Sommerfeld, M., Bourloutski, E. and Broder, D., Eu
ler/Lagrange calculations of bubbly flows with consid
eration of bubble coalescence., The Canadian Journal of
Chemical Engineering, 81: 508518, 2003.
Tomiyama, A., Matsuoka, T., Fukuda, T., and Sak
aguchi, T., A simple numerical method for solving
an incompressible twofluid model in a general curvi
linear coordinate system, In A. Seizawa, T. Fukano,
and J. Bataille, editors, Advances in Multiphase Flow,
pages 241252, Amsterdam, November 1995. Society of
Petroleum Engineers, Inc., Elsevier.
Vreman, A.W., An eddyviscosity subgridscale model
for turbulent shear flow: Algebraic theory and applica
tions., Physics of Fluids, Vol.16, 10, October 2004.
