Agglomeration Process Modeling Based on a PDE Approximation of

ARTICLE pubs.acs.org/IECR

Agglomeration Process Modeling Based on a PDE Approximation of the Safronov Agglomeration Equation Andrey V. Bekker and Iztok Livk* CSIRO Process Science and Engineering and Light Metals Flagship, Parker CRC for Integrated Hydrometallurgy Solutions, Waterford WA 6152, Australia ABSTRACT: A one-dimensional dynamic partial differential equation (PDE) agglomeration model is derived based on the continuous Safronov agglomeration equation. A regularized PDE agglomeration model, represented by a set of convectionreaction-diffusion PDEs, can be solved within a standard adaptive-mesh implicit numerical framework that does not require additional quadrature assumptions to evaluate the aggregation integral. The PDE agglomeration model is solved numerically using a general Newton’s-method-based implicit Galerkin finite-element algorithm. The applied algorithm uses an automatic Gear-type time step and nonuniform adaptive-mesh strategies, which aids solution convergence. A numerical solution of the model for an agglomeration degree of 99.9% closely matches an asymptotic analytic solution of the original Safronov equation, which confirms the accuracy of the numerical procedure used. It is also shown that the number density function predicted by the new PDE agglomeration model satisfactorily agrees with the analytic solution of the Smoluchowski agglomeration equation. For small particle sizes and first and zeroth full moments, the two models give similar solutions. However, for larger particle sizes and the second full moment, the difference between the two models increases with increasing degree of agglomeration. Industrially important gibbsite agglomeration is used as a case study to demonstrate the application of the new numerical approach for agglomeration modeling.

1. INTRODUCTION Agglomeration, aggregation, coagulation, and coalescence are some of the synonyms of a generic term for processes that enlarge the size of an entity. These processes have been shown to take place in a wide range of physical systems, such as cosmic object formation, aerosol evolution, particular and polymer reactions, and economic and sociological phenomena. Agglomeration of particles to form larger entities is also a common phenomenon observed in different crystallization systems. Numerical agglomeration models can be employed to describe the distribution of agglomerates in such systems by size and time. A mathematical model of agglomeration was first formulated in a discretized form by Smoluchowski.1 According to Dubovski,2 the main assumptions of the Smoluchowski agglomeration model are (i) a large number of chaotically moving particles, (ii) a spatially homogeneous and unlimited system, (iii) a rarefied system with pair particle interactions only, and (iv) a mean interparticle collision time (microscopic interaction time) that is substantially shorter than the time required for changes in particle size distribution. Despite the simplicity of the assumptions of the Smoluchowski agglomeration model, the resulting mathematical model is nonlinear, and a general analytic solution is not available.3 With the advent of new numerical methodologies and increased computing power, the complexity of population-balanceequation- (PBE-) based models, incorporating the Smoluchowski agglomeration equation, has been steadily growing. Such models can involve multiple internal and external coordinates,4 population balance coupled with computational fluid dynamics (CFD),5 and population balance equations with nonlinear stiff kinetics.6 Numerical solutions of such models are commonly attained by specialized algorithms that involve lengthy evaluations of the integral part of the Smoluchowski equation. Published 2011 by the American Chemical Society

Considerable progress in increasing the efficiency of the numerical algorithms used to solve the Smoluchowski agglomeration equation has been made recently. This equation can now be solved in its discrete form using a discretized population balance technique7 or in its continuous form using a finiteelement method (FEM).8 Both the discrete and continuous solution methods can involve explicit or implicit integration schemes with respect to time. In the second case, the Jacobian matrix of the Smoluchowski equation is of full rank. A semi-implicit discrete scheme for solving the Smoluchowski equation was presented by Jacobson and Turco.9 Their scheme results in a sparse band-diagonal Jacobian matrix. When using a fully implicit scheme for a nonlinear convolution integral, such as the Smoluchowski equation, it is not possible to construct a sparse band-diagonal Jacobian matrix, as required to use a commercially available implicit adaptive-mesh partial differential equation (PDE) solver.10 FE methodologies, developed previously using the Smoluchowski agglomeration equation, employ fixed-mesh solvers,6 which are less computationally efficient and require an additional quadrature assumption to evaluate the aggregation integral at each time step. To overcome these problems, a new PDE agglomeration model is formulated in this work based on the Safronov aggregation equation. The main advantage of this approach is that the aggregation integral is replaced by a PDE formulation that can be incorporated into the FEM solution framework. The use of the Safronov equation for modeling coagulation systems, and its relation to the Smoluchowski equation, was previously discussed by Dubovski.11 Received: September 20, 2010 Accepted: December 23, 2010 Revised: December 9, 2010 Published: February 16, 2011 3464

dx.doi.org/10.1021/ie101933r | Ind. Eng. Chem. Res. 2011, 50, 3464–3474

Industrial & Engineering Chemistry Research

ARTICLE

The objective of this work is to approximate the nonlinear integro-differential Safronov agglomeration model by a continuous reaction-convection-diffusion PDE operator, with the aim of solving the model with an implicit FEM solver using a general Newton’s-method-based Galerkin finite-element adaptive-mesh algorithm.10 A numerical solution of the newly developed PDE agglomeration model is compared to an asymptotic Safronov equation solution and a special-case analytic solution of the Smoluchowski equation. Gibbsite agglomeration is used as a case study to demonstrate the suitability of the new approach for an industrially important particle agglomeration system. The validation of the newly developed PDE agglomeration model against experimental gibbsite precipitation data is reported in a related article.12 Figure 1. Exponential number density functions (a = 0.1) for 0%, 99.9%, and 99.9999% agglomeration.

2. SMOLUCHOWSKI AGGLOMERATION EQUATION 2.1. Continuous Smoluchowski Agglomeration Equation.

A discrete agglomeration model, represented by a system of ordinary differential equations, was initially derived by Smoluchowski.1 Its continuous version, in the form of an integro-differential equation, is credited to M€uller13 and was solved for a special case analytically by Schumann.14 The Smoluchowski agglomeration equation describes the evolution of the number density function with time. Here, we focus on a continuous form of the Smoluchowski equation, which, for certain initial conditions, has an analytic solution and can serve as a benchmark for our numerical investigation. In the chemical engineering literature,4,15,16 the Smoluchowski equation is commonly referred to as general dynamic PBE for pure aggregation. For a one-dimensional internal coordinate system, a dimensionless transformation of the continuous Smoluchowski agglomeration equation with a constant kernel can be written as t∈½0, tend ,

x∈½xmin , xmax , y∈½xmin , xmax Z Z xmax ∂f ðx, tÞ K x ¼ f ðx - y, tÞ f ðx, tÞ dy - Kf ðx, tÞ f ðy, tÞ dy ∂t 2 xmin xmin ∂f ðx, tÞ ∂f ðx, tÞ ¼ ¼ 0 ∂x ∂x xmin

ð1Þ

xmax

f ðx, tÞjt ¼ 0 ¼ f0 ðxÞ ¼ expð - axÞ

where f(x,t) is the dimensionless number density of the particle size distribution, x and y represent the dimensionless particle size coordinate, and t is the dimensionless time. Parameters in the system are the dimensionless final time, tend; the minimum and maximum size of particles, xmin and xmax, respectively; the constant agglomeration kernel, K; and the width of the initial particle size distribution (PSD), a. The Smoluchowski equation with a constant agglomeration kernel describes a case where the probability of successful binary particle interactions is independent of particle size and process conditions. Boundary conditions represent impermeable barriers for the particle number density. An exponential initial PSD is assumed. 2.2. Analytical Solution of the Smoluchowski Agglomeration Equation. If xmin= 0 and xmax = ¥, the model of eq 1 has an analytic solution3,14 of the form 0 1 1 B ax C fSM ðx, tÞ ¼ f ðx, tÞ ¼ exp@Kt A Kt 2 1þ 1þ 2a 2a

To enable comparison of the agglomeration model results of this work with those of Hounslow,17 the agglomeration degree (Iagg) is defined as μ0, f ðtÞ ð3Þ Iagg ðtÞ ¼ 1 μ0, f ðt ¼ 0Þ were μ0,f(t) is the full zero moment of f(x,t) according to Z xmax μi, f ðtÞ ¼ xi f ðx, tÞ dx

ð4Þ

xmin

In practical terms, an agglomeration degree of 0.999, for example, means that, on average, an agglomerate consists of 1000 primary particles. A 99.9% agglomeration process thus reduces the initial number of particles by a factor of 1000, resulting in a 10-fold increase in the particle mean size. 2.3. Number Density and Cumulative Moments As Functions of the Agglomeration Degree. The agglomeration degree as a function of time can be analytically calculated from eqs 2 and 3 for the case of xmin= 0 and xmax = ¥ as 1 ð5Þ Iagg ðtÞ ¼ 1 Kt 1þ 2a Using the expression for the agglomeration degree in eq 3, the analytic Smoluchowski solution in eq 2 can be rewritten as a function of the agglomeration degree ð6Þ f ðx, Iagg Þ ¼ ð1 - Iagg Þ2 exp½- að1 - Iagg Þ The solution of eq 6 is presented in Figure 1 for a = 0.1 and two different values of Iagg, 0.999 and 0.999999. Compared to the value for the initial system, the number densities for Iagg = 0.999 and Iagg = 0.999999 are 6 and 12 orders of magnitude lower, respectively. At the same time, the mean particle size increases by 3 and 6 orders of magnitude, respectively, as can also be observed in Figure 1. The cumulative number of particles, μ0(x,t), and the cumulative volume of particles, μ1(x,t), are represented by the cumulative zeroth and first moments of the number density distribution. Cumulative moments are expressed here as functions of the agglomeration degree using eq 4, with variable upper limit, and eq 6 Z x f ðy, tÞ dy μ0 ðx, tÞ ¼ xmin ¼ 0

ð2Þ ¼ 3465

1 - Iagg ðtÞ ð1 - expf - ax½1 - Iagg ðtÞgÞ a

ð7Þ

dx.doi.org/10.1021/ie101933r |Ind. Eng. Chem. Res. 2011, 50, 3464–3474


ARTICLE

Figure 2. Initial and final cumulative zeroth, first, and second moments as functions of particle size for 99.9% agglomeration (a = 0.1).

Z μ1 ðx, tÞ ¼ ¼

x

xmin ¼ 0

yf ðy, tÞ dy

1 1 - f1 þ ax½1 - Iagg ðtÞg expf - ax½1 - Iagg ðtÞg a2 a2 ð8Þ Z x μ2 ðx, tÞ ¼ y2 f ðy, tÞ dy xmin ¼ 0

2 - fa2 x2 ½1 - Iagg ðtÞ2 þ 2ax½1 - Iagg ðtÞ þ 2g expf - ax½1 - Iagg ðtÞg ¼ a3 ½1 - Iagg ðtÞ

ð9Þ Note that, unlike the full moments defined previously, the cumulative moments are functions of x. Shown in Figure 2 are the first three cumulative moments (eq 7-9) for a = 0.1 and an agglomeration degree of 0.999. 2.4. Full Moments of the Smoluchowski Equation as Functions of the Agglomeration Degree. An expression for the full zeroth moment as a function of agglomeration degree can be obtained directly from the cumulative zeroth moment of eq 7 by setting x equal to ¥, to give μ0, f ðIagg Þ ¼

1 - Iagg a

ð10Þ

The zeroth moment decreases with increasing agglomeration degree. For an arbitrary initial PSD, it can be written as a function of time as14 μ0, f ðtÞ ¼

1 1 Kt þ μ0, f ðt ¼ 0Þ 2

ð11Þ

In a similar fashion, the first moment can be derived from eq 8 by setting x equal to ¥ to give μ1, f ðIagg Þ ¼ μ1, f ðtÞ ¼

1 ¼ constant a2

ð12Þ

This equation indicates that the first moment does not depend on the agglomeration degree and is preserved during agglomeration.

Figure 3. Normalized moments and maximum number densities of the Smoluchowski equation as functions of the agglomeration degree.

The second moment of the Smoluchowski equation as a function of time was estimated by Dubovski11 as Z ¥ μ2, f ðtÞ ¼ x2 f ðx, tÞ dx 0

¼ μ2, f ðt ¼ 0Þ þ tKμ1, f 2 ðtÞ

ð13Þ

For the statement in eq 1, with xmin= 0 and xmax = ¥, after analytic integration of eqs 2-9, the evolutions of the second moment with respect to time and agglomeration degree are given by 2 Kt ð14Þ μ2, f ðtÞ ¼ 3 þ 4 a a ! 2 1 μ2, f ðIagg Þ ¼ 3 ð15Þ a 1 - Iagg For the purpose of comparison, we define the full normalized moments as μi, f ðIagg Þ , i ¼ 0, 1, 2 ð16Þ μ0i, f ðIagg Þ ¼ μi, f ðIagg ¼ 0Þ The normalized moments are presented as a function of agglomeration degree in Figure 3. In addition, the evolution of the maximum number density is also shown. The maximum number density function is defined as f0 max = max{f/max[f(0)]} of eq 6. As expected, the first moment is constant, with the zeroth moment decreasing and the second moment increasing. As shown in Figure 3, the differences between the higher moments and the maximum of the number density function increase with increasing agglomeration degree. For an agglomeration degree of 99.9%, a 106 range of values needs to be taken into account for an adequate representation of the final number density and its zeroth and first moments. This means that, in this case, the relative numerical error should be less than 10-6. In the case of a 99.9999% agglomeration degree, numerical calculations can be challenging because the required relative error of 10-12 is getting close to the digital 64-bit double precision of 10-16. 2.5. Minimum and Maximum Particle Sizes. It is desirable to know the value of xmin > 0 for the statement in eq 1 that leads to a result that closely matches the analytic solution in eq 2. For the exponential number density case, xmin can be determined 3466



ARTICLE

using the relative error (ε) between the cumulative zeroth moment of eq 7 and the initial full zeroth moment ε ¼

μ0 ðx ¼ xmin , Iagg ¼ 0Þ ¼ 1 - expð- axmin Þ μ0 ðx ¼ ¥, Iagg ¼ 0Þ

nonlinear, but the Safronov equation does not contain a convolution integral. The Safronov agglomeration equation was originally stated as

ð17Þ

t∈ð0, tend ,

x∈ðxmin , xmax Þ, y∈ðxmin , xmax Þ Z x Z xmax ∂f ðx, tÞ ∂ ¼ - K ½f ðx, tÞ yf ðy, tÞ dy - Kf ðx, tÞ f ðy, tÞ dy ∂t ∂x xmin x ∂f ðx, tÞ ∂f ðx, tÞ ¼ ¼ 0 ∂x ∂x

where the minimum size is xmin ¼ - lnð1 - εÞ=a

ð18Þ

For ε = 0.001 and a = 0.1, xmin = 0.01. For higher moments and larger values of the agglomeration degree, the relative error will be less than ε. It can be seen from the number density and its moments, shown in Figures 1 and 2, that a minimum size value of xmin ¼ 0:001=a

ð19Þ

preserves all quantitative information about the number density and its moments. For a general PSD, a number-average particle size can be used for the parameter a, as defined by eq 35. The maximum particle size, xmax, can be determined based on the relative error between the first cumulative and full moments of the PSD. It is possible to calculate the first cumulative moment (eq 8) by increasing the upper integration limit, denoted by x = xmax. For example, an acceptable xmax value can be estimated to guarantee an error in the first moment of less than 0.1% ε ¼

μ1 ðx ¼ xmax , Iagg Þ < 0:001 μ1 ðx ¼ ¥, Iagg Þ

ð20Þ

The Smoluchowski equation (eq 1) has a parabolic property,11 which means that a perturbation of the solution extends to infinity along the particle size coordinate. For a fixed size range, this property of the Smoluchowski equation can lead to a violation of the conservation law for the full first moment and a strong deviation of higher moments during the agglomeration simulation. The other disadvantage of the Smoluchowski equation in numerical solution schemes comes from the discontinuity of the convolution integral, which results in an ill-posed problem so that a special regularization with recursive refinement of the solution is required.18 For this reason, we are seeking an approximation of eq 1 that exhibits mathematical properties better suited to a finite-element-based numerical solution of the agglomeration model.

3. ORIGINAL SAFRONOV AGGLOMERATION EQUATION 3.1. Safronov Agglomeration Model. While the Smoluchowski equation was originally derived to describe the aggregation of colloidal particles dispersed and stirred by Brownian motion under laminar conditions, another fundamental continuous agglomeration equation formulated by Safronov19 was used originally to describe preplanetary dispersed-phase aggregation in a turbulent gas environment. For crystallization applications, the Safronov equation is particularly suitable for agglomeration modeling in systems with broad particle size distributions under a turbulent flow regime. Dubovski11 showed that, in most cases, the Smoluchowski and Safronov aggregation equations give very similar results and, thus, can both be applied to the analysis of agglomeration in dispersed systems. The mathematical properties of the two equations, however, are quite different. The continuous Safronov equation exhibits hyperbolic properties, and a perturbation in the solution propagates with a finite speed.11 Both equations are

xmin

ð21Þ

xmax

f ðx, tÞjt ¼ 0 ¼ f0 ðxÞ ¼ expð - axÞ

Similarly to the Smoluchowski equation (eq 1), the agglomeration model in eq 21 also describes the evolution of the number density as a function of time during agglomeration with a sizeindependent agglomeration kernel, K. 3.2. Asymptotic Solution of the Safronov Equation. Although the Safronov aggregation equation exhibits some promising mathematical properties, no analytic solutions have yet been developed for the statement of eq 21. A comparison of numerical solutions of the discrete forms of the Safronov and Smoluchowski equations with K = 1, tend = 10, and a = 0.7 (final Iagg < 0.9) was carried out by Dubovski,2 who showed analytically that the evolutions of the full zeroth and first moments predicted by the Safronov and Smoluchowski equations are exactly the same. Whereas the general solution of eq 21 describes the evolution of the initial PSD, the asymptotic solution at infinite time, for a constant agglomeration kernel, has self-similar properties similar to those of the Smoluchowski equation.20 The self-similarity means that the finite-size initial PSD loses its detailed information with time and the solution asymptotically converges to a simple shape, as shown in Figure 4. As outlined by Bekker and Livk,21 we assume that the Safronov equation (eq 21) has a solution in the form of a traveling PSD wavefront with a rectangular form for the asymptotic time satisfying the condition t.

2 Kμ0, f ðt ¼ 0Þ

ð22Þ

Using eqs 10 and 12, we can solve the system 1 μ0, f ðtÞ ¼ fa ðtÞ xf ðtÞ ¼ 1 Kt þ μ0, f ðt ¼ 0Þ 2

μ1, f ðt ¼ 0Þ ¼ fa ðtÞ

xf 2 ðtÞ ¼ constant 2

ð23Þ

to obtain the height, fa(t), and width, xf(t), of a rectangular PSD as 2 32 fa ðtÞ ¼

6 6 6 2μ1, f ð0Þ4 1

7 7 7 1 Kt 5 þ μ0, f ð0Þ 2 1

" 1 Kt þ xf ðtÞ ¼ 2μ1, f ð0Þ μ0, f ð0Þ 2

# ð24Þ

For a general shape of the initial PSD, the asymptotic solution of the Safronov equation describing the number density function 3467



ARTICLE

Figure 4. Asymptotic and numerical agglomeration model results for a 99.9% agglomeration degree: (a) number density and (b) cumulative zeroth and first moments.

and its zeroth and first cumulative moments is 8 < fa ðtÞ, xexf ðtÞ f ðx, tÞ ¼ : 0, x > xf ðtÞ ( xfa ðtÞ, xexf ðtÞ μ0 ðx, tÞ ¼ xf ðtÞ fa ðtÞ, x > xf ðtÞ ( x2 fa ðtÞ=2, xexf ðtÞ μ1 ðx, tÞ ¼ xf 2 ðtÞ fa ðtÞ=2, x > xf ðtÞ

the solution of eq 24 is fSF ðx, tÞ ¼ f ðx, tÞ 8 0 12 > > > > < 1B 1 C , xex ðtÞ ¼ 2 1 þ Kt @ A f 2 1 þ Kt a 2a ¼ > > 2a > > : 0, x > xf ðtÞ ð25Þ

This solution is presented in Figure 4. The velocity of the front propagation, vx, of the Safronov equation solution is vx ¼

dxf ðtÞ ¼ Kμ1, f ðt ¼ 0Þ ¼ constant dt

ð26Þ

The full second moment of the Safronov equation solution can be obtained from eq 24 as μ2, f ðtÞ ¼ fa ðtÞ

Note that, in the above asymptotic solution, the amplitude is only one-half of that from the Smoluchowski equation solution at x = 0, as concluded by comparing eqs 2 and 29 and as illustrated in Figure 4a. The velocity of the front propagation, vx = dxf/dt = Kμ1,f(t = 0), in eqs 24 and 29 is constant if volume is taken as the internal coordinate for the PSD description. If the particle diameter, L, is chosen as the internal coordinate, then the front propagation velocity, vL, decreases with increasing particles size as

xf 3 ðtÞ 3

vL ¼

2 ¼ μ2, f ðt ¼ 0Þ þ t Kμ1, f 2 ðt ¼ 0Þ 3

ð27Þ

Compared with eq 13, the Safronov second moment evolves with two-thirds of the velocity of the Smoluchowski second moment. For the exponential initial PSD, the solution for the evolution of the full second moment from the Safronov equation (eq 21) is μ2, f ðtÞ ¼

2 2 Kt þ a3 3 a4

ð28Þ

A comparison of eq 27 with eq 13 provides evidence about an asymptotic behavior of the full second moment. For a general PSD shape, the Safronov equation solution's full second moment increases with time and asymptotically approaches a value that is two-thirds of the value of the Smoluchowski equation solution’s full second moment, as presented in eq 13. 3.3. Characteristics of the Safronov Asymptotic Solution. For the statement in eq 21 and an exponential initial PSD,

ð29Þ

dLðtÞ vx constant ¼ ¼ dt L2 3kV L2

ð30Þ

where kV is the volume shape factor, L = (x/kV)1/3. For example, if the maximum particle diameter during agglomeration increases 10-fold, then the PSD front propagation velocity decreases 102fold. The agglomeration front propagation behavior of the Safronov equation solution is somewhat similar to the Burgers equation22 front propagation. The steepness of the initially exponential PSD increases with time and asymptotically approaches a shock-wave quadratic form, as shown in Figure 7a. In the next section, the original Safronov equation statement (eq 21) is reformulated into a PDE form that is more amenable to numerical solution methods.

4. FORMULATION OF A PDE AGGLOMERATION MODEL FROM THE SAFRONOV EQUATION For both the Smoluchowski equation (eq 1) and the original Safronov equation (eq 21), a sparse band-diagonal Jacobian matrix (matrix of partial derivatives) cannot be constructed within a fully implicit numerical scheme. The problem arises 3468



ARTICLE

from the nonlinear integrals on the right-hand sides of the two equations. To overcome this problem with the fully implicit approximation of the nonlinear Safronov integral, eq 21 was reformulated as a system of PDEs. 4.1. One-Dimensional PDE Agglomeration Model. The PDE reformulation of eq 21 is straightforward through the use of a classical mathematical integral description. Two integrals in eq 21 are transformed through the use of new dependent variables, namely, the cumulative first and zeroth moments. The system of three dependent variables is closed by adding two moment equations in the form of partial derivatives of the cumulative zeroth and first moments t∈ð0, tend ,

x∈ðxmin , xmax Þ,

Table 1. Parameter Values Used in the Validation of the Developed PDE Agglomeration Model

∂f ðx, tÞ ∂ ¼ - K ½f ðx, tÞ μ1 ðx, tÞ þ Kf ðx, tÞ½μ0 ðx, tÞ ∂t ∂x - μ0 ðxmax , tÞ ∂μ0 ðx, tÞ ¼ f ðx, tÞ ∂x ∂μ1 ðx, tÞ ¼ xf ðx, tÞ ∂x

,

1.0 10-3

xmax tend

dimensionless maximum particle volume dimensionless process time

3.0 103 1

K

dimensionless agglomeration kernel

2000

a

dimensionless PSD constant

1

D

dimensionless artificial diffusion

10-9

xmax

xmax

f ðx, tÞjt ¼ 0 ¼ f0 ðyÞ ¼ expð- ayÞ, Z

xmax

¼

xf0 ðxÞ dx

x

μ0 ðx, tÞjt ¼ 0

f0 ðyÞ dy,

xmin

Z μ1 ðx, tÞjt ¼ 0 ¼

xmin

x

yf0 ðyÞ dy

ð33Þ

xmin

ð31Þ

xmin

The initial boundary value problem, as defined by eq 31, can be viewed as a convection-reaction system. To increase the stability of numerical solution for the case with K > 0, the left zero Neumann boundary condition for the first equation of the system in eq 31 should be replaced by the Dirichlet boundary condition. In eq 21, the lower integration limit is set to x = xmin, and the left boundary value for the number density is defined by the equation Z xmax dfxmin ðtÞ ¼ - Kfxmin ðtÞ f ðx, tÞ dx ð32Þ dt xmin Also, elliptic regularization is carried out to smooth the discontinuity between the values of the initial and boundary conditions in the convection equations by adding an artificial diffusion term.23 Based on the original Safronov equation, the regularized PDE agglomeration model is stated as x∈ðxmin , xmax Þ,

dimensionless minimum particle volume

xmax

f ðx, tÞjt ¼ 0 ¼ f0 ðyÞ ¼ expð - ayÞ, Z x f0 ðyÞ dy, μ0 ðx, tÞjt ¼ 0 ¼

t∈ð0, tend ,

xmin

¼ fxmin ðtÞ, μ0 ðx, tÞjxmin ¼ μ1 ðx, tÞjxmin ¼ 0 ∂f ðx, tÞ ∂μ0 ðx, tÞ ∂μ1 ðx, tÞ ¼ ¼ ¼ 0 ∂x ∂x ∂x

xmin

∂f ðx, tÞ ∂μ0 ðx, tÞ ∂μ1 ðx, tÞ ¼ ¼ ¼ 0 ∂x ∂x ∂x

μ1 ðx, tÞjt ¼ 0 ¼

value

dfxmin ðtÞ ¼ - Kfxmin ðtÞ μ0 ðxmax , tÞf ðx, tÞjxmin dt

∂f ðx, tÞ ¼ μ0 ðx, tÞ ¼ μ1 ðx, tÞ ¼ 0 ∂x

x

meaning

∂μ0 ðx, tÞ ∂ ∂μ ðx, tÞ ¼ D 0 þ f ðx, tÞ ∂x ∂x ∂x ∂μ1 ðx, tÞ ∂ ∂μ1 ðx, tÞ ¼ D þ xf ðx, tÞ ∂x ∂x ∂x

y∈ðxmin , xmax Þ

Z

parameter

y∈ðxmin , xmax Þ

∂f ðx, tÞ ∂ þ K ½f ðx, tÞ μ1 ðx, tÞ ∂t ∂x ∂ ∂f ðx, tÞ D ¼ þ Kf ðx, tÞ½μ0 ðx, tÞ - μ0 ðxmax , tÞ ∂x ∂x

To solve this PDE agglomeration model numerically, all terms of the system of equations in eq 33 were approximated by a fully implicit discrete numerical scheme with respect to time. For the implicit scheme, the resulting linearized system of algebraic equations has a sparse band-diagonal Jacobian coupling matrix, whereas the Jacobian matrices in the implicit cases of the original Safronov and Smoluchowski equations are full. The sparse banddiagonal Jacobian matrix of eq 33 is diagonally dominant, which allows stable numerical calculations. 4.2. Numerical Solution of the Agglomeration Model and Its Validation. A numerical solution of the regularized onedimensional PDE agglomeration operator, as defined by eq 33, was implemented within a commercially available PDE solver10 that used a general Newton’s-method-based implicit Galerkin finiteelement algorithm with adaptive nonuniform mesh and time-step strategies. A flowchart of the numerical algorithm can be found in Bekker and Livk.23 The numerical algorithm refines a finite-element mesh with a second-order polynomial basis until the relative error in any dependent variable is less than a specified tolerance, which was set to 10-7 for every cell of the mesh. The Gear method with a second-order fully implicit two-step backward-difference formula is used for variable time steps. The initial step is formed using one and two overlapping first-order steps, repeated with smaller time steps until the second-order error is less than the prescribed tolerance.10 The set of model parameters used in these calculations is presented in Table 1. 3469



ARTICLE

In accordance with eq 3, the final agglomeration degree for the parameter set given in Table 1 is Iagg ¼ 1 -

1 ¼ 0:999 Ktend 1þ 2a

ð34Þ

Numerical results for the agglomeration model of eq 33, obtained using the given set of parameters, are compared to the asymptotic Safronov solution given by eqs 25 and 29 in Figure 4a. In Figure 4b, the evolutions of the zeroth and first cumulative moments obtained numerically from the PDE agglomeration model are compared to the asymptotic solution of eq 25 for the Safronov equation. The initial exponential PSD and its moments, shown in Figure 7, were used as initial conditions. The very good agreement between the asymptotic and numerical solutions indicates that the implemented numerical solution of the PDE agglomeration model provides a close approximation of the original Safronov agglomeration equation. Shown in Figure 5 are comparisons between the asymptotic and numerical PDE agglomeration model solutions for the maximum value of the number density function and the full zeroth, first, and second moments. The corresponding analytic solutions for the zeroth, first, and second moments were obtained by eqs 11, 12, and 28, respectively. As demonstrated

Figure 5. Comparison between numerical solutions of the PDE agglomeration model (symbols) and analytic Safronov solutions (lines) of the maximum number density and zeroth, first, and second full moments.

in Figure 5, the asymptotic and numerical solutions are in very good agreement. Good agreement between the numerical solution of the PDE agglomeration model and the analytic solution of the Safronov agglomeration equation for agglomeration degrees of up to 99.9% has been demonstrated. The difference between the numerical number density and the asymptotic Safronov equation solution is ultimately controlled by the artificial diffusion and relative tolerance applied in the numerical algorithm, as discussed by Bekker and Livk.23 The maximum relative error between the numerical and analytic zeroth, first, and second moments does not exceed 0.6% for agglomeration degrees of up to 99.9%. Note that, as shown by Dubovski,2 the analytic expressions for the full zeroth and first moments are exactly the same for the Safronov and Smoluchowski equation solutions. However, the analytic expressions for the full second moment differ. 4.3. Comparison between the New PDE and Smoluchowski Agglomeration Models. Dubovski2 showed analytically that, for the Safronov equation, the first two moments, namely, the zeroth and first, are exactly the same as those obtained analytically for the Smoluchowski equation. Here, we explore relative errors by comparing the first two full moments, μ0,f and μ1,f, obtained numerically from the new PDE agglomeration model to the Smoluchowski analytic solution. The relative errors of the first two full moments for an agglomeration degree of 99.9% are presented in Figure 6, together with the evolution of the number-based mean particle size estimated

Figure 6. Evolution of the number-average particle size and relative errors of the zeroth and first full moments for the PDE agglomeration model.

Figure 7. Comparison between the numerical PDE agglomeration model solution (dashed line) and analytical Smoluchowski solution (solid line) for different agglomeration degrees: (a) number density function and (b) zeroth and first cumulative moments. 3470



ARTICLE

numerically from the PDE model and analytically from the Smoluchowski equation. The number mean particle size is defined as the ratio between the first and zeroth full moments as R¥ μ1, f ðtÞ xf ðx, tÞ dx ðtÞ ¼ R0 ¥ x0 ðtÞ ¼ ð35Þ μ0, f 0 f ðx, tÞ dx According to eqs 11 and 12, the theoretical value of x0 presented in Figure 6 for both the Smoluchowski and Safronov equations increases linearly with time as x0 ðtÞ ¼

μ1, f ð0Þ Kt 1 K þ μ1, f ð0Þ ¼ þ 2 t 2 a 2a μ0, f ð0Þ

ð36Þ

As the relative errors between the analytic Smoluchowski and PDE agglomeration solutions for the full zero and first moments do not exceed 0.6%, the match between the analytic and numerical x0 values is also very good. Shown in Figure 7a is a comparison between the analytic Smoluchowski and numerical PDE agglomeration model number density functions as predicted for agglomeration degrees of up to 99.9%. Independently of the agglomeration degree, the two number density functions match well in the small-particle-size region. With increasing particle size, number density value from the PDE agglomeration model remains lower than that from the Smoluchowski equation, until the two intercept at a larger particle size. For very large particle sizes, the Safronov solution drops almost vertically with increasing particle size, which clearly demonstrates a right-hand-side constraint of the Safronov agglomeration equation. Similarly, the first two cumulative moments from the PDE agglomeration model start to deviate from the Smoluchowski solution with increasing particle size, as shown in Figure 7b. The number density result presented here differs slightly from that of Dubovski,2 who concluded that, in comparison to the Smoluchowski solution, the Safronov solution gives lower values for small particle sizes, but nearly identical values for larger particle sizes. In his work, the discrete Safronov equation model was solved only for small agglomeration degrees using a crude mesh resolution. The relative difference between the Smoluchowski analytic solution and the Safronov number densities in any size interval [x1, x2] for particle sizes less than xf of eq 24 can be shown to be less than 70%, that is R R x2 x2 x1 fSM ðx, tÞ dx - x1 fSF ðx, tÞ dx hR i < 0:7 ð37Þ ε ¼ R x2 x2 max x1 fSM ðx, tÞ dx, x1 fSF ðx, tÞ dx These results confirm that the original Safronov agglomeration equation and the PDE approximation of the Safronov agglomeration model can both accurately resolve the agglomeration of small particles as a size-independent process. At the same time, modeling results for larger agglomerate sizes can be considered only an approximation of those predicted by the Smoluchowski equation. The difference observed between the Smoluchowski and Safronov solutions comes as a result of a size-dependent effective probability of binary particle interactions in the constant-kernel Safronov equation, which is different from the sizeindependent probability described by the Smoluchowski equation. In the Safronov equation, the large-size end of the PSD is automatically constrained, which is critical for preserving the

consistency of higher moments of the PSD, without the need to introduce any nonphysical model parameters. 4.4. Gibbsite Agglomeration Case Study. To illustrate the application of the new Safronov-based PDE agglomeration model, we apply the developed agglomeration model to describe the gibbsite agglomeration process taking place in an isothermal constant-supersaturation crystallizer with a solids retention system in place.24 The gibbsite agglomeration process has been studied extensively because of its importance in the manufacturing of alumina using the Bayer process.25,26 In the case investigated here, gibbsite seeds are charged into the crystallizer, where they agglomerate with time at constant supersaturation. A Safronov equation gibbsite agglomeration model can be stated based on eq 21 in dimensional form as T∈ð0, Tend ,

V ∈ðVmin , Vmax Þ, U∈ðVmin , Vmax Þ " # Z V ∂nðV , TÞ ∂ ¼ - β0 nðV , TÞ UnðU, TÞ dU ∂T ∂V Vmin Z Vmax - β0 nðV , TÞ nðU, TÞ dU ∂nðV , TÞ ∂V

V

Vmin

∂nðV , TÞ ¼ ∂V

nðV , TÞjt ¼ 0 ¼ n0 ðV Þ ¼

¼ 0 Vmax

N0 V exp V0 V0

ð38Þ

where V and U are the gibbsite particle volumes and T is the time. The dependent variable, n(V,T), is the number density function. The values of the model parameters used in the simulation are presented in Table 2. Equation 38 is related to eq 21 by the following normalization: f = n/n0, x = V/V0, y = U/V0, t = T/T0. If n0 = N0/V0, V0 = V0, and T0 = Tend, then the dimensional parameters of eq 37, presented in Table 2, are related to the dimensionless parameters of eq 21, presented in Table 1, as xmin = Vmin/V0, xmax = Vmax/V0, tend = Tend/T0, K = β0N0Tend, and a = V/V0. All dimensionless parameters are the same as those given in Table 1. The gibbsite mass in the system is estimated as ms= FsN0V0. The characteristic aggregation time for a 50% agglomeration degree is τa = 2/(N0β0). The numerical solution of the gibbsite agglomeration model (eq 38) was obtained using its dimensionless form of eq 21 and the corresponding regularized PDE agglomeration model of eq 33 with the dimensionless artificial diffusion coefficient of Table 1. The model-predicted number density functions at different agglomeration degrees are presented in Figure 8, which shows the number densities for 50% and 90% agglomeration degrees as a function of particle diameter. The initial number density function is also presented for comparison in each case. Comparison between the Safronov and Smoluchowski solutions shows that, for a moderate agglomeration degree, such as 50%, they are relatively closely matched. For higher agglomeration degrees, such as 90% and more, the difference between the two becomes larger, but the main trend in the evolution of the number density function is still closely preserved. For any agglomeration degree, the Safronov solution closely matches the Smoluchowski solution in the small-particle-size region, close to the particle size of 1.4 μm. Although not presented here, the two solutions match closely in terms of the first two full moments, as illustrated previously in Figures 5 and 6. 3471



ARTICLE

Table 2. Values of Gibbsite Agglomeration Parameters Used in the Simulation parameter

meaning

dmin Vmin = kVdmin

3

value

units

minimum particle diameter

1.4 10-6

m

minimum particle volume

1.44 10-18

m3

-6

dmax

maximum particle diameter

200 10

m

Vmax = kVdmax3

maximum particle volume

4.19 10-12

m3 s

tend

time for 99.9% agglomeration degree

50 3600

t90

time for 90% agglomeration degree

4 3600

s

β0

agglomeration kernel

2.78 10-15

m3 s-1

d0 V0 = kVd03

initial number-average particle diameter initial number-average particle volume

14 10-6 1.44 10-15

m m3

N0

initial total number of particles

4.0 1012

m-3

kV

volume shape factor

π/6

Fs

solids density

2350

kg m-3

ms

total solids mass

13.5

kg m-3

τa

time for 50% agglomeration degree

180

s

Figure 8. Numerical Safronov-based PDE agglomeration model and analytical Smoluchowski number density functions of gibbsite particles for (a) 50% and (b) 90% agglomeration degrees.

Shown in Figure 9 is the evolution of the volume median (dV50) and number mean (d0) gibbsite particle diameters predicted by the Safronov agglomeration model compared to those based on the Smoluchowski model results. The volume median diameter is defined as dV50 = LV50 = (V50/kV)1/3, where V50 is the volume median particle size, which gives the cumulative first moment equal to 50% of the first full moment R x50 ðtÞ μ1 ðx50 , tÞ xf ðx, tÞ dx ðtÞ ¼ 0R ¥ 0:5 ¼ μ1, f 0 xf ðx, tÞ dx R V50 ðtÞ VnðV , tÞ dV ð39Þ ¼ 0R ¥ 0 VnðV , tÞ dV The number mean gibbsite particle diameter, d0, is the diameter of a particle of the average volume, which can be estimated based on eqs 35 and 36 as x0 ðtÞV0 1=3 d0 ðt ¼ 0Þ ¼ ð40Þ d0 ðtÞ ¼ kV ½1 - Iagg ðtÞ1=3 A comparison between the numerical Safronov and analytic Smoluchowski predictions of the volume median and number

mean gibbsite particle diameters as functions of the agglomeration degree is shown in Figure 9. As shown in Figure 9a, the number-based mean particle diameter predicted by the Safronov-based PDE agglomeration model is the same as that predicted by the analytic Smoluchowski solution. This result is expected based on previous findings that the full zeroth and first moments of the Safronov and Smoluchowski agglomeration models are exactly the same. The two agglomeration models, however, predict slightly different volume median diameters, as shown in Figure 9b. The two volume-based median diameters start deviating from each other with increasing agglomeration degree, reaching a difference of around 4% for a 99.9% agglomeration degree. This deviation is a result of the difference between the Safronov- and Smoluchowski-predicted number density functions, as presented in Figure 8.

5. CONCLUSIONS A one-dimensional dynamic PDE agglomeration model, represented by a set of convection-reaction-diffusion partial differential equations (PDEs), was derived based on the continuous Safronov aggregation equation. A numerical solution of the developed PDE agglomeration model was implemented 3472



ARTICLE

Figure 9. Evolutions of the volume median and number mean gibbsite particle diameters as predicted by numerical Safronov-based PDE agglomeration model and analytical Smoluchowski solutions: (a) number-based mean diameter and (b) volume-based median diameter.

using a regularization of the new agglomeration model and a Gear-type Newton’s-method-based Galerkin finite-element algorithm that converged in terms of time and mesh resolution. Importantly, an implicit discretization of the new agglomeration model leads to a sparse band-diagonal Jacobian matrix, which represents a significant advantage over full Jacobian matrices resulting from the original Safronov and Smoluchowski agglomeration models. The new Safronov-based PDE agglomeration model is developed so as to fit into the adaptive-mesh-size FEM-based population balance solution framework, which does not require additional quadrature assumptions for the evaluation of the aggregation integral. The numerical solution of the Safronov-based PDE agglomeration model was compared to the analytic solution of the Smoluchowski equation and the asymptotic Safronov equation solution. For agglomeration degrees of up to 99.9%, the numerical results for the first two full moments are in excellent agreement with the asymptotic solution of the Safronov equation and the analytic solution of the Smoluchowski equation. Numerical solutions of the number density function show that, for small particle sizes, the two agglomeration models, Safronovbased and Smoluchowski, predicted very similar values. For larger particle sizes, however, the difference between the two models increased with increasing agglomeration degree. A comparison between the numerical Safronov-based PDE agglomeration solution and Smoluchowski analytic cumulative moments shows a trend similar to that observed for the number density function. The new Safronov-based PDE agglomeration model was applied for modeling the agglomeration of gibbsite crystals in an isothermal constant-supersaturation crystallizer. The results for this case showed that the number density predicted by the new Safronovbased agglomeration model approximately matched that predicted by the Smoluchowski model. The two models predicted consistent number mean diameters of gibbsite aggregates, but deviated by up to 4% in the prediction of the median volume diameter. The newly developed PDE agglomeration model will be incorporated into the crystallization model built within a commercial PDE solver to predict the behavior of a realworld gibbsite crystallization process. The presented approach, however, is not restricted to crystallization and can be used for modeling other systems that involve particle aggregation.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ REFERENCES (1) Smoluchowski, M. V. Drei vortr€age€uber diffusion, Brownische bewegung und koagulation von kolloidteilchen. Phys. Z. 1916, 17 (23), 557–585; Phys. Z. 1916, 17 (23), 585–599 (in German). (2) Dubovski, P. B. New discrete model of coagulation kinetics and properties of continuous analog. Math. Model. 2000, 12 (9), 3–15 (in Russian). (3) Agoshkov, V. I.; Dubovski, P. B.; Shutyaev, V. P. Methods for Solving Mathematical Physics Problems; Cambridge International Science Publishing: Cambridge, U.K., 2006. (4) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Processes in Engineering; Academic Press: New York, 2000. (5) Woo, X. Y.; Tan, R. B. H.; Chow, P. S.; Braatz, R. D. Simulation of Mixing Effects in Antisolvent Crystallization Using a Coupled CFD-PDF-PBE Approach. Cryst. Growth Des. 2006, 6, 1291–1303. (6) Alexopoulos, A. H.; Kiparissides, C. Part II: Dynamic evolution of the particle size distribution in particulate processes undergoing simultaneous particle nucleation, growth and aggregation. Chem. Eng. Sci. 2005, 60, 4157–4169. (7) Hounslow, M. J.; Ryall, R. L.; Marshall, V. R. A discretized population balance for nucleation, growth and aggregation. AIChE J. 1988, 34 (11), 1821–1832. (8) Gelbard, F.; Seinfeld, J. H. Numerical solution of the dynamic equation for particulate systems. J. Comput. Phys. 1978, 28, 357–375. (9) Jacobson, M. Z.; Turco, R. Modeling coagulation among particles of different composition and size. Atmos. Environ. 1994, 28, 1327–1338. (10) FlexPDE 5. Help, User Guide and Reference Manual; PDE Solutions Inc.: Spokane Valley, WA, 2005. (11) Dubovski, P. B. A 'triangle' of interconnected coagulation models. J. Phys. A 1999, 32, 781–793. (12) Bekker, A. V.; Li, T.; Livk, I. Comparison of FEM and DPB Numerical Methodologies for Dynamic Modeling of Isothermal Batch Gibbsite Crystallization. Ind. Eng. Chem. Res. 2011, DOI: 10.1021/ ie101942s. (13) M€uller, H. Zur allgemeinen Theorie der raschen Koagulation: Die Koagulation von St€abchen- und Bl€attchenkolloiden; die Theorie beliebig polydisperser Systeme und der Str€omungskoagulation. Kolloidchem. Beihefte 1928, 27, 223–250 (in German). (14) Schumann, T. E. W. Theoretical aspects of the size distribution of fog particles. Q. J. R. Meteorol. Soc. 1940, 66, 195–207. (15) Wright, H.; Ramkrishna, D. Solutions of inverse problems in population balances-I. Aggregation kinetics. Comput. Chem. Eng. 1992, 16, 1019–1038. 3473



ARTICLE

(16) Rigopoulos, S.; Jones, A. G. Finite-element scheme for solution of the dynamic population balance equation. AIChE J. 2003, 49, 1127– 1139. (17) Hounslow, M. J. A Discretized Population Balance for Nucleation, Growth and Aggregation. Ph.D. Thesis, University of Adelaide: Adelaide, Australia, 1990. (18) Luk’yanenko, V. A. Some algorithms for approximate solution of convolution integral equations. J. Math. Sci. 1996, 82, 3384–3390. (19) Safronov, V. S. Evolution of the Pre-Planetary Cloud and the Formation of the Earth and Planets; Nauka: Moscow, 1969 (in Russian) . (20) Laurencot, P. Convergence to self-similar solutions for a coagulation equation. Z. Angew. Math. Phys. 2005, 56, 398–411. (21) Bekker, A. V.; Livk, I. A PDE agglomeration model for Bayer precipitation. In CD-ROM Proceedings of the Chemeca 2009 Conference; ICMS Pty Ltd.: Victoria, Australia, 2009. (22) Burgers, J. M. Mathematical examples illustrating relations occurring in the theory of turbulent fluid motion. Trans. R. Neth. Acad. Sci. Amsterdam 1939, 17, 1–53. (23) Bekker, A. V.; Livk, I. An FEM Modeling of Gibbsite Crystallization with Constant and Nonlinear Kinetics. Ind. Eng. Chem. Res. 2011, DOI: 10.1021/ie1019326. (24) Ilievski, D. Development of a constant supersaturation, semibatch crystalliser and its application to investigating crystal agglomeration. J. Cryst. Growth 2001, 233, 846–862. (25) Ilievski, D.; Livk, I. An Agglomeration Efficiency Model for Gibbsite Crystallization in a Turbulently Stirred Vessel. Chem. Eng. Sci. 2006, 61, 2010–2022. (26) Livk, I.; Ilievski, D. A macroscopic agglomeration kernel model for gibbsite precipitation in turbulent and laminar flows. Chem. Eng. Sci. 2007, 62, 3787–3797.

3474


Agglomeration Process Modeling Based on a PDE Approximation of

Recommend Documents