7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Investigation of mass transfer in polydisperse gasliquid systems by using a
multivariate population balance
A. Buffo* D.L.MarchisioT M.Petitti* M.VanniT N.Mancinitand F.Podenzanit
Dip. Scienza Mat. Ing. Chimica, Politecnico di Torino, C.so Duca degli Abruzzi 24, 10129, Torino, Italy
t Division Refining & Marketing, Eni S.p.a., Via Felice Maritano 26, 20097, San Donato Milanese, Italy
antonio.buffo@polito.it and daniele.marchisio@polito.it
Keywords: population balance modelling, gasliquid stirred tanks, quadraturebased moment methods
Abstract
The investigation of mass transfer in gasliquid systems requires a formulation of a multivariate population balance
equation, in order to consider the effect of a inhomogeneous distribution of bubble size, velocity and composition. A
very efficient method to solve this equation is represented by the direct quadrature method of moments (DQMOM),
where the evolution of the properties of the bubbles is evaluated by tracking several moments of the Number
Density Function (NDF), with the advantage that the obtained equations can be easily implemented in commercial
Computational Fluid Dynamics (CFD) codes. In this work, a detailed model was developed, tested on a simplified
case study and coupled with a CFD code with the aim to provide a complete description of transport phenomena in a
real system. The airwater stirred tank investigated by Laakkonen et al. (2006) was here simulated and results were
compared with the experimental data.
Introduction
Agitated gasliquid reactors are widely used in the
chemical and biochemical processes. In these systems,
gas dispersion is very important since it strongly influ
ences gasliquid mass transfer. Mechanical agitation en
hances the homogeneity of the system improving the
overall performance of the reactor and flexibility of the
process.
The design and scaleup of agitated gasliquid tanks
is difficult, because chemical reactions usually have a
complex relationship with mass, momentum and heat
transfer; moreover, controlling phenomena often vary
with the reactor size. Traditionally, some semiempirical
correlations based on vast experimentation were used
for the scaleup of reactors. These correlations are ex
tremely useful for an initial analysis, but they consider
only volumeaveraged properties calculated over the en
tire system (i.e., global gas holdup) and they can be
applied only to specific vessel geometries. However,
large industrial stirred tank reactors are often strongly
inhomogeneous: ,igilii k.ii spatial inhomogeneities ex
ist even in laboratory scale tanks (Calderbank 1958)
and their effect on global properties must be consid
ered. Computational Fluid Dynamics (CFD) has become
a popular tool for the analysis of chemical reactors, be
cause it allows for the estimation of the fluid flow and
turbulent fields in almost any vessel size and geome
try. At the moment, multiphase CFD predictions are un
certain and require further development and validation
(Ranade 2002). CFD codes, together with the Popula
tion Balance Models (PBM) and rigorous mass transfer
models, offer a great potential for the detailed analysis
of local mass transfer in agitated gasliquid reactors.
As it is well known, gasliquid mass transfer is a com
mon ratelimiting step of reactor performance. Mass
transfer is often controlled by the liquid side through
the volumetric mass transfer coefficient (kLa). Corre
lations for this parameter are popular, because the mea
surement and modelling of the global mass transfer area
for unit volume is quite difficult (Calderbank 1958). Lo
cal interfacial area for gasliquid mass transfer depends
on the Bubble Size Distribution (BSD), which is known
to vary in agitated tanks depending on the operating con
ditions and for the same operating condition from point
to point in the vessel. For example, bubbles near the
impeller breakup due to the high shear experienced,
whereas bubbles accumulate and coalesce in stagnant
zones. Moreover, every single bubble has a different
concentration of chemical species; this is important in
order to determine the driving force for the local mass
transfer flux. Population balances are a generalized ap
proach for modelling the evolution of the local BSD, the
local distribution of concentration of the different chem
ical species and the local bubble velocity distribution.
Of course to solve the PBM one must describe bubble
breakage and coalescence by using phenomenological
models. The final result is the formulation of a multi
variate PBM capable of predicting both size and con
centration distributions.
In the present work such an approach is formulated,
discussed and implemented in a commercial CFD code.
The approach is then used to model bubble coalescence
and breakage as well as masstransfer. Model predic
tions are validated by comparison through experimental
data from Laakkonen et al. (2006, 2007).
PBM for gasliquid systems
A polydisperse gasliquid system can be successfully
described by using PBM. It is possible to consider a pop
ulation of bubbles characterized by their size L, their
velocity Ub and their chemical composition 4b (inter
nal coordinates), defining a Number Density Function
(NDF) so that the following quantity:
n(L, Ub, ?b; x, t) dL dUb d4b dx,
represents the expected number of bubbles, in the in
finitesimal volume dx dx 1 dx2 dx3 around the phys
ical point x (= (x, 2, X3) with size in the range be
tween L and L + dL, velocity in between Ub and Ub +
dUb and composition in between 4kb and 4b + d4b.
The evolution of this NDF is governed by the Gener
alized Population Balance Equation (GPBE) that reads
as follows:
On a a a
t + a (Ubin) + (Gn) + I (Abin) +
at 0dxi dL OUbi
+ f .) '9h(L, Ub, b;x,t), (1)
where UbV is the ith component of the bubble veloc
ity, Ub (with i E 1, 2, 3), G is the rate of continuous
change of bubble size usually related to molecular pro
cesses, where bubbles continuously grow (or shrink) be
cause of the addition (or depletion) of single molecules
(e.g., evaporation or condensation processes), Abi is the
ith component of the acceleration, Ab, acting on a bub
ble as a result of continuous forces, Ybj is the continu
ous rate of change of bubble composition with respect
to the chemical component j (with j E 1,... q) as a re
sult of chemical reactions and/or mass transfer between
phases. The term on the righthand side is the discon
tinuous event term, accounting in this case for instan
taneous changes of size, velocity, and composition due
to bubble collisions, coalescence and breakage. In this
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
work it is assumed that G and 4bj depend only on mass
transfer of the different chemical components. More
over, the acceleration acting on a bubble, Ab, is sim
ply considered as a result of buoyancy, gravity and drag
force; any other inferfacial force is in this context ne
glected, but could of course be added if necessary.
In the literature different methods are available to
solve the GPBE, but not all of them are able to treat
multivariate problems and not all of them result in
computationally efficient methods capable of describ
ing real systems with complex geometries. For example,
both Classes Methods (CMs) and Monte Carlo Methods
(Vanni 2000; Zucca et al. 2007) are generally not used
for these problems, due to their computational requests.
It was demonstrated that methods based on quadrature
approximations, such as the quadrature method of mo
ments (QMOM) (McGraw 1997; Marchisio et al. 2003),
and the direct quadrature method of moments (DQ
MOM) (Fan et al. 2004; Marchisio et al. 2005), could
be amenable for solving PBM coupled with CFD codes;
QMOM has been recently used for the simulation of gas
liquid systems (Petitti et al. 2010) however some limita
tions where in that case detected. First of all the fact
that, all the bubbles are moving with same velocity, does
not allow for accounting for some segregation phenom
ena that can be described only by tracking the distribu
tion of bubble velocity. Secondly the introduction of ad
dictional internal coordinates is tricky and not straight
forward. For this specific application is interesting to
explore the possibilities offered by DQMOM. DQMOM
involves the direct solution of the transport equations for
weights and abscissas of the quadrature approximation.
Moreover, each node of the quadrature approximation
can be treated as a distinct gas phase, with its own ve
locity, size and composition.
Applying the quadrature approximation, the NDF be
comes:
n(L, Ub, kb; x,t)
N
= E wa(L L,)65(Ub Ub,a)6(b 0b,a),
a=l
(2)
where Wa is the number density (i.e., number of bub
bles per unit of total volume) of the bubbles with size
equal to La, velocity equal to Ub,a and with composi
tion vector equal to 4b,. Each delta function is char
acterized by a different index a and it can be thought
of as a group of bubbles having the same size, veloc
ity, and composition, each group representing one of N
nodes of the quadrature approximation. In order to re
duce the number of equations, it is possible to consider
only two nodes (N = 2) with an acceptable lack of ac
curacy as explained by Marchisio et al. (2005); more
over, only two chemical components (resulting therefore
in only one chemical composition) are considered in this
work. Readers familiar with the method will recognize
that this is not a serious limitation since in DQMOM
internal coordinates can be easily added. For brevity in
what follows, space and time dependencies will be omit
ted.
By substituting Eq. (2) into Eq. (1), and applying a
moment transformation, it is possible to derive the fol
lowing transport equations:
a a
(wi)+ ,(w1Ubi,l)
at ore,
al, (3)
at (w2) + a (W2Ubi,2)= a2,
a a
(wiLi) + a (W1Ubi,lL1) b +
SwiGi(L1, Ub,, 1b,1),
a a
at (w2L2) + x (2Ubi,22) = b+
+ w2G2(L2, Ub,2, Ob,2),
a a
a (Wlb,1) + a (W1Ubi,lb,1) = di
+ Wlbl(L1, Ub,1, b,1),
S(w24b,2) + 2 (2bi,2 b,2) = da2
+W2b2 (L2, Ub,2, Pb,2),
a a
a W (WlUbj,1) (W 1 Ubi,1 Ubj,l)
a~t Oxzi
Cjl+
+ lAbj(L, Ub,1, b,1), (9)
(w2Ubj,2) + (W2Ubi,2Ubj,2) = Cj2+
at Oxi
+W2Abj(L2, Ub,2, b,2), (10)
where the subscripts 1 and 2 refer to the two nodes of
the quadrature approximation.
Since only mass transfer of a single component is con
sidered, Yb,i can be written as:
bi b(Li, Ubi, b,i) ( Hb,i), (11)
where kL is the mass transfer coefficient, 0c is the con
centration of the chemical specie in the continuous liq
uid phase and Ob,i is the concentration of the chemical
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
species in the gas bubble and H is the Henry constant.
Mass transfer affects also the size of the bubbles through
the bubble growth rate G. Through a simple mass bal
ance on a single bubble it is possible to write that:
i = G(Li, Ubi, Yb,i) 2kLM (0c Hb,i), (12)
Pb
where Mw is the molecular weight of the chemical
species considered. The acceleration acting on a single
bubble, Ab, can instead be written as follows:
ebAbi = bAb(Li, Ubi, Pb,i)
3 acCD,ic
4 L,
ap
0 + bS+
axi
Ubi (U Ub),
(13)
where g is the gravity acceleration, yc and eb are the
densities of the continuous (liquid) and disperse (gas)
(5) phases, ac is the liquid volume fraction and (Uc Ubi)
is the slip velocity, namely the velocity difference be
tween the continuous phase and the bubbles and CD,i
CD(Li) is the drag coefficient. There are several equa
tions available in the literature to estimate CD,i; in this
(6)
work the Haberman and Morton model (reported in Wal
lis 1969):
4 (Vc eb)Lg
CD,i 3 U,2 (14)
3 PcUt2
(7) is used. This model is based on the terminal velocity of
a bubble of diameter Li in a liquid at rest, Ut, where is
evaluated through the Clift correlation (reported in Clift
et al. 1978):
( (2.14a
Ut = + 0.505 gLi)
V ecL
Li > 1.3mm
(15)
where a is the surface tension. Moreover an adhoc cor
rection taking into account the effects of turbulence and
of the presence of other bubbles is also implemented in
this work (details are reported in Petitti et al. 2010). The
terminal velocity used to compute the drag force was
corrected with respect to the value obtained for an iso
lated air bubble rising in water at rest, in order to take
into account for the effect of turbulence. Although this
correction should be applied and calculated on the local
turbulence intensity and gas holdup, in this work a con
servative approach was used, considering a constant ter
minal velocity of 13 cm/s through out the entire stirred
tank.
The source terms of this set of equations (Eqs. 310)
are calculated by forcing the method to track some spe
cific moments of the NDF. The choice of the moments
is completely arbitrary and in this work only pure mo
ments are selected. This choice allows to decouple the
evolution of the different internal coordinates, as it will
become clearer below. It is examplificative to examine
how the source terms are calculated. Let us start with a1,
a*, b* and b*. These four quantities are calculated all to
gether by solving the following linear system, written in
matrix form as:
Aac d, (16)
where the coefficient matrix stands for:
1 1 0 0
2 2 (17)
A L2 2Li 2L2 (17)
2L 2L2 3Li 3LJ
and the vector of unknowns a is defined by:
S= [al, a2, b", b ]T (18)
where
2kg M
b = bi wi ( Hb,i). (19)
eb
The known right hand vector d can be written as follows:
d = [Sooo, S1oo, S2o, 30 oo]T (20)
where Skoo with k e 0, 1, 2, 3 corresponds to the source
term for the k* pure moment with respect to the size L.
By employing the quadrature approximation, this can be
written in the following closed form:
13 N N3
Skoo 2= W w (Li + L) /a(Li, L,)+
i=1 j=l
N N
Lw, 5 a(Li, Lj)wj+
i=1 j=l
N
+ b(L .(kO)
+ f = (1 '0
N
'i(Li)wi
i=l
where a(Li, Lj) is the coalescence kernel, b(Li) is the
breakage kernel and b() is the distributive daughter
function. For the moment these quantities are assumed
to depend only on bubble size. However, bubble compo
sition dependencies can be easily accounted for.
Repeating the same procedure for d* and d*, it is pos
sible to define:
A [1 1 ]2
A [21 212 '
a =[d*, d ]T,
d =[Solo, So20 + a al+ ]2a21]T,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where
6kL
d = di w H b,i), (25)
SLi
and Solo with I E 1, 2 corresponds to the source term for
the th pure moment with respect to the composition y:
Solo = 2 S wi a(L L )w J L + L 3 +
i=1 j=1 1 3
N N
wi a(Li, Lj)wj+
i=1 j 1
N
+ b(Lj)wjb t')
+i= 1
N
Slb(Li)wi.
The expression of bW') differs from b ko), as it will be
explained below.
The same procedure is applied to obtain the source
terms for the momentum equations c,1 and c ,2. These
terms account for the change in momentum of the two
bubble classes due to discontinuous events such as bub
ble collisions. For dilute systems, however, this effect
is neglegible with respect to drag and buoyancy forces,
therefore it will, in this contest, be neglected.
It is worth mentioning that this choice (k E
0, 1, 2, 3; 1 E 1, 2) allows us to track the first pure four
moments of the NDF with respect of the bubble size and
the moments of order one and two with respect to bubble
composition and velocity.
Before calculating Eqs. 21 and 26, it is necessary to
define the expressions for the distributive daughter func
tions, coalescence and breakage kernels. In this work,
the following breakage kernel will be used (Laakkonen
et al. 2006):
b(L) 6.0 e/3 x
x erfc 0.04 2/3L5/3 + 0.01 Cd13L4/
(27)
where p is the viscosity of the continuous phase and e
is the rate of energy dissipation.
To model the frequency of bubbles coalescence it is pos
sible to refer to kernel expression of Prince et al. (1990):
a(A, L) = 0.88 /3(A + L)2(A2/3 + L2/3)/2(A, L),
(22) (28)
where the coalescence efficiency can be written by using
the model of Coulaloglou et al. (1977):
(23)109 t .
(24) q(A,L) exp 610'A ) (29)
Concerning the distributive daughter function, it can be
realistic to assume binary breakage (as explained by An
dersson et al. 2006) resulting in:
b(k,) 3240L (30)
(k + 9)(k + 12)(k + 15)
When 1 = 0, Eq. 30 becomes:
b(ko) 3240L (31)
(k + 9)(k + 12)(k + 15)'
and when k = 0, Eq. 30 can be written:
b') = 2y. (32)
The mass transfer coefficient kL is estimated by using
the Lamont and Scott model (Lamont et al. 1970):
k= 0.4 D5 o ,25 (33)
where D is molecular diffusivity of the chemical species
in the liquid phase and v, is the kinematics viscosity.
Now all the transport equations constituting DQMOM
are closed and through their solution it is possible to re
construct the evolution of this polydisperse gasliquid
system. Coupling these equations to a commercial CFD
code, however, requires to rewrite them in terms of vol
ume fractions:
7 3
ai =Li wi, (34)
where ac is the volume fraction of the generic ith dis
perse phase. Thus, the final governing equations be
come:
a a
at (alb) + a (alebbi,l) 3
=e(b (2L L al, (35)
a a
S(29eb) + (a2Pebbi,2)
a a
at (albL1) + a (a1ebLUbi,) )
at oxi
/(2 ( ,
Qb 7rL3 b
a a
a (02bL2) + (02PbL2Ubi,2)
at 9xi
2 La) ,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
S(alebb,l) + 0, (aebbb,lUbi,1)
Qb (7 L1 + Yb,1 1b L a)) (39)
a a
(aQ2bb,2) + a (a2Pbb,2Ubi,2)
Pb ( 3 2 2+ b,2 ( 2 Lb2 L 22)) (40)
(6 2 2 (2 2 2 2 2L* 7 7
a a
a (cal bUbjI) + a (alpbUbi,lUbj,l)
a~t Oxzi
ap
S + alebgj+
Oxj
3 ocCDec
+ 1 IU Ubj,1 (Ucj
4 L,
Ubj,1), (41)
a (a2PbUbj,2) + a (2b bUbi,2Ubj,2)
2p
= a2 + a2gbgi+
axi
3 coCD c
+ 4 a D U Ubj,2 (Uj Ubj,2). (42)
4 L2
These equations must be solved together with the
equations for the continuous phase:
a a
at (ac&) + a (ccUci) = 0 (43)
aa a xi
S(accUcj) + (accUciUcj) =
hp
= c + acecgj+
+ ac(Pc + T,c) + )) +
2 OUck ,
3 xk
2 3 CD
4 ca k Uc Ub Icj bj,k)
k(=1
(44)
S(Ocecc) + (cecQcOci)
at 6xi
,?c (kL6a (4c Hfb,j) + kL 60(c
L, L2
Hb,2))
(45)
where PT,c is turbulence viscosity, modelled as:
2 L2 a2
Pb 2 L 3b
3
PT,c = 0.09 Pc
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Summary of gas and liquid properties used in
computations.
yc 998.2 kg m3
Pb 1.255 kg m3
fPc 1.0 103 Pas
Pb 1.78 105 Pas
a 0.07 N m1
The values of turbulent kinetic energy k and rate of tur
bulent energy dissipation c are obtained, in this work,
through k c Multiphase Mixture Model (Kataoka et al.
1989).
Test case, operative conditions and numerical
details
The gasliquid stirred tank showed in Fig. 1 is investi
gated in this work. The reactor is a 194 liter fourbaffled
reactor, agitated by a sixblade Rushton turbine with cir
cular metal porous sparger with diameter of 3.3 cm and
pores of 15 pm, located 10.5 cm below the impeller
and used for low gassing rates. This particular config
uration was already investigated by Laakkonen et al.
(2006, 2007) and detailed experimental measurements
are available for model validation. The continuous liq
uid phase is water and the disperse gaseous phase is air;
the properties used in the simulations are showed in Tab.
1.
The simulations were carried out with the Multiple
Reference Frame approach and were performed by con
sidering half of the geometry, since this is the smallest
portion of the reactor that satisfies the conditions of geo
metrical symmetry and periodicity. 230,000 hexahedri
cal cells were used to mesh the geometry, since both
the flow and the power number did not show, igMil Ik.iil
change with further grid refinement. The computational
grids were created by using GAMBIT, whereas simula
tions were performed with the commercial code FLU
ENT 12 and UserDefined Functions and Subroutines
were used to implement the DQMOM equations and the
drag force based on bubble terminal velocity. The trans
port equations for DQMOM were implemented by defin
ing additional scalars for the disperse gaseous phases.
CFD predictions were validated by comparison with ex
perimental data. The predicted BSDs were compared
with the experimental ones in the measurement points
reported in Fig. 1 for a wide range of operating condi
tions. The stirring rate ranged between 155 and 250 rpm
and the gas flow rate ranged from 0.018 to 0.093 volume
of gas per volume of reactor per minute (vvm), resulting
in global holdup values up to 1.5 %. Validation through
comparison with transient mass transfer experiments (by
Figure 1: Rappresentation of stirred tank simulated and
location of sample points in which BSD data was de
tected
using oxygen and nitrogen) was also carried out and only
preliminary results are reported here.
As far as the boundary conditions are concerned, the
section of the sparger from which the gas is introduced
in the stirred tank was defined as a velocity inlet; the
gas volume fraction was evaluated as portion of the ge
ometrical area available for gas flow and then the gas
velocity was determined by knowing the gas volumetric
flow rate. The liquid velocity components were set null.
The upper surface of the reactor was defined as a pres
sure outlet in order to let the gas exit from the system
and, if there is a backflow, it is formed by liquid only.
Setting the boundary condition for the inlet bubbles it
is a nontrivial problem. The bubbles at the metal porous
sparger were modelled as a lognormal distribution as
suming a standard deviation of 0.15, as usually observed
with metal porous spargers, and with a mean bubble size
at the inlet ( 1 . calculated with the correlation pro
posed by Kazakis et al. (2008):
/. 7, 35 We 1Re'1Fr1'8 dp171 (47)
where the dimensionless numbers of Reynols (Re), We
ber (We), Froude (Fr) are defined as follows:
ecU2 d
We 9c U
ScUgsds
and where U., is the gas superficial velocity based on
sparger area, d, is the sparger diameter and d, is the
mean pore size.
From this lognormal distribution the first mo
ments with respect to bubble size were calculated and
then through the ProductDifference algorithm (Gordon
1968) weights and abscissas were also determined. With
this procedure the inlet volume fractions for the two bub
ble classes (ai, a2) and their characteristic sizes (L1,
L2) were calculated. The inlet velocities were assumed
to be equal (Ub,i Ub,2 9.94 ms 1). The av
erage inlet composition, written in terms of the inlet
oxygen concentrations in the two bubble classes (4b,1,
Yb,2), was instead fixed to be identical to that of the ex
periments. For numerical and stability reasons also for
the composition a lognormal distribution with a stan
dard deviation of 0.15 and a mean concentration of 8.56
mol/m3 was considered.
Some tests were also conducted on a homogeneous
system in order to verify the model and to detect all pos
sible numerical problems. The set of equations was for
this simplified case solved in MATLAB.
Results and discussion
Firstly the tests carried out on the simple homogeneous
system will be presented, then the CFD simulations
and the comparison with experimental data will be dis
cussed.
a) PBM for homogeneous system. Several tests were
conducted using MATLAB for an hypothetical homoge
neous system constituted of a simple computational cell.
In all of them mass transfer was not considered, and the
attention was focused on bubble coalescence and break
up in order to evaluate the consistency of the multivariate
model.
An example of these tests is reported in Fig. 2, which
shows the time evolution of the quadrature nodes and
weights in the case of N = 2 and in the case of pure
coalescence and pure breakage. As it is seen for pure
coalescence, the number of bubbles per unit volume,
wl and W2, decreases with time, while the characteristic
size of these two groups of bubbles, L1 and L2. When
pure breakage was simulated, the number of bubbles per
unit volume (the sum of wl and w ) increases with time,
while the characteristic sizes of these two groups of bub
bles, L1 and L2, decrease.
In both cases the characteristic concentration of oxygen
in the two bubble classes, after a short period of time,
tends to the same value at about the average concentra
tion. When such a situation takes place, the sizes of the
two classes are very different from the initial ones and
it is likely that the transfer of bubbles between the two
classes has been so large that the memory of the differ
ent initial concentration has faded.
When simultaneous coalescence and breakage is simu
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
lated with realistic kernels and volumeaveraged turbu
lent dissipation rate (c) (Fig. 3) the process is dominated
by coalescence. All the weights and abscissas of quadra
ture approximation, wl, W2, L1, L2, Yb,i and Ob,2, as
sume correct trends in time and physically consistent fi
nal values. These tests showed the formal correctness
of the method that can therefore be implemented in the
CFD code.
Table 2: Comparison between experimental data and re
sults from simulation.
Distribution of d32 (mm)
155 rpm and 0,018 vvm
Sample point Experimental Data Model
R2 2.37 3.10
R4 2.48 2.56
R8 2.29 2.57
R9 1.65 2.63
R12 3.31 3.09
Distribution of d32 (mm)
220 rpm and 0,041 vvm
Sample point Experimental Data Model
R2 2.56 2.66
R4 3.34 3.04
R8 2.57 2.47
R9 1.76 2.50
R12 3.81 3.20
Distribution of d32 (mm)
250 rpm and 0,052 vvm
Sample point Experimental Data Model
R2 2.74 2.45
R4 2.93 3.31
R8 2.17 2.55
R9 2.01 2.65
R12 3.18 3.57
b) CFD simulations on the real geometry. The re
sults presented and discussed here mainly focus on the
prediction of BSD and their comparison with experi
mental data from the literature (Laakkonen et al. 2006,
2007); the predictions of gas distribution profiles and
global holdup for this configuration were the subject of
previous work (Petitti et al. 2010), and will be not dis
cussed here.
An example of a typical simulation is reported in Fig.
4. As it is possible to see, the quadrature approxima
tion, composed in this case of only two nodes (N = 2),
can be thought of as a description of the population of
bubbles in two classes. As the bubbles exit from the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10 x 3 15
?2 1 I
5 2 1 0.5
0 0 0 0
0 0.01 0.02 0.03 0 0.01 002 0.03 0 a01 002 0.03 0 0.01 0.02 0.03
1 t t t
0.02 0.03 1.5x10' 2 x107
0.02 I 1 2
j0.01 0,01 0.5 1 1
0 0 0  05
0 0.01 0.02 0.03 0 0.01 002 0.03 0 001 0.02 003 0 0.01 0.02 0.03
I t t t
10 10 10 10
95 95
S9 J9 > t 9 90
8.5 85
8 88 8
0 0.01 0.02 0.03 0.01 002 0.03 0 0.01 0.02 0.03 0 0.01 0.02 0.03
I t t t
Figure 2: Plots of weights and abscissas of NDF, calculated with DQMOM and kL 0, in the case of pure coalescence
(left) and pure breakage (right).
x104
5
10
0
0 0.01 0.02 0.03 0.04 0.05
x 10
6
4
C\N
2
0
0 0.01 0.02 0.03 0.04 0.05
x10
0.01 0.02 0.03 0.04 0.05
0 0.01 0.02 0.03 0.04 0.05
9
8
0 0.01 0.02 0.03 0.04 0.05
0.01 0.02 0.03 0.04 0.05
Figure 3: Plots of weights and abscissas of NDF, calculated with DQMOM and kL
coalescence and breakage.
0, in the case of simultaneous
9
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
i ,:. ,
S ,9.94et0D
9 54.OD
9 24e+0
S95e+0D
S8.65e0D
8.0 e+00
S' .... 7.45et00
S7.16e2OD
S86et00
6 194e+OD
96 1 e+00
5 "01 8 .5e O0
o, 0..0.537eO0
SO O 417e00
u 4 t 7o .5e+00
2.31 6et00
2 09e+O
S279e+D00
61.49et0+O
1 19e00 1 39300
4 95e01 795eDO
5.96e0O 4.6e 0
S00e00 017e+D00
99 e+0D 9.604e+00
743e00 8O 43e+O
8 40et00
S20eDO
6.620
528 e00, 1 337e+00
I .04e+00 I .0e+D00
4.57e O I 546eOO
4.32e+00 5 12e+00
3.8400 460e+DO
O+00 4.31e+00
S3600 4.80e+DO
3.12+00 3I7e+00
2 40eD 2 1.et00
2.16e,00 7O63e.00
1680e00 2.00e300
1.44 00 1.75e+00
1.20e+00 1.42e+DO
9.57.01 4.89e+00
477e01 5574.01
2.37201 269e01
2600e00 0 00e000
0.018 vvm.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
15.40
4.44
3.63
2.80
0.12
Operative condition
250 RPM 0.093 vvm
DfiuetLr (.am)
Figure 5: Comparison of Distribution of mean Sauter Diameter d32 (nun) and local Bubble Size Distribution.
sparger, they evolve due to coalescence, breakage and
mass transfer between liquid and gas, and their charac
teristic volume fractions (ai, a2), sizes (L1, L2), veloci
ties (Ub,1, Ub,2) and oxygen concentrations (0b,1, Yb,2)
change in tracking the evolution of the NDF as dictated
by the GPBE.
In Fig. 5, instead, the comparison between experi
mental data and CFD predictions is reported for one op
erative condition investigated. The comparison is here
limited to the BSD and the mean Sauter diameter. In
general, the model is able to describe qualitatively and
quantitatively the experimental trends also under other
operating conditions (see Tab. 2). The experimental
data show a clear trend for the Sauter diameter (d32): in
all the operative conditions it is possible to observe that
dR12 > dR4 > dR8 and dR2 > dRg, where the subscript
indicate the sampling points. Moreover, in general, d32
increases with the gas flow rate: in this case, in fact, the
bubbles collide and coalesce more frequently, whereas
the breakage frequency decreases. The values detected
in zone R2 and R8 are very similar, but for high flow
rates the d32 measured in R2 is higher than in R8; in
other words, for high flow rates, coalescence becomes
more important near the upper surface, where the bub
bles tend to accumulate.
The comparison with CFD results show that the agree
ment with the experimental data increases as the gas
flow rate and the rotational speed of the impeller are in
creased. At 155 RPM, all the predicted values are larger
and the trend for the experiments is not clearly observ
able.
At 220 RPM, there is good agreement with the experi
mental data; the trend for the Sauter diameter is clearly
observable: dR12 > dR4 > dR8 and dR2 > dRg, with
the increase of the distance between dR12 and dR9. For
this rotational speed and gas flow rate it is clear that the
bubbles with the smallest size stay on the impeller axis,
as the lower recirculation zone entraps the bubbles with
the largest size.
At 250 RPM, the general qualitative agreement is very
good. For 0.052 vvm, the plane near the impeller is not
well described: in fact slightly dR9 > dR2. Whereas for
0.093 vvm, the experimental trend is well represented
by the simulations, as confirmed by the progressive in
crease of coalescence near the upper surface.
Let us now discuss the results concerning mass trans
fer. In Fig. 6 the contour plot of the local volumetric
mass transfer coefficient, kL a, is reported over the whole
geometry. This variable is strongly inhomogeneous, in
fact the larger values are on the impeller plane and near
the buffles, where breakage caused by turbulence signif
icantly increases the number of bubbles. In these zones
there is the largest interfacial area per unit of volume a
and the higher mass transfer coefficient kL. On the con
trary, in the bottom part of reactor, where no bubbles are
present, the mass transfer coefficient nearly vanishes. It
is important to remind that most of these effects could
not be appreciated without using a population balance
model. Moreover with this approach not only the size
and the velocity distributions are correctly described, but
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
S2.15e+00
1A2e+00
0.71e01
Conclusions
0.00e+00
Figure 6: Spatial distribution of the local mass transfer
coefficient kLa (m s 1) for 250 rpm 0.093 vvm.
* Experimental Data [220 RPM 0,041 vm]
* Experimental Data [250 RPM 0,052 vm]
 CFD prediction [220 RPM 0,041 wm]
 CFD prediction [250 RPM 0,052 wm]
600
Time (s)
1200
Figure 7: Comparision between CDF predictions and
experimental data of oxygen concentration YL (mol
m 3) considering a monodisperse bubbles distribution
kLa (m s 1) for different operative conditions.
the composition distribution is also accounted for. In
fact, by solving Eqs. 39 and 40, the method is capable
of modelling the fact that smaller bubbles are charac
terized by higher mass transfer coefficients than bigger
ones.
Since the code still suffers of some numerical issues
that jeopardize its stability, in Fig. 7 only some prelim
inar results of two mass transfer tests for a monodis
perse bubble distribution are reported. These experi
ments were carried out fluxing into the stirred tank ni
trogen and then switching to oxygen. The time evolution
of species concentration in the liquid was tracked for a
single sample point (R4) (Laakkonen et al. 2006, 2007).
As it is seen there is indeed room for improvements.
Further steps of this work will focus on the solution of
all the numerical problems, related to the CFD imple
mentation.
In this work, a multivariate population balance for poly
disperse gasliquid systems was formulated. DQMOM
has been used to solve the PBM for a homogeneous sys
tem to consider only the effect of coalescence and break
age on the bubbles. Tests conducted on this simple ge
ometry show the consistency of the model in the case of
pure coalescence or breakage as well as combined coa
lescence and breakage problem.
Then the equations were rewritten in order to be cou
pled with a commercial CFD code to simulate the evolu
tion of the NDF for a real geometry, for which exper
imental data concerning local BSD and mass transfer
were available. The comparison between CFD predic
tions and experimental data for the real system indicates
generally good agreement, even using only two nodes
for the quadrature approximation. A first mass transfer
analysis was successful carried out but some numerical
issues did not allow a complete model validation. The
next steps of our work will focus on these issues and on
the final model validation.
References
Andersson R., Andersson B., On the breakup of fluid
particles in turbulent flows, AIChE Journal, Vol. 52, pp.
20202030, 2006
Calderbank P.H., Physical rate processes in industrial
fermentation, Part 1: The interfacial area in gasliquid
contacting with mechanical agitation, Trans. Inst. Chem.
Engrs., Vol. 36, pp. 443463, 1958
Clift R., Grace J.R., Weber M.E., Bubbles, Drops and
Particle, Academic Press, London, 1978
Coulaloglou C.A., Tavlarides L.L., Description of in
teraction process in agitated liquidliquid dispersion,
Chemical Engineering Science, Vol. 32, pp. 12891297,
1977
Fan R., Marchisio D.L., Fox R.O., Application of the di
rect quadrature method of moments to polydisperse gas
solid fluidized beds, Powder Technology, Vol. 139, pp.
720, 2004
Gordon R.G., Error bounds in equilibrium statistical me
chanics, Journal of Mathematical Physics, Vol. 9, pp.
655667, 1968
Kataoka I., Serizawa A., Basic equations of turbulence
in gasliquid twophase flow, International Journal of
Multiphase Flow, Vol. 15, pp. 843855, 1989
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Kazakis N.A., Mouza A.A., Paras S.V., Experimental
study of bubble formation at metal porous spargers: ef
fect of liquid properties and sparger characteristics on
the initial bubble size distribution, Chemical Engineer
ing Journal, Vol. 137, pp. 265281, 2008
Laakkonen M., Alopaeus V., Aittamaa J., Validation of
bubble breakage, coalescence and mass transfer mod
els for gasliquid dispersion in agitated vessel, Chemical
Engineering Science, Vol. 61, pp. 218228, 2006
Laakkonen M., Moilanen P., Alopaeus V., Aittamaa J.,
Modelling local bubble size distributions in agitated ves
sels, Chemical Engineering Science, Vol. 62, pp. 721
740, 2007
Lamont J.C., Scott D.S., An eddy cell model of mass
transfer into surface of a turbulent liquid, AIChE Jour
nal, Vol 16, pp. 513519, 1970
Marchisio D.L., Virgil R.D., Fox R.O., Quadrature
method of moments for aggregationbreakage process,
Journal of Colloidal and Interface Science, Vol. 258, pp.
322334, 2003
Marchisio D.L., Fox R.O., Solution of population bal
ance equations using the direct quadrature method of
moments, Journal of Aerosol Science, Vol. 36, pp. 43
73, 2005
McGraw R., Description of aerosol dynamics by the
quadrature methods of moments, Aerosol Science &
Technology, Vol. 27, pp. 255265, 1997
Petitti M., Nasuti A., Marchisio D.L., Vanni M., Baldi
G., Mancini N., Podenzani F., Bubble size distribution
in stirred gasliquid reactors with QMOM augmented by
a new correction algorithm, AIChE Journal, Vol. 56, pp.
3653, 2010
Prince M.J., Blanch H.W., Bubble coalescence and
breakup in airsparged bubble columns, AIChE Journal,
Vol. 36, pp. 14851499, 1990
Ranade V.V., Computational flow modelling for chem
ical reactor engineering, Academic Press, San Diego,
USA, 2002
Vanni M., Approximate population balance equations
for aggregationbreakage processes, Journal of Colloid
and Interface Science, Vol. 221, pp. 143160, 2000
Wallis G.B., One Dimensional Twophase Flow, Mc
Graw Hill, New York, 1969
Zucca A., Marchisio D.L., Vanni M., Barresi A.A., Val
idation of bivariate DQMOM for nanoparticle processes
simulation, AIChE Journal, Vol. 53, pp. 918931, 2007
