Master Curves for Aggregation and Gelation: Effects of Cluster

Feb 10, 2007 - The aggregation kinetics is described in terms of master curves characterizing the evolution of suitable dimensionless averages (measur...
0 downloads 0 Views 213KB Size
Download PDF
Ind. Eng. Chem. Res. 2007, 46, 1709-1720

1709

Master Curves for Aggregation and Gelation: Effects of Cluster Structure and Polydispersity Miroslav Soos,† Jan Sefcik,‡ and Massimo Morbidelli*,† Institute for Chemical and Bioengineering, Department of Chemistry and Applied Bioscience, ETH Zurich, 8093 Zurich, Switzerland, and Department of Chemical and Process Engineering, UniVersity of Strathclyde, James Weir Building, 75 Montrose Street, Glasgow G1 1XJ, Scotland, United Kingdom

A parametric study of the effects of cluster structure and polydispersity on the kinetics of aggregation and gelation is presented. The aggregation kinetics is described in terms of master curves characterizing the evolution of suitable dimensionless averages (measurable by light scattering) of the underlying cluster mass distribution (CMD), as a function of a suitable dimensionless time. Such master curves are shown to be dependent only on two dimensionless parameters: the cluster fractal dimension (df) and the ratio between Brownian and shear aggregation rates (κ*). Shear aggregation of fractal clusters leads to higher-order moments of the CMD diverging to infinity in a finite time, which is usually called runaway. It is determined that the parameter space is split into two distinct regions: either where aggregates gel because of space filling before runaway occurs or where gelation occurs as a consequence of runaway. It is also determined that cluster polydispersity and the fractal dimension significantly affect the runaway and gelation times. 1. Introduction Shear aggregation occurs when particles carried by a suspending fluid are brought into contact because of fluid velocity gradients (shear) and are able to stick together. The kinetics of shear aggregation in not-too-dense systems is second-order in particle concentration and the corresponding rate constant is proportional to the third power of the sum of the radii of the colliding particles. This strong dependency on particle dimension leads to self-accelerating aggregation dynamics, which typically results in wide particle distributions broadening with time until other processes (such as sedimentation, breakage, or gelation) take over and change the process dynamics. Here, we focus on the early behavior of shear aggregation processes governed by pure aggregation. These processes can be conveniently modeled by population balance equations (PBEs) with appropriate aggregation kernels. To obtain analytical solutions for shear aggregation kinetics, simplifying assumptions must be made. For example, if one assumes that the cluster distribution remains monodisperse during the aggregation process, it is possible to derive analytical relations for the time evolution of cluster size for pure Brownian aggregation, as well as for shear aggregation.1 These relations allow one to model the dynamic behavior of aggregating systems as a function of the initial solid volume fraction, the size of the primary particles, the shear rate, and the fractal dimension (df). However, similar analytical solutions are not available for the general case of polydisperse cluster mass distributions (CMDs) under shear conditions, and one must use numerical approaches to solve the corresponding population balance equation, accounting for more-realistic cluster distributions. Examples of the comparison of the result obtained for monodisperse distribution and numerical solution of general PBEs using a very coarse numerical grid were presented by Gonzalez et al.2 In this work, they * To whom correspondence should be addressed. Tel.: +41 44 6323033.Fax: +41446321082.E-mail: [email protected]. † Institute for Chemical and Bioengineering, Department of Chemistry and Applied Bioscience, ETH Zurich. ‡ Department of Chemical and Process Engineering, University of Strathclyde.

derived correlations between the two approaches; however, their results were dependent on the numerical grid that was used and, therefore, cannot be used for conditions outside the tested range. The objective of this work is to analyze the dynamics of shear aggregation in particulate systems, in terms of master curves that describe the evolution of suitable dimensionless variables (related to measurables from light scattering measurements) in dimensionless time. Such master curves are dependent on only a few dimensionless parameters that concisely describe the corresponding dynamics in a general way. In particular, we investigate the runaway behavior of these systems, which corresponds to the situation where sufficiently large aggregates are produced to lead to a very fast growth, or runaway, of the cluster mass distribution, which, in the limit, leads to the formation of a single cluster. This behavior is characterized by the fact that all moments of the distribution with an order equal to or larger than 2 diverge; that is, the mass-weighted average mass of the distribution grows to an infinite value. This phenomenon is similar to the formation of chemical gels in polymeric systems undergoing crosslinking3 (not to be confused with the purely physical gels produced by aggregation considered in this work), as well as to the thermal runaway in exothermic reacting systems.4 The runaway behavior of aggregating systems should not be confused with the physical gelation process. The latter corresponds to the situation where the aggregates, which grow increasingly open with size (because of their fractal nature), occupy all the available space. In other words, they form a network percolating through the entire system volume. Gelation does not necessarily imply runaway behavior, because the aggregates can fill the entire volume available without their average mass diverging to infinity. This has been confirmed in previous work by Sandku¨hler and co-workers,5-7 in the case of purely Brownian aggregating systems. It was shown that, in this case, runaway does not occur and the aggregate size grows continuously in time, always with a finite rate, until gelation occurs. In this work, we will specifically investigate the role that shear has on the aggregation behavior, with specific reference to gelation and runaway phenomena.

10.1021/ie061220+ CCC: $37.00 © 2007 American Chemical Society Published on Web 02/10/2007

1710

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007

2. Model Description 2.1. Population Balance Equation. The time evolution of the CMD due to aggregation of monodisperse primary particles can be described by the following population balance equation:

dn(ξ,t) 1 ) dt 2

∫1ξ-1 K(ξ - ξ′,ξ′)n(ξ - ξ′,t)n(ξ′,t) dξ′ ∞ n(ξ,t) ∫1 K(ξ,ξ′)n(ξ′,t) dξ′

(1)

where n(ξ,t) is the number concentration of clusters with dimensionless mass ξ (expressed in terms of the number of primary particles) at time t. During an aggregation event, a cluster with mass ξ - ξ′ collides with a cluster of mass ξ′, and if the collision is successful, an aggregate with mass ξ is formed. For relatively dilute conditions, where only binary collisions are significant,8 the matrix of rate constants (the aggregation kernel, K(ξ - ξ′,ξ′)) is independent of the particle concentration. 2.2. Aggregation Kernel. Aggregation processes can be generally driven by various physical mechanisms, such as fluid velocity gradients and fluctuations, Brownian motion, or sedimentation. Shear aggregation is driven by collisions of particles due to spatial fluid velocity gradients, and the corresponding rate constant can be expressed as follows:9 shear ) RG(Rc,ξ-ξ′ + Rc,ξ′)3 Kξ-ξ′,ξ′

(2)

where R is the aggregation rate prefactor, G is the shear rate, and Rc,ξ-ξ′ and Rc,ξ′ are the collision radii of the colliding aggregates. Under certain simplifying assumptions, it is possible to derive the value of the aggregation rate prefactor theoretically. For example, for a simple shear flow undisturbed by the presence of particles, it is equal to 4/3,10 whereas for turbulent flows, the value of 1.294 was derived by Saffman and Turner.11 However, the values of such an aggregation rate prefactor observed experimentally12-16 are typically one order of magnitude smaller than the values previously mentioned. Therefore, in the following, we treat the aggregation rate prefactor as an unknown parameter that is dependent on system properties, such as the primary particle size, particle interactions, and the type of fluid flow. Brownian aggregation is driven by collisions of particles due to their Brownian motion. This mechanism becomes more important for smaller (sub-micrometer) particles and the corresponding rate constant is given by9 Brownian ) Kξ-ξ′,ξ′

4π + Dξ′)(Rc,ξ-ξ′ + Rc,ξ′) (D W ξ-ξ′

(3)

where W is the stability ratio, which is dependent on the interaction forces between colliding aggregates, whereas Dξ-ξ′ and Dξ′ are their diffusion coefficients. In the case of negligible settling and inertial effects, the overall aggregation rate constant for fully destabilized suspended colloidal particles can be expressed as the sum of the two rate constants, corresponding to shear and Brownian aggregation:17 shear Brownian + Kξ-ξ′,ξ′ K(ξ - ξ′,ξ′) ) Kξ-ξ′,ξ′

(4)

We note that the population balance model considered here is restricted to conditions where multibody interactions can be neglected, which is obviously not the case beyond the gel point.

Figure 1. Typical values of the parameter κ*, as a function of the primary particle size in the case of fully destabilized conditions in water at 25 °C, for various values of the shear rate, for W ) 1 and R ) 0.1: (s) G ) 1 s-1, (- - -) G ) 10 s-1, (‚ ‚ ‚) G ) 100 s-1, and (- ‚ -) G ) 1000 s-1.

Furthermore, the model does not account for effects of the shear rate being nonuniform in time or space, which would be typically observed in realistic aggregating systems. For example, the growth of fractal clusters in aggregation suspensions leads to an increasing effective suspension viscosity, which might, in turn, influence local values of the shear rate as aggregation progresses. Such effects can be taken into account through a combination of population balance equations with a fluid dynamics simulation. Finally, the model neglects mechanisms such as sedimentation, breakage, or restructuring, which could affect aggregation and gelation kinetics in certain situations, especially if clusters grow sufficiently large so that inertial effects become significant. 2.3. Structure of AggregatessFractal Dimension. Because the PBE (eq 1) is a material balance, it is written in terms of mass, while the aggregation rate constants (K(ξ - ξ′,ξ′)) are expressed in terms of the sizes of the colliding clusters. Thus, it is necessary to introduce an appropriate relationship between the mass and size of the clusters. In the case of colloidal systems, it is well-established that aggregates constituted of larger number of primary particles form randomly shaped clusters that can be described by a fractal scaling relationship between the dimensionless mass ξ and the radius of gyration Rg,ξ:18

ξ ) kg

( ) Rg,ξ RP

df

(5)

where RP is the radius of the primary particles, kg is a scaling prefactor of the order of unity, and df is the fractal dimension. A similar relationship is also valid for the collision radius Rc,ξ, as well as for the hydrodynamic radius Rh,ξ.19 The latter is related to the diffusion coefficient Dξ in eq 3 through the StokesEinstein equation, so that

Rc,ξ ) RPξ1/df Dξ )

()

kBT 1 -1/df ξ 6πη RP

(6) (7)

In the scaling relationships for the collision and hydrodynamic radii, we use the scaling prefactors that are equal to unity, because possible variations of the prefactor value from unity can be included in the kernel prefactors W and R. 2.4. Dimensionless Form of the PBE: Master Curves. It is convenient to rewrite eq 1 into a dimensionless form using the dimensionless time,

τ)

3 RGφ0t 4π

(8)

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007 1711

and the dimensionless concentration,

x(ξ) )

n(ξ)VP φ0

(9)

where φ0 is the initial solid volume fraction and VP is the volume of primary particles (4/3πRP3). Note that φ0 ) n0VP, where n0 is the initial number concentration of primary particles. After substituting these dimensionless quantities into eq 1, we get

dx(ξ,τ) 1 ) dτ 2

∫1ξ-1 κ(ξ - ξ′,ξ′)x(ξ - ξ′,τ)x(ξ′,τ) dξ′ ∞ x(ξ,τ) ∫1 κ(ξ,ξ′)x(ξ′,τ) dξ′

M

x(ξ,t) ) (10)

where the dimensionless aggregation kernel is given by

κ(ξ - ξ′,ξ′) )

K(ξ - ξ′,ξ′) RGRP

3

) κ*((ξ - ξ′)-1/df + ξ′-1/df) ×

By considering the kernel expressions (eqs 2-4) together with the fractal scaling (eqs 5-7), the following expression for the dimensionless kinetic parameter κ* is obtained:

2kBT 3ηRWGRP3

(12)

This parameter represents the relative importance of the Brownian (2kBT/3ηW) aggregation mechanism, with respect to the shear (RGRP3) aggregation mechanism for primary particles, and it is controlled by the process operating conditions, such as temperature, fluid viscosity, primary particle size, stability ratio, and shear rate. From this equation, one can see that, when the system has a heterogeneous shear rate distribution or is approaching gelation (where the effective suspension viscosity increases and the local value of shear rate may vary in time), the value of κ* may vary in space or in time. Typical values of κ* for water at 25 °C are shown in Figure 1, as a function of the primary particle radius for various values of the shear rate. It is important to note that, similar to the static systems, where the value of the stability ratio W must be measured experimentally, the value of the aggregation rate prefactor (R) for shear aggregation also must be extracted from the appropriate experimental measurements. More details about the determination of R can be found in the work by Waldner et al.16 The aforementioned analysis shows that the dimensionless aggregation kernel is dependent only on two dimensionless parameters: the cluster fractal dimension df and the dimensionless kinetic parameter κ*, which, together, fully determine the time evolution of the dimensionless CMD. The same holds for any property of the distribution that can be expressed in terms of an arbitrary function of the dimensionless CMD, such as, for example, the various moments that can be determined from light scattering measurements.20 Therefore, we can conclude that the kinetics of aggregation can be represented through master curves (depending only on df and κ*) describing the evolution of appropriate dimensionless variables characterizing the CMD (for example those corresponding to measurables from light scattering) as a function of the dimensionless time τ.

Xi(τ)δ(ξ - ξi) ∑ i)1

(13)

and the transport equation for each class can be written as follows:

dXk(τ)

((ξ - ξ′)1/df + ξ′1/df) + ((ξ - ξ′)1/df + ξ′1/df)3 (11)

κ* )

2.5. Numerical Solution of the Population Balance Equation. Here, we consider two different methods for the numerical solution of the dimensionless PBE (eq 10), namely (i) the classes (or fixed-pivot) method with various grid choices,21-23 and (ii) the quadrature method of moments (QMOM).24,25 2.5.1. Classes Methods. Classes methods are based on discretizing the internal coordinate (here, the dimensionless mass) into a finite number of bins defined by pivots representing the discretized CMD. Each pivot has its characteristic mass ξi and all relevant quantities are based on this mass. The CMD is then approximated as follows:

)

dτ

1 k-1 k-1 2

∑ ∑ i)1 j)1

∞

κijXi(τ)Xj(τ) - Xk(τ)

κikXi(τ) ∑ i)1

(14)

where κij ) κ(ξi,ξj). To describe aggregation phenomena, it is often necessary to cover sizes ranging from a few nanometers (primary particles) to micrometers and even millimeters (final aggregates), so that one typically must use geometrically spaced pivots. However, because the performance of the classes methods is sensitive to the discretization grid used, we must test the effect of the computational grid (the number, as well as the spacing, of pivots) on the accuracy of the method. More details about the fixed-pivot method can be found in the work by Kumar and Ramkrishna.23 2.5.2. Method of Moments. The method of moments is based on solving the transport equations for the moments of the CMD, instead of solving the original population balance equations. The moment of order k of the CMD is defined as

µk(τ) )

∫0∞ x(ξ,τ)ξk dξ

(15)

Applying the moment transformation to the population balance equation, we obtain the transport equation for the kth-order moment, µk,:26 dµk(τ) dτ

)

1 2

∫

∞

1

ξk dξ

∫ ∫

ξ-1

1 ∞

1

κ(ξ - ξ′,ξ′)x(ξ - ξ′,τ)x(ξ′,τ) dξ′ -

ξkx(ξ,τ) dξ

∫

∞

1

K(ξ,ξ′)x(ξ′,τ) dξ′ (16)

To describe the evolution of the CMD through a closed set of equations, for its moments of order less than or equal to k, one must address the so-called “closure problem”.26 In the QMOM, this problem is solved using the quadrature approximation, as follows:24 M

∫1∞ x(ξ) f(ξ) dξ ≈ ∑ f(ξi) wi

(17)

i)1

where the quadrature approximation of order M is defined by M weights wi and M abscissas ξi. The corresponding moments of any order k are then calculated as follows: M

µk )

wiξik ∑ i)1

(18)

1712

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007

Finally, the quadrature weights and abscissas are expressed in terms of the corresponding moments (m0, ..., m2M-1), using the product-difference algorithm.27 The final transport equation for the kth-order moment µk is

dµk(τ) dτ

)

1

M

M

M

(19)

The important difference between the two previously discussed numerical methodssthe classes methods and QMOMsis the number of equations that ultimately must be solved. In the case of the classes method, the typical number of pivots is between 50 and 100 (and the same number of equations) to reach a good accuracy, whereas, in the case of QMOM, the number of nodes needed for comparable accuracy is typically 3-5 (which corresponds to 6-10 equations). Thus, using QMOM on one side significantly reduces the computational time; however, on the other side, we lose detailed information about the CMD, which is then characterized only by its moments. 2.6. Connection to Light Scattering. To validate population balance modeling by comparison with experiments, we need to choose suitable quantities that can be measured and calculated from the CMD predicted by the PBE. For this, it is required to have a reliable model for quantifying the behavior of a population of aggregates, with respect to such a selected quantity. In this work, as a form of example, we consider static light scattering (SLS) measurements, because sound models for optical properties of fine particles and fractal aggregates are available within the range of validity of the RDG theory.18 Using appropriate optical models, we can calculate the same average quantities for a given population of aggregates as those obtained from SLS measurements;5,20 that is, for example, the mean radius of gyration (Rg) and the mass-weighted-average mass, which is proportional to the zero angle intensity (I(0)). These quantities are, in fact, moments of the underlying cluster mass distribution of the population of aggregates and describe certain well-defined average properties of this population. The dimensionless mean radius of gyration is given by

x

∫ξ x(ξ)ξ2Fg,ξ2 dξ ∫ξ x(ξ)ξ2 dξ

(20)

where Fg,ξ is the dimensionless radius of gyration of the aggregate of mass ξ, which is defined as Fg,ξ ) Rg,ξ/RP and x(ξ) is the dimensionless number density of aggregates of dimensionless mass ξ. Based on the definition (eq 20), one can see that 〈Fg〉 is not a simple moment of the CMD and its order is dependent on the value of the fractal dimension df using eq 5. The zero angle intensity is given (up to a multiplicative constant) by

I(0) )

∫ξ x(ξ)ξ2 dξ ∫ξ x(ξ)ξ dξ

(21)

so that the dimensionless zero angle intensity (I(0)norm, scaled by the initial zero angle intensity) reduces to the second-order moment of the CMD:

I(0)norm )

I(0) ) I(0)t)0

φnorm )

M

(ξi + ξj)kκijwiwj - ∑ ∑ ξikκijwiwj ∑ ∑ 2 i)1 j)1 i)1 j)1

〈Rg2〉1/2 ) 〈Fg〉 ) RP

We will also consider the occupied volume fraction (φocc) and the corresponding normalized value (φnorm):

∫ξ x(ξ)ξ2 dξ

(22)

φocc ) φ0

V

∫ξ x(ξ)VPξ dξ ) ∫ξ x(ξ)ξ3/d dξ f

(23)

which is, again, not a simple moment of the CMD, but it is dependent on the value of the fractal dimension df. The occupied volume fraction is a quantity relevant for assessing the conditions that lead to the formation of colloidal gels,5,6 which occurs when a certain critical value of φocc is reached. 3. Results and Discussion 3.1. Analysis of the Monodisperse Systems. We consider first the approach presented by Bremer et al.,1 assuming that the cluster distribution can be approximated as monodisperse at all times. In this case, a closed-form expression of the time evolution of the cluster size has been derived for shear aggregation, which can be expressed in dimensionless form as follows:

τ)

df 1 (1 - Fdf-3) × 4 3 - df

(for df < 3)

(24)

where F ) R/Rp is the dimensionless aggregate radius. In the presence of shear aggregation alone, i.e., κ* ) 0, we see from eq 24 that (for df < 3) there exists a finite time τ* where F diverges to infinity:

τ* ) lim

Ff∞

df df 1 1 × (1 - Fdf-3) ) × 4 3 - df 4 3 - df (for df < 3) (25)

This will be referred to as the runaway time in the particular case of a monodisperse distribution of aggregates in the presence of purely shear aggregation, i.e., κ* ) 0. If also Brownian aggregation is present in the system, the time evolution of the cluster size can be determined numerically from the relation derived by Bremer et al.,1 using, again, the monodispersity assumption, which, in dimensionless form, reduces to

τ)

df 2

df-1

∫1F(τ) 2FF3 + κ* dF

(26)

Because of the fact that the dependency of the dimensionless cluster size F on the dimensionless time τ is determined only by the fractal dimension df and the kinetic parameter κ*, eq 26 can be regarded as the master curve that describes the aggregation kinetics of monodisperse distributions in the presence of both Brownian and shear mechanisms. By taking the limit for F f ∞, from eq 26, we can compute the runaway time τ*, as a function of the kinetic parameter κ*. The obtained results are shown in Figure 2 for various values of the fractal dimension df. As one can see, when the aggregation process is controlled by the shear mechanism (κ* < 1), the runaway time is only a function of df, according to eq 25. On the other hand, when the Brownian aggregation is gradually turned on by increasing the dimensionless kinetic parameter κ* (e.g., by reducing the shear rate G), the runaway time (in dimensionless form) becomes progressively shorter and is dependent on both df and κ*. After a certain transition period, this dependence of the τ* vs κ* plot in a log-log scale reaches a limiting slope equal to (df - 3)/3, as illustrated in Figure 3. The same scaling was also found by analytically solving eq 26

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007 1713

Figure 2. Comparison of runaway time (τ*), as a function of the kinetic parameter κ* for various fractal dimension (df) values calculated for a monodisperse cluster mass distribution (CMD) by integration of eq 26. Increasing the value of κ* increases the relative importance of the Brownian mechanism, compared to the shear aggregation mechanism. The difference in the df values between the lines is equal to 0.2 ((s) df ) 1.8, (- - -) df ) 2.0, (- ‚ -) df ) 2.2, (- ‚ ‚ -) df ) 2.4, and (- - -) df ) 2.6).

Figure 3. Slope of the τ* vs κ* from the region of κ* > 104 shown in Figure 2, as a function of the fractal dimension df. The points are wellapproximated by a linear relationship, with a slope equal to (df - 3)/3.

in the limiting case when the dimensionless kinetic parameter κ* approaches infinity, as shown in the Appendix. Let us now consider two values of κ*: 1 , κ*1 < κ*2, corresponding to two different values of the shear rate (G1 > G2). Figures 2 and 3 show that, while the dimensionless runaway time τ* decreases with increasing κ*, the dimensional runaway time t* ≈ τ*/G (see eq 8) increases by a factor of t*2/t*1 (t*2/t*1 ) G1τ*2/G2τ*1 ) κ*2τ*2/κ*1τ*1). The opposite trend that is observed for the dimensional and dimensionless runaway times, respectively, is due to the fact that the dimensionless runaway time τ* decreases slower with κ* than the inverse of κ*, so that the product of κ* and τ* actually decreases with increasing κ*. We note that, in the limit of a pure Brownian aggregation mechanism, where κ* becomes infinitely large (aggregation proceeds under static conditions), there is no runaway, because the dimensional runaway time becomes infinitely large. Therefore, in the limit of pure Brownian aggregation, it is appropriate to use another choice of the dimensionless time, namely τκ*, as was previously shown by Sandku¨hler and co-workers.5-7 3.2. Analysis of Polydisperse Systems. 3.2.1. Validation of Numerical Methods. According to the expression of the runaway time given by eq 25 for monodisperse distributions, we expect that, similarly, for a polydisperse system, the moments of the CMD of an order of higher than 1 also diverge at a finite time, which becomes progressively shorter for decreasing cluster df. This faster growth rate of the aggregates with lower df makes the numerical integration of the PBEs difficult, which is investigated in the following discussion using the classes method with various grids23 and the quadrature method of moments.24,25 Fractal dimensions ranging from df ) 1.8 (typical for diffusionlimited cluster aggregation in static conditions) up to ) 2.6 are considered, to cover the typical df values that are observed in experiments, for both static and shear conditions.5-7,19,28

Figure 4. Effect of the maximum size of the geometrical computational grid (sectional method) on the calculated time evolution of two moments of the CMD for df ) 1.8 and pure shear kernel: (0) 19 grid points, (O) 21 grid points, (4) 23 grid points, and (3) 25 grid points.

The effect of the length of the discretization grid on the time evolution of two moments of the CMD, i.e., 〈Fg〉 and I(0)/I(0)t)0, computed through the classes method, is illustrated in Figure 4a-d for df ) 1.8. Here, we have used a classical version of the classes method with a simple geometrical binary grid,21 which is also known as the sectional method,22 where the dimensionless cluster mass is discretized using pivots with masses equal to 2i-1, where i is the corresponding pivot number. From Figures 4c and 4d, we can see that, in the early stages of the aggregation process, the same values for both moments are obtained with the different lengths of the discretization grid. However, as time proceeds, the aggregates grow and we see that the computed values of the two moments change with the length of the computational grid. In particular, the time when the moments start to diverge is different for the different grid lengths, as well as for the two moments. The time of numerical divergence becomes shorter as the maximum pivot mass increases (i.e., the grid becomes longer) as is shown in Figures 4a and 4b. This seemingly counter-intuitive effect is, in fact, a numerical artifact that is due to the delay of the runaway behavior by stopping the self-accelerating broadening of the distribution earlier for shorter grids. To investigate the role of the grid size further, in the following discussion, we use a series of progressively finer computational grids, where we use linear scaling for the first n pivots, followed by a geometrical grid with a grid factor equal to (n + 1)/n. We present results for n ) 3 and n ) 16, together with those corresponding to the previous case of binary grid (sectional method) with n ) 1. To assess the effect of the grid spacing only, we chose the number of pivots such that the last pivot size was (approximately) the same for all tested cases. The obtained results are summarized in Figure 5 for the highest moment of the CMD considered in this work, i.e., the dimensionless mean radius of gyration 〈Fg〉. It is observed that, as the grid becomes finer (larger n values), the time when the solution diverges becomes larger and the CMD reaches the end of the grid later. From this analysis, we conclude that, because the computational grid would need to be extremely fine to obtain reliable results, the classes method is not suitable for computing the critical time where the CMD moments runaway. A better alternative is the QMOM,24 which is obviously not restricted to a finite discretization grid. On the other hand, using this method, we lose detailed information about the CMD. A comparison of the classes method with different discretization grids with the QMOM with different number of nodes is

1714

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007

Figure 5. Effect of the computational grid structure on the time evolution of the dimensionless radius of gyration 〈Fg〉 calculated for df ) 1.8 and pure shear kernel: (0) pure geometric grid; (O) first 3 linear grid points, followed by a geometric grid with a grid factor of 1.3333; and (4) first 16 linear grid points, followed by a geometric grid with a grid factor of 1.0625. Figure 7. Comparison of the time evolution of normalized zero angle intensity I(0)/I(0)t)0 and dimensionless radius of gyration 〈Fg〉 calculated for pure shear aggregation (κ* ) 0) for various df values: (s) df ) 1.8, (- - -) df ) 2.0, (- ‚ -) df ) 2.2, (- ‚ ‚ -) df ) 2.4, and (- - -) df ) 2.6. Arrows indicate the runaway time τ* corresponding to the vertical asymptote in the time evolution of the considered moment (illustrated for df ) 2.2).

Figure 6. Comparison of the time evolution of the normalized zero angle intensity and the dimensionless radius of gyration computed for a pure shear kernel with df ) 2.2, using the fixed-pivot method with different computational grids ((0) pure geometric grid; (O) first 3 linear grid points, followed by a geometric grid with a grid factor of 1.3333; and (4) first 16 linear grid points, followed by a geometric grid with a grid factor of 1.0625) and the quadrature method of moments (QMOM) with different numbers of nodes ((- - -) N ) 2, (- ‚ -) N ) 3, and (s) N ) 5).

presented in Figure 6 for the two moments considered previously in Figure 4a-d and fractal dimension df ) 2.2. It is observed that both methods converge to the same solution as the discretization becomes finer and that the QMOM provides reliable results when the number of nodes is g5, which is a similar conclusion to that presented by Diemer and Ehrman,29 obtained for a simulation of flame-synthesized aerosol particles. This trend was observed for all fractal dimensions considered in this work (i.e., df ) 1.8-2.6). In the following discussion,we compute the CMD moments with QMOM using five nodes and, when needed, the detailed CMDs with the classes method with appropriately fine grid, as validated against the QMOM. 3.2.2. Gelation and Runaway in Pure Shear Aggregation. The time evolutions of the normalized zero angle intensity and the dimensionless mean radius of gyration calculated for pure shear aggregation (κ* ) 0) are shown in Figure 7 for various df values. It is observed that, as the fractal dimension increases, the dimensionless runaway time τ* (indicated by the arrow for the case of df ) 2.2) increases and runaway occurs at larger aggregate sizes. The fact that both quantities considered here diverge at the same time τ* value (for a given value of df) can be explained from their definition (see eqs 20 and 21). As one can see, both are calculated as the ratio between a higher moment and a lower moment; however, in both cases, the order of the moment in the numerator is g2 and, therefore, as runaway

Figure 8. Comparison of the time evolution of the normalized occupied volume fraction calculated for pure shear aggregation (κ* ) 0) for various values of df: (s) df ) 1.8, (- - -) df ) 2.0, (- ‚ -) df ) 2.2, (- ‚ ‚ -) df ) 2.4, and (- - -) df ) 2.6. Arrows indicate the gelation time τG (using φ0 ) 0.1) corresponding to the point where φocc ) 0.5 and runaway time τ* (illustrated for the case where df ) 2.2).

occurs, they diverge to infinity. To investigate the occurrence of gelation in Figure 8, we investigate the time evolution of the normalized occupied volume fraction, given by eq 23, for the same df values considered in Figure 7. From eq 23, one can observe that the occupied volume fraction is also a moment of the CMD, and, therefore, it also diverges as runaway occurs, as shown in Figure 8 and indicated by the arrow in the case of df ) 2.2. From these results, we can compute the gelation time (τG), following the work of Sandku¨hler et al.5,6 for purely Brownian systems, as the time where the occupied volume fraction (φocc) is equal to 0.5. An example of the gelation time τG for df ) 2.2 and φ0 ) 0.1 is also indicated in Figure 8. The so-obtained values, which are clearly a function of the initial volume fraction φ0, are summarized in Figure 9 for various df values. It can be observed that, for sufficiently low values of the initial solid volume fraction, gelation is actually driven by runaway in the sense that the time when φocc reaches the value of 0.5 practically coincide with the runaway time, i.e., τG = τ*, when divergence of all higher moments occurs. We note that τ* for various df values is independent of the initial solid volume fraction and its value is equal to the plateau value for φ0 f 0, which increases with increasing df (see Figure 7). On the other hand, for φ0 g 0.1, the gelation time is significantly shorter than the runaway time (τG , τ*) and, as df increases, this difference progressively increases. We can then conclude that, for large values of the initial solid volume fraction (g0.1), which is the condition of typical industrial interest, the concept

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007 1715

Figure 9. Values of the dimensionless time τG when the occupied volume fraction reaches a value of φocc ) 0.5 (typical for volume filling under stagnant conditions), as a function of the initial solid volume fraction φ0 calculated for pure shear aggregation (κ* ) 0) for various df values: (b) df ) 1.8, (4) df ) 2.2, and (9) df ) 2.6. Table 1. Comparison of the Runaway Time τ* Calculated with eq 25, Assuming a Monodisperse Distribution of Aggregates, with the Results Obtained Using the Full Population Balance Equation (PBE) for Various Fractal Dimension Values

Figure 10. Comparison of the time evolution of the normalized zero angle intensity (I(0)/I(0)t)0), as a function of dimensionless time for various values of κ* and two different df values: (- - -) df ) 1.8 and (s) df ) 2.2.

Runaway Time, τ* fractal dimension, df

monodisperse model, eq 25

full PBE

1.65 1.8 2.0 2.2 2.4 2.6

0.306 0.375 0.500 0.688 1.000 1.625

0.082 0.133 0.226 0.369 0.606 1.109

of volume filling by the aggregates, characterized by φocc ≈ 0.5 as in static systems,5,6 is an appropriate indicator of the occurrence of gelation. In fact, after gelation occurs, it is obvious that the system is no longer correctly described by the population balance equation (eq 1) and, consequently, the runaway behavior that they predict at later times is not realistic. To define the effect of the CMD polydispersity on the runaway, the dimensionless runaway time obtained from the full PBE solution are compared in Table 1 with those computed from the monodisperse model (eq 25), in the case of pure shear aggregation. We can see that the runaway time predicted by the model, assuming monodisperse distribution, is larger by a factor of 2-3, compared to the numerical solution of the PBE, accounting for the polydispersity of the CMD. This polydispersity effect becomes smaller with increasing df. A similar trend was found also by Gonzalez et al.,2 although their numerical results should be taken with care, because they are dependent on the adopted computational grid (see section 3.1). 3.2.3. Role of the Aggregation Mechanism. Let us now consider aggregating systems where both Brownian and shear mechanisms are present and their relative importance is determined by the value of the parameter κ*. For this, we first analyze the behavior of I(0)norm, i.e., the second-order moment of the CMD, as a function of time for various values of the parameter κ* for two values of fractal dimension (df ) 1.8 and 2.2), as shown in Figure 10. It is observed that, in all cases, I(0)norm diverges at finite time values, which indicates the occurrence of runaway. We see that runaway occurs faster for lower df values, as well as for larger contributions of the Brownian aggregation mechanism (i.e., for increasing values of κ*), with respect to dimensionless time τ, although the absolute value of the time when runaway occurs increases as κ* increases (∼33 times as we increase the value of κ* from 5.6 to 560 by reducing the shear rate; see section 3.1). In addition, with increasing contribution of the Brownian aggregation mechanism, we can clearly see the change in the shape of the time evolution of I(0)norm at the early stage of the process. For example, for κ* ) 56 and 560, we see that aggregation is

Figure 11. Time evolution of the normalized CMD calculated for two df values ((a) df ) 1.8, (b) df ) 2.2) corresponding to the simulations examined in Figure 10. From top to down, the relative importance of Brownian aggregation, with respect to shear aggregation, increases: (a1, b1) κ* ) 0, (a2, b2) κ* ) 5.6, (a3, b3) κ* ) 56, and (a4, b4) κ* ) 560.

initially controlled by a Brownian mechanism and, therefore, the aggregates grow with rates that decrease over time. However, after the aggregates have grown sufficiently large, the pure shear mechanism takes over and the shape of the time evolution of I(0)norm becomes similar to that for κ* ) 0. To better understand the effect of the aggregation mechanism, in Figure 11, we plot the CMDs at various times, corresponding to the transitions to runaway shown in Figure 10. The two columns (a1-a4) and (b1-b4) correspond to two different fractal dimension values (df ) 1.8 and 2.2, respectively), while going down in the column, a larger contribution of the Brownian aggregation mechanism is introduced, i.e., increasing the κ* value. In each figure, the CMDs are shown for five equally spaced normalized dimensionless times: τ/τ* ) 0.14, 0.28, 0.42, 0.56, and 0.7. It is worth mentioning that, in all these simulations (that is, for time values of τ/τ* e 0.7, the results for 〈Fg〉 and

1716

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007

Figure 12. Time evolution of the normalized CMD calculated for two df values ((a) df ) 1.8 and (b) df ) 2.2). Top panel shows combined shear and Brownian aggregation (κ* ) 560), and the bottom panel shows pure Brownian aggregation (κ* ) ∞).

I(0)norm calculated with both numerical methods mentioned previously coincide. In addition, to compare their relative shape, all the CMDs have been normalized by the concentration of primary particles at the corresponding time, i.e., x(i,τ)/x(1,τ). As expected, the distribution broadens as time increases. By comparing the relative shapes of the CMDs, it is observed that, in the case of df ) 1.8, the normalized CMD does not reach any limiting shape before the moments start diverging, or at least up to τ/τ* ) 0.7. If we continue in the simulation at longer times, the right-hand side of the CMD exhibits runaway, i.e., the concentration of larger aggregates starts to grow very fast. The relative importance of Brownian aggregation increases as κ* increases, and a maximum appears in the CMD. This is a characteristic feature of distributions obtained for fully destabilized conditions in stagnant systems (diffusion-limited clustercluster aggregation, DLCA), where Brownian aggregation is the only aggregation mechanism. Under these conditions, aggregation between small and large aggregates prevails over aggregation between large ones, and, therefore, a maximum develops in the CMD. The situation is different in the case of df ) 2.2, where the distribution attains, at least in the intermediate aggregate mass range, a certain limiting shape for a certain period of time (before undergoing runaway), which exhibits a linear behavior in the log-log plot with a slope of ∼2. Similar to the case of df ) 1.8, also here, with increasing importance of Brownian aggregation, with respect to the shear mechanism, a maximum in the CMD appears, but the limiting shape feature in the intermediate aggregate range remains. This means that, in this case (df > 2), the two sides of the CMD are controlled by different mechanisms: in the left part, Brownian aggregation prevails, leading to its characteristic maximum, whereas in the right part, shear aggregation dominates. The different behavior observed for the two df values can be explained by considering that, for a given mass, aggregates with lower df values are larger and, therefore, aggregate faster. Accordingly, the onset of runaway occurs at shorter dimensionless times (or smaller aggregate masses), as is clearly apparent in Figure 10. In Figure 12, we compare the CMDs calculated for pure Brownian aggregation (bottom) with those calculated for combined shear and Brownian aggregation (top), for two values of the fractal dimension: df ) 1.8 (a1, a2) and df ) 2.2 (b1, b2). We see that, as time proceeds, the maximum in the CMD becomes more pronounced with an increasing importance of

Brownian aggregation, with respect to shear aggregation, and the df value becomes larger. It is worth reiterating that, under pure Brownian aggregation conditions, none of the CMD moments diverge, and the distribution remains rather narrow and grows continuously, over time achieving a self-similar shape, typical of DLCA under stagnant conditions.5-7 However, when shear aggregation is present, the larger clusters are affected more by shear aggregation than by Brownian aggregation, and, therefore, when the clusters reach a sufficiently large size, the shear mechanism takes over and the distribution rapidly broadens and increases its polydispersity and eventually runaway occurs. On the other hand, the existence of a maximum in the CMD indicates that, also in this case, the aggregation of small particles is controlled by the Brownian mechanism. To complete the picture about the role of aggregation mechanism and fractal dimension on the dimensionless runaway time, in Figure 13, we present the comparison of the dimensionless runaway times calculated for combination of Brownian and shear aggregation mechanisms obtained by numerical solution of the full PBE for various values of the kinetic parameter κ* (eq 12). It is observed that the qualitative trend of the calculated curves is similar to that calculated by approximating the distribution of aggregates to be monodisperse (shown in Figure 2). When the aggregation process is controlled by the shear mechanism (κ* < 1), the dimensionless runaway time is only a function of fractal dimension df, which is the same conclusion reached in the case of monodisperse distributions (eq 25). On the other hand, when Brownian aggregation is gradually turned on by increasing the dimensionless kinetic parameter κ*, the dimensionless runaway time becomes progressively shorter and is dependent on both df and κ*. We note that, also for polydisperse CMD, in the limit when the process is controlled by Brownian mechanism, the absolute runaway time increases, because of the dimensionless scaling used, as was discussed in section 3.1. Similar to the case of monodisperse distribution, also here, after a certain transition period, this decay reaches a limit and the slope of the τ* vs κ* in the log-log scale becomes equal to (df - 3)/3, as shown in Figure 3 and in the Appendix. Comparing now the results presented in Figures 2 and 13 quantitatively, it is observed that the dimensionless runaway time is much longer in the case of monodisperse distributions than for polydisperse CMDs, which is consistent with our previous observations. As was shown in Figure 9, for relatively high values of the initial solid volume fraction, the volume filling (gelation) can occur before runaway, i.e., τG < τ*. Therefore, in Figure 14, we plot the dimensionless time when the absolute value of the occupied volume fraction is equal to 0.5, τG indicated by symbols, as a function of κ* for three selected fractal dimension values (df ) 1.8, 2.2, 2.6) and assuming initial solid volume fractions equal to 0.01 (top) and 0.1 (bottom). The lines represent the runaway time in dimensionless form (τ*), shown in Figure 13. It is observed that, for small values of κ*, τ* = τG, which means that runaway starts before gelation and the consequent very fast increase of the aggregate sizes leads almost immediately to gelation. Accordingly, the gelation time is not a function of κ* but rather only of df. As the relative importance of the Brownian aggregation mechanism, compared to the shear aggregation mechanism, increases (in this particular case, by decreasing the value of the shear rate G), the gelation time becomes significantly shorter than the runaway time. This occurs at lower κ* values for the larger initial solid volume fraction (φ0 ) 0.1) than for the smaller one (φ0 ) 0.01). It is also

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007 1717

Figure 13. Comparison of runaway time, as a function of the kinetic parameter κ*, for various df values calculated with the full population balance equation (PBE) model. Increasing the value of the kinetic parameter κ* increases the relative importance of the Brownian aggregation mechanism, with respect to the shear mechanism. The difference in the df values between the lines is equal to 0.2 ((0) df ) 1.8, (1) df ) 2.0, (O) df ) 2.2, ([) df ) 2.4, and (4) df ) 2.6).

Figure 15. Critical values, κ*c, represented by the arrow, separating the region of operating conditions where gelation occurs as a consequence of runaway (represented by the solid line) from that where it occurs before runaway (represented by the dashed line); df ) 2.2 and φ0 ) 0.01.

Figure 16. Phase plane κ*c vs φ0 for various df values: (b) df ) 1.8, (4) df ) 2.2, and (1) df ) 2.6. The region below the lines indicates operating conditions where gelation occurs as a consequence of runaway, whereas the region above the lines indicates conditions where gelation occurs before runaway.

Figure 14. Comparison of the runaway time τ*, shown in Figure 13 (lines), with the gelation time τG, as a function of the kinetic parameter κ* for various df values calculated with the full PBE model (df ) 2.6 (triangles), df ) 2.2 (circles), df ) 1.8 (squares)) and two initial solid volume fraction (φ0) values.

observed that this deviation appears sooner for lower df values. After a certain transition period, one can observe that decay in τG, with increasing κ*, starts to be linearly proportional (in a log-log scale) to κ* for all investigated df values, as shown in Figure 14. Because the increase in the relative importance of Brownian aggregation, compared to shear aggregation, was, in this particular case, realized by decreasing the value of shear rate (other parameters in eq 12 remained unchanged), such a linear decay is caused by the fact that the absolute value of the time when we reach the volume filling conditions is constant and no contribution of shear is necessary to form a gel. Consequently, increasing the value of κ*, in this case, by reducing the shear rate, leads to the proportional reduction of the dimensionless time τ (which is proportional to the shear rate G). Based on the aforementioned analysis, it is clear that, depending on the value of the solid volume fraction, two distinct regions exist in the κ* values: one where aggregates fill the

space before runaway occurs and another where volume filling occurs as a consequence of runaway. In Figure 15, we illustrate this behavior, where the critical κ*C value separating these two regions is indicated by an arrow for a particular fractal dimension and initial solid volume fraction (i.e., df ) 2.2 and φ0 ) 0.01). The solid line represents the runaway time τ* and the dashed line represents the gelation time τG, where the available volume was filled. Such critical values of κ*C are shown in Figure 16 , as a function of the initial solid volume fraction (φ0) for three fractal dimension values. The region below the lines represents conditions where the volume filling occurs as a consequence of runaway, whereas the area above the lines represents conditions where gelation occurs before runaway. As one can see, with decreasing the fractal dimension or increasing the initial solid volume fraction, the boundary between runaway and volume filling region moves toward smaller values of κ*. Because we have already discussed and compared runaway times determined from monodisperse and full CMD models in Table 1, we can now compare the master curves as a function of the time scaled by the runaway time (both in dimensionless form) τ/τ*, calculated from the full PBE model and the monodisperse model (see eqs 25 and 26). Such a comparison is illustrated in Figure 17, in terms of the normalized radius of gyration 〈Fg〉 for a fractal dimension df ) 2.2 and two kinetic parameter values: κ* ) 0 and 560. We note that, although the analytical solution, assuming monodisperse CMD, differs by a factor of 2-3 in the dimensionless runaway time (compared to the numerical PBE solution), the overall shape of the scaled master curves is similar, although, as expected, the evolution of both moments of the CMD is faster for polydisperse distributions, compared to monodisperse distributions. The same behavior has been observed for other df values. 3.2.4. Time Evolution of the Moment Ratio. Until now, we have discussed our results in terms of time evolution of

1718

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007

Figure 17. Comparison of the time evolution of the normalized radius of gyration 〈Fg〉 calculated using the analytical equations for monodisperse distributions ((- - -) eqs 25 and 26 and (s) the full PBE model) for df ) 2.2: (a) pure shear aggregation (κ* ) 0) and (b) combined shear and Brownian aggregation (κ* ) 560).

Figure 18. Time evolution of the moment ratio I(0)norm/〈Fg〉 for various κ* values and df ) 1.8: (s) κ* ) 0, (- - -) κ* ) 5.6, (- ‚ -) κ* ) 56, and (- ‚ ‚ -) κ* ) 560.

Figure 19. Time evolution of the moment ratio I(0)norm/〈Fg〉 for various κ* values and df ) 2.2: (s) κ* ) 0, (- - -) κ* ) 5.6, (- ‚ -) κ* ) 56, and (- ‚ ‚ -) κ* ) 560.

individual moments of the CMD, whose qualitative shapes are rather similar to each other. Some more insight into the evolution of the CMD during shear aggregation can be achieved by considering a ratio of two different moments of the distribution. Such a ratio provides information about the broadness and shape of the CMD. Because both I(0)norm and 〈Fg〉 can be obtained from experiments (light scattering measurement), let us consider the ratio of these two moments. As was shown previously in the context of Figure 11, the relative importance between shear and Brownian aggregation significantly affects the shape of the CMD and, consequently, also the runaway time. Therefore, in Figures 18 and 19, we show the ratio I(0)norm/〈Fg〉 as a function of the dimensionless time for various values of kinetic parameter κ*, for two values of the fractal dimension (df ) 1.8 and 2.2, respectively), corresponding to the same CMDs illustrated in Figure 11. The

analysis of the moment ratio (see eqs 21 and 20) shows that, at longer times, where ξi f ∞, the moment ratio diverges to infinity. As can be seen from Figure 18, the time evolution of the moment ratio for df ) 1.8 at the beginning of the process is strongly dependent on the kinetic parameter κ* and either monotonically decreases in the case of pure shear aggregation (κ* ) 0), remains rather constant in the case of relatively weak contributions of the Brownian mechanism (κ* ) 5.6), or increases in the case of significant contributions of the Brownian mechanism (κ* ) 56 and 560). In the case of df ) 2.2, the moment ratio monotonically increases for all values of κ*, as shown in Figure 19. We note that, for df ) 1.8, no limiting shape of the CMD was attained for all values of κ* before divergence, whereas in the case of df ) 2.2, the limiting shape of the CMD was achieved for all values of κ*, as it can also be observed in Figure 11. In the case of pure shear aggregation and low df values, the distribution broadens substantially, because aggregation between large aggregates is favored over that between small and large aggregates. Because the mean radius of gyration (eq 20) is strongly affected by the presence of large aggregates, compared to the zero angle intensity (see eq 21), which characterizes the average mass of the aggregates, even small concentrations of large aggregates will significantly increase the value of 〈Fg〉. On the other hand, because of their low concentration, their contribution to the increase of I(0) is negligible, which finally leads to the low values of the moment ratio. In the case of df ) 2.2, the slower pace of the broadening of the CMD will lead to a situation where large aggregates contribute not only to the increase of 〈Fg〉 but their concentration is large enough to influence also the value of I(0), which leads to the appearance of a certain limiting shape of the CMD. Increasing the relative importance of the Brownian aggregation, which favors aggregation between large and small aggregates, leads to a rather narrow CMD. Therefore, this can affect I(0) more then 〈Fg〉, leading to an increase of the value of the moment ratio. After the largest aggregates reach the size which can be affected by the shear aggregation mechanism, a scenario similar to that described previously appears and the CMD starts to broaden and, therefore, the moment ratio decreases. Note that the above-mentioned observations can be used to estimate the initial df value in aggregating dispersions of primary particles.16,30 For example, if we observe experimentally that the value of the moment ratio I(0)norm/〈Fg〉 monotonically decreases in the early stages of aggregation, we can conclude that the Brownian mechanism does not have a dominant role and, at the same time, the initial df value of the aggregates is