- Email: [email protected]
Estimation of nucleation kinetics from crystal size distribution transients of a continuous crystallizer
Chemical
Engineerme
Printed in Great
Snencc
Vol.
46. No.
3. pp. 807
818. 1991.
KG?-2509191 %3X0 + 0.00 8 1991 Pergamon Press pk
Britam.
ESTIMATION OF NUCLEATION KINETICS FROM CRYSTAL SIZE DISTRIBUTION TRANSIENTS OF A CONTINUOUS CRYSTALLIZER JOHAN Laboratory
JAGER,
SJOERD
for Process Equipment,
DE
WOLF,’
HERMAN
J. M. KRAMER*
and ESSO J. DE JONG Delft University of Technology, Netherlands
(First received 30 August
1989; accepted
Leeghwaterstraat
in revised form 8 January
44, 2628 CA Delft,
1990)
Abstract-Four techniques which can be used for the estimation of nucleation kinetics under unsteadystate conditions are applied to transients of a 20-l evaporative ammonium sulfate crystallizer. A moment and a Laplace transformation technique appear not to be feasible because they rely heavily on the measurement of the number of small crystals. The present crystal size distribution (CSD) measurement technique is inaccurate in this size range and therefore both methods suffer severely from measurement noise. The third technique, using the method of characteristics, indicates the appearance of so-called secondary nucleation kinetics. The technique is unsuitable for estimating the kinetic parameters quantitatively, again as a result of the measurement noise. Therefore a fourth technique, which uses a non-linear parameter estimation routine, has been developed. This technique is very successful in estimating the parameters in a secondary nucleation model. With the estimated parameters the second and third moments of the CSD can be predicted. If additional moments are to be predicted it is shown that other phenomena which affect the crystal size distribution, such as crystal attrition, must be included in the dynamic model of the process.
to describe the nucleation rate (B) by an empirical power law. This power law comprises the supersaturation (S), which is considered to be the main driving force for nucleation, the crystal surface (k,m,), crystal mass (k,p,m,) or a higher moment mj, which reflects the idea that the presence of crystals plays a dominant role (secondary nucleation), and the stirrer speed (N) which governs fracture processes (Daudey, 1987):
INTRODUCTION
unit operation in chemical engineering practice. It is employed in the production of bulk chemicals such as ammonium sulfate, sodium chloride and sucrose, and the production of fine chemicals such as aspartame. Usually the process produces a solid product which has a wide crystal size distribution (CSD). This CSD dictates the behaviour of the product in succeeding operations, such as filtration, drying, transport and storage, and is also used to define customer specifications. It is, therefore, desirable to control the CSD which is produced by the process. The first requirement is an on-line CSD measurement technique. Recent developments have made this possible (Jager et al., 1989). Another requirement for process control is an accurate dynamic model. This mode1 can be obtained by estimating the parameters of “a black box” model or, alternatively, by determining a model based on physical laws. This paper proceeds with the second approach. Crystallization can be described in terms of dynamic heat and mass balances, combined with a dynamic population balance which describes the CSD. The population balance contains specific population events such as crystal growth and crystal birth. The latter is especially a complicated process involving fracture mechanics, two-phase flow hydrodynamics and solid state physics. Hence, the usual approach is
Crystallization
from
solution
is an
important
B = kbS”(?nj~N”.
(1)
If the supersaturation is difficult to measure, it may be exchanged by the growth rate G, which is a function of the supersaturation: G = k,S’
(2)
B = k,Gi(mj)l’N”.
(3)
The usual approach is to perform steady-state experiments under various operating conditions. From the resulting CSDs, the growth and nucleation rates can be estimated, often by assuming the distribution to be log-normal (Randolph and Larson, 1988). By curve fitting, the empirical constants can be found for different moments (Shah and Garside, 1980). There are two objectives in extending this analysis to unsteady-state conditions: (1) Under steady-state conditions the moments for an idealized MSMPR (mixed suspension mixed product research) crystallizer, defined in the Appendix, at a residence time z, correlate according to
‘Laboratory for Measurement and Control, Delft University of Technology, Mekelweg 2, 2628 CD Delft, Netherlands. ‘Author to whom correspondence should be addressed.
mi = jmj807
,rG.
(4)
808
(2)
JOHAN
JAGER
Therefore, it is not possible to determine the constants in the nucleation model from steadystate data. The verification of secondary nucleation kinetics under transient conditions. Recently, experimental evidence has been presented which indicates that the steady-state constants do not describe correctly the nucleation process for the unsteady state (Heffels, 1986). The occurrence of a so-called nuclei stock, which contains slow or non-growing crystals, and which can be activated at momentarily higher supersaturations, was postulated (Heffels, 1986).
However, the growth and nucleation rates, and consequently the parameters in eq. (3), cannot easily be evaluated under unsteady-state conditions. Four techniques which could serve this purpose are explored in this paper, using experimental transients produced in a 20-l continuous evaporative crystallizer. A moment and a Laplace technique were applied to CSD transients from batch and continuous crystallizers by Tavare and Garside (1982). However, the experimental scatter in their results was pronounced. The method of characteristics was first reported by Daudey and Jong (1984), but with loo few results to judge its merits. The fourth technique, which uses a general non-linear parameter estimation technique, has been developed in order to estimate the parameters of the nucleation kinetics. Each of these techniques requires a dynamic model which wil1 be discussed first. THEORY
Dynamic modeling The crystallizer contents of an MSMPR crystallizer are fully defined by the heat balance, the total mass balance, the solute balance and the population balance (Wolf et al., 1989). The differential equations which represent the dynamic heat, mass and solute balances can be simplified by applying two constraints for a constant-volume, isothermal and ClassII crystallizer with size-independent crystal growth (Appendix). The Class-II assumption implies that the ammonium sulfate system is a fast growing system in which the crystallizer concentration is constant with time and approaches the equilibrium concentration. The assumption of size-independent crystal growth can be implied, for the intermediate size range, from a log-normal population plot. The heat and mass balance simplify into a relation for the product flowrate (Q,), as a function of the feed flow (Q,), the feed concentration (C,), and the heat input (P): Qp(tl
=
Cki - kzCj(t)lQf(t) - p(t)
et al.
G(t)
(5)
The solute balance reduces with the Class-II assumption to a relation for the growth rate(G), as a function of the crystallizer volume (V), the crystallizer concentration (C), and the second moment (m2):
Q,(t)C
-
(6)
- C) ’
3k, W(t)(p,
At a constant feed flow, eq. (5) predicts that the product flow will decrease upon an increase in the heat input. Consequently, the growth rate will immediately increase upon an increase in the heat input according to eq. (6), and then decrease, because of the subsequent increase in the second moment. The dynamic population balance describes the CSD transients. If crystal breakage, attrition and agglomeration are neglected, crystals are formed at very small size, and all crystals have the same growth rate (McCabe, 1929), the well-known population balance (Randolph and Larson, 1988) for the population density, n, at size L results:
anw,t) + n(L, t) Q,(t) = 0. ad5 t) ~ + G(t)~ V at 3L The boundary differential
condition
equation
for this hyperbolic
(7)
partial
is given by no(t)
=
B(f) G(f)
which changes,
by substitution
n,(r) The jth moment
of eq. (3), into
= k,G’-‘(t)mS(t)N”. of the distribution
(9) is defined
by
cc n(L, 2)LjdL.
mi(t) =
(10)
I0 Techniques Three of the techniques which are discussed, denoted the method of moments (Tavare, 1986), the Laplace method (Tavare, 1986; Tavare and Garside, 1982), and the method of characteristics (Daudey and Jong, 1984), calculate the growth and nucleation rates appertaining during an experiment. The parameters in eq. (3) can then be estimated using a linear least-squares routine. The fourth method, which combines two mathematical techniques, the method of lines and non-linear parameter estimation, estimates the unknown parameters without the intermediate step of calculating the growth and nucleation rate. The fourth method will be denoted non-linear parameter estimation. The method ofmoments. This method uses the moment transformation to transform the population balance into a set of ordinary differential equations describing the moments of the distribution:
mo(t)Q,(t)
dmdt) =dt
K3
QJt)c/(t)
=
dmj(t) -= dt
V
+ B(t)
- mj(t)Qp(t) . +Jmj_,(t)G(t) V
(11) for j z 1.
(12) By replacing proximations,
the derivatives with finite-difference apthe nucleation rate can be calculated
Estimation of nucleation kinetics from crystal six distribution from the transient zeroeth growth rate can be calculated and higher moments: Am,(t
moment whereas from the transient
+ fAt, + mo(t + fAt)~,(t) At
V
the first
= B(T + +At) (13)
Amj(t + &At) + mj(t + )A~)Q~(c) At V = G(t + $At) jfi,_ ,(r + +At) (14) where Amj(t + + At) = mj(r) - mj(t + At)
m,(t +
fAt,
=
CmAt)+ mAt + “)I
J
2
In this paper the growth rate has been calculated from eq. (6) rather than from eq. (14) in order to avoid numerical differentiation, which is very sensitive to noise on the experimental data. Unfortunately this cannot be avoided for the nucleation rate. The calculated nucleation rate will depend on the measured zeroth moment which is dominated by the small crystals. An average rather than an effective nucleation rate, therefore, results. The Laplace method. The result of the Laplace transformation of the population balance is a firstorder ordinary differential equation in the transformed population density, n’(s, t), using the Laplace variable s: n’(s, t) =
CI) n(L,t)exp(--sL)dL. s0
dn’(s, t) + n’(s, t)Q,(t) dt
V
(15)
= - G(t)sn’(s,
t) + B(t).
(16)
The transformed population density can be calculated from the measured population density. The derivative in eq. (16) can be calculated using a finite-difference approximation. Over a range of s values a linear relation results, with B(t) and G(t) as unknown parameters The slope of the straight line Y(t) = B(t) - G(t)X(t)
transients
for CSD measurement techniques which are most accurate in the intermediate size ranges. If negative s values are used, the larger size ranges are emphasized. At first sight this might prove to be favourable: however, it also puts a very strong emphasis on the “tail” of the distribution. The weight function [exp ( - sL)] approaches infinity as L tends to infinity, where the measurement accuracy is nil. The s values to be taken cannot be estimated in advance but must be evaluated by experiment, which introduces uncertainties. This second method also calculates an average nucleation rate. Method of characteristics. The method of characteristics seeks special curves in the L-t plane, called characteristic curves, along which the solution of the partial differential equation is reduced to the integration of an ordinary differential equation. Within the present context, the partial differential equation is transformed into an ordinary differential equation along curves where eq. (18), for the crystal size of a specific class of crystals (L*) holds: (18) By eq. (18) a collection of curves, denoted characteristics, is defined. Along these characteristics the population density in time for this specific class of crystals (n*) is given by (19) These two equations are easily solved. The growth rate G can be calculated from eq. (6). Using this growth rate, the birth time (tbirlh) of a particular size class attained at a later time (t,,,,l.) can be estimated using flami.lc L*( Llnple 1 -
G(t)dt
= 0.
s fbirth
(20)
The population density at the time of birth can be calculated, using the product flowrate from eq. (5), by integrating eq. (19):
(17)
nO*(tbirth)
=
Lmpk - tmJQp(tl n*(Lampk)exp V
where
(21) dn’(s, r) + n’(s, t)QJt) y(t) =dt x(t)
V-~ _______
= sn’(s, t)
represents the growth rate while the abscissa will give the nucleation rate. Closer examination of the transformation reveals that n’(s, t) results from weighting of the measured population densities by [exp ( - sL)], and integrating over the size range L. This gives the opportunity to weight the data and reduce the measurement noise. For positive s values this results in a weight function which decreases with size, which is unfavourable
and
B*(tbirth)= 4 (tbirth)G(tbirth)+
(22)
For a measuring device with 32 size classes, the nucleation rate can be calculated for each of these classes. The birth-rates for the different size classes do not necessarily coincide and interpolated values therefore must be compared. These are identical if McCabe’s AL law holds. However, if growth dispersion occurs the nucleation rate will differ for different size classes. Growth dispersion states that crystals have different growth rates when exposed to the same supersaturation (Janse and Jong, 1976; Berglund
810
JOHANJAGER~~~I.
and Larson, 1984). If this phenomenon is restricted to the smaller size range, the nucleation rates calculated from size classes exceeding this range can be interpreted as an effective nucleation rate. The effective nucleation rate accounts for the number of nuclei, which grow to a for the sizing technique detectable size only (Jancic and Grootscholten, 1984). At larger size ranges the nucleation rate will be affected increasingly by measurement noise. The difference. between the sample time and the time of birth which appears in the exponential of eq. (21) is larger for larger sizes. This tends to promote an error, due to the time of birth term in the calculated nucleation rate. For small si7e ranges the measured population density is itself of limited accuracy. The intermediate size ranges are, therefore, to be preferred.
Non-linear parameter estimation. Rather than calculate the growth and the nucleation rates separately and then fit them to a model, the first stage can be avoided by estimating the unknown parameters in eq. (9) directly. A linear, least-squares routine cannot be used for this purpose, since the parameters occur in the boundary condition of a hyperbolic partial differential equation. An alternative technique must there-
fore be used. Using starting values for the unknown nucleation parameters one can simulate the overall model consisting of the population balance [eq. (7)], the nucleation model [eq. (9)], the growth rate constraint [eq. (6)] and the constraint for the product flowrate [eq. (511. The kinetic parameters can be used to adjust the simulated to the experimental responses. The adjustment is made by a non-linear optimization technique. The procedure is terminated when a chosen criterion is met. Three aspects of this technique will be discussed in detail. Firstly, the technique which is used to simulate the responses will be discussed. Subsequently, the number of parameters to be estimated is, by using steady-state information, reduced to one. This is a convenient way to adjust the solution to the initial and final values of the step responses which were investigated. The use of steady-state information does not compromise the first objective of this paper, because eq. (4) is certainly violated immediately after a step in heat input [see eqs (6) and (12)]. The criterion which was used in the optimization will finally be discussed. The partial differential equation [eq. (7)] cannot be solved numerically using standard techniques and it must therefore be transformed. It was decided to choose a technique which is not restricted to simple crystallizer configurations, and the moment transformation was therefore not considered (Marchal et a(., 1988). An alternative technique to transform eq. (7) into a form which lends itself to simulation is the method of lines (Cooper and Clough, 1985; Tsuruoka and Randolph, 1987; Wolf er al., 1989). In this procedure the size axis is discretized and the length
derivative is substituted by a finite-difference approximation. Eventually, a number of first-order differential equations for n(L) in discrete gridpoints replaces the original partial differential equation. Several finite-difference approximations can be chosen, depending on the required accuracy. Throughout this paper the size axis is discretized into x + 1 points and a fourth-order central difference, based on an equidistant grid, is used to approximate the derivative with respect to crystal size:
WL,
t)
=~r(n(L,-,,r)-sn(L,~,,t)
dL
+ WL+,,
t) - fi(L,+,,
t)i + o(~L4).
(23)
This central difference cannot he used at the boundaries of the grid L,, L,, Lx_ 1 and L,, since n(L_,, I), t) and n(L,+*, t) miss. Therefore n(L-I, t). n(L,*,, eq. (23) is replaced by two forward differences on the left boundary [eqs (24a) and (24b)] and two backward differences on the right boundary [eqs (25a) and (25b)] (Wolf et al., 1989):
hdt)
&
-z
a_c
[-25no(t)
- 36n(L,,
+ 48n(L,,
t) + 16n(L,,
t)
t) - h(L,,
t)]
+ 0(AL4)
anw,. -=
(24a)
t) &
aI5
[-3n,(tI
+ 18n(L,,
-
lOn(L,,
t)
t)
t) - 6n(L,,
+ n(.J+, t)] + O(AL4) WL,-
1,
JL
t)
= &
c -n(L,-.,
- 18n(L,_,,
Wb)
t) + 6n(I+
J, t)
t) + IOn(L,_ I, I)
+ 3n(L,, t)] + O(AL“)
ML, -= aL
t)
&
C+3$&.,
2) - 16n(L,-,,
+ 36n(L,_ 2, t) - 48n(L,_ + 25n(L,,
(25a)
t)] + O(AL’).
t)
I, t) (25b)
Replacing eq. (24a) by the boundary condition for no, and substituting eqs (24b), (25a) and (25h) into eq. (7) then results in a set of ordinary differential equations at discrete gridpoints L,. In the overall model, the population density must be integrated numerically in the calculation of m2 in eq. (6) and m, in eq. (9). The size axis was divided into 8, 31, 101 and 201 gridpoints, lying between zero size and a maximum of 2 mm, respectively. The changes in the simulated responses where negligible upon an increase from 31 to 101 gridpoints. Thus, the condition of grid independence is satisfied for both the responses and the estimated parameters at the current number of 201 gridpoints. The set of differential equations is solved using a Runge-Kutta routine.
Estimation of nucleation kinetics from crystal size distribution At a constant stirrer speed, the term N” in eq. (9) can be included in the nucleation constant k,. Then, for likely values of j equals 2 or 3, three parameters (k,, i, k) must be estimated in the boundary condition [eq. (9)]. At steady state, n, can be calculated from the jth and the (j - l)th moment: .j+lmj’:
n, =
_I
/+
If two steady-state conditions (1) and (2) are known correctly from experimental data, two degrees of freedom are removed by two constraints:
logn,,
I = logk,
+ (i -
1)log G1 + klog(m,),
(27)
log n,, 2 = log k, + (i - 1) log G, + k log @I~)~.
(28)
If k is the parameter to be optimized, i and k, are determined by these constraints according to Klog (mJ)l - log (m&l i= log no, I - log no, 2 (log G, - log G,) + 1 (29) logk,
= log no*, -(i-
l)logG,
-klog(m,),.
(30)
The optimization of k can then be achieved with any non-linear line search technique. A technique based upon the Fibonacci series (Gill et al., 1981) was chosen because of its ease of implementation. Either one moment or a combination of moments can be optimized. Typically, simulated datapoints are subtracted from the experimental ones, squared and summed. This sum is minimized. Weighting factors can then be used to emphasize the moments which are considered important: Sum = i
i:
x=l j=I
wjCmj..h(~,) - mj,ex*(tx)12 (31)
where n = number
of datapoints, z = number of moments, wj = weighting factor for jth moment. EXPERIMENTAL
STRATEGY
Throughout this investigation a 20-l continuous evaporative crystallizer which is operated at about 100 mbar and a stirrer speed of 700 rpm was used. A detailed description of the experimental set-up is given by Daudey (1987). A micro-computer which controls the temperature, input flow, and heat load was added to the original configuration in order to reduce undesired disturbances caused by random variations in the process conditions. This crystallizer
transients
was operated at a temperature of 50°C and manually sampled every 5 min. The resulting samples were subsequently filtered, washed and dried. The total solids content was evaluated using the sample volume and the amount of solids after the filtration procedure. This stage proved to be very important for the crystal size analysis. Opaque samples give different results compared to completely transparent ones. It is expected that the main cause of the measurement noise is due to this phenomenon. The samples were analyzed 3 times using a Malvern 26OOc laser diffraction particle sizer and the results were averaged. Differences between these three measurements were small compared to differences between samples. In these investigations, the Malvern analyzer was fitted with a lOOO-mm lens which is appropriate for a measurement range of 0.519OOpm, which is divided into 32 size classes (Jager et al., 1990). In order to match the ammonium sulfate system, the software of the instrument was adapted (Boxman and Scarlett, 1988). The crystals were analyzed in the dry state using a dry-powder feeder. In order to avoid classification during this procedure, care was taken to measure the complete sample. An on-line CSD measurement technique was available on a 970-l pilot plant crystallizer. It could not be applied on the 20-I crystallizer because the volume contained by the dilution unit is also approximately 20 1. A typical experimental procedure was as follows. First, a start-up of the crystallizer at the required process conditions. After about eight residence times one or two stepwise changes in heat input or input flowrate were made, and then the experiment was continued for about eight further residence times. Over 15 such experiments were conducted. In this paper three of these experiments will be discussed, all of which consisted of changes in heat input. The experimental conditions are listed in Table 1. The changes in heat input effect changes in the vapour flow which, in turn, result in variations in the residence time, since the input flow was kept constant. The feed of the crystallizer is saturated at 46°C. The population densities and the moments of the distributions were calculated from the measured volume fractions (AV), the experimental magma density (MS!), the average size in each interval (L) and the interval widths (AL) given elsewhere (Jager et al., 1990): A V(L,, t) Msl( t) n(L,
t) =
p,k,L:AL Msl(t)
mj(t) = ~
32 AV(L,, c
X=l
p
Table I, Experimental conditions for three different experiments Experiment steps at (s) Heat input (kW) Residence time (s)
1 17,700 3.4-2.0 2450-2270
811
2 16,200-19,200 1.9-5.3-l .9 24OG293~24Cil
3 16.500 2.C3.3 2270-2430
L,3
t)L; (33)
812
JOHAN
JAGER
et
al.
RESULTS
1
The results of the first experiment were typical for most of the experiments and, unless otherwise stated, have been used to produce the graphs in this paper. From eqs (32) and (33) it is clear that the total amount of solids (Msl) is used to calculate the population densities and the moments. Therefore the measurement error in Msl contributes to that in the population densities and the moments. The amount of solids, Msl, was determined experimentally. However, it can also be calculated. Substituting eq. (6) into eq. (12) for j = 3 yields an ordinary differential equation in the third moment:
m,(t) + CfQJt)
dm,(t) -=-dt The solution is given by
- CQJt)
k,(p, - C)
r(t)
v
.
The total amount
CQ, I - C,Q,s t m%--C)V .
71-1-
CQ,., Up,
15
20
2s
-
(34) change
in heat
input
30
3s
[setI ~101
tlme
LO
vs predicted slurry density
for a step of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
hhl
-?)I
x [I--exp(
(9)~ = 72
10
5
Fig. 1. Experimental
+ C(m& -
(m3)1=
kt.,
0 0
of eq. (34) for a step in heat input (1 + 2)
m,(t - tstep) = (m&
2o t
- C_rQ,,2 - C) v
(35a)
800
25
T
'tstcp
(35’4 (35c)
of solids is given by
Msl = pek,ms.
(36)
Both the experimental and the theoretical values are shown in Fig. 1 and show a very good agreement. In order to reduce the measurement noise in the other moments and the population densities, they have been calculated from eqs (32) and (33) using the predicted slurry density (Msl) [eqs (35) and (36)] rather than the experimental values. The growth rate, combined with the second moment from which it was calculated according to eq. (6), is depicted in Fig. 2. A very high growth rate existed during the beginning which decreased gradually towards the first steady state. On the step in heat input, the growth rate decreased sharply, after which a new steady-state value was eventually attained. Subsequently, the nucleation rate was calculated using the different calculation methods.
Method of moments Figure 3 shows the results using the method of moments. Hardly any trends can be discerned. The zeroth moment from which the nucleation rate is calculated is too distorted by the measurement noise. Even a fourth-order finite-difference approximation for the derivative did not improve the results. Data filtering, which was subsequently used to reduce the measurement noise, was not effective, without substantial loss of information, for the number and spacing of data points.
200 100
1
0 'I
0
, 0
s
10
15
tstrp 20
2s -
30
35
J.0
time [set I xl03
Fig. 2. Growth rate, and the second moment growth
L&ace
from which the rate is calculated, for a step change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
transformation
Figure 4 shows the Laplacian nucleation rates, which show the same tendency. They are calculated with 100 s values, lying between zero and a value of SL 50 = 1 as proposed by Tavare (1986). No improvement results if other values are used. This results from the weighting of the data. Table 2 shows the population density of a steady-state sample with no = 7.7 x lOI and Gr = 1.6 x 10m4, the corresponding weighting factors, and their cumulative contribution to the transformed population density for a value of SL 50 = 1. The cumulative contribution is given by
(37)
The measurement noise is greater for the smaller size ranges and the weighting by the Laplace method is,
Estimation of nucleation kinetics fmm crystal size distribution
transients
813
5 -20
,trtq 0
5
10
15
20 -
25 30 35 tint [sec.]x10)
Fig. 3. Nucleation rate calculated with the method of moments for a stcD change in heat inout of 3.4 kW to 2.0 kW at 17,700 s in expe%ment 1.
trt.,
Fig. 4. Nucleation for a step change
0
60
0
5
10
15
20
25
-
time
30 [sec.]
35
LO
~103
Fig. 5. Nucleation rate calculated with the method of characteristics from the sixth (62.5-75 gm), the twentieth (6C@-7OOpm) and the thirtieth size class (160&17OQ~m), respectively, for a step change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
I
20
25
-
tim
30 lsrc.1
35
1
-103
t*t.r
rate calculated with the Laplace method in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
5
10
15
20
25
30
35
Fig. 6. Nucleation rate calculated with the method characteristics from the nineteenth (500-600 bm), Table 2. Cumulative contribution of different size classes to the overall transformed population density for a theoretical crystal size distribution
n&J Class 1
6 11 16 21 26 31
(No./m’) 7.4 5.0 2.8 1.0 6.9 3.0 1.3
x x x x x x x
10’3 1OL3 lot’ IOL’ lOL1 lOL0 IO9
exp( - s&A 0.99 0.90 0.76 0.60 0.31 0.14 0.06
“bi”’ 0
9.0 44.4 74.7 93.7 99.8 loo.0 loo.0
therefore, ineffective for the measurement method employed throughout these investigations. The calculated nucleation rates vary strongly and even have negative values. Moreover, the time of the step cannot be recovered.
of the
twentieth (600-700 pm) and the twenty-first size class (700-800 pm), respectively, for a step change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
Method
c~~characteristics
Figure 5 shows the results for the sixth, twentieth and thirtieth size class, which range from 62.5 to 75, 600-700 and 16OC-1700 pm, respectively. The scatter for the sixth and thirtieth size classes is pronounced, which is in accordance with the prediction earlier made. Figure 5 also shows an additional feature of this technique. For times greater than 13,218 s, nucleation rates are not available for the thirtieth size class. These nucleation rates should have been caIculated from samples which were beyond the duration of the experiment. The reproducibility of the twentieth size class appears to be best. Figure 6 gives the nucleation rate for the nineteenth, twentieth and twenty-first size class, which range from 500 to 600, 600-700 and
JOHAN JAGER et al
Fig. 7. Nucleation rate calculated with the method of characteristics from the twentieth size class (6OG7OU pm) for a step change in heat input of 1.9 kW to 5.3 kW at 16,200 s, a step change in heat input of 5.0 kW to 1.9 kW at 19,200 s in experiment 2, and a step change in heat input of 2.0 kW to 3.3 kW at 16,500 s in experiment 3, respectively.
Fig. 8. Steady-state population densities for a sample immediately before a step in heat input (34 CSD) and a sample at the end of experiment I (100 CSD).
700-800 km, respectively. It shows the reproducibility for the birth-rate of three succeeding classes. The results indicate that a steady state has not been attained by time t = 17,700 s. This conclusion is not supported by the experimental moments. The artefact probably results from the variations in the growth rate which cause the calculated nucleation rate to have a large value at the beginning of the experiment. It also illustrates how the technique is affected by the measurement noise. This effect does not appear in Fig. 7 which shows the results of the twentieth size class for the other two experiments. Both Fig. 6 and Fig. 7 indicate the agreement between the birth-rates which are calculated using the In (n) vs L technique and those using the method of characteristics. They both agree with the expected tendencies. The experiments start with a high nucleation rate which levels off, and is then succeeded by a sharp reaction on the heat input. However, it does not seem to be justified to use these values for quantitative correlations. The decrease in the nucleation rate from the sixth to the twentieth size class in Fig. 5 suggests the occurrence of growth dispersion within the small crystal sizes. This is evident from the curvature of the steadystate CSDs for small sizes as shown in Fig. 8. Non-linear parameter estimation The third moment is determined by the process conditions only and does not, as indicated by eq. (35), depend on the kinetic parameters. The second moment was chosen as the additional moment to be optimized since it is measured accurately by the Malvern particle sizer. The third moment is used in the calculation of the nucleation rate, as indicated in eq. (38). The results for optimization of the parameter k are shown in Fig. 9, in which a good fit to the second moment is seen. From these results the nucleation
Fig. 9. Experimental vs simulated second moment for an optimized parameter combination using the third moment in the nucleation rate, in which the transient of the second moment is optimized. The results refer to a change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
model can be fitted by the following
expression:
B = 7.3 x 10” G2.5m;.8.
(38)
If the power of the growth rate (i) is estimated the parameter values are identical. A similar fit is obtained using the second moment as the correct moment in the nucleation kinetics. This shows that it is virtually impossible to discriminate between the second and third moment from these CSD transients. This fact brings into question the importance of which moment should be used in the cakulation of the nucleation rate. Simulated results for three other parameter combinations, which are calculated according to eqs (29) and (30) indicate that, although the nucleation rate change significantly, the influence of the nucleation parameters on the second moment is limited.
Estimation of nucleation
L "
kinetics
from crystal
550 f -
450 t.t.p
400 0
20
815
I 25
30 -
Fig. 10. Simulated vs experimental
hme
35 (set)
response
40 ~103
for the mass-
based average size (m4/m3), in which an optimized parameter combination derived from the second moment is used in the simulation.
I
transients
From the steady-state population densities shown in Fig. 8 the MSMPR deviations are interpreted in the small size range probably as growth dispersion and in the larger size range probably as crystal attrition. Figure 8 also illustrates an enhanced curvature in the population density plot prior to the step. This suggests an increase in attrition rates at higher slurry densities, which is also apparent from the sieve analysis. Unless both these phenomena are added, the simulated model will fail to predict the responses for more than two moments correctly. Conversely, adding these phenomena will change the parameter values found in this study.
600
500
size distribution
I ts3.p
Fig. 11.
Simulated vs experimental response for the zerceth in which an optimized parameter combination derived from the second moment is used in the simulation.
moment,
Subsequently the simulated result for the massbased average size (m4/m3) in the first experiment is shown in Fig. 10. This is indicative of the behaviour of the fourth moment since the third moment is independent of the nucleation parameters. It shows a displacement of the steady-state value, and, consequently, a displacement in the unsteady-state region. When the average size itself was chosen to be optimized, using the steady-state nucleation rates calculated for j = 4 in eq. (26), a good agreement was found, albeit at the cost of a displacement of the value of the second moment. The zeroeth moment is even further displaced as indicated in Fig. 11. This clearly indicates that the MSMPR assumptions do not hold, otherwise at least the steady-state moments would be predicted correctly for a combination of two moments. This appears from eq. (4), in which for a given residence time only one degree of freedom is left. Consequently, all the moments can be calculated from any combination of two known moments.
DISCUSSION
Four techniques have been tested. It appears that two of the techniques are not feasible in our application. Their dependence upon the smaller size range, where the measurement noise is strongest, reduces their usefulness. They both suffer from the sensitivity of numerical differentiation to the measurement noise. This is a serious disadvantage of these techniques. A more positive result was expected from the Laplace technique. However, the weighting accomplished is unsuccessful for a volume-based sizing technique. Moreover, a total nucleation rate is calculated rather than an effective one. Therefore the calculated nucleation rate fails to predict the behaviour in the larger crystal size range, where only the effective nuclei will eventually appear. The method of characteristics results in qualitatively correct results, although the outcome does depend on the size class chosen to recalculate the nucleation rate. From the larger size ranges an effective nucleation rate is calculated whereas from the smaller size range a total nucleation will result. The results for the nucleation rates given in Figs 6 and 7 show a strong relationship with the heat input, or consequently the growth rate, the third moment or maybe the second moment. This suggests that the nucleation kinetics can be represented by a secondary nucleation model. The necessity of assuming a nuclei stock does not appear from these results for ammonium sulfate. Although the method predicts the expected tendencies it does not seem justifiable to use the results for parameter estimation since the scatter is pronounced. Therefore a fourth technique which estimates the kinetic parameters directly, rather than through the calculation of growth and nucleation rates, has been developed. This technique is applied successfully to estimating kinetic parameters. From these results the influence of the kinetic parameters appears to be limited. This can be mathematically interpreted using the moment equations. The nucleation rate enters into the zeroeth-moment equation. Only after three integrations will a change in nucleation rate affect the third moment. This indicates that changes in the nucleation rate are dampened before they end up in the volume distribution. The necessity to add other population events such as
JOHAN JAGER et al.
816
crystal attrition and growth dispersion appears. The fourth technique can be extended to the estimation of additional parameters involved in these phenomena, if they can be modelled adequately. It can easily be applied to Class-I systems too, by adding the appropriate growth kinetics, which probably contain unknown parameters also. This involves the estimation of more parameters and requires, consequently, a more sophisticated optimization algorithm. This analysis will be presented in a subsequent paper.
Acknowledgements-The authors wish to thank the Netherlands Technical Foundation (STW), AKZO, DOW, DSM, Dupont, Rhone Poulenc, Suiker Unie and Unilever for their financial support of the UNIAK programme. Ir. M. G. van Kooten, Ir. L. H. H. Hiemstra and Ir. A. N. Vlug contributed to this paper during their final year MSc study on this topic. NOTATION
a B
c
cs
CONCLUSIONS
Four techniques have been tested for their ability to calculate growth and nucleation rates under unsteady-state conditions. -The moments technique and the Laplace transformation appear not to be feasible while they use numerical differentiation which is very sensitive to measurement noise. The moments technique relies upon measurement of the zeroeth moment which is difficult to assess with the present sizing technique. This also handicaps the Laplace technique as datapoints at decreasing sizes are increasingly important. This does not comply with the accuracy of the measuring technique. -The method of characteristics produces sensible results. The outcome depends on the size class from which the nucleation rate is recalculated as a result of the curved population density plots. A strong qualitative correlation with the heat input appears. Quantitative correlation appears to be unrealistic as a result of the scatter. Therefore a fourth technique has been developed. -The fourth technique appears to be extremely effective in calculating the nucleation parameters required to predict the response of the second moment or the average size separately, additional to the response of the third moment which does not depend on the nucleation kinetics. It appears that both the average size and the second moment can not be reliably predicted so long as crystal attrition and or growth dispersion are excluded from the model. The technique is flexible enough to incorporate these phenomena, if the appropriate models are available. With respect to the ammonium study it appears that:
sulfate system
under
CP G i, i’, k
k, k,, kb
k, k k,, 4, k, L MS1
mj N
n0 4-c t)
d”,(X) n’(s, t) P
Q
Rv S
T V
WI x X Y
Y Greek letters d. W,, ALx & P
-Secondary nucleation kinetics prove to be effective in describing nucleation kinetics, Surprisingly, the system itself reacts only slightly to changes in the kinetic parameters. Whether the second or the third moment is used proves to have little influence on the final results. The necessity to add a nuclei stock is not apparent. The influence of other population events is not negligible.
empirical constant nucleation rate, No./m3 s concentration, kg/m3 clear liquor concentration at saturation, kg/m3 clear liquor specific heat, kJ/kg K growth rate m/s empirical constants used to denote the moment j shape factor ( = 3.14) nucleation constant, mass transfer coefficient, m/s shape factor ( = 0.26) constants crystal size, m slurry density, kg/m3 slurry jth moment of the distribution stirrer speed, rps population density at zero size, No./m4 population density at size L, No./m4 cumulative contribution of the xth size class to the Laplace transformed population density Laplace transformed population density, No./m3 heat input, W flowrate, m3/s heat of evaporation, kJ/kg supersaturation [(C - C,)/C,] Laplace variable, l/m time, s temperature, K crystallizer volume, m3 weighing factor index x-co-ordinate index function value
r
t)
volume fraction in the xth size class width of a size interval x, m void fraction density, kg/m3 residence time, s
Indices
5 P V no index
denotes denotes denotes denotes denotes denotes
crystal feed liquid product vapour crystallizer
contents
Estimation
of nucleation
kinetics
from crystal -no -no
REFERENCES Berglund, K. A. and Larson, M. A., 1984, Modeling of growth rate dispersion of citric acid monohydrate in continuous crystallizers. A.I.Ch.E. J. 30. 280. Boxman, A. and Scarlett, B., 1988, On-line measurement of crystal size and shape using combined optical techniques, in SPIE Vol. 952, Proceedings of "Laser Technologies in Industry” (Edited by 0. D. D. Soares), p. 520. Cooper, D. J. and Clough, C. E., 1985, Real-time estimation of dynamic particle-size distributions in a fluidized bed: theoretical foundation. 4.I.Ch.E. J. 31, 1202. Daudey, P. J., 1987, Crystallization of ammonium sulfate, secondary nucleation and growth kinetics in suspension. PhD thesis, Delft University of Technology. Daudey, P. J. and Jong, E. J. de, 1984, The dynamic behaviour of NaCl in a 91 liter MSMPR-crystallizer, in Industrial Crystallization 84 (Edited by E. J. de Jong and S. J. Jancic), pp. 447-450. Elsevier, Amsterdam. Gill, P. E., Murray, W. and Wright, M. H., 1981, Practical Optimization. Academic Press, London. Heffels, S. H., 1986, Product size distributions in continous PhD thesis, Delft and batch sucrose crystallizers. University of Technology. lager, J., Wolf, S. de, Klapwijk, W. K. and Jong, E. J. de, 1989, A new design for on-line CSD-measurements in product slurries. Presented at the 10th Symposium on Industrial Crystallization, Bechyne Castle. Jager, J., Wolf, S. de, Kramer, H. J. M. and Jong, E. J. de, 1990, On-line particle size measurement in dense slurries. in press. Jancic, S. J. and Grootscholten, P. A. M., 1984, Industrial Crystallizarion. Delft University Press, Dordrecht. lanse, A. H. and Jong, E. J. de, 1976, The occurrence of growth dispersion and its consequences, in Industrial Crystallization (Edited by J. W. Mullin). Plenum Press, New York. Marchal, P., David, R., Klein, J. P. and Villermaux, J., 1988, Crystallization and precipitation engineering-I. An efficient method for solving population balance in crystallization with agglomeration, Chem. Engng Sci. 43, 59. McCabe, W. L., 1929, Crystal growth in aqueous solutions. lnd. Engng Chem. 21, 30. Randolph, A. D. and Larson, M. A., 1988, Theory of Particulate Processes. Academic Press, London. Shah, M. B. and Garside, J., 1980, Crystallization kinetics from MSMPR crystallizers. Ind. Engng Chem. Process Des. Deu. 10, 509. Tavare, N. S., 1986, Crystallization kinetics from transients of an MSMPR-crystallizer. Can. J. them. Engng 64, 752. Tavare, N. S. and Garside, J., 1982, Simultaneous estimation of crystal nucleation and growth kinetics from batch experiments. Chem. Engng Res. Des. 64, 109. Tsuruoka, C. and Randolph, A. D., 1987, State space rep resentation of the dynamic crystallizer population balance: application to CSD controller design A.J.Ch.E. Symp. Ser. No. 253, 83, 104. Wolf, S. de, Jager, J., Kramer, H. I. M., Eek, R. and Bosgra, 0. H., 1989, Derivation of state-space models of continuous crystallizers. Preprints of IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, 21-23 August 1989, Maastricht.
APPENDIX
A: DERIVATION AND GROWTH
OF PRODUCT
FLOWRATE
--constant volume, constant temperature, -neglectable heat of crystallization, --crystal-free feed flow. -salt-free condensate.
transients
817
growth dispersion or size-dependent growth agglomeration, crystal breakage or attrition,
rate,
the salt balance, the total mass balance and the heat balance can be simplified into a constraint for the growth rate and a constraint for the product flowrate. In this derivation the time variable will be omitted in the notation. The salt balance is given by $ t
CCs+ p,(l - 41 VI = C,Q, - CCE+ PAI - 41 Q, (Al)
where the void fraction
is defined
by
E = 1 - k,m,
(A21
and thus ds = Pk”d!$
(A3)
dt The moment
equation
for i = 3 follows from eq. (12):
dm3 dt
_=
--+33Gm, r
Combining eqs (Al)-(A4) and using the Class-II dC/dt = 0 results in a growth rate constraint:
G=
is given by
+ AU - 41Qp - P,Q,, W4
p,Q, - Cw Inserting
assumption
Q&,--Q& % Vm,( pc - Cl
The total mass balance
eqs (A2)-(A4)
in eq. (A6) gives VGm,.
The heat balance
(A7)
is given by
- Ccc,p,~+ P,c,,(~ - 41 TQ, - p,Q,(R, + cpu77 + P. WI The heat balance can be simplified into a constraint vapour flow by inserting eqs (A2)-(A4) in eq. (A8):
*(P,c, - ~01.
+ P - 3k,% Combining lC
+
RATE CONSTRAINT
In this Appendix it will be shown that, evaporative MSMPR crystallizer with:
size distribution
this constraint
C(P, - Pr)
7% R, + c,,T
- P&J
+ cpv7-HPC - Cl 1
ML
(A9)
with eqs (A5) and (A7) results
P,(P. - C)
*UP&
for the
in
CP,(R, + cpu *)I Q,
for a Class-II
*) - TJC,,,P,l -
-
*cP&
- P&l)
(P, - C)
1C, >
Q, - P
s(R. r
+ e,,T)
(AI’3
818
JOHAN JAGER et al.
Or
Qp=
k,
k, = P,& + cpv~) - T,c,,/P, k _; 2
1 "c'(R,+ QT) c
k,=
(k, - k,C,)Q, - P
C( PC - PI) -~
Tc,,
(All) (A121
+ ~Chc, - w,,) b,(R, + cpoZ-11. PdR” + C&s”WP, - q 1
_ T(pcc, -p&!) (PC--_
I-
[
6413)
PdP, - Cl
R,+c,,T 6414)
Engineerme
Printed in Great
Snencc
Vol.
46. No.
3. pp. 807
818. 1991.
KG?-2509191 %3X0 + 0.00 8 1991 Pergamon Press pk
Britam.
ESTIMATION OF NUCLEATION KINETICS FROM CRYSTAL SIZE DISTRIBUTION TRANSIENTS OF A CONTINUOUS CRYSTALLIZER JOHAN Laboratory
JAGER,
SJOERD
for Process Equipment,
DE
WOLF,’
HERMAN
J. M. KRAMER*
and ESSO J. DE JONG Delft University of Technology, Netherlands
(First received 30 August
1989; accepted
Leeghwaterstraat
in revised form 8 January
44, 2628 CA Delft,
1990)
Abstract-Four techniques which can be used for the estimation of nucleation kinetics under unsteadystate conditions are applied to transients of a 20-l evaporative ammonium sulfate crystallizer. A moment and a Laplace transformation technique appear not to be feasible because they rely heavily on the measurement of the number of small crystals. The present crystal size distribution (CSD) measurement technique is inaccurate in this size range and therefore both methods suffer severely from measurement noise. The third technique, using the method of characteristics, indicates the appearance of so-called secondary nucleation kinetics. The technique is unsuitable for estimating the kinetic parameters quantitatively, again as a result of the measurement noise. Therefore a fourth technique, which uses a non-linear parameter estimation routine, has been developed. This technique is very successful in estimating the parameters in a secondary nucleation model. With the estimated parameters the second and third moments of the CSD can be predicted. If additional moments are to be predicted it is shown that other phenomena which affect the crystal size distribution, such as crystal attrition, must be included in the dynamic model of the process.
to describe the nucleation rate (B) by an empirical power law. This power law comprises the supersaturation (S), which is considered to be the main driving force for nucleation, the crystal surface (k,m,), crystal mass (k,p,m,) or a higher moment mj, which reflects the idea that the presence of crystals plays a dominant role (secondary nucleation), and the stirrer speed (N) which governs fracture processes (Daudey, 1987):
INTRODUCTION
unit operation in chemical engineering practice. It is employed in the production of bulk chemicals such as ammonium sulfate, sodium chloride and sucrose, and the production of fine chemicals such as aspartame. Usually the process produces a solid product which has a wide crystal size distribution (CSD). This CSD dictates the behaviour of the product in succeeding operations, such as filtration, drying, transport and storage, and is also used to define customer specifications. It is, therefore, desirable to control the CSD which is produced by the process. The first requirement is an on-line CSD measurement technique. Recent developments have made this possible (Jager et al., 1989). Another requirement for process control is an accurate dynamic model. This mode1 can be obtained by estimating the parameters of “a black box” model or, alternatively, by determining a model based on physical laws. This paper proceeds with the second approach. Crystallization can be described in terms of dynamic heat and mass balances, combined with a dynamic population balance which describes the CSD. The population balance contains specific population events such as crystal growth and crystal birth. The latter is especially a complicated process involving fracture mechanics, two-phase flow hydrodynamics and solid state physics. Hence, the usual approach is
Crystallization
from
solution
is an
important
B = kbS”(?nj~N”.
(1)
If the supersaturation is difficult to measure, it may be exchanged by the growth rate G, which is a function of the supersaturation: G = k,S’
(2)
B = k,Gi(mj)l’N”.
(3)
The usual approach is to perform steady-state experiments under various operating conditions. From the resulting CSDs, the growth and nucleation rates can be estimated, often by assuming the distribution to be log-normal (Randolph and Larson, 1988). By curve fitting, the empirical constants can be found for different moments (Shah and Garside, 1980). There are two objectives in extending this analysis to unsteady-state conditions: (1) Under steady-state conditions the moments for an idealized MSMPR (mixed suspension mixed product research) crystallizer, defined in the Appendix, at a residence time z, correlate according to
‘Laboratory for Measurement and Control, Delft University of Technology, Mekelweg 2, 2628 CD Delft, Netherlands. ‘Author to whom correspondence should be addressed.
mi = jmj807
,rG.
(4)
808
(2)
JOHAN
JAGER
Therefore, it is not possible to determine the constants in the nucleation model from steadystate data. The verification of secondary nucleation kinetics under transient conditions. Recently, experimental evidence has been presented which indicates that the steady-state constants do not describe correctly the nucleation process for the unsteady state (Heffels, 1986). The occurrence of a so-called nuclei stock, which contains slow or non-growing crystals, and which can be activated at momentarily higher supersaturations, was postulated (Heffels, 1986).
However, the growth and nucleation rates, and consequently the parameters in eq. (3), cannot easily be evaluated under unsteady-state conditions. Four techniques which could serve this purpose are explored in this paper, using experimental transients produced in a 20-l continuous evaporative crystallizer. A moment and a Laplace technique were applied to CSD transients from batch and continuous crystallizers by Tavare and Garside (1982). However, the experimental scatter in their results was pronounced. The method of characteristics was first reported by Daudey and Jong (1984), but with loo few results to judge its merits. The fourth technique, which uses a general non-linear parameter estimation technique, has been developed in order to estimate the parameters of the nucleation kinetics. Each of these techniques requires a dynamic model which wil1 be discussed first. THEORY
Dynamic modeling The crystallizer contents of an MSMPR crystallizer are fully defined by the heat balance, the total mass balance, the solute balance and the population balance (Wolf et al., 1989). The differential equations which represent the dynamic heat, mass and solute balances can be simplified by applying two constraints for a constant-volume, isothermal and ClassII crystallizer with size-independent crystal growth (Appendix). The Class-II assumption implies that the ammonium sulfate system is a fast growing system in which the crystallizer concentration is constant with time and approaches the equilibrium concentration. The assumption of size-independent crystal growth can be implied, for the intermediate size range, from a log-normal population plot. The heat and mass balance simplify into a relation for the product flowrate (Q,), as a function of the feed flow (Q,), the feed concentration (C,), and the heat input (P): Qp(tl
=
Cki - kzCj(t)lQf(t) - p(t)
et al.
G(t)
(5)
The solute balance reduces with the Class-II assumption to a relation for the growth rate(G), as a function of the crystallizer volume (V), the crystallizer concentration (C), and the second moment (m2):
Q,(t)C
-
(6)
- C) ’
3k, W(t)(p,
At a constant feed flow, eq. (5) predicts that the product flow will decrease upon an increase in the heat input. Consequently, the growth rate will immediately increase upon an increase in the heat input according to eq. (6), and then decrease, because of the subsequent increase in the second moment. The dynamic population balance describes the CSD transients. If crystal breakage, attrition and agglomeration are neglected, crystals are formed at very small size, and all crystals have the same growth rate (McCabe, 1929), the well-known population balance (Randolph and Larson, 1988) for the population density, n, at size L results:
anw,t) + n(L, t) Q,(t) = 0. ad5 t) ~ + G(t)~ V at 3L The boundary differential
condition
equation
for this hyperbolic
(7)
partial
is given by no(t)
=
B(f) G(f)
which changes,
by substitution
n,(r) The jth moment
of eq. (3), into
= k,G’-‘(t)mS(t)N”. of the distribution
(9) is defined
by
cc n(L, 2)LjdL.
mi(t) =
(10)
I0 Techniques Three of the techniques which are discussed, denoted the method of moments (Tavare, 1986), the Laplace method (Tavare, 1986; Tavare and Garside, 1982), and the method of characteristics (Daudey and Jong, 1984), calculate the growth and nucleation rates appertaining during an experiment. The parameters in eq. (3) can then be estimated using a linear least-squares routine. The fourth method, which combines two mathematical techniques, the method of lines and non-linear parameter estimation, estimates the unknown parameters without the intermediate step of calculating the growth and nucleation rate. The fourth method will be denoted non-linear parameter estimation. The method ofmoments. This method uses the moment transformation to transform the population balance into a set of ordinary differential equations describing the moments of the distribution:
mo(t)Q,(t)
dmdt) =dt
K3
QJt)c/(t)
=
dmj(t) -= dt
V
+ B(t)
- mj(t)Qp(t) . +Jmj_,(t)G(t) V
(11) for j z 1.
(12) By replacing proximations,
the derivatives with finite-difference apthe nucleation rate can be calculated
Estimation of nucleation kinetics from crystal six distribution from the transient zeroeth growth rate can be calculated and higher moments: Am,(t
moment whereas from the transient
+ fAt, + mo(t + fAt)~,(t) At
V
the first
= B(T + +At) (13)
Amj(t + &At) + mj(t + )A~)Q~(c) At V = G(t + $At) jfi,_ ,(r + +At) (14) where Amj(t + + At) = mj(r) - mj(t + At)
m,(t +
fAt,
=
CmAt)+ mAt + “)I
J
2
In this paper the growth rate has been calculated from eq. (6) rather than from eq. (14) in order to avoid numerical differentiation, which is very sensitive to noise on the experimental data. Unfortunately this cannot be avoided for the nucleation rate. The calculated nucleation rate will depend on the measured zeroth moment which is dominated by the small crystals. An average rather than an effective nucleation rate, therefore, results. The Laplace method. The result of the Laplace transformation of the population balance is a firstorder ordinary differential equation in the transformed population density, n’(s, t), using the Laplace variable s: n’(s, t) =
CI) n(L,t)exp(--sL)dL. s0
dn’(s, t) + n’(s, t)Q,(t) dt
V
(15)
= - G(t)sn’(s,
t) + B(t).
(16)
The transformed population density can be calculated from the measured population density. The derivative in eq. (16) can be calculated using a finite-difference approximation. Over a range of s values a linear relation results, with B(t) and G(t) as unknown parameters The slope of the straight line Y(t) = B(t) - G(t)X(t)
transients
for CSD measurement techniques which are most accurate in the intermediate size ranges. If negative s values are used, the larger size ranges are emphasized. At first sight this might prove to be favourable: however, it also puts a very strong emphasis on the “tail” of the distribution. The weight function [exp ( - sL)] approaches infinity as L tends to infinity, where the measurement accuracy is nil. The s values to be taken cannot be estimated in advance but must be evaluated by experiment, which introduces uncertainties. This second method also calculates an average nucleation rate. Method of characteristics. The method of characteristics seeks special curves in the L-t plane, called characteristic curves, along which the solution of the partial differential equation is reduced to the integration of an ordinary differential equation. Within the present context, the partial differential equation is transformed into an ordinary differential equation along curves where eq. (18), for the crystal size of a specific class of crystals (L*) holds: (18) By eq. (18) a collection of curves, denoted characteristics, is defined. Along these characteristics the population density in time for this specific class of crystals (n*) is given by (19) These two equations are easily solved. The growth rate G can be calculated from eq. (6). Using this growth rate, the birth time (tbirlh) of a particular size class attained at a later time (t,,,,l.) can be estimated using flami.lc L*( Llnple 1 -
G(t)dt
= 0.
s fbirth
(20)
The population density at the time of birth can be calculated, using the product flowrate from eq. (5), by integrating eq. (19):
(17)
nO*(tbirth)
=
Lmpk - tmJQp(tl n*(Lampk)exp V
where
(21) dn’(s, r) + n’(s, t)QJt) y(t) =dt x(t)
V-~ _______
= sn’(s, t)
represents the growth rate while the abscissa will give the nucleation rate. Closer examination of the transformation reveals that n’(s, t) results from weighting of the measured population densities by [exp ( - sL)], and integrating over the size range L. This gives the opportunity to weight the data and reduce the measurement noise. For positive s values this results in a weight function which decreases with size, which is unfavourable
and
B*(tbirth)= 4 (tbirth)G(tbirth)+
(22)
For a measuring device with 32 size classes, the nucleation rate can be calculated for each of these classes. The birth-rates for the different size classes do not necessarily coincide and interpolated values therefore must be compared. These are identical if McCabe’s AL law holds. However, if growth dispersion occurs the nucleation rate will differ for different size classes. Growth dispersion states that crystals have different growth rates when exposed to the same supersaturation (Janse and Jong, 1976; Berglund
810
JOHANJAGER~~~I.
and Larson, 1984). If this phenomenon is restricted to the smaller size range, the nucleation rates calculated from size classes exceeding this range can be interpreted as an effective nucleation rate. The effective nucleation rate accounts for the number of nuclei, which grow to a for the sizing technique detectable size only (Jancic and Grootscholten, 1984). At larger size ranges the nucleation rate will be affected increasingly by measurement noise. The difference. between the sample time and the time of birth which appears in the exponential of eq. (21) is larger for larger sizes. This tends to promote an error, due to the time of birth term in the calculated nucleation rate. For small si7e ranges the measured population density is itself of limited accuracy. The intermediate size ranges are, therefore, to be preferred.
Non-linear parameter estimation. Rather than calculate the growth and the nucleation rates separately and then fit them to a model, the first stage can be avoided by estimating the unknown parameters in eq. (9) directly. A linear, least-squares routine cannot be used for this purpose, since the parameters occur in the boundary condition of a hyperbolic partial differential equation. An alternative technique must there-
fore be used. Using starting values for the unknown nucleation parameters one can simulate the overall model consisting of the population balance [eq. (7)], the nucleation model [eq. (9)], the growth rate constraint [eq. (6)] and the constraint for the product flowrate [eq. (511. The kinetic parameters can be used to adjust the simulated to the experimental responses. The adjustment is made by a non-linear optimization technique. The procedure is terminated when a chosen criterion is met. Three aspects of this technique will be discussed in detail. Firstly, the technique which is used to simulate the responses will be discussed. Subsequently, the number of parameters to be estimated is, by using steady-state information, reduced to one. This is a convenient way to adjust the solution to the initial and final values of the step responses which were investigated. The use of steady-state information does not compromise the first objective of this paper, because eq. (4) is certainly violated immediately after a step in heat input [see eqs (6) and (12)]. The criterion which was used in the optimization will finally be discussed. The partial differential equation [eq. (7)] cannot be solved numerically using standard techniques and it must therefore be transformed. It was decided to choose a technique which is not restricted to simple crystallizer configurations, and the moment transformation was therefore not considered (Marchal et a(., 1988). An alternative technique to transform eq. (7) into a form which lends itself to simulation is the method of lines (Cooper and Clough, 1985; Tsuruoka and Randolph, 1987; Wolf er al., 1989). In this procedure the size axis is discretized and the length
derivative is substituted by a finite-difference approximation. Eventually, a number of first-order differential equations for n(L) in discrete gridpoints replaces the original partial differential equation. Several finite-difference approximations can be chosen, depending on the required accuracy. Throughout this paper the size axis is discretized into x + 1 points and a fourth-order central difference, based on an equidistant grid, is used to approximate the derivative with respect to crystal size:
WL,
t)
=~r(n(L,-,,r)-sn(L,~,,t)
dL
+ WL+,,
t) - fi(L,+,,
t)i + o(~L4).
(23)
This central difference cannot he used at the boundaries of the grid L,, L,, Lx_ 1 and L,, since n(L_,, I), t) and n(L,+*, t) miss. Therefore n(L-I, t). n(L,*,, eq. (23) is replaced by two forward differences on the left boundary [eqs (24a) and (24b)] and two backward differences on the right boundary [eqs (25a) and (25b)] (Wolf et al., 1989):
hdt)
&
-z
a_c
[-25no(t)
- 36n(L,,
+ 48n(L,,
t) + 16n(L,,
t)
t) - h(L,,
t)]
+ 0(AL4)
anw,. -=
(24a)
t) &
aI5
[-3n,(tI
+ 18n(L,,
-
lOn(L,,
t)
t)
t) - 6n(L,,
+ n(.J+, t)] + O(AL4) WL,-
1,
JL
t)
= &
c -n(L,-.,
- 18n(L,_,,
Wb)
t) + 6n(I+
J, t)
t) + IOn(L,_ I, I)
+ 3n(L,, t)] + O(AL“)
ML, -= aL
t)
&
C+3$&.,
2) - 16n(L,-,,
+ 36n(L,_ 2, t) - 48n(L,_ + 25n(L,,
(25a)
t)] + O(AL’).
t)
I, t) (25b)
Replacing eq. (24a) by the boundary condition for no, and substituting eqs (24b), (25a) and (25h) into eq. (7) then results in a set of ordinary differential equations at discrete gridpoints L,. In the overall model, the population density must be integrated numerically in the calculation of m2 in eq. (6) and m, in eq. (9). The size axis was divided into 8, 31, 101 and 201 gridpoints, lying between zero size and a maximum of 2 mm, respectively. The changes in the simulated responses where negligible upon an increase from 31 to 101 gridpoints. Thus, the condition of grid independence is satisfied for both the responses and the estimated parameters at the current number of 201 gridpoints. The set of differential equations is solved using a Runge-Kutta routine.
Estimation of nucleation kinetics from crystal size distribution At a constant stirrer speed, the term N” in eq. (9) can be included in the nucleation constant k,. Then, for likely values of j equals 2 or 3, three parameters (k,, i, k) must be estimated in the boundary condition [eq. (9)]. At steady state, n, can be calculated from the jth and the (j - l)th moment: .j+lmj’:
n, =
_I
/+
If two steady-state conditions (1) and (2) are known correctly from experimental data, two degrees of freedom are removed by two constraints:
logn,,
I = logk,
+ (i -
1)log G1 + klog(m,),
(27)
log n,, 2 = log k, + (i - 1) log G, + k log @I~)~.
(28)
If k is the parameter to be optimized, i and k, are determined by these constraints according to Klog (mJ)l - log (m&l i= log no, I - log no, 2 (log G, - log G,) + 1 (29) logk,
= log no*, -(i-
l)logG,
-klog(m,),.
(30)
The optimization of k can then be achieved with any non-linear line search technique. A technique based upon the Fibonacci series (Gill et al., 1981) was chosen because of its ease of implementation. Either one moment or a combination of moments can be optimized. Typically, simulated datapoints are subtracted from the experimental ones, squared and summed. This sum is minimized. Weighting factors can then be used to emphasize the moments which are considered important: Sum = i
i:
x=l j=I
wjCmj..h(~,) - mj,ex*(tx)12 (31)
where n = number
of datapoints, z = number of moments, wj = weighting factor for jth moment. EXPERIMENTAL
STRATEGY
Throughout this investigation a 20-l continuous evaporative crystallizer which is operated at about 100 mbar and a stirrer speed of 700 rpm was used. A detailed description of the experimental set-up is given by Daudey (1987). A micro-computer which controls the temperature, input flow, and heat load was added to the original configuration in order to reduce undesired disturbances caused by random variations in the process conditions. This crystallizer
transients
was operated at a temperature of 50°C and manually sampled every 5 min. The resulting samples were subsequently filtered, washed and dried. The total solids content was evaluated using the sample volume and the amount of solids after the filtration procedure. This stage proved to be very important for the crystal size analysis. Opaque samples give different results compared to completely transparent ones. It is expected that the main cause of the measurement noise is due to this phenomenon. The samples were analyzed 3 times using a Malvern 26OOc laser diffraction particle sizer and the results were averaged. Differences between these three measurements were small compared to differences between samples. In these investigations, the Malvern analyzer was fitted with a lOOO-mm lens which is appropriate for a measurement range of 0.519OOpm, which is divided into 32 size classes (Jager et al., 1990). In order to match the ammonium sulfate system, the software of the instrument was adapted (Boxman and Scarlett, 1988). The crystals were analyzed in the dry state using a dry-powder feeder. In order to avoid classification during this procedure, care was taken to measure the complete sample. An on-line CSD measurement technique was available on a 970-l pilot plant crystallizer. It could not be applied on the 20-I crystallizer because the volume contained by the dilution unit is also approximately 20 1. A typical experimental procedure was as follows. First, a start-up of the crystallizer at the required process conditions. After about eight residence times one or two stepwise changes in heat input or input flowrate were made, and then the experiment was continued for about eight further residence times. Over 15 such experiments were conducted. In this paper three of these experiments will be discussed, all of which consisted of changes in heat input. The experimental conditions are listed in Table 1. The changes in heat input effect changes in the vapour flow which, in turn, result in variations in the residence time, since the input flow was kept constant. The feed of the crystallizer is saturated at 46°C. The population densities and the moments of the distributions were calculated from the measured volume fractions (AV), the experimental magma density (MS!), the average size in each interval (L) and the interval widths (AL) given elsewhere (Jager et al., 1990): A V(L,, t) Msl( t) n(L,
t) =
p,k,L:AL Msl(t)
mj(t) = ~
32 AV(L,, c
X=l
p
Table I, Experimental conditions for three different experiments Experiment steps at (s) Heat input (kW) Residence time (s)
1 17,700 3.4-2.0 2450-2270
811
2 16,200-19,200 1.9-5.3-l .9 24OG293~24Cil
3 16.500 2.C3.3 2270-2430
L,3
t)L; (33)
812
JOHAN
JAGER
et
al.
RESULTS
1
The results of the first experiment were typical for most of the experiments and, unless otherwise stated, have been used to produce the graphs in this paper. From eqs (32) and (33) it is clear that the total amount of solids (Msl) is used to calculate the population densities and the moments. Therefore the measurement error in Msl contributes to that in the population densities and the moments. The amount of solids, Msl, was determined experimentally. However, it can also be calculated. Substituting eq. (6) into eq. (12) for j = 3 yields an ordinary differential equation in the third moment:
m,(t) + CfQJt)
dm,(t) -=-dt The solution is given by
- CQJt)
k,(p, - C)
r(t)
v
.
The total amount
CQ, I - C,Q,s t m%--C)V .
71-1-
CQ,., Up,
15
20
2s
-
(34) change
in heat
input
30
3s
[setI ~101
tlme
LO
vs predicted slurry density
for a step of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
hhl
-?)I
x [I--exp(
(9)~ = 72
10
5
Fig. 1. Experimental
+ C(m& -
(m3)1=
kt.,
0 0
of eq. (34) for a step in heat input (1 + 2)
m,(t - tstep) = (m&
2o t
- C_rQ,,2 - C) v
(35a)
800
25
T
'tstcp
(35’4 (35c)
of solids is given by
Msl = pek,ms.
(36)
Both the experimental and the theoretical values are shown in Fig. 1 and show a very good agreement. In order to reduce the measurement noise in the other moments and the population densities, they have been calculated from eqs (32) and (33) using the predicted slurry density (Msl) [eqs (35) and (36)] rather than the experimental values. The growth rate, combined with the second moment from which it was calculated according to eq. (6), is depicted in Fig. 2. A very high growth rate existed during the beginning which decreased gradually towards the first steady state. On the step in heat input, the growth rate decreased sharply, after which a new steady-state value was eventually attained. Subsequently, the nucleation rate was calculated using the different calculation methods.
Method of moments Figure 3 shows the results using the method of moments. Hardly any trends can be discerned. The zeroth moment from which the nucleation rate is calculated is too distorted by the measurement noise. Even a fourth-order finite-difference approximation for the derivative did not improve the results. Data filtering, which was subsequently used to reduce the measurement noise, was not effective, without substantial loss of information, for the number and spacing of data points.
200 100
1
0 'I
0
, 0
s
10
15
tstrp 20
2s -
30
35
J.0
time [set I xl03
Fig. 2. Growth rate, and the second moment growth
L&ace
from which the rate is calculated, for a step change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
transformation
Figure 4 shows the Laplacian nucleation rates, which show the same tendency. They are calculated with 100 s values, lying between zero and a value of SL 50 = 1 as proposed by Tavare (1986). No improvement results if other values are used. This results from the weighting of the data. Table 2 shows the population density of a steady-state sample with no = 7.7 x lOI and Gr = 1.6 x 10m4, the corresponding weighting factors, and their cumulative contribution to the transformed population density for a value of SL 50 = 1. The cumulative contribution is given by
(37)
The measurement noise is greater for the smaller size ranges and the weighting by the Laplace method is,
Estimation of nucleation kinetics fmm crystal size distribution
transients
813
5 -20
,trtq 0
5
10
15
20 -
25 30 35 tint [sec.]x10)
Fig. 3. Nucleation rate calculated with the method of moments for a stcD change in heat inout of 3.4 kW to 2.0 kW at 17,700 s in expe%ment 1.
trt.,
Fig. 4. Nucleation for a step change
0
60
0
5
10
15
20
25
-
time
30 [sec.]
35
LO
~103
Fig. 5. Nucleation rate calculated with the method of characteristics from the sixth (62.5-75 gm), the twentieth (6C@-7OOpm) and the thirtieth size class (160&17OQ~m), respectively, for a step change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
I
20
25
-
tim
30 lsrc.1
35
1
-103
t*t.r
rate calculated with the Laplace method in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
5
10
15
20
25
30
35
Fig. 6. Nucleation rate calculated with the method characteristics from the nineteenth (500-600 bm), Table 2. Cumulative contribution of different size classes to the overall transformed population density for a theoretical crystal size distribution
n&J Class 1
6 11 16 21 26 31
(No./m’) 7.4 5.0 2.8 1.0 6.9 3.0 1.3
x x x x x x x
10’3 1OL3 lot’ IOL’ lOL1 lOL0 IO9
exp( - s&A 0.99 0.90 0.76 0.60 0.31 0.14 0.06
“bi”’ 0
9.0 44.4 74.7 93.7 99.8 loo.0 loo.0
therefore, ineffective for the measurement method employed throughout these investigations. The calculated nucleation rates vary strongly and even have negative values. Moreover, the time of the step cannot be recovered.
of the
twentieth (600-700 pm) and the twenty-first size class (700-800 pm), respectively, for a step change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
Method
c~~characteristics
Figure 5 shows the results for the sixth, twentieth and thirtieth size class, which range from 62.5 to 75, 600-700 and 16OC-1700 pm, respectively. The scatter for the sixth and thirtieth size classes is pronounced, which is in accordance with the prediction earlier made. Figure 5 also shows an additional feature of this technique. For times greater than 13,218 s, nucleation rates are not available for the thirtieth size class. These nucleation rates should have been caIculated from samples which were beyond the duration of the experiment. The reproducibility of the twentieth size class appears to be best. Figure 6 gives the nucleation rate for the nineteenth, twentieth and twenty-first size class, which range from 500 to 600, 600-700 and
JOHAN JAGER et al
Fig. 7. Nucleation rate calculated with the method of characteristics from the twentieth size class (6OG7OU pm) for a step change in heat input of 1.9 kW to 5.3 kW at 16,200 s, a step change in heat input of 5.0 kW to 1.9 kW at 19,200 s in experiment 2, and a step change in heat input of 2.0 kW to 3.3 kW at 16,500 s in experiment 3, respectively.
Fig. 8. Steady-state population densities for a sample immediately before a step in heat input (34 CSD) and a sample at the end of experiment I (100 CSD).
700-800 km, respectively. It shows the reproducibility for the birth-rate of three succeeding classes. The results indicate that a steady state has not been attained by time t = 17,700 s. This conclusion is not supported by the experimental moments. The artefact probably results from the variations in the growth rate which cause the calculated nucleation rate to have a large value at the beginning of the experiment. It also illustrates how the technique is affected by the measurement noise. This effect does not appear in Fig. 7 which shows the results of the twentieth size class for the other two experiments. Both Fig. 6 and Fig. 7 indicate the agreement between the birth-rates which are calculated using the In (n) vs L technique and those using the method of characteristics. They both agree with the expected tendencies. The experiments start with a high nucleation rate which levels off, and is then succeeded by a sharp reaction on the heat input. However, it does not seem to be justified to use these values for quantitative correlations. The decrease in the nucleation rate from the sixth to the twentieth size class in Fig. 5 suggests the occurrence of growth dispersion within the small crystal sizes. This is evident from the curvature of the steadystate CSDs for small sizes as shown in Fig. 8. Non-linear parameter estimation The third moment is determined by the process conditions only and does not, as indicated by eq. (35), depend on the kinetic parameters. The second moment was chosen as the additional moment to be optimized since it is measured accurately by the Malvern particle sizer. The third moment is used in the calculation of the nucleation rate, as indicated in eq. (38). The results for optimization of the parameter k are shown in Fig. 9, in which a good fit to the second moment is seen. From these results the nucleation
Fig. 9. Experimental vs simulated second moment for an optimized parameter combination using the third moment in the nucleation rate, in which the transient of the second moment is optimized. The results refer to a change in heat input of 3.4 kW to 2.0 kW at 17,700 s in experiment 1.
model can be fitted by the following
expression:
B = 7.3 x 10” G2.5m;.8.
(38)
If the power of the growth rate (i) is estimated the parameter values are identical. A similar fit is obtained using the second moment as the correct moment in the nucleation kinetics. This shows that it is virtually impossible to discriminate between the second and third moment from these CSD transients. This fact brings into question the importance of which moment should be used in the cakulation of the nucleation rate. Simulated results for three other parameter combinations, which are calculated according to eqs (29) and (30) indicate that, although the nucleation rate change significantly, the influence of the nucleation parameters on the second moment is limited.
Estimation of nucleation
L "
kinetics
from crystal
550 f -
450 t.t.p
400 0
20
815
I 25
30 -
Fig. 10. Simulated vs experimental
hme
35 (set)
response
40 ~103
for the mass-
based average size (m4/m3), in which an optimized parameter combination derived from the second moment is used in the simulation.
I
transients
From the steady-state population densities shown in Fig. 8 the MSMPR deviations are interpreted in the small size range probably as growth dispersion and in the larger size range probably as crystal attrition. Figure 8 also illustrates an enhanced curvature in the population density plot prior to the step. This suggests an increase in attrition rates at higher slurry densities, which is also apparent from the sieve analysis. Unless both these phenomena are added, the simulated model will fail to predict the responses for more than two moments correctly. Conversely, adding these phenomena will change the parameter values found in this study.
600
500
size distribution
I ts3.p
Fig. 11.
Simulated vs experimental response for the zerceth in which an optimized parameter combination derived from the second moment is used in the simulation.
moment,
Subsequently the simulated result for the massbased average size (m4/m3) in the first experiment is shown in Fig. 10. This is indicative of the behaviour of the fourth moment since the third moment is independent of the nucleation parameters. It shows a displacement of the steady-state value, and, consequently, a displacement in the unsteady-state region. When the average size itself was chosen to be optimized, using the steady-state nucleation rates calculated for j = 4 in eq. (26), a good agreement was found, albeit at the cost of a displacement of the value of the second moment. The zeroeth moment is even further displaced as indicated in Fig. 11. This clearly indicates that the MSMPR assumptions do not hold, otherwise at least the steady-state moments would be predicted correctly for a combination of two moments. This appears from eq. (4), in which for a given residence time only one degree of freedom is left. Consequently, all the moments can be calculated from any combination of two known moments.
DISCUSSION
Four techniques have been tested. It appears that two of the techniques are not feasible in our application. Their dependence upon the smaller size range, where the measurement noise is strongest, reduces their usefulness. They both suffer from the sensitivity of numerical differentiation to the measurement noise. This is a serious disadvantage of these techniques. A more positive result was expected from the Laplace technique. However, the weighting accomplished is unsuccessful for a volume-based sizing technique. Moreover, a total nucleation rate is calculated rather than an effective one. Therefore the calculated nucleation rate fails to predict the behaviour in the larger crystal size range, where only the effective nuclei will eventually appear. The method of characteristics results in qualitatively correct results, although the outcome does depend on the size class chosen to recalculate the nucleation rate. From the larger size ranges an effective nucleation rate is calculated whereas from the smaller size range a total nucleation will result. The results for the nucleation rates given in Figs 6 and 7 show a strong relationship with the heat input, or consequently the growth rate, the third moment or maybe the second moment. This suggests that the nucleation kinetics can be represented by a secondary nucleation model. The necessity of assuming a nuclei stock does not appear from these results for ammonium sulfate. Although the method predicts the expected tendencies it does not seem justifiable to use the results for parameter estimation since the scatter is pronounced. Therefore a fourth technique which estimates the kinetic parameters directly, rather than through the calculation of growth and nucleation rates, has been developed. This technique is applied successfully to estimating kinetic parameters. From these results the influence of the kinetic parameters appears to be limited. This can be mathematically interpreted using the moment equations. The nucleation rate enters into the zeroeth-moment equation. Only after three integrations will a change in nucleation rate affect the third moment. This indicates that changes in the nucleation rate are dampened before they end up in the volume distribution. The necessity to add other population events such as
JOHAN JAGER et al.
816
crystal attrition and growth dispersion appears. The fourth technique can be extended to the estimation of additional parameters involved in these phenomena, if they can be modelled adequately. It can easily be applied to Class-I systems too, by adding the appropriate growth kinetics, which probably contain unknown parameters also. This involves the estimation of more parameters and requires, consequently, a more sophisticated optimization algorithm. This analysis will be presented in a subsequent paper.
Acknowledgements-The authors wish to thank the Netherlands Technical Foundation (STW), AKZO, DOW, DSM, Dupont, Rhone Poulenc, Suiker Unie and Unilever for their financial support of the UNIAK programme. Ir. M. G. van Kooten, Ir. L. H. H. Hiemstra and Ir. A. N. Vlug contributed to this paper during their final year MSc study on this topic. NOTATION
a B
c
cs
CONCLUSIONS
Four techniques have been tested for their ability to calculate growth and nucleation rates under unsteady-state conditions. -The moments technique and the Laplace transformation appear not to be feasible while they use numerical differentiation which is very sensitive to measurement noise. The moments technique relies upon measurement of the zeroeth moment which is difficult to assess with the present sizing technique. This also handicaps the Laplace technique as datapoints at decreasing sizes are increasingly important. This does not comply with the accuracy of the measuring technique. -The method of characteristics produces sensible results. The outcome depends on the size class from which the nucleation rate is recalculated as a result of the curved population density plots. A strong qualitative correlation with the heat input appears. Quantitative correlation appears to be unrealistic as a result of the scatter. Therefore a fourth technique has been developed. -The fourth technique appears to be extremely effective in calculating the nucleation parameters required to predict the response of the second moment or the average size separately, additional to the response of the third moment which does not depend on the nucleation kinetics. It appears that both the average size and the second moment can not be reliably predicted so long as crystal attrition and or growth dispersion are excluded from the model. The technique is flexible enough to incorporate these phenomena, if the appropriate models are available. With respect to the ammonium study it appears that:
sulfate system
under
CP G i, i’, k
k, k,, kb
k, k k,, 4, k, L MS1
mj N
n0 4-c t)
d”,(X) n’(s, t) P
Q
Rv S
T V
WI x X Y
Y Greek letters d. W,, ALx & P
-Secondary nucleation kinetics prove to be effective in describing nucleation kinetics, Surprisingly, the system itself reacts only slightly to changes in the kinetic parameters. Whether the second or the third moment is used proves to have little influence on the final results. The necessity to add a nuclei stock is not apparent. The influence of other population events is not negligible.
empirical constant nucleation rate, No./m3 s concentration, kg/m3 clear liquor concentration at saturation, kg/m3 clear liquor specific heat, kJ/kg K growth rate m/s empirical constants used to denote the moment j shape factor ( = 3.14) nucleation constant, mass transfer coefficient, m/s shape factor ( = 0.26) constants crystal size, m slurry density, kg/m3 slurry jth moment of the distribution stirrer speed, rps population density at zero size, No./m4 population density at size L, No./m4 cumulative contribution of the xth size class to the Laplace transformed population density Laplace transformed population density, No./m3 heat input, W flowrate, m3/s heat of evaporation, kJ/kg supersaturation [(C - C,)/C,] Laplace variable, l/m time, s temperature, K crystallizer volume, m3 weighing factor index x-co-ordinate index function value
r
t)
volume fraction in the xth size class width of a size interval x, m void fraction density, kg/m3 residence time, s
Indices
5 P V no index
denotes denotes denotes denotes denotes denotes
crystal feed liquid product vapour crystallizer
contents
Estimation
of nucleation
kinetics
from crystal -no -no
REFERENCES Berglund, K. A. and Larson, M. A., 1984, Modeling of growth rate dispersion of citric acid monohydrate in continuous crystallizers. A.I.Ch.E. J. 30. 280. Boxman, A. and Scarlett, B., 1988, On-line measurement of crystal size and shape using combined optical techniques, in SPIE Vol. 952, Proceedings of "Laser Technologies in Industry” (Edited by 0. D. D. Soares), p. 520. Cooper, D. J. and Clough, C. E., 1985, Real-time estimation of dynamic particle-size distributions in a fluidized bed: theoretical foundation. 4.I.Ch.E. J. 31, 1202. Daudey, P. J., 1987, Crystallization of ammonium sulfate, secondary nucleation and growth kinetics in suspension. PhD thesis, Delft University of Technology. Daudey, P. J. and Jong, E. J. de, 1984, The dynamic behaviour of NaCl in a 91 liter MSMPR-crystallizer, in Industrial Crystallization 84 (Edited by E. J. de Jong and S. J. Jancic), pp. 447-450. Elsevier, Amsterdam. Gill, P. E., Murray, W. and Wright, M. H., 1981, Practical Optimization. Academic Press, London. Heffels, S. H., 1986, Product size distributions in continous PhD thesis, Delft and batch sucrose crystallizers. University of Technology. lager, J., Wolf, S. de, Klapwijk, W. K. and Jong, E. J. de, 1989, A new design for on-line CSD-measurements in product slurries. Presented at the 10th Symposium on Industrial Crystallization, Bechyne Castle. Jager, J., Wolf, S. de, Kramer, H. J. M. and Jong, E. J. de, 1990, On-line particle size measurement in dense slurries. in press. Jancic, S. J. and Grootscholten, P. A. M., 1984, Industrial Crystallizarion. Delft University Press, Dordrecht. lanse, A. H. and Jong, E. J. de, 1976, The occurrence of growth dispersion and its consequences, in Industrial Crystallization (Edited by J. W. Mullin). Plenum Press, New York. Marchal, P., David, R., Klein, J. P. and Villermaux, J., 1988, Crystallization and precipitation engineering-I. An efficient method for solving population balance in crystallization with agglomeration, Chem. Engng Sci. 43, 59. McCabe, W. L., 1929, Crystal growth in aqueous solutions. lnd. Engng Chem. 21, 30. Randolph, A. D. and Larson, M. A., 1988, Theory of Particulate Processes. Academic Press, London. Shah, M. B. and Garside, J., 1980, Crystallization kinetics from MSMPR crystallizers. Ind. Engng Chem. Process Des. Deu. 10, 509. Tavare, N. S., 1986, Crystallization kinetics from transients of an MSMPR-crystallizer. Can. J. them. Engng 64, 752. Tavare, N. S. and Garside, J., 1982, Simultaneous estimation of crystal nucleation and growth kinetics from batch experiments. Chem. Engng Res. Des. 64, 109. Tsuruoka, C. and Randolph, A. D., 1987, State space rep resentation of the dynamic crystallizer population balance: application to CSD controller design A.J.Ch.E. Symp. Ser. No. 253, 83, 104. Wolf, S. de, Jager, J., Kramer, H. I. M., Eek, R. and Bosgra, 0. H., 1989, Derivation of state-space models of continuous crystallizers. Preprints of IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, 21-23 August 1989, Maastricht.
APPENDIX
A: DERIVATION AND GROWTH
OF PRODUCT
FLOWRATE
--constant volume, constant temperature, -neglectable heat of crystallization, --crystal-free feed flow. -salt-free condensate.
transients
817
growth dispersion or size-dependent growth agglomeration, crystal breakage or attrition,
rate,
the salt balance, the total mass balance and the heat balance can be simplified into a constraint for the growth rate and a constraint for the product flowrate. In this derivation the time variable will be omitted in the notation. The salt balance is given by $ t
CCs+ p,(l - 41 VI = C,Q, - CCE+ PAI - 41 Q, (Al)
where the void fraction
is defined
by
E = 1 - k,m,
(A21
and thus ds = Pk”d!$
(A3)
dt The moment
equation
for i = 3 follows from eq. (12):
dm3 dt
_=
--+33Gm, r
Combining eqs (Al)-(A4) and using the Class-II dC/dt = 0 results in a growth rate constraint:
G=
is given by
+ AU - 41Qp - P,Q,, W4
p,Q, - Cw Inserting
assumption
Q&,--Q& % Vm,( pc - Cl
The total mass balance
eqs (A2)-(A4)
in eq. (A6) gives VGm,.
The heat balance
(A7)
is given by
- Ccc,p,~+ P,c,,(~ - 41 TQ, - p,Q,(R, + cpu77 + P. WI The heat balance can be simplified into a constraint vapour flow by inserting eqs (A2)-(A4) in eq. (A8):
*(P,c, - ~01.
+ P - 3k,% Combining lC
+
RATE CONSTRAINT
In this Appendix it will be shown that, evaporative MSMPR crystallizer with:
size distribution
this constraint
C(P, - Pr)
7% R, + c,,T
- P&J
+ cpv7-HPC - Cl 1
ML
(A9)
with eqs (A5) and (A7) results
P,(P. - C)
*UP&
for the
in
CP,(R, + cpu *)I Q,
for a Class-II
*) - TJC,,,P,l -
-
*cP&
- P&l)
(P, - C)
1C, >
Q, - P
s(R. r
+ e,,T)
(AI’3
818
JOHAN JAGER et al.
Or
Qp=
k,
k, = P,& + cpv~) - T,c,,/P, k _; 2
1 "c'(R,+ QT) c
k,=
(k, - k,C,)Q, - P
C( PC - PI) -~
Tc,,
(All) (A121
+ ~Chc, - w,,) b,(R, + cpoZ-11. PdR” + C&s”WP, - q 1
_ T(pcc, -p&!) (PC--_
I-
[
6413)
PdP, - Cl
R,+c,,T 6414)