- Email: [email protected]
Toward a unified framework for the derivation of breakage functions based on the statistical theory of turbulence
Chemical Engineering Science 60 (2005) 6584 – 6595 www.elsevier.com/locate/ces
Toward a unified framework for the derivation of breakage functions based on the statistical theory of turbulence M. Kostogloua, b,∗ , A.J. Karabelasb, c a Department of Chemical Technology, School of Chemistry, Aristotle University, University Box 116, 541 24 Thessaloniki, Greece b Chemical Process Engineering Research Institute, P.O. Box 1517, 54006 Thessaloniki, Greece c Department of Chemical Engineering, University Box 455, Aristotle University of Thessaloniki, GR 54124 Thessaloniki, Greece
Received 31 March 2005; received in revised form 19 May 2005; accepted 24 May 2005 Available online 11 July 2005
Abstract The evolution of droplet or bubble size distribution in turbulent flow is of great significance in a variety of technological fields. Modeling this evolution by employing a population balance approach requires knowledge of the so-called breakage functions (rate and kernel). Over the years, a large number of phenomenological breakage functions with various degrees of sophistication have been proposed in the literature. Among them, those based on the statistical theory of turbulence are of particular interest, in that they attempt to take into account the structure of the flow field responsible for breakage. The purpose of the present work is to present a unified framework for developing this type of breakage functions and to show how existing models can be derived in a systematic and consistent way. The key parameters in this modeling approach are identified, which have to be determined by comparison with experimental data. It is also shown that the breakage functions, obtained within the framework presented here, lead to predictions of a droplet size evolution whose main features are consistent with experimental observations. It is suggested that this framework is an important step toward the development of a standard approach for modeling droplet size evolution in turbulent flow. 䉷 2005 Elsevier Ltd. All rights reserved. Keywords: Bubble size distribution; Mathematical modeling; Breakage process; Statistical approach to turbulence
1. Introduction Predicting the evolution of size distribution of fluid particles (bubble/droplets) in a turbulent flow field is one of the key problems of chemical engineering science. Knowledge of the size distribution is critical for the design of a large variety of equipment, including two phase-flow pipelines, bubble columns (for chemical, biochemical or other physical processes, e.g. Dudukovic et al., 1999), liquid–liquid dispersion in columns (e.g extraction; Attarakih et al., 2004a) or agitated vessels (e.g. emulsion polymerization; Chatzi et al., 1989). The size distribution of fluid particles is generally the outcome of two competing mechanisms, i.e., those ∗ Corresponding author. Department of Chemical Technology, School of Chemistry, Aristotle University, University Box 116, 541 24 Thessaloniki, Greece. Tel.: +30 2310997767. E-mail address: [email protected] (M. Kostoglou).
0009-2509/ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2005.05.051
of coalescence and breakage. In many cases coalescence is not important either naturally (e.g. for small concentration of the dispersed phase) or by design (e.g. if coalescence inhibition is applied). In this study the case of dilute dispersions is considered where only breakage is significant. The first attempts to model breakage focused on obtaining a relationship of the maximum particle size (capable of surviving in a turbulent flow field) with the physicochemical and hydrodynamic parameters of the system (Hinze, 1955). Later the dynamic nature of the breakage phenomenon was dealt with through the solution of an appropriate kinetic equation, the so-called breakage population balance (Valentas and Amudson, 1966). However, to solve such an equation for the evolution of the size distribution, the breakage functions (rate and kernel) must be known. The literature on efforts to determine breakage functions, over the past 40 years, is very large. Among the several ways proposed for obtaining
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
the breakage functions, i.e., direct observation of breakage events, direct numerical simulation, “inverse” problem solution (by employing experimental data), purely algorithmic functions and phenomenological models, the latter seems to be the most promising one. The main advantage of phenomenological models is their explicit dependence on the system physicochemical parameters. On the contrary, the frequently used algorithmic models do not include an explicit dependence on the system state but only fitting parameters. Furthermore, regarding the “inverse” breakage problem, it was recently shown (Kostoglou and Karabelas, 2005) that it is an ill-posed one (because of multiplicity of solutions) requiring a restricted space of candidate breakage functions in order to obtain a unique solution. This space may be defined by means of a phenomenological model. There are a large number of phenomenological models for breakage functions in the literature. However, great uncertainties exist regarding the reliability of such models; for instance, similar physical arguments used by various authors lead to quite different forms of the breakage functions. Additionally, breakage functions (phenomenological and algorithmic) of quite different form often seem to achieve an equally good fit of experimental size distributions (e.g. Luo et al., 2004; Chen et al., 1998; Ruiz and Padilla, 2004; Hesketh et al., 1991; Martinez-Bazan et al., 1999a,b). In the authors’ opinion, despite significant progress made in recent years, the state of the art in the development of breakage functions is unsatisfactory. There is little doubt that the best approach to follow is the formulation of a generalized (phenomenological) framework for developing breakage functions, which should be tuned and validated by carefully obtained (new and already available) experimental size distributions over a broad range of parameter values. An important step in this direction has been made by Luo and Svendsen (1996), who related the breakage functions with statistical quantities of homogeneous turbulence in a systematic and detailed way. Despite their frequent use, the final forms of the breakage functions obtained by these authors have some drawbacks; thus, several attempts to improve them, based on the same framework, have recently appeared in the literature (Lehr and Mewes, 2001; Lehr et al., 2002; Hagesaether et al., 2002a; Wang et al., 2003). However, it will be noted that among the above works there is no detailed study of the implications of the developed breakage functions on the maximum particle size dmax and on the evolving particle size distribution as well as on their relation with the system parameters. These significant issues are addressed in this paper. In the present work, a unified framework is proposed for the derivation of breakage functions based on the statistical description of isotropic turbulence. Models of this class reported in the literature are reviewed and critically evaluated. Additionally, the implications of developing particular types of breakage functions for the particle size distribution are presented and discussed. The structure of the paper is as follows: In the next section, the breakage population bal-
6585
ance is outlined. A unified framework for the derivation of breakage functions is presented in detail next. The breakage equation is then discretized in a specific way to conform to the particular form of the breakage functions. Finally, fairly extensive results, for the breakage rate and evolution of particle size distribution, are presented and discussed. 2. Problem formulation 2.1. The breakage equation The breakage process can be described in general by the following linear partial integrodifferential equation: ∞ jf (x, t) p(x, y)b(y)f (y, t) dy − b(x)f (x, t), (1) = jt x where t is the time, x the particle volume, f (x, t) the particle number density distribution, b(x) the breakage frequency (the term rate is also customarily used) and p(x, y) the probability distribution of particles of volume x resulting from the breakup of a particle of volume y. For a given initial particle size distribution f0 (x) = f (x, 0), Eq. (1) can be solved for f (x, t). The function p(x, y) should satisfy the following conditions (Ziff and McGrady, 1987): (a) Conservation of mass: y xp(x, y) dx = y. (2) 0
This equation stipulates that the total volume of particles resulting from the breakup of a particle of volume y must be equal to y. (b) z y y xp(x, y) dx (y − x)p(x, y) dx where z < 2 0 y−z (3) For binary fragmentation (two fragments per parent particle), the above restriction is simplified, being equivalent to a symmetric kernel in the sense p(x, y) = p(y − x, y). The number of particles resulting from breakage of a single particle of volume y is given as y v(y) = p(x, y) dx. (4) 0
To be able to use the fluid particle diameter d as the characteristic size, instead of particle volume, the following relations are introduced (assuming spherical particles):
3 2 f (x, t) = 2 fd (d, t). (5) d 6 d Further, a dimensionless fragment size probability function is employed: d 3 d 3 d 3 P (fv , d) = (6) p fv , with fv = x/y. 6 6 6
x=
6586
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
2.2. Particle breakage functions based on statistical approaches to isotropic turbulence 2.2.1. General framework assumptions According to the classical statistical approach to isotropic turbulence (see Baldyga and Bourne, 1999, for a comprehensive review), the flow is represented by a population of eddies exhibiting certain dynamics. The latter may be assumed very fast in comparison with the droplet dynamics, so the equilibrium properties of the eddy population may be considered. Each eddy is characterized by its size (diameter) and its kinetic energy e. The generalized (with respect to both eddy size and energy) number density of eddies is n Pe () where 12 n d is the concentration of eddies of e size between 1 and 2 , and e12 Pe de represents the fraction of eddies of size having energies between e1 and e2 . ¯ )) exp(−e/e( ¯ )) where e( ¯ ) is the Typically, Pe = (1/e( mean kinetic energy of eddies of size , depending on flow conditions. The eddies of size are considered to “collide” with the droplets of size d, with a frequency which can be determined by using concepts from the kinetic theory of gases; i.e.,
(d) = Ccs u n .
(7)
This frequency is the product of the following quantities: • Particle–eddy collision cross section, Ccs . The classical relation Ccs = (/4)(d + )2 is used here but alternatives exist; e.g., Ccs = d is used by Lehr and Mewes (2001). • Relative particle–eddy velocity. The particle is assumed “not moving”, so the relative velocity is equal to the mean velocity of eddies of size , u ; estimates of the latter are obtained from the theory of isotropic turbulence. Some authors (Hagesaether et al., 2002b) have added to this velocity the bubble velocity due to buoyancy. Such bubble motion leads to scavenging of eddies and increased breakage rate, but its importance has not been assessed quantitatively. • Number density (with respect to eddy size) n of eddies of size . Only the eddies in the inertial sub-range are considered to be effective for particle breakage. Expressions for n , u , Pe () will not be given here since they can be found elsewhere (e.g. Luo and Svendsen, 1996; Wang et al., 2003). The collision frequency does not depend on the eddy energy but the outcome of the collision does. Eddies having a sufficiently high energy level can cause breakage of droplets upon collision. The fraction of such collisions leading to breakage is called collision efficiency. Before proceeding in the quantitative analysis, the following four assumptions of the model should be stated: (i) The presence of droplets does not influence the equilibrium properties of the eddy population. A formalism for the relaxation of this assumption has been developed by Jairazbhoy and Tavlarides (2000).
(ii) Successful collisions lead to two daughter particles (i.e., binary breakage). This is a major assumption, with experimental support from detailed observations of gas bubble breakage events. The situation may be different for liquid droplets for which the experimental distributions are usually fitted assuming breakage in more than two fragments (e.g. Chatzi et al., 1989). The extension of the present approach to multiple breakage is not a trivial task, necessitating the introduction of additional parameters (e.g. number of fragments). (iii) The droplets always remain in the same state. In particular, it is assumed that the energy, absorbed from eddies incapable of causing breakage, is dissipated instantaneously (before the next collision event) so that there are no synergistic effects between successive collision events. The alternative type of transient particle behavior may be considered in the case of breakage of water droplets in high-velocity jets and its implementation requires dynamic modeling at the particle level, rendering inadequate the mean field theory examined here. (iv) As the viscosity of the dispersed phase increases, the influence of the viscous dissipation of the eddy energy on the breakage rate will increase as well. Although the effect of viscosity was taken into account by older correlations for dmax and by phenomenological models for breakage rate (see Attarakih et al., 2004b), it is not considered by the more recent models based on statistical theory of turbulence. The latter models are restricted to cases of negligible dispersed phase viscosity. Obviously, an extension of the theory to include very viscous droplets is needed. The above major assumptions are cited by all the authors using the statistical approach to turbulent breakage in fluid–fluid dispersions. The following additional working assumptions are flexible enough, providing a lot of opportunities to the developer of a particular model; these are essentially restrictions regarding the size of the eddy colliding with the parent droplet and the size of the (smaller of the 1/3 two) daughter droplet df = fv d: (a) A droplet–eddy collision can be effective only if the eddy size is smaller than the droplet size d. In general, this criterion has some physical basis (see Wang et al., 2003), but the precise upper limit set to the active eddy size is rather arbitrary. It will be pointed out, in this respect, that the breakage rate is very sensitive to the upper limit of active eddy size since large eddies contain large amounts of energy and are very efficient in breaking up droplets. Hagesaether et al. (2002a) ignored completely this criterion assuming that the breakage efficiency depends only on the energy and not on the size of the eddy. Generalizing this approach, the upper active eddy size may be considered as a parameter (max ) taking values from d to infinity and requiring further research for its more precise specification. (b) Criterion setting an upper limit to the smaller of the two daughter particle sizes (fv max ) (surface energy criterion): The energy content of the colliding eddy must be greater than that corresponding to new surface created by
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
6587
the breakage event. The form of this criterion is unique and there are no alternative interpretations. The maximum size fv max of the smaller daughter particle, resulting from the collision of a parent particle of size d with an eddy of energy e, is obtained from the solution of the equation fv2/3 max + (1 − fv
2/3 max )
= min
e emax , , d 2 d 2
(8)
where is the interfacial tension between the two fluids. The value emax = (21/3 − 1)d 2 is the maximum amount of energy which can be absorbed by the fluid particle during a breakage event and corresponds to two equal fragments. Thus, any eddy with e > emax will lead to two equal fragments. (c) Criterion setting a lower limit to daughter particle size (fv min ). Several forms of this criterion appear in the literature. According to the capillary pressure formulation (Wang et al., 2003), the dynamic pressure of the turbulent eddy c u2 /2 must be larger than the capillary pressure of the smaller fragment of breakup. The value used for the capillary pressure by Wang et al. (2003) is /df where df is the diameter of the smaller fragment. These authors state that the capillary pressure is /r, where r is a characteristic radius of curvature but it is not clear how the df used in their derivation is related to a radius of curvature. A more obvious value of the capillary pressure 2/df is suggested by Lehr et al. (2002) who assume that the fragment has a nearly cylindrical shape at the moment of breakup. A different formulation of the criterion is developed by Hagesaether et al. (2002a), who consider that the energy density of the turbulent eddy must be larger than that of the smaller fragment. Quantitatively, this means that c u2 /2 6/df . Generalizing, one may consider that the smaller fragment size is determined by the requirement c u2 /2 /df , where is a numerical parameter; this inequality leads to fv
min
=
3 6ed
3 .
(9)
In addition to this criterion for the smaller fragment size, Lehr et al. (2002) state that df must not be greater than the colliding eddy (i.e., df ). 2.2.2. Computation of breakage functions Fig. 1 is useful for depicting the above criteria. In this graph the curves fv max and fv min are shown versus energy e for a particular eddy and parent particle size and d, respectively. The value ecrit can be obtained from the solution of equation fv min (ecrit ) = fv max (ecrit ) and clearly denotes the minimum eddy energy which can lead to a breakage event. Finally, e1 = max(fv−1min (fv1 ), fv−1max (fv1 )) is the smallest eddy energy leading to a fragment of fractional volume fv1 . From the above analysis it is clear that the breakage efficiency Pb (, d) of collisions between particles of diameter d
Fig. 1. Curves fv max (e), fv min (e) and location of points (e1 , fv1 ), (ecrit , fv crit ), (emax , 21 ).
and eddies of size is ∞ ecrit (, d) Pb (, d) = Pe () de = exp − . e( ¯ ) ecrit (,d)
(10)
The total breakage rate is obtained by integrating the product of the collision rate and breakage efficiency for all eddy sizes ∞ b(d) = Pb (, d) (d) d 0 max ( + d)2 = 0.923(1 − )1/3 11/3 min ecrit (, d) × exp − (11) d , e( ¯ ) where is the volume fraction of the dispersed phase and the turbulent energy dissipation rate. The lower limit of the active eddy size is taken to be proportional to the Kolmogorov microscale of turbulence with a proportionality constant between 11.4 and 31.4 (Wang et al., 2003) in order to maintain compatibility with the used eddy size distribution which corresponds to the inertial subrange of turbulence. It is noted that, unlike the case of max , the precise value of min is of reduced significance for the computed breakage rate since small eddies have very small breakage efficiencies. In the case where only the surface energy criterion is considered, ecrit = 0; i.e., all the eddies are effective in breaking up droplets. Very small eddies give one very small fragment and another close to the parent particle. In the case where only the capillary pressure criterion is considered ecrit is obtained from the solution of fv min (ecrit )= 21 . The computation of the function P (fv , d) is much more involved than that of the rate b(d). The fundamental quantity for the computation of the general fragment size distribution is the probability distribution D(fv , e, , d) representing the probability of having a fragment of size fv from the successful collision of an eddy (with properties , e) with
6588
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
a particle of size d. The function D is non-zero only for fv min < fv < fv max and obeys the following normalization condition: fv max D(fv , e, , d) dfv = 1. (12) fv
min
It is convenient to write the function D in the following form: D(fv , e, , d) = D1 (fv ; e, , d)(U (fv − U (fv min )),
max )
(13)
where U is the Heaviside step function. The advantage of the new function D1 is that only explicit influences of e, , d are included in the argument, whereas the implicit ones through fv min and fv max are omitted. In general, the shape of D1 can differ depending on the values of e, , d, but up to now the usual assumption is that the shape of D1 with respect to the variable s = (fv − fv min )/(fv max − fv min ) is constant (i.e., D1 (fv ; e, , d) = D1 (s)). A notable exception can be found in Hagesaether et al. (2002a); these authors consider that the probability for a particular breakage event is larger as the differences of the eddy energy and energy density from those needed for the particular event are larger; i.e., D(fv , e, , d) ∝ [e − (fv2/3 + (1 − fv )2/3 − 1)
6e 6 2 × d ] . − 3 df 1/3 v
(fv − F ) = (e − F
Any shape of D1 (fv ) is possible but the usual choices are the uniform one (equal probability for all permitted fragment sizes) (Wang et al., 2003) and the Dirac delta function corresponding to the capillarity criterion fv min (Lehr et al., 2002). Having chosen the form of the function D, the fragment distribution function Pf (fv , , d) for a collision event between an eddy of size and a particle of size d is developed through the integration of D with respect to energy spectrum of the eddies: ∞ Pf (fv , , d) = D(fv , e, , d)Pe () de. (14) 0
The rate of production of fragments of size fv is given as max Pf (fv , , d) (d) d. (15) bf (fv , d) = min
The fragment size distribution can be obtained from P (fv , d) =
bf (fv , d) . b(d)
−1
dF (fv )) de
−1 ,
(18)
where the notation F −1 (fv ) stands for the solution of the equation fv = F (e) ⇒ e = F −1 (fv ). Substituting in Eq. (14) leads to Pf (fv , , d) = PF −1 (fv )
dF de
−1 e=F −1 (fv )
.
(19a)
This relation cannot be used in the case of F independent of e because the slope of the F (e) curve becomes zero; a typical example is the case F = fmax where F = 21 for e > emax . Eq. (19a) is then transformed to Pf (fv , , d) = PF −1 (fv )
dF de
−1
+ (fv − 1/2)
e=F −1 (fv ) ∞ emax
Pe () de.
(19b)
Case (ii): Uniform probability distribution (i.e., D1 = 1/(fv max − fv min )). In this case Eq. (14) can be written as
(16)
The total breakage rate b(d) can be obtained as an alternative to Eq. (11) by integration of the rate bf over all the fragment sizes, i.e., 1/2 bf (fv , d) dfv . (17) b(d) = 0
Although the use of Eq. (17) could be used for the computation of the denominator in Eq. (16) (i.e., a normalized P (fv , d) is needed), it is not recommended for the computation of the breakage rate b(d), which alternatively can be made through Eq. (11). The reason is that Eq. (17) requires a complicated computational procedure, including the use of function D, which is susceptible to computational error. However, this source of error can be avoided (bypassed through Eq. (11)) since b(d) does not really depend on D. The function P (fv , d) results from a double integral which in general is computed numerically. This computational load can be significantly reduced for some simple forms of the function D1 used in the literature; a summary of such cases follows. Case (i): D1 (fv ) = (fv − F ), where F is a given combination of fv max and fv min . For example, F = fv max means that the particle breaks into two fragments obeying exactly the surface energy criterion; F = fv min means that the smaller fragment obeys exactly the capillarity criterion; F = (fv min + fv max )/2 means that the size of fragment is the average of fv min and fv max . F is a function of energy e through the e-dependence of fv max and fv min . From the properties of Dirac delta functions one obtains
Pf (fv , , d) =
∞ max[fv−1min (fv ),fv−1max (fv )]
× Pe () de.
fv
1 − fv max
min
(20)
It is observed that the integrand does not depend on fv . The dependence on fv is only through the lower integration limit, so a sequential computation of the function Pf at the points fvi = i/(2M) for i = 1, 2, 3, . . . , M is possible, as follows:
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
For fvi > fv crit Pf (fvi , , d) = Pf (fvi+1 , , d) +
bf (fv , d) =
fv−1max (fvi+1 )
fv−1max (fvi ) × (fv max − fv min )−1 Pe () de
(21a)
with Pf (fvM , , d) =
∞ emax
1 − fv 2
−1 min
Pe () de.
(21b)
For fvi < fv crit Pf (fvi , , d) = Pf (fvi−1 , , d) +
fv−1min (fvi )
fv−1min (fvi−1 ) × (fv max − fv min )−1 Pe () de
with Pf (fv1 , , d) =
∞ fv−1min (fv1 )
(fv
× Pe () de.
max
− fv
min )
(22a)
−1
(22b)
The integrals in (21b), (22b) include an exponential term in Pe () and with a simple coordinate transformation can be cast in a form appropriate for the use of the Gauss–Laguerre quadrature. The rest of the integrals are computed using the orthogonal rule which is more accurate than the trapezoidal one (Pozrikidis, 1998). This procedure, proposed by Wang et al. (2004), greatly facilitates the computation of the breakage kernel since it reduces the order of integration from two to one. It is noted that Wang et al. (2003, 2004) have assumed that no breakage occurs for fv max − fv min < , where is a parameter having a small value. It is believed that this addition has been made for computational purposes but it tends to reduce the breakage rate since ecrit is now obtained from the solution of fv max (ecrit ) − fv min (ecrit ) = . Providing a unified framework based on the statistical theory of turbulence, which will encompass breakage models already proposed in the literature, would greatly facilitate further developments and the correct derivation of additional models. However, as a counter example of incorrect application of this approach, one may refer to the breakage functions derived in the pioneering work of Luo and Svendsen (1996). The key assumption in their derivation is that only the surface energy criterion is obeyed and all the fragment sizes have a uniform probability between 0 and fv max . According to the present analysis, making use of this assumption leads to the following breakage function: max ( + d)2 b(d) = 0.923(1 − )1/3 d 11/3 min 0.923(1 − )1/3 −2/3 = [2(min − −2/3 max ) 3 −5/3 −8/3 2 −8/3 + 10d(min − −5/3 max ) + 8d (min − max )], (23a)
max min
∞ 0
1 fv
max
Pe () de (d) d.
6589
(23b)
The above rate is the maximum one and corresponds to breakage efficiency equal to 1. The derivation of Luo and Svendsen (1996) is incorrectly based on the concept of conditional breakage probability for a given fragment size fv , overlooking that eddies incapable of producing fragments of size fv can be efficient in producing fragments of smaller size. This inconsistency is recognized by Wang et al. (2003) only as a drawback of the model and not as an error in the derivation. Due to this erroneous derivation the model of Luo and Svendsen (1996) (having collision efficiency 1) leads to smaller breakage rate compared to that of Wang et al. (2003) which has a collision efficiency smaller than 1. Although the physical assumptions of the model of Luo and Svendsen (1996) have been extensively criticized and relaxed, the aforementioned error has been essentially overlooked; thus even today the model (with no amendment) is repeatedly used for simulation and comparison with experiments (e.g. Chen et al., 2005; Dhanasekharan et al., 2005; Yeoh and Tu, 2005). 2.3. Numerical solution of the breakage equation Unlike the case of solid fragmentation, the extent of breakage in fluid–fluid dispersions is not very large; thus, the usual finite difference methods (in the particle diameter domain) are extensively used in the literature for the discretization of Eq. (1) (Ribeiro et al., 1995). Recently, for the solution of Eq. (1), Hagesaether et al. (2002b) developed an alternative method. The latter has the drawback that it permits only a geometrical grid with a ratio 2 with respect to particle volume. The need for a finer and arbitrary grid renders the sectional approach of Kumar and Ramkrishna (1996) the preferred choice, even for fluid–fluid applications. According to Kumar and Ramkrishna (1996) the PSD can be represented by a sum of discrete particle sizes xi with number concentration Ni (i = 1, 2, 3, . . . , L), or mathematically f (x, t) =
L
Ni (x − xi ).
(24)
i=1
Consequently, the following system of equations must be solved for the particle concentration Ni : L
dNi nik b(xk )Nk − b(xi )Ni , = dt
(25)
k=i
ni,k = (1 − ik )Z(xi , xi+1 , xk ) − (1 − 1k ) × Z(xi , xi−1 , xk ), a2 a2 − y Z(a1 , a2 , a3 ) = p(y, a3 ) dy. a 2 − a1 a1
(26) (27)
6590
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
The moments are computed as follows: Mk =
L i=1
Ni xik .
(28)
The method is used in bubble breakup applications with geometrical grid of ratio R = 2 (Chen et al., 2005), R = 1.7 and R = 1.45 (Wang et al., 2005). This choice is appropriate for the case of combined breakage–coalescence in which the particle size distribution does not evolve very much. However, in the case of pure breakage treated here, the size distribution continuously evolves becoming narrower, as will be shown in the next section. Thus, after some time, the geometric grid is not sufficiently dense to accurately describe the particle size distribution. The following grid construction procedure is proposed, combining high accuracy and computational efficiency: First, a particle volume x0 is determined for which the particles undergo insignificant breakage during the required simulation time (i.e., b(x0 )tfin = 0.001). Then, the region between the size of the largest initial particles xmax and x0 is discretized in L2 points using a geometrical grid with each interval equal to R times the previous one. No breakage is considered for particles of size smaller than x0 . The region from 0 to x0 is discretized using L1 equidistant points. This choice leads to a controlled resolution of the particle size distribution for t < tfin and has the computational advantage that the matrix Pi,j is not needed for j < L1 (L = L1 + L2 ). Also with this type of discretization the well-known finite domain error is of no concern. It is noted that in the method proposed by Hagesaether et al. (2002b) the resolution with respect to fv of the function P (fv , d) is dictated by the discretization of Eq. (1). However, in general the resolution required for P (fv , d) is much finer (and on different grid points) than that required for the sectional method; thus, the two stages of computation must be decomposed. Having chosen the grid points xi (i = 1, 2, . . . , L), the matrix Pij = P (fvi , dj ) (i = 1, 2, . . . , M, j =L1 , 2, . . . , L) is constructed. The integrals in Eq. (27) are determined using the trapezoidal rule at K points and the integrand is computed at arbitrary values of fv by linearly interpolating between the values of the matrix Pij . In the present work the values K =70, L1 =10, L2 =50, R = 1/0.85 are selected (after several trials) because they provide satisfactory accuracy for the cases treated here, and lead to practically grid-independent results. It is noted that the combination of a very coarse grid (e.g. the nine classes in Dhanasekharan et al., 2005) with very steep kernels (e.g. that of Luo and Svendsen, 1996), and coarse discretization of integrals in (27) (e.g. the five Gaussian points in Attarakih et al., 2004b) can lead to questionable results.
3. Results and discussion In this section the form of the breakage rate and the evolution of the droplet size distribution will be obtained for
Fig. 2. The function Pb (fv ) in the neighborhood of its singularity for d =6 mm, =4.2 mm, =1 m2 /s3 , =72 m N/m and for several numbers M of discretization points fvi .
several sets of parameters and will be discussed in some detail. The analysis will be restricted to the standard model presented by Wang et al. (2003); i.e., uniform breakage probability for fv min < fv < fv max and fv max − fv min > . The function Pb (fv , , d) is singular going to infinity at fv crit = fv min (ecrit ). A non-zero value of can remove the singularity and facilitate the computations but reduces significantly the breakage rate without having any physical meaning. Thus, cannot be considered to be a simple parameter for computations, and has no place in the model. Of course, the numerical computation of the singular function Pb is still a problem but it must be overcome using numerical means, without altering the physics of the model. Two important points will be stressed here. First, the breakage rate can be computed directly from Eq. (11) avoiding any relation with Pb ; second, the singularity of Pb with respect to fv is integrable (i.e., the integral of Pb in any finite interval around the singularity is finite). The only numerical requirement is the use of a very fine grid (large M value) with respect to fv in order to capture the shape of Pb in the region of singularity. To make this point explicit, the function Pb is shown in the region of singularity for the case d = 6 mm, = 4.2 mm, = 1 m2 /s3 and several values of M, in Fig. 2. A plot of Pb for the same conditions but with = 0.001 can be found in Wang et al. (2004). It is obvious that all the grids give the same result, except at the points close to the singularity. As the grid becomes more dense, on one hand the correct shape of Pb is followed closer, and on the other any discrepancy is restricted to a smaller domain, leading to a smaller overall error. In the present case the total breakage rate is obtained with three digits of accuracy for M = 1000. In order to have an accurate description of the breakage kernel, at least M = 200 points must be used for fv . The seemingly obvious suggestion of choosing a coarse grid and making it denser in the vicinity of fv crit is not really good since fv crit depends on ; so multiple grids must be used making the situation even more complicated. The significant effect of the artificial parameter (used for singularity masking) on the breakage rate is shown in
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
Fig. 3. Breakage rate b ( = 1 m2 /s3 , = 72 m N/m) computed by Eq. (17) versus parameter , for several values of particle diameter and discretization points M.
Fig. 3 ( = 1 m2 /s3 , = 72 m N/m). The breakage rate has been computed using the wrong approach (Eq. (17)) to show the effect of discretization. At values small enough with respect to fv discretization, there is no influence on the breakage rate; however, above a crossover value 0 there is a significant effect on b(d). Indeed, increasing leads to decreasing b(d), as expected. Using a finer grid, the crossover value 0 decreases. In the ideal limit of infinitely dense grid (or alternatively using the correct way, Eq. (11), for b(d)) 0 = 0. One can notice inFig. 3 that the breakage rate for d =5 mm and =0 is not the same for M =200 and M =400. As already mentioned, a value M = 1000 is needed for a converged breakage rate value. Evaluation of the breakage rate using Eq. (11) involves relatively small computational load. The trapezoidal rule is used for the numerical evaluation of the integral and the bisection method for ecrit computation from the transcendental equation fv min (e) = fv max (e). The breakage rate as a function of the droplet diameter for surface tension =25 m N/m and three values of is shown in Fig. 4. In the b(d) curves one can discern two parts: i.e., a rather smooth variation with d for large d values, and a very strong dependence on d for relatively small d values. The sharp part of the curve becomes even sharper with increasing . It is remarkable that reduction of the particle diameter by a few percentage points leads to reduction of the breakage rate by several orders of magnitude. An attempt to fit the b(d) curves with algebraic relations was unsuccessful. Also, the general form of the breakage rate given by Lehr et al. (2002) but with unknown parameter values was incapable of fitting the curvesshown in Fig. 4. This confirms the assertion that the much simpler expression suggested by Lehr et al. (2002) cannot be used in place of the approach examined here (see Wang et al., 2005). In Fig. 5, the breakage rate is depicted versus the interfacial tension of the dispersion for several sets of
6591
Fig. 4. Breakage rate b versus particle diameter d for = 25 m N/m and several values of turbulent energy dissipation rate .
Fig. 5. Breakage rate b versus interfacial tension for several values of particle diameter d and turbulent energy dissipation rate .
parameters and d. As increases the breakage rate decreases, as one would have expected. The relative reduction of b for the same increase of seems to be larger as the absolute value of b is smaller. Results for the breakage kernel will not be presented here since everything needed can be found in the recent literature (see Wang et al., 2003). It will be pointed out, though, that for high breakage rates the kernel is of “M” shape with the two peaks very close to fv = 0. As the breakage rate is reduced, the peaks are displaced inward and finally collapse to a single peak at fv = 21 . The simple kernel of Lehr et al. (2002) exhibits the same behavior. Extensive comparisons between the kernels based on the statistical approach to turbulence are presented by Wang et al. (2003). It appears that the key feature determining whether a kernel exhibits a “U” shape or a “U–M” transition (with d decreasing) is the existence or not of fragments of very small size(existence of minimum fragment size criterion).
6592
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
Fig. 6. Evolution of diameter dm (corresponding to mean particle volume) for = 50 m N/m and several values of and initial monodisperse size d0 (small time scale).
An issue of paramount importance is the effect of breakage functions, described above, on the evolution of droplet size distribution in a turbulent flow field. In the following simulations, aimed at assessing such effects, the initial droplet distribution is monodisperse with diameter d0 . The diameter dm (called hereafter mean diameter) is the one corresponding to the mean volume of the droplets xm = M1 /M0 (i.e., dm = (6xm /)1/3 ). The evolution of dm for = 50 m N/m and three values of is shown in Fig. 6; for each value of two different initial diameters are considered. It is very interesting that after some time the system “forgets” its initial condition and the mean size evolution curves collapse to one. This is shown in Fig. 6 for only two different initial conditions but it can be generalized for any initial condition. The time of collapse of the evolution curves into a single one depends on the initial condition. This means that the breakage rate exhibits a “self-regulating” tendency, leading to a long time (seemingly) asymptotic evolution curve, independent of the initial size distribution. By inspecting Fig. 6 one would be tempted to conclude that dm reaches a final constant value and that the approach to this “steady state” is faster for larger values of . This trend would be in accordance with the classical view of existence of a maximum stable particle size dmax . However, the picture obtained from Fig. 6 is quite deceptive. In Fig. 7 the same evolution curves are shown in a long time scale; obviously, what might be considered as “steady state” in Fig. 6 corresponds to a very slow “sliding” of dm to smaller values. During the past 15 years, some experimental studies on the size evolution of dispersions in turbulent flow have cast doubt on the existence of a dmax or equivalently of a true steady state droplet size distribution. Kostoglou and Karabelas (2001) tried to explain these results by arguing that dmax may or may not exist but, from a practical point of view, it does not matter since it is experimentally inaccessible due to the form of the breakage rate function. According to an explanation offered by Baldyga and Bourne (1999) dmax is a sliding function of time due to the intermittent nature of turbulence; i.e., as time passes the chance of a drop
Fig. 7. Evolution of diameter dm (corresponding to mean particle volume) for = 50 m N/m and several values of (large time scale).
Fig. 8. Evolution of diameter dm (corresponding to mean particle volume) for = 1 m2 /s3 and several values of .
to encounter a more violent turbulent burst increases leading to reduction of dmax . The phenomenological approach to breakage presented here appears to be compatible with both the above interpretations of experimental observations regarding a non-existent steady state. According to the present approach, although a true steady state may not really exist, the PSD evolution is very slow rendering its existence of no practical concern, as Kostoglou and Karabelas (2001) have already suggested. Furthermore, from the results of Figs. 6 and 7 one may consider that a pseudo-steady state exists depending on the time scale of observation. This seems to agree with the concept of dmax “sliding” with time (Baldyga and Bourne, 1999). The evolution of dm for = 1 m2 /s3 and several values of is shown in Fig. 8. As already pointed out, the curves for large values of time are independent of d0 ; also, as expected, the droplet size decreases with decreasing interfacial tension. An important question is what happens to the evolving particle size distribution; i.e., is there a long time asymptote for the size distribution regardless of the initial distribution, or is this the case only for the mean size dm but with a different distribution shape? The present modeling effort suggests that an asymptotic size distribution does exist. Size distributions of droplets independent of the initial
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
Fig. 9. Long time (t = 10 000 s) normalized particle size distributions fd (d)/M0 for = 1 m2 /s3 and several values of .
Fig. 10. Evolution of the normalized particle size distributions fd (d)/M0 for = 0.1 m2 /s3 , = 50 m N/m and d0 = 4 mm.
condition are shown in Fig. 9 for t =10 000 s, =1 m2 /s3 and three values of . The distributions hereafter are normalized with the corresponding total droplet number concentration M0 inorder to make their presentation feasible in the same graph. It is pointed out that the “hat-type” shape and the slight bimodality of the distributions are not artificial due to the grid used; similar results are obtained using finer grids. Additionally, the size distribution becomes narrower as the surface tension () decreases. The evolution of the droplet size distribution for the case d0 = 4 mm, = 50 m N/m and = 0.1 m2 /s3 is shown in Fig. 10. Initially a tri-modal PSD appears consisting of the unbroken initial particles and the two modes of fragments (U-shape breakage kernel for large particles). At later times, the PSD becomes unimodal. Surprisingly enough, at very large time values, a weak bimodality appears and tends to evolve. The origin of this bimodality is not obvious but it is believed to be due to the combination of the shape (bell shape toward equal size breakage) of the breakage kernel for small droplet sizes with the extremely large slopes of the b(d) curve (e.g. Fig. 4). This creates discrete accumulation modes in the droplet size domain, having extremely different breakage rates. This difference in breakage rate may restrict the number of modes to two. The smaller mode can-
6593
Fig. 11. Evolution of the normalized particle size distributions fd (d)/M0 for = 1 m2 /s3 , = 50 m N/m and d0 = 4 mm.
Fig. 12. The maximum particle diameter at large time (t =10 000 s) versus interfacial tension , for two values of .
not suffer significant breakage (and generate a third mode) before the extinction of the larger one. The corresponding PSD evolution, for the case = 1 m2 /s3 , is shown inFig. 11. At the initial stages of breakage the size distribution is bimodal; a discrete third mode corresponding to the large size fragments does not appear here. Otherwise, the evolution is quite similar to that of the previous cases, with a bimodality appearing at large time values. The maximum droplet size at the time of 10 000 s is shown in Fig. 12 versus surface tension for two values of . It is very interesting that both curves can be fitted with a power law having exponent 0.56. The latter is close to the 0.6 exponent in the well-known relation of Hinze (1955) for dmax . Similarly, the maximum droplet size at time 10 000 s is shown in Fig. 13 versus mean dissipation rate for = 50 and 70 m N/m. Here again, there is a perfect fit of both curves with an exponent −0.43 of a power function, as compared to the −0.40 exponent of the expression of Hinze (1955). It will be noted that these dmax dependencies in the Hinze expression have been confirmed by Karabelas (1978) for turbulent pipe flow of dilute liquid–liquid dispersions at small flow times (10–30 s) and mean values in therange 0.4–8.8 m2 /s3 . Thus, according to the present theory, whereas a true steady state value dmax may not really exist, the maximum droplet size at a specified (sufficiently long) time is independent of the initial size distribution and
6594
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
d D(fv , e, , d)
Fig. 13. The maximum particle diameter at large time (t =10 000 s) versus mean dissipation rate for two values of .
exhibits a parametric dependence very similar to that assigned to dmax by the classical approach. As an example to demonstrate the importance of the generalized framework presented here one may refer to the recent attempt by Chen et al. (2005) to simulate bubble columns using the original Luo and Svendsen (1996) breakage model; these authors state that a 10-fold increase of computed breakage rate is required in order to achieve agreement with the experimental data. However, as pointed out already, the correct computation of the Luo and Svendsen rate, with the use of a much larger than d effective eddy size max , leads to much larger breakage rates without abandoning the general framework.
4. Concluding remarks A generalized framework for the derivation of breakage functions based on the statistical theory of turbulence is presented here. It is shown how several recently proposed breakage functions emerge from this framework. Identified as key parameters for the derivation of breakage functions are the maximum eddy size active for breakage (max ), the minimum size of fragment resulting from a breakage event (i.e., parameter in Eq. (9)) and the function D1 (fv , e, , d) defining the fragment size distribution resulting from a particle–eddy collision. The evolution of particle size distribution based on a particular member of the breakage functions from the generalized framework (i.e., the one proposed by Wang et al., 2003) is examined in detail to show that the functions of the framework are compatible with experimental observations and with their theoretical interpretation about the bubble/droplet evolution in turbulent dispersions.
Ccs
Greek symbols
Notation b(x) bf (fv , d)
parent particle diameter dimensionless size probability function for fragments resulting from the collision of an eddy with size and energy e with a particle of size d D1 (fv ; e, , d) modified form of D (see Eq. (13)) e eddy energy ecrit the lowest energy of an eddy of size capable of breaking a particle of size d emax required energy for breakage in two equal fragments of a particle of size d e( ¯ ) mean kinetic energy of eddies of size f (x, t) particle number density distribution (volume based) fd (d, t) particle number density distribution (diameter based) fv ratio of fragment to parent particle volume fv max maximum fractional volume of the smaller of two fragments (see Eq. (8)) fv min minimum fractional volume of the fragment (due to capillarity criterion) fv crit fragment volume fraction corresponding to ecrit K number of discretization point for numerical evaluation of integrals in Eq. (27) L total number of classes for the numerical solution of the breakage equation M number of discretization points for fv in the interval [0, 0.5] Mk kth moment of the particle size distribution n number density of eddies with respect to eddy size p(x, y) fragment size probability function P (fv , d) dimensionless fragment size probability function Pe () probability density function of energy of eddies with size Pb (, d) efficiency of collisions between particles of diameter d and eddies of size Pf (fv , , d) dimensionless size probability function for fragments resulting from the collision of an eddy with size with a particle of size d u relative continuous phase particle velocity x, y particle volume
breakage frequency for particle of volume x frequency of production of fragments with volume fv d from breakage of particles of size d particle–eddy collision cross section
min , max c
volume fraction of the dispersed phase generalization parameter for the capillarity criterion turbulent energy dissipation rate eddy diameter minimum and maximum eddy sizes effective for particle breakage density of thecontinuous phase
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
(d)
interfacial tension between particles and continuous phase density with respect to of collision frequency of particles of size d with eddies of size
References Attarakih, M.M., Bart, H.J., Faqir, N.M., 2004a. Solution of the droplet breakage equation for interacting liquid–liquid dispersions: a conservative discretization approach. Chemical Engineering Science 59, 2547–2565. Attarakih, M.M., Bart, H.J., Faqir, N.M., 2004b. Numerical solution of the spatially distributed population balance equation describing the hydrodynamics of interacting liquid–liquid dispersions. Chemical Engineering Science 59, 2567–2592. Baldyga, J., Bourne, J.R., 1999. Turbulent Mixing and Chemical Reactions. Wiley, New York. Chatzi, E.G., Gavrielides, A.D., Kiparissides, C., 1989. Generalized model for prediction of the steady state drop size distributions in batch stirred vessels. Industrial Engineering Chemistry Research 28, 1704–1711. Chen, Z., Pruss, J., Warnecke, H.J., 1998. A population balance model for disperse systems: drop size distribution in emulsion. Chemical Engineering Science 53, 1059–1066. Chen, P., Sanyal, J., Dudukovic, M.P., 2005. Numerical simulation of bubble columns flows: effect of different breakup and coalescence closures. Chemical Engineering Science 60, 1085–1101. Dhanasekharan, K.M., Sanyal, J., Jain, A., Haidari, A., 2005. A generalized approach to model oxygen transfer in bioreactors using population balances and computational fluid dynamics. Chemical Engineering Science 60, 213–218. Dudukovic, M.P., Larachi, F., Mills, P.L., 1999. Multiphase reactors—revisited. Chemical Engineering Science 54, 1975–1995. Hagesaether, L., Jakobsen, H.A., Svendsen, H.F., 2002a. A model for turbulent binary breakup of dispersed fluid particles. Chemical Engineering Science 57, 3251–3267. Hagesaether, L., Jakobsen, H.A., Svendsen, H.F., 2002b. Modeling of the dispersed-phase size distribution in bubble columns. Industrial Engineering Chemistry Research 41, 2560–2570. Hesketh, R.P., Etchells, A.W., Russell, T.W.F., 1991. Bubble breakage in pipeline flow. Chemical Engineering Science 46, 1–9. Hinze, J.O., 1955. Fundamentals of the hydrodynamic mechanics of splitting in dispersion processes. A.I.Ch.E. Journal 1, 289–295. Jairazbhoy, V., Tavlarides, L.L., 2000. A numerical technique for the solution of integrodifferential equations arising from balances over populations of drops in turbulent flows. Computers and Chemical Engineering 23, 1725–1735. Karabelas, A.J., 1978. Droplet size spectra generated in turbulent pipe-flow of dilute liquid–liquid dispersions. A.I.Ch.E. Journal 24, 170–180.
6595
Kostoglou, M., Karabelas, A.J., 2001. A contribution towards predicting the evolution of droplet size distribution in flowing dilute liquid/liquid dispersions. Chemical Engineering Science 56, 4283–4292. Kostoglou M., Karabelas A.J., 2005. On the-self similar solution of fragmentation equation: Numerical evaluation with implication for the inverse problem. Journal of Colloid and Interface Science 284, 571–581. Kumar, S., Ramkrishna, D., 1996. On the solution of population balance equations by discretization—I. A fixed pivot technique. Chemical Engineering Science 51, 1311–1332. Lehr, F., Mewes, D., 2001. A transport equation for the interfacial area density applied to bubble columns. Chemical Engineering Science 56, 1159–1166. Lehr, F., Mewes, D., Millies, M., 2002. Bubble-size distributions and flow fields in bubble columns. A.I.Ch.E. Journal 48, 2426–2443. Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. A.I.Ch.E. Journal 42, 1225–1233. Luo, G.S., Li, H.B., Tang, X.J., Wang, J.D., 2004. Drop breakage in a coalescence–dispersion pulsed-sieve plate extraction column. Chemical Engineering Journal 102, 185–191. Martinez-Bazan, C., Montanes, J.L., Lasheras, J.C., 1999a. On the breakup of an air bubble injected into a fully developed turbulent flow. Part 1. Breakup frequency. Journal of Fluid Mechanics 401, 157–182. Martinez-Bazan, C., Montaés, J.L., Lasheras, J.C., 1999b. On the breakup of an air bubble injected into a fully developed flow. Part 2. Size PDF of the resulting daughter bubbles. Journal of Fluid Mechanics 401, 183–207. Pozrikidis, C., 1998. Numerical Computation in Science and Engineering. Oxford University Press, New York. Ribeiro, L.M., Regueiras, P.F.R., Guimaraes, M.M.L., Madureira, C.M.N., Cruz-Pinto, J.J.C., 1995. The dynamic behaviour of liquid–liquid agitated dispersions—I. The hydrodynamics. Computers and Chemical Engineering 19, 333–343. Ruiz, M.C., Padilla, R., 2004. Analysis of breakage functions for liquid–liquid dispersions. Hydrometallurgy 72, 245–258. Valentas, K.J., Amudson, N.R., 1966. Breakage and coalescence in dispersed phase systems. Industrial Engineering Chemistry Fundamentals 5, 533–542. Wang, T., Wang, J., Jin, J., 2003. A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow. Chemical Engineering Science 58, 4629–4637. Wang, T., Wang, J., Jin, J., 2004. An efficient numerical algorithm for a novel theoretical breakup kernel function of bubble/droplet in a turbulent flow. Chemical Engineering Science 59, 2593–2595. Wang, T., Wang, J., Jin, J., 2005. The population balance model for gas liquid flows—influence of bubble coalescence and break-up models. Industrial Engineering Chemistry Research, in press. Yeoh, G.H., Tu, J.Y., 2005. Thermal-hydrodynamic modelling of bubbly flows with heat and mass transfer. A.I.Ch.E. Journal 51, 8–27. Ziff, R.M., McGrady, E.D., 1987. Shattering transition in fragmentation. Physical Review Letters 58, 892–895.
Toward a unified framework for the derivation of breakage functions based on the statistical theory of turbulence M. Kostogloua, b,∗ , A.J. Karabelasb, c a Department of Chemical Technology, School of Chemistry, Aristotle University, University Box 116, 541 24 Thessaloniki, Greece b Chemical Process Engineering Research Institute, P.O. Box 1517, 54006 Thessaloniki, Greece c Department of Chemical Engineering, University Box 455, Aristotle University of Thessaloniki, GR 54124 Thessaloniki, Greece
Received 31 March 2005; received in revised form 19 May 2005; accepted 24 May 2005 Available online 11 July 2005
Abstract The evolution of droplet or bubble size distribution in turbulent flow is of great significance in a variety of technological fields. Modeling this evolution by employing a population balance approach requires knowledge of the so-called breakage functions (rate and kernel). Over the years, a large number of phenomenological breakage functions with various degrees of sophistication have been proposed in the literature. Among them, those based on the statistical theory of turbulence are of particular interest, in that they attempt to take into account the structure of the flow field responsible for breakage. The purpose of the present work is to present a unified framework for developing this type of breakage functions and to show how existing models can be derived in a systematic and consistent way. The key parameters in this modeling approach are identified, which have to be determined by comparison with experimental data. It is also shown that the breakage functions, obtained within the framework presented here, lead to predictions of a droplet size evolution whose main features are consistent with experimental observations. It is suggested that this framework is an important step toward the development of a standard approach for modeling droplet size evolution in turbulent flow. 䉷 2005 Elsevier Ltd. All rights reserved. Keywords: Bubble size distribution; Mathematical modeling; Breakage process; Statistical approach to turbulence
1. Introduction Predicting the evolution of size distribution of fluid particles (bubble/droplets) in a turbulent flow field is one of the key problems of chemical engineering science. Knowledge of the size distribution is critical for the design of a large variety of equipment, including two phase-flow pipelines, bubble columns (for chemical, biochemical or other physical processes, e.g. Dudukovic et al., 1999), liquid–liquid dispersion in columns (e.g extraction; Attarakih et al., 2004a) or agitated vessels (e.g. emulsion polymerization; Chatzi et al., 1989). The size distribution of fluid particles is generally the outcome of two competing mechanisms, i.e., those ∗ Corresponding author. Department of Chemical Technology, School of Chemistry, Aristotle University, University Box 116, 541 24 Thessaloniki, Greece. Tel.: +30 2310997767. E-mail address: [email protected] (M. Kostoglou).
0009-2509/ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2005.05.051
of coalescence and breakage. In many cases coalescence is not important either naturally (e.g. for small concentration of the dispersed phase) or by design (e.g. if coalescence inhibition is applied). In this study the case of dilute dispersions is considered where only breakage is significant. The first attempts to model breakage focused on obtaining a relationship of the maximum particle size (capable of surviving in a turbulent flow field) with the physicochemical and hydrodynamic parameters of the system (Hinze, 1955). Later the dynamic nature of the breakage phenomenon was dealt with through the solution of an appropriate kinetic equation, the so-called breakage population balance (Valentas and Amudson, 1966). However, to solve such an equation for the evolution of the size distribution, the breakage functions (rate and kernel) must be known. The literature on efforts to determine breakage functions, over the past 40 years, is very large. Among the several ways proposed for obtaining
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
the breakage functions, i.e., direct observation of breakage events, direct numerical simulation, “inverse” problem solution (by employing experimental data), purely algorithmic functions and phenomenological models, the latter seems to be the most promising one. The main advantage of phenomenological models is their explicit dependence on the system physicochemical parameters. On the contrary, the frequently used algorithmic models do not include an explicit dependence on the system state but only fitting parameters. Furthermore, regarding the “inverse” breakage problem, it was recently shown (Kostoglou and Karabelas, 2005) that it is an ill-posed one (because of multiplicity of solutions) requiring a restricted space of candidate breakage functions in order to obtain a unique solution. This space may be defined by means of a phenomenological model. There are a large number of phenomenological models for breakage functions in the literature. However, great uncertainties exist regarding the reliability of such models; for instance, similar physical arguments used by various authors lead to quite different forms of the breakage functions. Additionally, breakage functions (phenomenological and algorithmic) of quite different form often seem to achieve an equally good fit of experimental size distributions (e.g. Luo et al., 2004; Chen et al., 1998; Ruiz and Padilla, 2004; Hesketh et al., 1991; Martinez-Bazan et al., 1999a,b). In the authors’ opinion, despite significant progress made in recent years, the state of the art in the development of breakage functions is unsatisfactory. There is little doubt that the best approach to follow is the formulation of a generalized (phenomenological) framework for developing breakage functions, which should be tuned and validated by carefully obtained (new and already available) experimental size distributions over a broad range of parameter values. An important step in this direction has been made by Luo and Svendsen (1996), who related the breakage functions with statistical quantities of homogeneous turbulence in a systematic and detailed way. Despite their frequent use, the final forms of the breakage functions obtained by these authors have some drawbacks; thus, several attempts to improve them, based on the same framework, have recently appeared in the literature (Lehr and Mewes, 2001; Lehr et al., 2002; Hagesaether et al., 2002a; Wang et al., 2003). However, it will be noted that among the above works there is no detailed study of the implications of the developed breakage functions on the maximum particle size dmax and on the evolving particle size distribution as well as on their relation with the system parameters. These significant issues are addressed in this paper. In the present work, a unified framework is proposed for the derivation of breakage functions based on the statistical description of isotropic turbulence. Models of this class reported in the literature are reviewed and critically evaluated. Additionally, the implications of developing particular types of breakage functions for the particle size distribution are presented and discussed. The structure of the paper is as follows: In the next section, the breakage population bal-
6585
ance is outlined. A unified framework for the derivation of breakage functions is presented in detail next. The breakage equation is then discretized in a specific way to conform to the particular form of the breakage functions. Finally, fairly extensive results, for the breakage rate and evolution of particle size distribution, are presented and discussed. 2. Problem formulation 2.1. The breakage equation The breakage process can be described in general by the following linear partial integrodifferential equation: ∞ jf (x, t) p(x, y)b(y)f (y, t) dy − b(x)f (x, t), (1) = jt x where t is the time, x the particle volume, f (x, t) the particle number density distribution, b(x) the breakage frequency (the term rate is also customarily used) and p(x, y) the probability distribution of particles of volume x resulting from the breakup of a particle of volume y. For a given initial particle size distribution f0 (x) = f (x, 0), Eq. (1) can be solved for f (x, t). The function p(x, y) should satisfy the following conditions (Ziff and McGrady, 1987): (a) Conservation of mass: y xp(x, y) dx = y. (2) 0
This equation stipulates that the total volume of particles resulting from the breakup of a particle of volume y must be equal to y. (b) z y y xp(x, y) dx (y − x)p(x, y) dx where z < 2 0 y−z (3) For binary fragmentation (two fragments per parent particle), the above restriction is simplified, being equivalent to a symmetric kernel in the sense p(x, y) = p(y − x, y). The number of particles resulting from breakage of a single particle of volume y is given as y v(y) = p(x, y) dx. (4) 0
To be able to use the fluid particle diameter d as the characteristic size, instead of particle volume, the following relations are introduced (assuming spherical particles):
3 2 f (x, t) = 2 fd (d, t). (5) d 6 d Further, a dimensionless fragment size probability function is employed: d 3 d 3 d 3 P (fv , d) = (6) p fv , with fv = x/y. 6 6 6
x=
6586
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
2.2. Particle breakage functions based on statistical approaches to isotropic turbulence 2.2.1. General framework assumptions According to the classical statistical approach to isotropic turbulence (see Baldyga and Bourne, 1999, for a comprehensive review), the flow is represented by a population of eddies exhibiting certain dynamics. The latter may be assumed very fast in comparison with the droplet dynamics, so the equilibrium properties of the eddy population may be considered. Each eddy is characterized by its size (diameter) and its kinetic energy e. The generalized (with respect to both eddy size and energy) number density of eddies is n Pe () where 12 n d is the concentration of eddies of e size between 1 and 2 , and e12 Pe de represents the fraction of eddies of size having energies between e1 and e2 . ¯ )) exp(−e/e( ¯ )) where e( ¯ ) is the Typically, Pe = (1/e( mean kinetic energy of eddies of size , depending on flow conditions. The eddies of size are considered to “collide” with the droplets of size d, with a frequency which can be determined by using concepts from the kinetic theory of gases; i.e.,
(d) = Ccs u n .
(7)
This frequency is the product of the following quantities: • Particle–eddy collision cross section, Ccs . The classical relation Ccs = (/4)(d + )2 is used here but alternatives exist; e.g., Ccs = d is used by Lehr and Mewes (2001). • Relative particle–eddy velocity. The particle is assumed “not moving”, so the relative velocity is equal to the mean velocity of eddies of size , u ; estimates of the latter are obtained from the theory of isotropic turbulence. Some authors (Hagesaether et al., 2002b) have added to this velocity the bubble velocity due to buoyancy. Such bubble motion leads to scavenging of eddies and increased breakage rate, but its importance has not been assessed quantitatively. • Number density (with respect to eddy size) n of eddies of size . Only the eddies in the inertial sub-range are considered to be effective for particle breakage. Expressions for n , u , Pe () will not be given here since they can be found elsewhere (e.g. Luo and Svendsen, 1996; Wang et al., 2003). The collision frequency does not depend on the eddy energy but the outcome of the collision does. Eddies having a sufficiently high energy level can cause breakage of droplets upon collision. The fraction of such collisions leading to breakage is called collision efficiency. Before proceeding in the quantitative analysis, the following four assumptions of the model should be stated: (i) The presence of droplets does not influence the equilibrium properties of the eddy population. A formalism for the relaxation of this assumption has been developed by Jairazbhoy and Tavlarides (2000).
(ii) Successful collisions lead to two daughter particles (i.e., binary breakage). This is a major assumption, with experimental support from detailed observations of gas bubble breakage events. The situation may be different for liquid droplets for which the experimental distributions are usually fitted assuming breakage in more than two fragments (e.g. Chatzi et al., 1989). The extension of the present approach to multiple breakage is not a trivial task, necessitating the introduction of additional parameters (e.g. number of fragments). (iii) The droplets always remain in the same state. In particular, it is assumed that the energy, absorbed from eddies incapable of causing breakage, is dissipated instantaneously (before the next collision event) so that there are no synergistic effects between successive collision events. The alternative type of transient particle behavior may be considered in the case of breakage of water droplets in high-velocity jets and its implementation requires dynamic modeling at the particle level, rendering inadequate the mean field theory examined here. (iv) As the viscosity of the dispersed phase increases, the influence of the viscous dissipation of the eddy energy on the breakage rate will increase as well. Although the effect of viscosity was taken into account by older correlations for dmax and by phenomenological models for breakage rate (see Attarakih et al., 2004b), it is not considered by the more recent models based on statistical theory of turbulence. The latter models are restricted to cases of negligible dispersed phase viscosity. Obviously, an extension of the theory to include very viscous droplets is needed. The above major assumptions are cited by all the authors using the statistical approach to turbulent breakage in fluid–fluid dispersions. The following additional working assumptions are flexible enough, providing a lot of opportunities to the developer of a particular model; these are essentially restrictions regarding the size of the eddy colliding with the parent droplet and the size of the (smaller of the 1/3 two) daughter droplet df = fv d: (a) A droplet–eddy collision can be effective only if the eddy size is smaller than the droplet size d. In general, this criterion has some physical basis (see Wang et al., 2003), but the precise upper limit set to the active eddy size is rather arbitrary. It will be pointed out, in this respect, that the breakage rate is very sensitive to the upper limit of active eddy size since large eddies contain large amounts of energy and are very efficient in breaking up droplets. Hagesaether et al. (2002a) ignored completely this criterion assuming that the breakage efficiency depends only on the energy and not on the size of the eddy. Generalizing this approach, the upper active eddy size may be considered as a parameter (max ) taking values from d to infinity and requiring further research for its more precise specification. (b) Criterion setting an upper limit to the smaller of the two daughter particle sizes (fv max ) (surface energy criterion): The energy content of the colliding eddy must be greater than that corresponding to new surface created by
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
6587
the breakage event. The form of this criterion is unique and there are no alternative interpretations. The maximum size fv max of the smaller daughter particle, resulting from the collision of a parent particle of size d with an eddy of energy e, is obtained from the solution of the equation fv2/3 max + (1 − fv
2/3 max )
= min
e emax , , d 2 d 2
(8)
where is the interfacial tension between the two fluids. The value emax = (21/3 − 1)d 2 is the maximum amount of energy which can be absorbed by the fluid particle during a breakage event and corresponds to two equal fragments. Thus, any eddy with e > emax will lead to two equal fragments. (c) Criterion setting a lower limit to daughter particle size (fv min ). Several forms of this criterion appear in the literature. According to the capillary pressure formulation (Wang et al., 2003), the dynamic pressure of the turbulent eddy c u2 /2 must be larger than the capillary pressure of the smaller fragment of breakup. The value used for the capillary pressure by Wang et al. (2003) is /df where df is the diameter of the smaller fragment. These authors state that the capillary pressure is /r, where r is a characteristic radius of curvature but it is not clear how the df used in their derivation is related to a radius of curvature. A more obvious value of the capillary pressure 2/df is suggested by Lehr et al. (2002) who assume that the fragment has a nearly cylindrical shape at the moment of breakup. A different formulation of the criterion is developed by Hagesaether et al. (2002a), who consider that the energy density of the turbulent eddy must be larger than that of the smaller fragment. Quantitatively, this means that c u2 /2 6/df . Generalizing, one may consider that the smaller fragment size is determined by the requirement c u2 /2 /df , where is a numerical parameter; this inequality leads to fv
min
=
3 6ed
3 .
(9)
In addition to this criterion for the smaller fragment size, Lehr et al. (2002) state that df must not be greater than the colliding eddy (i.e., df ). 2.2.2. Computation of breakage functions Fig. 1 is useful for depicting the above criteria. In this graph the curves fv max and fv min are shown versus energy e for a particular eddy and parent particle size and d, respectively. The value ecrit can be obtained from the solution of equation fv min (ecrit ) = fv max (ecrit ) and clearly denotes the minimum eddy energy which can lead to a breakage event. Finally, e1 = max(fv−1min (fv1 ), fv−1max (fv1 )) is the smallest eddy energy leading to a fragment of fractional volume fv1 . From the above analysis it is clear that the breakage efficiency Pb (, d) of collisions between particles of diameter d
Fig. 1. Curves fv max (e), fv min (e) and location of points (e1 , fv1 ), (ecrit , fv crit ), (emax , 21 ).
and eddies of size is ∞ ecrit (, d) Pb (, d) = Pe () de = exp − . e( ¯ ) ecrit (,d)
(10)
The total breakage rate is obtained by integrating the product of the collision rate and breakage efficiency for all eddy sizes ∞ b(d) = Pb (, d) (d) d 0 max ( + d)2 = 0.923(1 − )1/3 11/3 min ecrit (, d) × exp − (11) d , e( ¯ ) where is the volume fraction of the dispersed phase and the turbulent energy dissipation rate. The lower limit of the active eddy size is taken to be proportional to the Kolmogorov microscale of turbulence with a proportionality constant between 11.4 and 31.4 (Wang et al., 2003) in order to maintain compatibility with the used eddy size distribution which corresponds to the inertial subrange of turbulence. It is noted that, unlike the case of max , the precise value of min is of reduced significance for the computed breakage rate since small eddies have very small breakage efficiencies. In the case where only the surface energy criterion is considered, ecrit = 0; i.e., all the eddies are effective in breaking up droplets. Very small eddies give one very small fragment and another close to the parent particle. In the case where only the capillary pressure criterion is considered ecrit is obtained from the solution of fv min (ecrit )= 21 . The computation of the function P (fv , d) is much more involved than that of the rate b(d). The fundamental quantity for the computation of the general fragment size distribution is the probability distribution D(fv , e, , d) representing the probability of having a fragment of size fv from the successful collision of an eddy (with properties , e) with
6588
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
a particle of size d. The function D is non-zero only for fv min < fv < fv max and obeys the following normalization condition: fv max D(fv , e, , d) dfv = 1. (12) fv
min
It is convenient to write the function D in the following form: D(fv , e, , d) = D1 (fv ; e, , d)(U (fv − U (fv min )),
max )
(13)
where U is the Heaviside step function. The advantage of the new function D1 is that only explicit influences of e, , d are included in the argument, whereas the implicit ones through fv min and fv max are omitted. In general, the shape of D1 can differ depending on the values of e, , d, but up to now the usual assumption is that the shape of D1 with respect to the variable s = (fv − fv min )/(fv max − fv min ) is constant (i.e., D1 (fv ; e, , d) = D1 (s)). A notable exception can be found in Hagesaether et al. (2002a); these authors consider that the probability for a particular breakage event is larger as the differences of the eddy energy and energy density from those needed for the particular event are larger; i.e., D(fv , e, , d) ∝ [e − (fv2/3 + (1 − fv )2/3 − 1)
6e 6 2 × d ] . − 3 df 1/3 v
(fv − F ) = (e − F
Any shape of D1 (fv ) is possible but the usual choices are the uniform one (equal probability for all permitted fragment sizes) (Wang et al., 2003) and the Dirac delta function corresponding to the capillarity criterion fv min (Lehr et al., 2002). Having chosen the form of the function D, the fragment distribution function Pf (fv , , d) for a collision event between an eddy of size and a particle of size d is developed through the integration of D with respect to energy spectrum of the eddies: ∞ Pf (fv , , d) = D(fv , e, , d)Pe () de. (14) 0
The rate of production of fragments of size fv is given as max Pf (fv , , d) (d) d. (15) bf (fv , d) = min
The fragment size distribution can be obtained from P (fv , d) =
bf (fv , d) . b(d)
−1
dF (fv )) de
−1 ,
(18)
where the notation F −1 (fv ) stands for the solution of the equation fv = F (e) ⇒ e = F −1 (fv ). Substituting in Eq. (14) leads to Pf (fv , , d) = PF −1 (fv )
dF de
−1 e=F −1 (fv )
.
(19a)
This relation cannot be used in the case of F independent of e because the slope of the F (e) curve becomes zero; a typical example is the case F = fmax where F = 21 for e > emax . Eq. (19a) is then transformed to Pf (fv , , d) = PF −1 (fv )
dF de
−1
+ (fv − 1/2)
e=F −1 (fv ) ∞ emax
Pe () de.
(19b)
Case (ii): Uniform probability distribution (i.e., D1 = 1/(fv max − fv min )). In this case Eq. (14) can be written as
(16)
The total breakage rate b(d) can be obtained as an alternative to Eq. (11) by integration of the rate bf over all the fragment sizes, i.e., 1/2 bf (fv , d) dfv . (17) b(d) = 0
Although the use of Eq. (17) could be used for the computation of the denominator in Eq. (16) (i.e., a normalized P (fv , d) is needed), it is not recommended for the computation of the breakage rate b(d), which alternatively can be made through Eq. (11). The reason is that Eq. (17) requires a complicated computational procedure, including the use of function D, which is susceptible to computational error. However, this source of error can be avoided (bypassed through Eq. (11)) since b(d) does not really depend on D. The function P (fv , d) results from a double integral which in general is computed numerically. This computational load can be significantly reduced for some simple forms of the function D1 used in the literature; a summary of such cases follows. Case (i): D1 (fv ) = (fv − F ), where F is a given combination of fv max and fv min . For example, F = fv max means that the particle breaks into two fragments obeying exactly the surface energy criterion; F = fv min means that the smaller fragment obeys exactly the capillarity criterion; F = (fv min + fv max )/2 means that the size of fragment is the average of fv min and fv max . F is a function of energy e through the e-dependence of fv max and fv min . From the properties of Dirac delta functions one obtains
Pf (fv , , d) =
∞ max[fv−1min (fv ),fv−1max (fv )]
× Pe () de.
fv
1 − fv max
min
(20)
It is observed that the integrand does not depend on fv . The dependence on fv is only through the lower integration limit, so a sequential computation of the function Pf at the points fvi = i/(2M) for i = 1, 2, 3, . . . , M is possible, as follows:
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
For fvi > fv crit Pf (fvi , , d) = Pf (fvi+1 , , d) +
bf (fv , d) =
fv−1max (fvi+1 )
fv−1max (fvi ) × (fv max − fv min )−1 Pe () de
(21a)
with Pf (fvM , , d) =
∞ emax
1 − fv 2
−1 min
Pe () de.
(21b)
For fvi < fv crit Pf (fvi , , d) = Pf (fvi−1 , , d) +
fv−1min (fvi )
fv−1min (fvi−1 ) × (fv max − fv min )−1 Pe () de
with Pf (fv1 , , d) =
∞ fv−1min (fv1 )
(fv
× Pe () de.
max
− fv
min )
(22a)
−1
(22b)
The integrals in (21b), (22b) include an exponential term in Pe () and with a simple coordinate transformation can be cast in a form appropriate for the use of the Gauss–Laguerre quadrature. The rest of the integrals are computed using the orthogonal rule which is more accurate than the trapezoidal one (Pozrikidis, 1998). This procedure, proposed by Wang et al. (2004), greatly facilitates the computation of the breakage kernel since it reduces the order of integration from two to one. It is noted that Wang et al. (2003, 2004) have assumed that no breakage occurs for fv max − fv min < , where is a parameter having a small value. It is believed that this addition has been made for computational purposes but it tends to reduce the breakage rate since ecrit is now obtained from the solution of fv max (ecrit ) − fv min (ecrit ) = . Providing a unified framework based on the statistical theory of turbulence, which will encompass breakage models already proposed in the literature, would greatly facilitate further developments and the correct derivation of additional models. However, as a counter example of incorrect application of this approach, one may refer to the breakage functions derived in the pioneering work of Luo and Svendsen (1996). The key assumption in their derivation is that only the surface energy criterion is obeyed and all the fragment sizes have a uniform probability between 0 and fv max . According to the present analysis, making use of this assumption leads to the following breakage function: max ( + d)2 b(d) = 0.923(1 − )1/3 d 11/3 min 0.923(1 − )1/3 −2/3 = [2(min − −2/3 max ) 3 −5/3 −8/3 2 −8/3 + 10d(min − −5/3 max ) + 8d (min − max )], (23a)
max min
∞ 0
1 fv
max
Pe () de (d) d.
6589
(23b)
The above rate is the maximum one and corresponds to breakage efficiency equal to 1. The derivation of Luo and Svendsen (1996) is incorrectly based on the concept of conditional breakage probability for a given fragment size fv , overlooking that eddies incapable of producing fragments of size fv can be efficient in producing fragments of smaller size. This inconsistency is recognized by Wang et al. (2003) only as a drawback of the model and not as an error in the derivation. Due to this erroneous derivation the model of Luo and Svendsen (1996) (having collision efficiency 1) leads to smaller breakage rate compared to that of Wang et al. (2003) which has a collision efficiency smaller than 1. Although the physical assumptions of the model of Luo and Svendsen (1996) have been extensively criticized and relaxed, the aforementioned error has been essentially overlooked; thus even today the model (with no amendment) is repeatedly used for simulation and comparison with experiments (e.g. Chen et al., 2005; Dhanasekharan et al., 2005; Yeoh and Tu, 2005). 2.3. Numerical solution of the breakage equation Unlike the case of solid fragmentation, the extent of breakage in fluid–fluid dispersions is not very large; thus, the usual finite difference methods (in the particle diameter domain) are extensively used in the literature for the discretization of Eq. (1) (Ribeiro et al., 1995). Recently, for the solution of Eq. (1), Hagesaether et al. (2002b) developed an alternative method. The latter has the drawback that it permits only a geometrical grid with a ratio 2 with respect to particle volume. The need for a finer and arbitrary grid renders the sectional approach of Kumar and Ramkrishna (1996) the preferred choice, even for fluid–fluid applications. According to Kumar and Ramkrishna (1996) the PSD can be represented by a sum of discrete particle sizes xi with number concentration Ni (i = 1, 2, 3, . . . , L), or mathematically f (x, t) =
L
Ni (x − xi ).
(24)
i=1
Consequently, the following system of equations must be solved for the particle concentration Ni : L
dNi nik b(xk )Nk − b(xi )Ni , = dt
(25)
k=i
ni,k = (1 − ik )Z(xi , xi+1 , xk ) − (1 − 1k ) × Z(xi , xi−1 , xk ), a2 a2 − y Z(a1 , a2 , a3 ) = p(y, a3 ) dy. a 2 − a1 a1
(26) (27)
6590
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
The moments are computed as follows: Mk =
L i=1
Ni xik .
(28)
The method is used in bubble breakup applications with geometrical grid of ratio R = 2 (Chen et al., 2005), R = 1.7 and R = 1.45 (Wang et al., 2005). This choice is appropriate for the case of combined breakage–coalescence in which the particle size distribution does not evolve very much. However, in the case of pure breakage treated here, the size distribution continuously evolves becoming narrower, as will be shown in the next section. Thus, after some time, the geometric grid is not sufficiently dense to accurately describe the particle size distribution. The following grid construction procedure is proposed, combining high accuracy and computational efficiency: First, a particle volume x0 is determined for which the particles undergo insignificant breakage during the required simulation time (i.e., b(x0 )tfin = 0.001). Then, the region between the size of the largest initial particles xmax and x0 is discretized in L2 points using a geometrical grid with each interval equal to R times the previous one. No breakage is considered for particles of size smaller than x0 . The region from 0 to x0 is discretized using L1 equidistant points. This choice leads to a controlled resolution of the particle size distribution for t < tfin and has the computational advantage that the matrix Pi,j is not needed for j < L1 (L = L1 + L2 ). Also with this type of discretization the well-known finite domain error is of no concern. It is noted that in the method proposed by Hagesaether et al. (2002b) the resolution with respect to fv of the function P (fv , d) is dictated by the discretization of Eq. (1). However, in general the resolution required for P (fv , d) is much finer (and on different grid points) than that required for the sectional method; thus, the two stages of computation must be decomposed. Having chosen the grid points xi (i = 1, 2, . . . , L), the matrix Pij = P (fvi , dj ) (i = 1, 2, . . . , M, j =L1 , 2, . . . , L) is constructed. The integrals in Eq. (27) are determined using the trapezoidal rule at K points and the integrand is computed at arbitrary values of fv by linearly interpolating between the values of the matrix Pij . In the present work the values K =70, L1 =10, L2 =50, R = 1/0.85 are selected (after several trials) because they provide satisfactory accuracy for the cases treated here, and lead to practically grid-independent results. It is noted that the combination of a very coarse grid (e.g. the nine classes in Dhanasekharan et al., 2005) with very steep kernels (e.g. that of Luo and Svendsen, 1996), and coarse discretization of integrals in (27) (e.g. the five Gaussian points in Attarakih et al., 2004b) can lead to questionable results.
3. Results and discussion In this section the form of the breakage rate and the evolution of the droplet size distribution will be obtained for
Fig. 2. The function Pb (fv ) in the neighborhood of its singularity for d =6 mm, =4.2 mm, =1 m2 /s3 , =72 m N/m and for several numbers M of discretization points fvi .
several sets of parameters and will be discussed in some detail. The analysis will be restricted to the standard model presented by Wang et al. (2003); i.e., uniform breakage probability for fv min < fv < fv max and fv max − fv min > . The function Pb (fv , , d) is singular going to infinity at fv crit = fv min (ecrit ). A non-zero value of can remove the singularity and facilitate the computations but reduces significantly the breakage rate without having any physical meaning. Thus, cannot be considered to be a simple parameter for computations, and has no place in the model. Of course, the numerical computation of the singular function Pb is still a problem but it must be overcome using numerical means, without altering the physics of the model. Two important points will be stressed here. First, the breakage rate can be computed directly from Eq. (11) avoiding any relation with Pb ; second, the singularity of Pb with respect to fv is integrable (i.e., the integral of Pb in any finite interval around the singularity is finite). The only numerical requirement is the use of a very fine grid (large M value) with respect to fv in order to capture the shape of Pb in the region of singularity. To make this point explicit, the function Pb is shown in the region of singularity for the case d = 6 mm, = 4.2 mm, = 1 m2 /s3 and several values of M, in Fig. 2. A plot of Pb for the same conditions but with = 0.001 can be found in Wang et al. (2004). It is obvious that all the grids give the same result, except at the points close to the singularity. As the grid becomes more dense, on one hand the correct shape of Pb is followed closer, and on the other any discrepancy is restricted to a smaller domain, leading to a smaller overall error. In the present case the total breakage rate is obtained with three digits of accuracy for M = 1000. In order to have an accurate description of the breakage kernel, at least M = 200 points must be used for fv . The seemingly obvious suggestion of choosing a coarse grid and making it denser in the vicinity of fv crit is not really good since fv crit depends on ; so multiple grids must be used making the situation even more complicated. The significant effect of the artificial parameter (used for singularity masking) on the breakage rate is shown in
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
Fig. 3. Breakage rate b ( = 1 m2 /s3 , = 72 m N/m) computed by Eq. (17) versus parameter , for several values of particle diameter and discretization points M.
Fig. 3 ( = 1 m2 /s3 , = 72 m N/m). The breakage rate has been computed using the wrong approach (Eq. (17)) to show the effect of discretization. At values small enough with respect to fv discretization, there is no influence on the breakage rate; however, above a crossover value 0 there is a significant effect on b(d). Indeed, increasing leads to decreasing b(d), as expected. Using a finer grid, the crossover value 0 decreases. In the ideal limit of infinitely dense grid (or alternatively using the correct way, Eq. (11), for b(d)) 0 = 0. One can notice inFig. 3 that the breakage rate for d =5 mm and =0 is not the same for M =200 and M =400. As already mentioned, a value M = 1000 is needed for a converged breakage rate value. Evaluation of the breakage rate using Eq. (11) involves relatively small computational load. The trapezoidal rule is used for the numerical evaluation of the integral and the bisection method for ecrit computation from the transcendental equation fv min (e) = fv max (e). The breakage rate as a function of the droplet diameter for surface tension =25 m N/m and three values of is shown in Fig. 4. In the b(d) curves one can discern two parts: i.e., a rather smooth variation with d for large d values, and a very strong dependence on d for relatively small d values. The sharp part of the curve becomes even sharper with increasing . It is remarkable that reduction of the particle diameter by a few percentage points leads to reduction of the breakage rate by several orders of magnitude. An attempt to fit the b(d) curves with algebraic relations was unsuccessful. Also, the general form of the breakage rate given by Lehr et al. (2002) but with unknown parameter values was incapable of fitting the curvesshown in Fig. 4. This confirms the assertion that the much simpler expression suggested by Lehr et al. (2002) cannot be used in place of the approach examined here (see Wang et al., 2005). In Fig. 5, the breakage rate is depicted versus the interfacial tension of the dispersion for several sets of
6591
Fig. 4. Breakage rate b versus particle diameter d for = 25 m N/m and several values of turbulent energy dissipation rate .
Fig. 5. Breakage rate b versus interfacial tension for several values of particle diameter d and turbulent energy dissipation rate .
parameters and d. As increases the breakage rate decreases, as one would have expected. The relative reduction of b for the same increase of seems to be larger as the absolute value of b is smaller. Results for the breakage kernel will not be presented here since everything needed can be found in the recent literature (see Wang et al., 2003). It will be pointed out, though, that for high breakage rates the kernel is of “M” shape with the two peaks very close to fv = 0. As the breakage rate is reduced, the peaks are displaced inward and finally collapse to a single peak at fv = 21 . The simple kernel of Lehr et al. (2002) exhibits the same behavior. Extensive comparisons between the kernels based on the statistical approach to turbulence are presented by Wang et al. (2003). It appears that the key feature determining whether a kernel exhibits a “U” shape or a “U–M” transition (with d decreasing) is the existence or not of fragments of very small size(existence of minimum fragment size criterion).
6592
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
Fig. 6. Evolution of diameter dm (corresponding to mean particle volume) for = 50 m N/m and several values of and initial monodisperse size d0 (small time scale).
An issue of paramount importance is the effect of breakage functions, described above, on the evolution of droplet size distribution in a turbulent flow field. In the following simulations, aimed at assessing such effects, the initial droplet distribution is monodisperse with diameter d0 . The diameter dm (called hereafter mean diameter) is the one corresponding to the mean volume of the droplets xm = M1 /M0 (i.e., dm = (6xm /)1/3 ). The evolution of dm for = 50 m N/m and three values of is shown in Fig. 6; for each value of two different initial diameters are considered. It is very interesting that after some time the system “forgets” its initial condition and the mean size evolution curves collapse to one. This is shown in Fig. 6 for only two different initial conditions but it can be generalized for any initial condition. The time of collapse of the evolution curves into a single one depends on the initial condition. This means that the breakage rate exhibits a “self-regulating” tendency, leading to a long time (seemingly) asymptotic evolution curve, independent of the initial size distribution. By inspecting Fig. 6 one would be tempted to conclude that dm reaches a final constant value and that the approach to this “steady state” is faster for larger values of . This trend would be in accordance with the classical view of existence of a maximum stable particle size dmax . However, the picture obtained from Fig. 6 is quite deceptive. In Fig. 7 the same evolution curves are shown in a long time scale; obviously, what might be considered as “steady state” in Fig. 6 corresponds to a very slow “sliding” of dm to smaller values. During the past 15 years, some experimental studies on the size evolution of dispersions in turbulent flow have cast doubt on the existence of a dmax or equivalently of a true steady state droplet size distribution. Kostoglou and Karabelas (2001) tried to explain these results by arguing that dmax may or may not exist but, from a practical point of view, it does not matter since it is experimentally inaccessible due to the form of the breakage rate function. According to an explanation offered by Baldyga and Bourne (1999) dmax is a sliding function of time due to the intermittent nature of turbulence; i.e., as time passes the chance of a drop
Fig. 7. Evolution of diameter dm (corresponding to mean particle volume) for = 50 m N/m and several values of (large time scale).
Fig. 8. Evolution of diameter dm (corresponding to mean particle volume) for = 1 m2 /s3 and several values of .
to encounter a more violent turbulent burst increases leading to reduction of dmax . The phenomenological approach to breakage presented here appears to be compatible with both the above interpretations of experimental observations regarding a non-existent steady state. According to the present approach, although a true steady state may not really exist, the PSD evolution is very slow rendering its existence of no practical concern, as Kostoglou and Karabelas (2001) have already suggested. Furthermore, from the results of Figs. 6 and 7 one may consider that a pseudo-steady state exists depending on the time scale of observation. This seems to agree with the concept of dmax “sliding” with time (Baldyga and Bourne, 1999). The evolution of dm for = 1 m2 /s3 and several values of is shown in Fig. 8. As already pointed out, the curves for large values of time are independent of d0 ; also, as expected, the droplet size decreases with decreasing interfacial tension. An important question is what happens to the evolving particle size distribution; i.e., is there a long time asymptote for the size distribution regardless of the initial distribution, or is this the case only for the mean size dm but with a different distribution shape? The present modeling effort suggests that an asymptotic size distribution does exist. Size distributions of droplets independent of the initial
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
Fig. 9. Long time (t = 10 000 s) normalized particle size distributions fd (d)/M0 for = 1 m2 /s3 and several values of .
Fig. 10. Evolution of the normalized particle size distributions fd (d)/M0 for = 0.1 m2 /s3 , = 50 m N/m and d0 = 4 mm.
condition are shown in Fig. 9 for t =10 000 s, =1 m2 /s3 and three values of . The distributions hereafter are normalized with the corresponding total droplet number concentration M0 inorder to make their presentation feasible in the same graph. It is pointed out that the “hat-type” shape and the slight bimodality of the distributions are not artificial due to the grid used; similar results are obtained using finer grids. Additionally, the size distribution becomes narrower as the surface tension () decreases. The evolution of the droplet size distribution for the case d0 = 4 mm, = 50 m N/m and = 0.1 m2 /s3 is shown in Fig. 10. Initially a tri-modal PSD appears consisting of the unbroken initial particles and the two modes of fragments (U-shape breakage kernel for large particles). At later times, the PSD becomes unimodal. Surprisingly enough, at very large time values, a weak bimodality appears and tends to evolve. The origin of this bimodality is not obvious but it is believed to be due to the combination of the shape (bell shape toward equal size breakage) of the breakage kernel for small droplet sizes with the extremely large slopes of the b(d) curve (e.g. Fig. 4). This creates discrete accumulation modes in the droplet size domain, having extremely different breakage rates. This difference in breakage rate may restrict the number of modes to two. The smaller mode can-
6593
Fig. 11. Evolution of the normalized particle size distributions fd (d)/M0 for = 1 m2 /s3 , = 50 m N/m and d0 = 4 mm.
Fig. 12. The maximum particle diameter at large time (t =10 000 s) versus interfacial tension , for two values of .
not suffer significant breakage (and generate a third mode) before the extinction of the larger one. The corresponding PSD evolution, for the case = 1 m2 /s3 , is shown inFig. 11. At the initial stages of breakage the size distribution is bimodal; a discrete third mode corresponding to the large size fragments does not appear here. Otherwise, the evolution is quite similar to that of the previous cases, with a bimodality appearing at large time values. The maximum droplet size at the time of 10 000 s is shown in Fig. 12 versus surface tension for two values of . It is very interesting that both curves can be fitted with a power law having exponent 0.56. The latter is close to the 0.6 exponent in the well-known relation of Hinze (1955) for dmax . Similarly, the maximum droplet size at time 10 000 s is shown in Fig. 13 versus mean dissipation rate for = 50 and 70 m N/m. Here again, there is a perfect fit of both curves with an exponent −0.43 of a power function, as compared to the −0.40 exponent of the expression of Hinze (1955). It will be noted that these dmax dependencies in the Hinze expression have been confirmed by Karabelas (1978) for turbulent pipe flow of dilute liquid–liquid dispersions at small flow times (10–30 s) and mean values in therange 0.4–8.8 m2 /s3 . Thus, according to the present theory, whereas a true steady state value dmax may not really exist, the maximum droplet size at a specified (sufficiently long) time is independent of the initial size distribution and
6594
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
d D(fv , e, , d)
Fig. 13. The maximum particle diameter at large time (t =10 000 s) versus mean dissipation rate for two values of .
exhibits a parametric dependence very similar to that assigned to dmax by the classical approach. As an example to demonstrate the importance of the generalized framework presented here one may refer to the recent attempt by Chen et al. (2005) to simulate bubble columns using the original Luo and Svendsen (1996) breakage model; these authors state that a 10-fold increase of computed breakage rate is required in order to achieve agreement with the experimental data. However, as pointed out already, the correct computation of the Luo and Svendsen rate, with the use of a much larger than d effective eddy size max , leads to much larger breakage rates without abandoning the general framework.
4. Concluding remarks A generalized framework for the derivation of breakage functions based on the statistical theory of turbulence is presented here. It is shown how several recently proposed breakage functions emerge from this framework. Identified as key parameters for the derivation of breakage functions are the maximum eddy size active for breakage (max ), the minimum size of fragment resulting from a breakage event (i.e., parameter in Eq. (9)) and the function D1 (fv , e, , d) defining the fragment size distribution resulting from a particle–eddy collision. The evolution of particle size distribution based on a particular member of the breakage functions from the generalized framework (i.e., the one proposed by Wang et al., 2003) is examined in detail to show that the functions of the framework are compatible with experimental observations and with their theoretical interpretation about the bubble/droplet evolution in turbulent dispersions.
Ccs
Greek symbols
Notation b(x) bf (fv , d)
parent particle diameter dimensionless size probability function for fragments resulting from the collision of an eddy with size and energy e with a particle of size d D1 (fv ; e, , d) modified form of D (see Eq. (13)) e eddy energy ecrit the lowest energy of an eddy of size capable of breaking a particle of size d emax required energy for breakage in two equal fragments of a particle of size d e( ¯ ) mean kinetic energy of eddies of size f (x, t) particle number density distribution (volume based) fd (d, t) particle number density distribution (diameter based) fv ratio of fragment to parent particle volume fv max maximum fractional volume of the smaller of two fragments (see Eq. (8)) fv min minimum fractional volume of the fragment (due to capillarity criterion) fv crit fragment volume fraction corresponding to ecrit K number of discretization point for numerical evaluation of integrals in Eq. (27) L total number of classes for the numerical solution of the breakage equation M number of discretization points for fv in the interval [0, 0.5] Mk kth moment of the particle size distribution n number density of eddies with respect to eddy size p(x, y) fragment size probability function P (fv , d) dimensionless fragment size probability function Pe () probability density function of energy of eddies with size Pb (, d) efficiency of collisions between particles of diameter d and eddies of size Pf (fv , , d) dimensionless size probability function for fragments resulting from the collision of an eddy with size with a particle of size d u relative continuous phase particle velocity x, y particle volume
breakage frequency for particle of volume x frequency of production of fragments with volume fv d from breakage of particles of size d particle–eddy collision cross section
min , max c
volume fraction of the dispersed phase generalization parameter for the capillarity criterion turbulent energy dissipation rate eddy diameter minimum and maximum eddy sizes effective for particle breakage density of thecontinuous phase
M. Kostoglou, A.J. Karabelas / Chemical Engineering Science 60 (2005) 6584 – 6595
(d)
interfacial tension between particles and continuous phase density with respect to of collision frequency of particles of size d with eddies of size
References Attarakih, M.M., Bart, H.J., Faqir, N.M., 2004a. Solution of the droplet breakage equation for interacting liquid–liquid dispersions: a conservative discretization approach. Chemical Engineering Science 59, 2547–2565. Attarakih, M.M., Bart, H.J., Faqir, N.M., 2004b. Numerical solution of the spatially distributed population balance equation describing the hydrodynamics of interacting liquid–liquid dispersions. Chemical Engineering Science 59, 2567–2592. Baldyga, J., Bourne, J.R., 1999. Turbulent Mixing and Chemical Reactions. Wiley, New York. Chatzi, E.G., Gavrielides, A.D., Kiparissides, C., 1989. Generalized model for prediction of the steady state drop size distributions in batch stirred vessels. Industrial Engineering Chemistry Research 28, 1704–1711. Chen, Z., Pruss, J., Warnecke, H.J., 1998. A population balance model for disperse systems: drop size distribution in emulsion. Chemical Engineering Science 53, 1059–1066. Chen, P., Sanyal, J., Dudukovic, M.P., 2005. Numerical simulation of bubble columns flows: effect of different breakup and coalescence closures. Chemical Engineering Science 60, 1085–1101. Dhanasekharan, K.M., Sanyal, J., Jain, A., Haidari, A., 2005. A generalized approach to model oxygen transfer in bioreactors using population balances and computational fluid dynamics. Chemical Engineering Science 60, 213–218. Dudukovic, M.P., Larachi, F., Mills, P.L., 1999. Multiphase reactors—revisited. Chemical Engineering Science 54, 1975–1995. Hagesaether, L., Jakobsen, H.A., Svendsen, H.F., 2002a. A model for turbulent binary breakup of dispersed fluid particles. Chemical Engineering Science 57, 3251–3267. Hagesaether, L., Jakobsen, H.A., Svendsen, H.F., 2002b. Modeling of the dispersed-phase size distribution in bubble columns. Industrial Engineering Chemistry Research 41, 2560–2570. Hesketh, R.P., Etchells, A.W., Russell, T.W.F., 1991. Bubble breakage in pipeline flow. Chemical Engineering Science 46, 1–9. Hinze, J.O., 1955. Fundamentals of the hydrodynamic mechanics of splitting in dispersion processes. A.I.Ch.E. Journal 1, 289–295. Jairazbhoy, V., Tavlarides, L.L., 2000. A numerical technique for the solution of integrodifferential equations arising from balances over populations of drops in turbulent flows. Computers and Chemical Engineering 23, 1725–1735. Karabelas, A.J., 1978. Droplet size spectra generated in turbulent pipe-flow of dilute liquid–liquid dispersions. A.I.Ch.E. Journal 24, 170–180.
6595
Kostoglou, M., Karabelas, A.J., 2001. A contribution towards predicting the evolution of droplet size distribution in flowing dilute liquid/liquid dispersions. Chemical Engineering Science 56, 4283–4292. Kostoglou M., Karabelas A.J., 2005. On the-self similar solution of fragmentation equation: Numerical evaluation with implication for the inverse problem. Journal of Colloid and Interface Science 284, 571–581. Kumar, S., Ramkrishna, D., 1996. On the solution of population balance equations by discretization—I. A fixed pivot technique. Chemical Engineering Science 51, 1311–1332. Lehr, F., Mewes, D., 2001. A transport equation for the interfacial area density applied to bubble columns. Chemical Engineering Science 56, 1159–1166. Lehr, F., Mewes, D., Millies, M., 2002. Bubble-size distributions and flow fields in bubble columns. A.I.Ch.E. Journal 48, 2426–2443. Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. A.I.Ch.E. Journal 42, 1225–1233. Luo, G.S., Li, H.B., Tang, X.J., Wang, J.D., 2004. Drop breakage in a coalescence–dispersion pulsed-sieve plate extraction column. Chemical Engineering Journal 102, 185–191. Martinez-Bazan, C., Montanes, J.L., Lasheras, J.C., 1999a. On the breakup of an air bubble injected into a fully developed turbulent flow. Part 1. Breakup frequency. Journal of Fluid Mechanics 401, 157–182. Martinez-Bazan, C., Montaés, J.L., Lasheras, J.C., 1999b. On the breakup of an air bubble injected into a fully developed flow. Part 2. Size PDF of the resulting daughter bubbles. Journal of Fluid Mechanics 401, 183–207. Pozrikidis, C., 1998. Numerical Computation in Science and Engineering. Oxford University Press, New York. Ribeiro, L.M., Regueiras, P.F.R., Guimaraes, M.M.L., Madureira, C.M.N., Cruz-Pinto, J.J.C., 1995. The dynamic behaviour of liquid–liquid agitated dispersions—I. The hydrodynamics. Computers and Chemical Engineering 19, 333–343. Ruiz, M.C., Padilla, R., 2004. Analysis of breakage functions for liquid–liquid dispersions. Hydrometallurgy 72, 245–258. Valentas, K.J., Amudson, N.R., 1966. Breakage and coalescence in dispersed phase systems. Industrial Engineering Chemistry Fundamentals 5, 533–542. Wang, T., Wang, J., Jin, J., 2003. A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow. Chemical Engineering Science 58, 4629–4637. Wang, T., Wang, J., Jin, J., 2004. An efficient numerical algorithm for a novel theoretical breakup kernel function of bubble/droplet in a turbulent flow. Chemical Engineering Science 59, 2593–2595. Wang, T., Wang, J., Jin, J., 2005. The population balance model for gas liquid flows—influence of bubble coalescence and break-up models. Industrial Engineering Chemistry Research, in press. Yeoh, G.H., Tu, J.Y., 2005. Thermal-hydrodynamic modelling of bubbly flows with heat and mass transfer. A.I.Ch.E. Journal 51, 8–27. Ziff, R.M., McGrady, E.D., 1987. Shattering transition in fragmentation. Physical Review Letters 58, 892–895.