Influence of the Estimation Procedure on the Accuracy and Precision

Ind. Eng. Chem. Res. 2001, 40, 5005-5013

5005

Influence of the Estimation Procedure on the Accuracy and Precision of Aluminum Trihydroxide Crystallization Kinetics from Dynamic Data Tian S. Li, Iztok Livk, and Dean Ilievski* AJ Parker CRC for Hydrometallurgy, CSIRO Minerals, P.O. Box 90, Bentley, WA 6982, Australia

Differential and integral parameter estimation techniques were used in this work to estimate aluminum trihydroxide crystallization kinetics (i.e., growth, agglomeration, and source term rates) in batch and semibatch, constant-composition crystallization systems. Of the four differential techniques tested, the methods of Bramley et al. (J. Colloid Interface Sci. 1996, 183, 155-165) and Livk et al. (AIChE J. 1999, 45 (7), 1593-1596) were the most accurate. Comparisons of experimental and predicted crystal size distributions showed that the kinetics estimates obtained by the integral method lead to better agreement with the experimental data than those obtained by the differential methods. The integral and differential techniques both give similar estimates of the crystal growth rate but deviate more in the case of agglomeration and source term rates. Overall, the uncertainties of the kinetics estimates obtained by the integral and differential methods were of comparable magnitudes. 1. Introduction The crystal size distribution (CSD) is a critical quality parameter of the crystalline aluminum trihydroxide produced in the commercial Bayer process. The mechanisms that primarily control the product CSD are agglomeration, growth, and nucleation. However, understanding of these mechanisms and of how they interact with the process variables is still incomplete. One of the reasons for this is the difficulty encountered in obtaining reliable estimates for the kinetics of these mechanisms. They occur simultaneously, making it difficult to identify the individual contributions to the overall CSD. This is a common problem in the study of industrial crystallization systems, which is illustrated by the, all too common, occurrence of different values for crystallization kinetic parameters being reported by different authors studying the same system.1-3 Taking the aluminum trihydroxide system, for example, White and Bateman4 and Veesler and Boistelle5 tabulate the activation energies for the growth rate reported by different laboratories and show that they vary from 53 to 120 kJ mol-1. In a further example of this problem, Li et al.6 found that predictions of the aluminum trihydroxide CSD obtained using available literature correlations for the crystallization kinetics in this system7-9 were in poor agreement with their experimental CSD data. Poor reproducibility of crystallization kinetics estimates may be attributed to the following: (1) small random and systematic variations in the local conditions10 that affect the measurements from which the kinetics are estimated; (2) differences in (or poor) experimental design that affect the form and information content of the experimental data; (3) large uncertainties associated with the parameter estimation procedures. The first issue can be considered a property of the laboratory, its crystallizers, operating procedures, sampling methods, and analysis methods. Ideally, a * To whom correspondence should be addressed. Tel.: 61-893348050. Fax: 61-893348001. E-mail: dean.ilievski@ minerals.csiro.au.

laboratory would have conducted sufficient repeats to quantify this uncertainty. Melikhov et al.10 proposed models to describe the variability distribution of certain crystal properties due to fluctuations in certain process conditions. The importance of crystallization experimental design to kinetic parameter estimation and parameter confidence limits is treated in detail by Rawlings et al.3 The effect of different crystallizer configurations on the accuracy and precision of the estimated aluminum trihydroxide kinetics was investigated by Li et al.11 Their experiments showed that, even at the same supersaturation and operating conditions, there were differences in the kinetics estimates. They attributed the differences to uncertainties in the parameter estimates, arising from interplay between the different functional forms of the output CSD data from the different configurations with the uncertainty in the measurement of CSD and uncertainties associated with the parameter estimation procedure. Using Monte Carlo methods combined with dynamic crystallizer simulations, Li et al.11 investigated the influence of experimental design on the accuracy and precision of the aluminum trihydroxide crystallization kinetics estimates. They report that the crystallizer configuration (e.g., batch, semibatch, etc.) can be important under certain conditions (e.g., at lower supersaturation, where rates are low), but operating conditions in the crystallizer have a greater impact on the parameter estimate uncertainty. Rawlings et al.3 and Tavare12 report that different parameter estimation techniques used on the same experimental data can result in different estimates for the kinetics. This observation is the motivation for the current paper. The main objective is to evaluate the accuracy and precision of a variety of crystallization kinetics estimation techniques. Aluminum trihydroxide crystallization, a complex and commercially important process, was selected because experimental studies (Ilievski and White,7 Misra and White,8 and Halfon and Kaliaguine13,14) have shown that the crystal growth and agglomeration kinetics are

10.1021/ie0102145 CCC: $20.00 Published 2001 by the American Chemical Society Published on Web 10/04/2001

5006

Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001

best described by size-independent models. This eliminates the problem of having to incorporate the size dependency into the parameter estimation procedures, but it does limit the conclusions to size-independent kinetic models of growth and agglomeration. However, size-independent crystal growth, i.e., the McCabe ∆L law (Mullin15), and size-independent agglomeration (Mumtaz et al.16) are widely observed for other systems. 2. Estimation of Crystallization Kinetics Different methods for estimating the kinetics of agglomeration, growth, and nucleation simultaneously from experimental crystallization data are reviewed in the work of Tavare.12 More recently, Bramley et al.1 and Livk et al.2 have presented techniques which take advantage of a discretized population balance (DPB) approach to determine crystallization kinetics. Methods for the estimation of crystallization kinetics from dynamic data are generally based on the population balance (PB) theory developed by Randolph and Larson17 and can be divided into two main categories. The first category, the so-called differential methods, uses the differential form of the PB model (e.g., eq 1) together with time derivatives of the crystallization data. The second category, the so-called integral methods, uses the integrated form of the PB model with the experimental crystallization data. The underlying principle of the kinetics estimation procedures is to identify the values of the crystallization parameters, i.e., G(t), β(t), and Bu(t) in eq 1, that make the PB model of the crystallizer (e.g., eqs 1 or 5 or the moments formulation) consistent with the dynamic experimental data. 2.1. PB Model. The PB-based model of a crystallization system usually takes the form of an integropartial differential equation; e.g., for a nonstationary, wellmixed, constant-volume system, the model is

∂n ∂n +G ) Buδ(L - L1) + B - D + ∂t ∂L 1 ( Qinnin V

∑

∑Qoutnout)

(1)

where n is the number density function, t is time, G is the size-independent linear growth rate8 and is defined as dL/dt, L is the characteristic size of the crystals, Q is the volumetric flow into and out of the system, and V is the volume of a system. The source function, Bu, is the number rate of crystals entering the first size interval, L1, of the particle size analyzer and is a measure of the nucleation rate. δ is the Dirac delta function. B and D are the birth and death rates due to agglomeration. For size-independent agglomeration,7,18 they are given by

B(L) )

βL2 2

∫0L

n(t,[L3-λ3]1/3) n(t,λ) (L3 - λ3)2/3

dλ

(2)

and

D(L) ) βn(t,L)

∫0∞n(t,λ) dλ

(3)

where β is the aggregation kernel, which is a measure of the frequency of collisions between crystals of sizes L and λ that are successful in producing an agglomerate.

In some cases it is convenient to transform the PB equation into a set of ordinary differential equations, known as the moment equations,17 where the jth moment of the CSD is defined as

µj )

∫0∞Ljn dL

(4)

The moments can be related to properties of the crystallization product. For example, the first four moments are proportional to the total number, length, surface area, and volume of crystals per unit volume of suspension, respectively. 2.2. DPB Model. The PB models used in this work have no analytical solution. The CSD was simulated using a discretized version of the PB.19,20 For the ith size interval, eq 1 can be written as the following DPB:

dNi 1 ) G[n(L h i-1) - n(L h i)] + Bu,i + Ra,i + dt V ( QinNin - QoutNout) (5)

∑

∑

where the first term on the right-hand side of eq 5 accounts for the crystals growing into and out of the ith size interval. The discretized number density function, n, is defined as

n(L h i) )

Ni Li+1 - Li

(6)

where L h i is the mean particle size of the ith size interval. Bu,i represents the rate of the crystals entering into the field of view of a particle size analyzer and is nonzero only for i ) 1. The term Ra,i accounts for the change of the number of crystals in the ith size interval due to a binary aggregation. A discretized form of Ra,i, proposed by Litster et al.21 and modified by Wynn,22 was used in the simulation for the geometric discretization Li+1/Li ) 21/3q, where q is an integer greater than zero.

Ra,i ) Ra,i(βi,j,Ni,Nj)

(7)

In eq 7, the aggregation kernel βi,j is a measure of the h j collide and frequency at which crystals of size L h i and L form an agglomerate; q ) 6 was selected for the simulations, except where stated otherwise. The DPB model was integrated together with the initial conditions using a fourth-order Runge-Kutta algorithm. 2.3. Kinetics Estimation. Differential Methods. Most applications of the differential methods reported in the literature are based on the principle of matching moments.12,23 In the work of Livk et al.,2 the authors demonstrate that, for a general case, the differential method can also be posed as a nonlinear optimization problem. In this work, four different differential methods are considered and are briefly described below. 2.3.1. Simultaneous Differential Method. Bramley et al.1 The method of Bramley et al.1 uses the DPB developed by Hounslow et al.,20 which they transform to derive equations for the zeroth and third moments (i.e., j ) 0 and 3 in eq 4) and for the number of particles in the first size interval, N1. Application of their method to the general crystallizer configuration described by eq 5 gives the following system of linear

Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 5007

equations, which is solved simultaneously at each time instant.

dµ0 1 ) Bu + Φ0β + ( dt V

(8)

dµ3 1 h 13Bu + ( ) Φ3G + L dt V

(9)

∑Qinµ0,in - ∑Qoutµ0,out)


dN1 ) Φ2G + Bu + Φ1β + dt 1 ( QinN1,in V

∑

∑QoutN1,out)

(10)

The coefficients Φ0, Φ1, Φ2, and Φ3 can be calculated from the experimental data as described in Bramley et al.1 The jth moment of the experimental CSD can be estimated from ms

µj )

L h ijNi ∑ i)1

(11)

where Ni is the number of crystals in the ith size h i is interval, ms is the number of size intervals, and L the size interval midpoint. The time derivative terms in eqs 8-10 were determined by numerical differentiation from the experimental data. 2.3.2. Sequential Differential Method. Livk et al.2 In an attempt to make the kinetic estimation more robust, Livk et al.2 posed the crystallization kinetics estimation problem as a nonlinear optimization, with the objective function formulated in terms of time derivatives of measured crystallization properties. They show that if nucleation does not contribute significantly to the increase of the third moment, the optimization problem can be reduced to three lower-dimensional objective functions, which are optimized sequentially with respect to the kinetics parameters. In the special case of the size-independent growth rate and agglomeration kernel, the authors report that the original nonlinear optimization problem can be reduced to a sequence of three algebraic equations that can be solved explicitly for G, β, and Bu at each time instant. Using the notation from the previous section, these equations can be stated as

G)

β)

[

1 dµ3 1 - ( Φ3 dt V

1 Φ0 - Φ1

{

∑

[

V Bu )

ms dN

∑ i)2 dt

1

]

(12)

}

(13)


i

(Qin

ms

∑ i)2

dµ0 1 - Φ0β - ( dt V

(

Φ3

-

L h3

)

- Φ2 G -

Ni,in) -

∑

ms

Ni,out)] ∑ i)2

(Qout


(14)

When stated as above, this technique becomes very similar to the method of Bramley et al.1 if the size discretization scheme with Li+1 ) 21/3 × Li is used. Livk et al.2 showed that for aluminum trihydroxide crystallization their approach and that of Bramley and coworkers led to almost identical results.

Table 1. Comparisons of Agglomeration Kernel, Source Term Rate, and Growth Rate Estimates Determined from the Different Differential Methods against Selected Kinetics Values selected kinetics simultaneous differential sequential differential moments matching SMM

β (m3 s-1)

Bu (m-3 s-1)

G (m s-1)

3.11 × 10-16 3.11 × 10-16 3.11 × 10-16 3.64 × 10-16 3.14 × 10-16

1.67 × 107 1.67 × 107 1.67 × 107 2.46 × 107 1.69 × 107

1.39 × 10-9 1.39 × 10-9 1.40 × 10-9 1.40 × 10-9 1.40 × 10-9

2.3.3. Moments Matching Method. This method differs from the previous differential methods in that one of the three equations necessary for closure involves the 5th- and 6th-order moments. The use of these higher moments is a consequence of the complexity of the moment transformation on the agglomeration term; for the simplest case of size-independent agglomeration, only expressions whose moments are multiples of 3 can be evaluated analytically.20 The resultant equations are

dµ0 1 1 ) Bu - βµ02 + ( dt 2 V


(15)


(16)

dµ3 1 ) 3Gµ2 + ( dt V

dµ6 1 ) 6Gµ5 + βµ32 + ( dt V


(17)

The values of G, β, and Bu are found by reconciling the moments of the experimental CSD, i.e., eq 11, with eqs 15-17. 2.3.4. SMM Method. This method is a subset of the moments matching method and is only applicable for dynamic crystallization systems with constant kinetics. Under this assumption, Bu and β can be estimated from eq 15 with linear regression analysis, where the moments are estimated from the measured CSD data using eq 11. Similarly, G can be determined from eq 16 by linear regression. Application of this method to experimental crystallization data is reported by Li et al.24 2.3.5. Accuracy of the Different Differential Techniques. The accuracy of each of the above differential techniques was assessed by applying them to simulated CSD data, generated using selected values of G, β, and Bu. The selected kinetics rates were realistic estimates for the aluminum trihydroxide system (see Table 1). A semibatch, constant-composition crystallizer (i.e., constant kinetics) system with clear liquor overflow (i.e., full solids retention and constant volume) was selected for the simulation so that all of the differential methods, including the simplified moments matching (SMM), could be compared. The CSDs were generated using the DPB formulation of this model system with resolution q ) 2. The moments were estimated from the simulated CSD data using eq 11. A three-point numerical differentiation scheme was used to estimate the time derivative terms in the differential methods. Comparisons between the kinetics estimates obtained from different differential techniques and the selected kinetics are presented in Table 1. Table 1 shows that the simultaneous and sequential differential methods provide the most accurate estimates of β, Bu, and G estimates, followed by the SMM method. This result is expected, because both differential methods employed a DPB model similar to that used to simulate the CSD.

5008


The relatively poor quality of the estimates obtained with the moments matching method is most likely a consequence of the method’s dependence on higher moments (i.e., µ5 and µ6) and in particular the numerical differentiation of µ6. The higher moments weight toward the tail of the CSD, where the information content is relatively low given the lower particle numbers and broader size intervals. However, a relatively accurate estimate of G is obtained, suggesting that the errors introduced by numerical differentiation of µ3 are not significant. The reasonable estimates obtained by the SMM method also support the conclusion that it is the reliance on higher moments that compromises the kinetics estimates from the moments matching method. 2.4. Kinetics Estimation. Integral Method. The so-called integral approach has been widely used for the estimation of crystallization kinetics.7,25-27 The integral method can be viewed as a nonlinear parameter optimization. The kinetics are estimated by minimizing an objective function, which is a measure of the deviation between the experimental data and model predictions. In this work, a least-squares objective function, Ψ, of the following form was used: ms m

min Ψ(θ) ) θ

ωi[Nk,i exp - Nk,i(θ)]2 ∑ ∑ i)1k)1

(18)

where Nk,i and Nk,iexp represent the model-predicted and experimental number of crystals in the ith size interval taken at the kth sampling time, respectively. The set of kinetic parameters, θ, consists of the vectors of unknown values of growth, G, agglomeration, β, and source term rates, Bu, at certain time instants during the crystallization run

θ(t)T ) [G(t), β(t), Bu(t)]

(19)

The calculation procedure for the integral method is the same as the procedure used in the optimization of dynamic systems.28 The kinetics were optimized at each sampling time, with their values at times between the two sampling points estimated by linear interpolation. In each step of the optimization procedure, first the model equations (eqs 5-7) are integrated using the set of kinetics parameter values (these estimates are either initial estimates or feedback from the previous optimization step) and initial conditions. The objective function is then evaluated and a new set of values of kinetics is calculated. The updated set of kinetics was estimated using the sequential quadratic programming algorithm, available in the Matlab software package.29 The crystallization process model was solved using the DPB with q ) 6. In the formulation of the criterion in eq 18, it is assumed that the errors are homoscedastic (i.e., of similar magnitude) and that unity weighting is reasonable. 3. Accuracy of the Parameter Estimation Methods for Aluminum Trihydroxide Crystallization The accuracy of the parameter estimates from each of the different techniques was checked by comparing the simulated CSDs, obtained using the estimated kinetics, with the experimental CSDs. Data obtained in a laboratory aluminum trihydroxide crystallization system were used.

3.1. Aluminum Trihydroxide Crystallization Data. Aluminum trihydroxide [Al(OH)3] crystallization from sodium aluminate solutions was performed in two different crystallizer configurations: a batch system and a semibatch, constant-composition system with solids retention (CCSR). Detailed descriptions of the Al(OH)3sodium aluminate system, the crystallizer configurations, and the operating procedures are given in the work of Ilievski18 and Li et al.11 A brief description is given below. A 4-L laboratory stainless steel crystallizer, of the same design and fitted with a draft tube, was used in both crystallizer configurations. The crystallization experiments were carried out at 80 °C and agitated at 600 rpm with an A310 impeller. The synthetic caustic aluminate solution was prepared by digesting technicalgrade Al(OH)3 in an aqueous sodium hydroxide solution. A seed charge of 15 g L-1 of Al(OH)3 crystals was added initially. The constant-composition system was operated in a semibatch mode with clear liquor overflow and solids retained in the vessel. The composition in the crystallizer was maintained constant (i.e., a supersaturation, A/A*, of 1.69) by regulating the feed liquor rate using a simple feedback control loop, based on an inline conductivity measurement. The suspension density, liquor composition, and CSD were determined from samples taken periodically. Considerable effort was made to ensure representative particulate sampling. The CSD was determined using a Coulter counter multi-sizer II. The liquor analysis was performed by an automated Watts-Utley-type titration.30 All data were checked for mass balance consistency. Figure 1 shows the evolution of the zeroth moment, crystal numbers in the first size interval, and the third moment during both the batch and CCSR crystallization experiments. 3.2. Estimates of Crystallization Kinetics by the Different Techniques. 3.2.1. Estimation of Batch Kinetics. The kinetics estimates obtained by the differential and integral methods are plotted in Figure 2 as functions of the supersaturation ratio, A/A*, where A is the aluminum species concentration in solution and the asterisk denotes equilibrium. High supersaturation ratio values correspond to the early stages of the batch crystallization process. Figure 2 shows that the simultaneous and sequential differential methods and the integral method give comparable estimates of the growth, G, and source term, Bu, kinetics. The agglomeration kinetics estimates, β, are also similar at higher supersaturations, where agglomeration rates are higher. The moments matching method, however, only gives reasonable estimates for the growth rate and gives negative values for both β and Bu. The unrealistic estimates of β and Bu are attributed to the dependence of this method on the 5th and 6th moments of the experimental CSD, which show a much larger scatter than the 1st, 2nd and 3rd moments. 3.2.2. Constant Composition Kinetics Estimates. Averaged kinetics estimates obtained by the differential and integral methods from the constant-composition crystallization experiment are shown in Table 2. Again all methods give comparable estimates of the growth kinetics. The moments matching method gives unrealistic estimates for β and Bu. All of the other methods estimate similar values of β. The values of Bu estimated


Figure 1. Transient batch and constant-composition experimental data: (a) zeroth moment; (b) crystal numbers in the first size interval; (c) third moment.

from the differential techniques are comparable but differ from the value obtained by the integral technique. 3.3. Evaluation of the Estimation Methods by CSD Simulation. The agreement between the measured and the model-predicted CSD for the different sets of kinetics estimates, reported above, is shown in Figure 3. Note that Figure 3a pertains to the batch crystallization experiment and Figure 3b to a constantcomposition experiment (CCSR) at 165 min. Nk,i was computed for each set of kinetics estimates using the DPB model described previously (q ) 6); for the batch crystallizer model, the flow terms were removed. Figure 3 shows that simulating the CSD with the kinetics obtained by the moments matching method results in a very poor agreement with the experimental

Figure 2. Crystallization kinetics determined by the differential and integral methods from the batch experiment crystallization data: (a) agglomeration kernel; (b) source term rate; (c) growth rate. Table 2. Averaged Agglomeration Kernel, Source Term Rate, and Growth Rate Estimates Estimated by the Different Estimation Techniques from the Constant Composition Crystallization Data estimation method

β (m3 s-1)

Bu (m-3 s-1)

integral 1.89 × 10-16 1.19 × 107 simultaneous differential 1.89 × 10-16 1.56 × 107 sequential differential 1.89 × 10-16 1.56 × 107 moments matching -1.19 × 10-16 -6.11 × 107 SMM 1.92 × 10-16 1.56 × 107

G (m s-1) 1.00 × 10-9 8.33 × 10-10 8.40 × 10-10 1.08 × 10-9 1.00 × 10-9

CSD. The estimates from the integral method give the best agreement regardless of the crystallizer configuration. The deviation between the experimental and

5010


purpose of simulating the CSD but do not give any indication of the precision (i.e., the uncertainty in the parameter estimates). This uncertainty arises from errors in the data being propagated by the data reduction method. The estimated kinetic parameters can be considered random variables, and the probability distributions for each of the estimates, from which the uncertainties can be determined, could be generated if sufficient replicate experiments were available. In this work, using the experimental uncertainties described below, 1000 computer-generated replicate CSD data sets, N, were obtained by loading the original CSD data, N h , with random normally distributed errors

N)N h + zσNi

Figure 3. Comparisons between the experimental and predicted CSDs using the kinetics estimates from different parameter estimation techniques for two experimental configurations: (a) batch and (b) CCSR, at 165 min. Table 3. Sum of Squares of the Differences, Ψ, between the Experimental and Predicted CSD for Batch and CCSR Data estimation method

Ψ (batch)

Ψ (CCSR)

integral simultaneous differential sequential differential moments matching SMM

3.2 × 1020 7.1 × 1020 7.1 × 1020 5.2 × 1021 N/A

1.6 × 1020 4.8 × 1020 4.8 × 1020 5.0 × 1022 2.2 × 1020

predicted CSDs was quantified as a sum of squares of the differences between the predicted and experimental Nk,i, i.e., Ψ defined by eq 18, with ωi set to 1. Table 3 shows that Ψ is lowest when the integral method’s kinetics estimates are used. The results at other sampling times for both batch and CCSR experiments also showed that Ψ was lowest for the integral method’s kinetics estimates. This is not surprising because Ψ is used by the integral method as the objective function for the estimation of crystallization kinetics. As expected, the highest Ψ was associated with the moments matching method. The estimates from both differential methods lead to Ψ values of 2-3 times those from the integral method. It is interesting to note that, for the constant-composition system, the estimates determined from the SMM result in lower Ψ values than those from the other differential methods. 4. Uncertainties in the Kinetics Estimates The above results indicate that the kinetics estimates from the integral method are more accurate for the

(20)

where z is a 1000 × ms × m array of normally distributed random numbers with zero mean and variance 1. More details on the application of Monte Carlo simulation to aluminum trihydroxide crystallization studies can be found in Li et al.11 4.1. Error in the Experimental Data. The errors in the solids concentration and liquor composition are relatively small. The main sources of the error in experimental CSD data are sampling error, sample preparation error, and particle size analysis instrument error (Coulter counter multisizer). Poor counting statistics will affect errors in some of the larger size intervals where particle numbers are small. Li et al.31 report on a reproducibility study based on six repeats of a batch aluminum trihydroxide crystallization experiment conducted at conditions similar to those being studied here. The key findings from this study were as follows: (1) The desupersaturation and CSD curves show good reproducibility. (2) The uncertainty in the measured crystal numbers in a size interval varies greatly and depends on the numbers in the size interval, and hence on the CSD. (3) If the tails of the CSD distributions are ignored, then the normalized standard deviation of the numbers in a size interval, σˆ Ni, which is the ratio of the standard deviation of the counts in a size interval (σNi) and the average of the counts (N h i), ranges from 5% to 15%. (4) The normalized standard deviations in the zeroth and third moments, σˆ µ0 and σˆ µ3, are around 5%. (5) The normalized standard deviations in the crystal number in the first size interval, σˆ N1, are between 20 and 30%. 4.2. Uncertainties in Kinetics Estimates from the Differential Methods. 4.2.1. Uncertainties in Kinetics Estimates from the Simultaneous and Sequential Differential Methods. CSD from the batch experiment reported in section 3 was replicated using eq 20 with σNi ) 0.1. Both different differential methods were then applied to the replicate CSD data sets, generating for each method a distribution of the respective kinetics estimates. The normalized standard deviations of these estimates were obtained by dividing the absolute value of the standard deviation by the mean of the estimate. The normalized standard deviations obtained for the simultaneous and sequential differential methods agree within 1% and are shown in Figure 4 as functions of the crystallization time. The values in Figure 4 have been smoothed to highlight the overall trends during the precipitation run. Note that the numerical values shown in Figure 4 are specific to the described experiment and may differ at other conditions.


Figure 4. Normalized standard deviations in the kinetics estimates obtained from a Monte Carlo simulation of a batch crystallization experiment using the simultaneous and sequential differential methods.

Figure 4 shows that the uncertainties in the agglomeration kinetics estimates, σβ, are low initially but increase rapidly with the crystallization time (i.e., decreasing supersaturation). The uncertainty in the growth rate parameter, σG, shows a similar behavior, but the sensitivity to changes in supersaturation is lower. On the other hand, σBu shows a more complex behavior; the uncertainty only moderately increases with time. Two factors that have a major effect on the relative uncertainties in the kinetics estimates can be identified: (i) rates of change of the total particle numbers (µ0), numbers in the first size interval (N1), and numbers in the crystal volume (µ3), and (ii) the relative sensitivity of the agglomeration, growth, and source term estimates. The first factor controls the error introduced by numerical differentiation, which can be considerable. The second factor is concerned with the relationship between the magnitudes of the absolute error and the kinetic parameters, which can lead to large relative errors when the parameter values are low, e.g., at low supersaturation. A sensitivity analysis of the equations for the simultaneous method (eqs 8-10) reveals that the uncertainty in the G parameter estimate is primarily controlled by the uncertainty in the derivative of the third moment, µ˘ 3. The contribution of the nucleation term in eq 9 is negligibly small relative to the growth term. The uncertainties in the estimates of β and Bu are influenced by the uncertainties in the time derivatives of both µ0 and N1 and through the growth kinetics also by the uncertainties in µ˘ 3. A detailed picture is complex because the coefficients in the simultaneous set of equations are also functions of the experimental CSD. Batch crystallization of aluminum trihydroxide under Bayer process conditions is typified by rapid changes initially (i.e., first 0.5 to 1 h) in µ0, N1, and µ3, as illustrated by Figure 1. This should be the period where the sampling frequency is greatest. Also, as shown in Figure 2 and observed by other authors,7,8,13,14 the magnitudes of β and G are significantly lower at lower supersaturations (longer precipitation times). These factors are consistent with the trends for the relative uncertainties in β and G shown in Figure 4. The trend for the relative uncertainty in Bu is more difficult to reconcile. However, unlike β and G, in this case, the

Figure 5. Normalized standard deviations of the kinetics estimates obtained by a Monte Carlo analysis of a batch experiment using the integral method.

absolute value of Bu does not change very much over the precipitation time; cf. Figure 2b. 4.2.2. Uncertainties in the Estimates from the SMM Method for Constant Composition Experiment. The uncertainties associated with the kinetics parameter estimates from the constant-composition experiment obtained using the SMM technique were estimated using the statistics from the linear regression analysis. The estimated normalized standard deviations σˆ β, σˆ G, and σˆ Bu were 10%, 16%, and 40%, respectively. They are of the same order of magnitude as those obtained by the Monte Carlo analysis, reported above, on the differential methods; however, the absolute values should not be compared because they pertain to different experiments. 4.3. Uncertainties in the Estimates from the Integral Method. A Monte Carlo simulation method was also used to assess the uncertainties of the kinetics parameter estimates obtained from the integral method. The batch crystallization example and the same error structure were used as described above in section 4.2. However, because the procedure was very computationally intensive, a smaller number (50) of computergenerated replicate data sets was considered. The relative standard deviations in the batch kinetic parameter estimates are presented in Figure 5. The trends are similar to those observed in Figure 4 for the analysis of the differential techniques. The relative uncertainties in the kinetic parameters are lower at the beginning of the crystallization run, where supersaturation and the parameter values are higher, and increase toward the end of a run. For the integral method, this observation can be explained by the greater sensitivity of the objective function to the values of initial kinetics, an effect that is well documented in studies investigating problems in optimal control (e.g., Teo and Jennings32). It was found that the method of interpolation between two sampling times affects the accuracy and uncertainty of parameter estimates. For example, it was observed that using piecewise constant kinetics instead of a linear interpolation over an interval between two sampling points reduces both the uncertainty and, unfortunately, the accuracy of the kinetics estimates. 4.3.1. Short-Cut for Estimating Parameter Uncertainties from the Integral Method. The parameter uncertainties associated with the integral method were also estimated using statistics from the nonlinear optimization procedure. Using approximate statistics,33 information on the standard deviations of estimated

5012


Figure 6. Uncertainties in the kinetics estimates obtained from the approximate statistics arising from application of the integral method on the batch experimental data.

kinetics can be obtained from the diagonal elements of the covariance matrix m

Vθ ) σr2(

∑ BkTBk)-1 k)1

(21)

where the variance of data, σr2, can be estimated from the residuals and Bk is a matrix defined by

Bk,iR ) -

∂ek,i ∂θR

5. Conclusions

(22)

and the residual ek,i is defined as

ek,i ) Nk,iexp - Nk,i

(23)

The standard deviations of the individual estimates can then be estimated from

σ ) xdiag(Vθ)

This somehow unexpected result can be explained by two main reasons. First, the uncertainties associated with the differential method are lessened because of the fact that kinetics estimation is based on the cumulative properties of the size distribution, i.e., its moments. This consequently reduces the impact of a random measurement error on the corresponding kinetics estimates. In the case of Bu, however, the estimation of kinetics heavily relies on the time derivatives of the number of particles in the first size interval. This leads to higher uncertainties associated with the estimation of source term rates in the case of the differential method. In the integral method, numbers in individual size intervals are directly used in the estimation. Although a disadvantage from the point of view of uncertainties, this leads to a better agreement between the experimental and model-predicted size distributions. At the same time, this may be the main reason for higher uncertainties associated with β observed in the case of integral method. On the other hand, however, the estimation of Bu was proven to be much more robust, resulting in lower uncertainties. This can largely be credited to a feedback nature of the integral method, where a small change in the source term leads to a completely different CSD evolution, rather than affecting only the first size interval of a particle distribution.

(24)

The standard deviations in the parameter estimates computed following this procedure will be a rough approximation, correct to within an order of magnitude. Figure 6 shows the parameter uncertainties estimated following this procedure. The results are comparable to those obtained by the Monte Carlo simulation, with only the uncertainties for Bu notably overestimated. Tests using a number of different data sets found that this approach gives reasonable estimates of uncertainties provided that the uncertainties themselves are small but significantly overestimates uncertainties once they become larger. This observation is consistent with the approximation of the objective function value, implicit in this analysis, which is of acceptable accuracy only in the vicinity of a minimum. 4.4. Comparison of Uncertainties from the Differential and Integral Methods. From the presented results, it is evident that the kinetics estimates from both methods, differential and integral, are associated with considerable uncertainties. In addition to experimental errors, they mainly originate from either the estimation of the time derivatives or from the nature of the optimization problem. Comparison of the uncertainties from both methods shows that they are of similar magnitude, particularly for estimates of growth kinetics. However, the results also show that the integral method leads to lower uncertainties in Bu, whereas the differential method tends to give lower uncertainties for β.

A number of differential and integral kinetics estimation techniques were evaluated in terms of the accuracy and precision of the estimates from dynamic crystallization data. The techniques were applied to the estimation of the agglomeration, growth, and nucleation kinetics from laboratory aluminum trihydroxide crystallization experiments. A comparison of the estimates shows that in the case of the crystal growth rate the methods are generally in good agreement but differ more in the case of agglomeration and source term rates. The integral method was found to lead to the best agreement between the model and experiment. While the simultaneous and sequential differential techniques were the most accurate among the differential methods, the moments matching method was found to be inadequate. The SMM, however, was accurate but only suitable for the special case of constant kinetics. Uncertainties in the kinetics estimates from the various methods were estimated by Monte Carlo analysis. It was found that the uncertainties of kinetics estimates are determined by the interplay of the following factors: the rate of change of an observation, the relative contribution of the kinetics, and, in the case of the integral method, the formulation of the optimization problem. The uncertainties in the kinetics estimates were found to be largest at the end of a batch crystallization run, where supersaturation is depleted and rates are low. The uncertainties in the parameter estimates were comparable for the three most accurate techniques, namely, the integral, simultaneous differential, and sequential differential methods. Acknowledgment The authors gratefully acknowledge support from the Australian Governments Cooperative Research Centre (CRC) program, through the AJ Parker CRC for Hydrometallurgy. The authors express their gratitude and acknowledge their debt to Prof. E. T. White from University of Queensland for seeding interest in this


“dry” but overlooked and important subject and for the many stimulating discussions on industrial crystallization problems. Nomenclature A ) alumina concentration, g L-1 Al2O3 A* ) equilibrium alumina concentration, g L-1 Al2O3 B ) residual matrix B ) birth rate due to agglomeration, µm-1 m-3 s-1 Bo ) nucleation rate, m-3 s-1 Bu ) source term, m-3 s-1 CSD ) crystal size distribution CCSR ) constant-composition solids retention D ) death rate due to agglomeration, µm-1 m-3 s-1 DPB ) discretized population balance e ) residual error G ) linear growth rate, m s-1 L ) crystal size, µm L h 1 ) mean size of the first size interval, µm L h ) array of mean sizes in eq 13 m ) number of sampling points ms ) number of size intervals n ) number density function, µm-1 m-3 N ) number of crystals in each size interval, m-3 N1 ) number of crystals in the first size interval, m-3 PB ) population balance q ) adjustable parameter determining the size discretization resolution Q ) volumetric flow rate, m3 s-1 Ra ) discretized aggregation rate SMM ) simplified moments matching t ) time, min or h V ) precipitator volume, m3 Vθ ) covariance matrix Greek Letters β ) agglomeration kernel, m3 s-1 λ ) crystal size, µm θ ) vector containing kinetic parameters σ ) standard deviation σr ) standard deviation from residuals, eq 21 µj ) jth moment of the CSD, µmj Ψ ) objective function ωi ) weight for the ith measurement in eq 18 Subscript R ) element of vector θ

Literature Cited (1) Bramley, A. S.; Hounslow, M. J.; Ryall, R. L. Aggregation during Precipitation from Solution: A Method for Extracting Rates from Experimental Data. J. Colloid Interface Sci. 1996, 183, 155165. (2) Livk, I.; Ilievski, D.; Pohar, C. Estimation of Batch Precipitation Kinetics Using a Simplified Differential Method. AIChE J. 1999, 45 (7), 1593-1596. (3) Rawlings, J. B.; Miller, S. M.; Witkowski, W. R. Model Identification and Control of Solution Crystallisation Processes: A review. Ind. Eng. Chem. Res. 1993, 32, 1275-1296. (4) White, E. T.; Bateman, S. H. Effect of Caustic Concentration on the Growth Rate of Al(OH)3 Particles. Light Met. 1988, 257-162. (5) Veesler, S.; Boistelle, R. Growth Kinetics of Hydrargillite Al(OH)3 from Caustic Soda Solutions. J. Cryst. Growth 1994, 142, 177-183. (6) Li, T. S.; Rohl, A. L.; Ilievski, D. Modelling a Constant Supersaturation Precipitator. In Mixing and Crystallisation; Sen Gupta, B., Ibrahim, S., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2000; pp 227-241. (7) Ilievski, D.; White, E. T. Agglomeration during Precipitation: Agglomeration Mechanism Identification for Al(OH)3 Crystals in Stirred Caustic Aluminate Solutions. Chem. Eng. Sci. 1994, 49, 3227-3239.

(8) Misra, C.; White, E. T. Kinetics of Crystallisation of Aluminium Trihydroxide from Seeded Caustic Aluminate Solutions. Chem. Eng. Prog. Symp. Ser. 1971, 67 (110), 53-65. (9) Misra, C.; White, E. T. Crystallisation of Bayer Aluminium Trihydroxide. J. Cryst. Growth 1971, 6, 172-178. (10) Melikhov, I. V.; Podkopov, V. M.; Ilyin, B. A.; Kozlovskaya, E. D. Periodical Crystallisation Variability. Chem. Eng. Sci. 1996, 51 (5), 671-682. (11) Li, T. S.; Livk, I.; Ilievski, D. The Influence of Crystalliser Configuration on the Accuracy and Precision of Gibbsite Crystallisation Kinetics Estimates. Chem. Eng. Sci. 2001, 56, 2511-2519. (12) Tavare, N. S. Industrial Crystallisation: Process Simulation Analysis and Design; Plenum Press: New York, 1995. (13) Halfon, A.; Kaliaguine, S. Alumina Trihydrate Crystallization: Part 1. Secondary Nucleation and Growth Kinetics. Can. J. Chem. Eng. 1976, 54, 160-167. (14) Halfon, A.; Kaliaguine, S. Alumina Trihydrate Crystallization: Part 2. A Model of Agglomeration. Can. J. Chem. Eng. 1976, 54, 168-172. (15) Mullin, J. W. Crystallization; Butterworth-Heinemann: Oxford, U.K., 1993. (16) Mumtaz, H. S.; Hounslow, M. J.; Seaton, N. A.; Paterson, W. R. Orthokinetic Aggregation during Precipitation: A Computational Model for Calcium Oxalate Monohydrate. Chem. Eng. Res. Des. 1996, 75 (A2), 152-159. (17) Randolph, A. D.; Larson, M. Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallisation, 2nd ed.; Academic Press: Toronto, 1988. (18) Ilievski, D. Development of a Constant Supersaturation Semi-Batch Crystalliser and Its Application to Investigating Agglomeration. J. Cryst. Growth 2001, in press. (19) Marchal, P.; David, R.; Klein, J. P.; Villermaux, J. Crystallisation and Precipitation Engineerings1: An efficient Method for Solving Population Balance in Crystallisation with Agglomeration. Chem. Eng. Sci. 1988, 43, 59-67. (20) Hounslow, M. J.; Ryall, R. L.; Marshall, V. R. A Discretized Population Balance for Nucleation, Growth, and Agglomeration. AIChE J. 1988, 34, 1821-1832. (21) Litster, J. D.; Smit, D. J.; Hounslow, M. J. Adjustable Discretized Population Balance for Growth and Agglomeration. AIChE J. 1995, 41, 591-603. (22) Wynn, E. J. W. Improved Accuracy and Convergence of Discretized Population Balance of Litster et al. AIChE J. 1996, 42, 2084-2086. (23) Tavare, N. S.; Garside, J. Simultaneous Estimation of Crystal Nucleation and Growth Kinetics from Batch Experiments. Chem. Eng. Res. Des. 1986, 64, 109-118. (24) Li, T.; Rohl, A.; Ilievski, D. Modelling Non-Stationary Precipitation Systems: Sources of Error and Their Propagation. Chem. Eng. Sci. 2000, 55 (24), 6037-6047. (25) Farrell, R. J.; Tsai, Y. C. Modelling, Simulation and Kinetic Parameter Estimation in Batch Crystallisation Processes. AIChE J. 1994, 40 (4), 586-593. (26) Miller, S. M.; Rawlings, J. B. Model Identification and Control Strategies for Batch Cooling Crystallizers. AIChE J. 1994, 40 (8), 1312-1327. (27) Hostomsky, J.; Jones, A. G. Calcium Carbonate Crystallisation, Agglomeration and Form during Continuous Precipitation from Solution. J. Phys. D: Appl. Phys. 1991, 24, 165-170. (28) Edgar, T. F.; Himmelblau, D. M. Optimization of Chemical Processes; McGraw-Hill: New York, 1988. (29) The Mathworks Inc. Optimization Toolbox User’s Guide; The Mathworks Inc.: Natick, MA, 1997. (30) Watts, H. L.; Utley, D. W. Sodium Gluconate as a Complexing Agent in the Volumetric Analysis of Aluminium Compounds. Anal. Chem. 1956, 28, 1731-1735. (31) Li, T. S.; Rohl, A. L.; Ilievski, D. An Investigation on the Uncertainty in Batch Gibbsite Precipitation Experimental Data. 28th Australasian Chemical Engineering Conference, Perth, Australia, 2000; Paper 151. (32) Teo, K. L.; Jennings, L. S. Optimal Control with a Cost on Changing Control. J. Optim. Theory Appl. 1991, 68 (2), 335-357. (33) Bard, Y. Nonlinear Parameter Estimation; Academic Press: New York, 1974.

Received for review March 5, 2001 Revised manuscript received August 15, 2001 Accepted August 24, 2001 IE0102145

Influence of the Estimation Procedure on the Accuracy and Precision

Recommend Documents