Estimation of Nucleation and Growth Kinetics of Ammonium Sulfate

Abdolsamad Tadayon, Sohrab Rohani,* and Michael K. Bennett. Department of Chemical and Biochemical Engineering, The University of Western Ontario,...
5 downloads 0 Views 294KB Size
Ind. Eng. Chem. Res. 2002, 41, 6181-6193

6181

Estimation of Nucleation and Growth Kinetics of Ammonium Sulfate from Transients of a Cooling Batch Seeded Crystallizer Abdolsamad Tadayon, Sohrab Rohani,* and Michael K. Bennett Department of Chemical and Biochemical Engineering, The University of Western Ontario, London, Ontario, Canada N6A 5B9

Nucleation and growth kinetics of ammonium sulfate are determined using a 1.5-L cooling batch seeded crystallizer, a rigorous mathematical model of the process, and a nonlinear optimization program. Supersaturation and crystal size distribution are employed as the major outputs of the model in the optimization objective function. The former is estimated by an online density meter, and the latter is measured by a Malvern Mastersizer. The size-independent growth kinetic is shown to lack the necessary parameters to accurately fit the experimental data. The information gathered from a large number of measurements in this study indicates that ammonium sulfate possesses a size-dependent growth kinetic. It is shown that a second-order size-dependent growth model reasonably describes the experimental measurements, while a combination of a second-order size dependence for small and medium crystals and a slightly lower dependence (order of 1.85) for large crystals results in the least error between experimental and simulated results. Introduction Solution crystallization is an important separation process for the production of bulk chemicals such as potash and sucrose as well as pharmaceuticals and photographic materials. Purity and crystal size distribution (CSD) are the two main quality indices of crystals. Depending on the type of material, small or large crystals, usually with uniform CSD, may be desired. Pharmaceuticals for inhalation or injection purposes are often required to have a small size, while for bulk chemicals, large crystals are desirable in order to promote flowability and to avoid caking during storage. CSD is the result of the competition between nucleation and growth kinetics in the solution.1 Smaller crystals tend to form when the operating conditions and the types of materials favor nucleation over growth. Crystallization kinetics explain how much of the supersaturation is consumed by the nucleation and growth and how much remains in the solution as operating conditions change. The kinetics are required for process simulation, optimization, and advanced crystallization control. Of the numerous techniques for crystallization kinetics estimation, the method devised by Randolph and Larson1 using a mixed suspension mixed product removal (MSMPR) crystallizer is perhaps the most popular. However, the technique is laborious and timeconsuming and needs an elaborate crystallization arrangement.2 Frequent plugging of the tubes containing a supersaturated solution makes the operation tedious and imposes major disturbances on the system. In addition, the analysis does not directly give the kinetic parameters, and further computations for extracting these parameters are required.3 Moreover, experimental evidence3 indicates that the steady-state kinetic parameters do not correctly describe unsteady-state operations. The problem is partly due to the underlying assumptions such as the size-independent growth rates * To whom correspondence should be addressed. E-mail: [email protected].

in the construction of the simple steady-state CSD equation: n ) n0 exp(-L/Gτ). A problem with this technique is that a small error between the experimental data and the fitted line results in an error in the estimated G and n0. When the kinetic correlations are used in first principle models, they often lead to large errors for the estimated solution concentration and the supersaturation. This is more pronounced when the semilog plot of n vs L appears to be nonlinear. An alternative approach for the estimation of kinetics is employment of unsteady-state data from batch experiments. Temperature, CSD, concentration, and magma density are all changing with time, providing information over a wide range of operating conditions in one experiment. Kinetic estimation using this approach requires fewer experiments than the MSMPR approach. Also, the batch crystallizer is easier to operate, and the experiments are less time-consuming. However, the flexibility and the ease of batch experiments are partially offset by their computational complexity.2 Steady-state crystallizers can be described by simple ordinary differential equations (ODEs), while the variables in the batch crystallizers are changing with time and crystal size, leading to hyperbolic partial differential equations (PDEs). These equations either need to be solved by transformation to an ODE or by use of an appropriate numerical technique. In this study, the crystallization kinetics of ammonium sulfate, an important inorganic salt, in a batch cooling crystallizer are investigated. Background Several studies have been concerned with the determination of kinetics of the ammonium sulfate system. Mullin et al.4 used a simple method devised by Nyvlt and measured the nucleation and growth rates of ammonium sulfate. The supersaturation was measured based on the measurement of the saturation temperature. Nucleation was determined with only two crystals present in the vessel, and the growth rate of two major

10.1021/ie020163r CCC: $22.00 © 2002 American Chemical Society Published on Web 10/29/2002

6182

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002

faces of a single crystal was determined with a microscope. The method did not take into account the presence of seed crystals in a quantitative manner. Tavare5 employed the method of initial derivatives for the estimation of the kinetics of ammonium sulfate. The original technique, developed by Garside et al.,6 was based on the calculation of time derivatives from the initial portion of a single seeded isothermal desupersaturated curve. It was later extended to batch cooling crystallizers with variable solution temperature.5 The technique neglects mass deposition due to nucleation and therefore is applicable only to growth rate calculation. The advantage of the method is its simplicity, while its sensitivity to experimental errors, full reliance on concentration measurements (absence of CSD in the calculations), and negligible nucleation remain major disadvantages. The moment method2 converts the population balance to a set of ODEs:

mj(t) )

∫0∞n(L,t) Lj dL

(1)

The zeroth and subsequent moments for a batch crystallizers are

dm0(t) ) B(t) dt dmj(t) ) jmj-1(t) G(t) dt

(2) jg1

(3)

These equations can be easily used to obtain the nucleation and growth rates of a crystallization system. The method is based on the measurement of the population density alone, and the direct contribution of the concentration measurement to the computation of the nucleation and growth rates is absent. Using methods that employ both CSD and concentration increases the confidence in the estimated kinetics. The growth and nucleation rates obtained from eqs 2 and 3 and the experimental data are used to estimate the kinetic constants from the following equations:

G ) kg∆Cg

(4)

B ) kb Gi M

(5)

In these correlations, supersaturation, ∆C, needs to be measured experimentally while, magma density, M, may be measured experimentally or computed from the drop in the original concentration. Jones and Mullin7 distinguished between the seed crystals (S crystals) and the new crystals (N crystals) in a batch crystallizer and derived the following equation for supersaturation:

-d∆C dWS dWN 3WS0LSG 3FckvANG + ) + ) dt dt dt ka L 3S

(6)

S0

Equations 2, 3, and 6 are the major balance equations that describe a seeded batch crystallizer. Using information provided by these equations and the experimental measurement of the CSD and supersaturation in an optimization objective function gives more reliable kinetic parameters than by just using eqs 2 and 3 and the measured CSD.

The most accurate approach for the estimation of kinetic parameters is accomplished by using a nonlinear optimization algorithm, a full mathematical description of a batch crystallizer, and reliable experimental data. A major advantage of this method is that the estimation of growth and nucleation rates and their kinetic parameters in the mathematical method is performed simultaneously. The kinetic constants in the mathematical model are adjusted by the optimization algorithm to minimize an objective function. The procedure is terminated when kinetic parameters in the model result in a minimum error between the estimated and experimental data. Both the moment and the s-plane methods require the assumption of a constant nucleation and growth rate during a sampling time. Quite often the sampling time for the CSD samples is long and the above assumptions are not valid. These unrealistic assumptions can be avoided by using more rigorous numerical methods such as finite difference methods instead of moment methods. In finite difference schemes such as Lax-Wendroff,8 Crank-Nicholson,9 or the combined Lax-Wendroff and Crank-Nicholson,10 all state variables including the population density and concentration are computed at every computational time instant. The time interval is usually short, on the order of 1 s. The seed crystals constitute part of the initial conditions of the model, and therefore the seeds and new crystals are not required to be distinguished from each other. Variables such as the crystal mean size, population density, concentration, and magma density may be used in the optimization objective function. Reducing sampling times for the measurement of each of the above variables increases the number of measurements employed in the objective function and therefore decreases the uncertainty and increases the reliability of the estimated parameters. It is noted that the optimization algorithm does not require a constant sampling time for all or individual variables in the objective function. The variables in the objective function should be selected carefully. If both growth and nucleation kinetics are to be estimated, an indication of the size distribution (population, length, area, or mass) and concentration should be included in the objective function. Including only concentration, magma density, or supersaturation in the performance index may result in accurate prediction of these variables but will result in an error in the prediction of CSD. Such observations have been made by Miller and Rawlings.11 The reason is that supersaturation in the crystallizer can be consumed by both the nucleation and growth processes, while an objective function with no CSD measurement cannot distinguish the difference between supersaturation consumed by these two phenomena. Population, area, or mass density distribution alone may be used in the objective function, because their units are a combination of the CSD and solvent mass, and therefore they can be utilized to compute the solution concentration. If a precise population density is predicted by the model, an accurate simulated concentration or supersaturation level is also obtained. Because of measurement and noise errors in the population density, the model prediction of this variable is usually not very accurate. This inaccuracy, especially when it occurs for large particles, results in an error in the prediction of supersaturation. This, in turn, propagates the error into growth and nucleation rate calculations. To avoid these problems, both the

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002 6183

population density and concentration need to be included in the objective function. Instruments such as a Malvern Mastersizer report CSD as percent volume, area, length, or number of particles in certain size intervals. As a result, the reported magnitudes do not contain indications of the mass or volume of the solvent. Therefore, using this information alone without further treatment to convert to the population density does not provide sufficient information for calculation of the concentration or supersaturation. Theory Model Development. Major equations describing batch crystallization dynamics include solute, solvent, energy, and population balances. Simultaneous solution of these equations gives the crystal magma density, solution concentration, temperature, and population density. If the crystallizer temperature is measured experimentally, the energy balance becomes idle and can be removed from the set of equations. The population balance for a seeded batch crystallizer in which crystals are born at size zero, crystal agglomeration and breakage are absent, and crystal shapes are uniform is

∂n(L,t) ∂(G(L) n(L,t)) + )0 ∂t ∂L

(7)

In this study both size-independent and -dependent growth rates are studied. The initial condition for eq 7 is the population density of seed crystals. The boundary condition is

B0 n(t,0) ) G|L)0

(8)

The solute balance for the crystallizer is

dm3(t) dC(t) dM(t) ) -Fkv )dt dt dt

(9)

where m3(t) is the third moment and is defined in eq 1. The auxiliary equations include nucleation, growth, and solubility information:

G ) kg∆Cg exp(-Eg/RT)

(10)

B ) kbGiMj exp(-Eb/RT)

(11)

The stirrer speed is maintained constant in all experiments. Solution of the Population Balance Equation. Equation 7 falls in the general class of hyperbolic PDEs known as flux conservation equations. An analytical solution for these equations does not exist. Several finite difference schemes for the solution of this class of equations are available. The Lax-Wendroff8 method has been used for this purpose. Expending the second term in eq 7 and discretization using this technique12 lead to

) npq - (Gλ1 + G˙ lnpq)∆t + np+1 q ∆t2 (12) 2

(G2λ2 + 2G˙ lGλ1 + G˙ l2npq)

p p p where λ1 ) (nq+1 - nq-1 )/2∆x and λ2 ) (nq+1 - 2npq + p nq-1)/∆x2 and G˙ l is the derivative of the growth rate

with respect to size. The approach is explicit and gives an error order of accuracy of [(∆L)2, ∆t)]. Bennett13 argued that implicit methods are unconditionally stable and therefore are preferred over the explicit methods. He employed the Crank-Nicholson method:9 p-1 p+1 p p λ3nq-1 + λ4np+1 - λ3nq+1 ) -λ3nq-1 + λ5npq + λ3nq+1 q (13)

where λ3 ) -G∆t/4∆L, λ4 ) 1 + G˙ l∆t/2, and λ5 ) 1 G˙ l∆t/2. This results in a system of equations, forming a tridiagonal matrix with an error order of accuracy of [(∆L)2, (∆t)2]. The tridiagonal can be solved using the Thomas algorithm.14 Bennett and Rohani10 have shown that while both the Lax-Wendroff and Crank-Nicholson algorithms work well for crystallization systems with smooth population density functions, they fail when seeds or evolving population densities contain a discontinuity. The discontinuity in the crystal population density is inevitable because seed crystals are usually much larger than the newly formed crystals and the population density for the size between seeds and the nuclei is zero, resulting in a discontinuous population density. To solve this problem, these investigators combined the Lax-Wendroff and Crank-Nicholson methods: p+1 p+1 + (1 + 2λ6 - λ8)np+1 - (λ6 + λ7)nq+1 ) (-λ6 + λ7)nq-1 q p + (1 - 2λ6 + λ8)npq + (λ6 + λ7) (14) (λ6 - λ7)nq-1

where λ6 ) G2∆t2/4∆L2, λ7 ) (-G∆t + GG˙ l∆t2)/4∆L, and λ8 ) (-2G˙ l∆t + G ¨ l2∆t2)/2. This system of equations also results in a tridiagonal matrix and can be solved with the Thomas algorithm. Optimization Algorithm. The Nelder-Mead15 Simplex approach was used for the minimization of the objective function. A simplex is a geometrical figure formed by r + 1 points in an r-dimensional space. In two-dimensional space, a simplex is a triangle; in threedimensional space, it is a pyramid. The algorithm is based on evaluating the objective function at the vertexes of a simplex, finding the worst vertex of the simplex, forming its symmetrical image through the center of the opposite face, and generating a new simplex based on the new vertex. This shrinks the volume of the simplex if the newly found point is better than the older point. This step is repeated until the volume of the simplex is smaller than some user-defined preset value. Any commercially available software may be used for this purpose. Function fminsearch in the optimization tool-box of Matlab (Mathwork, MA) uses the Nelder-Mead Simplex approach and can be employed for the minimization of the objective function. Experimental Section. In this study, the solution concentration (C) of ammonium sulfate in the crystallizer was estimated by measuring the solution density (F) and temperature (T). Before the kinetic runs were done, some experiments were carried out to develop a mathematical relationship between C, F, and T. The experimental apparatus is shown in Figure 1. The experiments were performed in a 500-mL doublejacketed glass vessel (Bellco Glass, Vineland, NJ) placed on a magnetic stirrer. The solution temperature in the vessel was varied by manipulating the temperature of the flowing water through the jacket. Water manipulation was carried out by a variable-setpoint water bath

6184

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002

Figure 1. Experimental apparatus for developing a correlation between F, C, and T.

system (RTE 220; Neslab instruments, Inc., Portsmouth, NH). Samples were continuously withdrawn from the solution by a peristaltic pump at 6.2 mL/min, passed through an online density meter (MPDS2000; Anton Paar, Graz, Austria), and then recycled to the solution vessel. The density meter had an accuracy of 10-5 g/mL and was set to measure the density as well as the temperature of the solution passing through its sensing zone. A LabView (National Instruments, Austin, TX) hardware/software system was employed for data acquisition. This system was set to receive the density and the temperature data from the density meter and display these variables on two charts online. Solutions with known concentrations (69.452, 71.407, 73.393, and 75.842 g/100 g of water) were prepared by adding an appropriate mass of analytical-grade ammonium sulfate (BDH Inc., Toronto, Canada) and deionized distilled water to the experimental vessel, covering the top and mixing until all crystals were dissolved. The temperature of the solution was decreased from 30 °C to the nucleation points approximately in 10 steps. At each step, the solution temperature was maintained constant until a constant density and temperature by the density meter were recorded. The concentration and density meter readings were recorded for all experimental runs. This information was fitted in a deterministic relation between C, F, and T, and the parameters were estimated. This equation along with F and T measured by the density meter was used to estimate the concentration of the solution in the vessel during the kinetic experiments (eq 15). Solubility Measurement. The same experimental setup (Figure 1) was utilized for the measurement of the saturation concentration of ammonium sulfate. Three hundred grams of deionized distilled water and an excess amount of solids were added to the vessel. The temperature of the vessel was fixed at 0 °C, and the magma was stirred until no change in the density meter reading was observed for 1 h. The flow rate of the solution sample in the density meter was adjusted to allow only a clear solution to be withdrawn from the vessel. This low flow rate imposed a transportation lag of around 90 s in reading the density of the solution in the vessel. This delay time was later considered in the

analysis of the results. After density meter readings were recorded at a temperature, the solution temperature was increased to a new value. This procedure was repeated for seven temperatures between 0 and 31 °C. The solubility concentrations were estimated from these measurements using the equation relating C, F, and T. The solution temperatures and their associated estimated concentrations were fitted to a third-order polynomial equation, and the correlation parameters were estimated (eq 16). Seeds Preparation. To prepare the seeds for experiments, an appropriate proportion of crystals and deionized distilled water were added to the crystallizer to give a saturated solution at 35 °C. To speed up the dissolution process, the temperature was first increased to 40 °C and then immediately cooled to 23 °C by flowing water at 22 °C through the crystallizer jacket. This resulted in the formation of fine crystals through primary homogeneous nucleation. In these runs, it was intended to produce seeds with sizes smaller than 600 µm. Slow cooling rates provided newly generated crystals with more time to grow and produced crystals larger than 600 µm. Fast cooling of solutions saturated at temperatures higher than 33 °C resulted in the formation of undesirable nail-shaped crystals. After each run, the entire contents of the crystallizer was filtered using a 20-µm filter paper. The crystals were then washed with a solution of 2-propanol (EM Science, Merck KgaA, Darmstadt, Germany) saturated with ammonium sulfate. The crystals were collected and dried at room temperature and then in an oven at 60 °C. This procedure was repeated until around 300 g of seeds was prepared. Crystals larger than 525 µm and smaller than 20 µm were removed using a Gilsonic autosiever (Gilson Autosiever, Gilsonic Co. Inc., Worthington, OH). The crystals were divided into 10 parts using the coning and quartering method. Approximately 5-g samples were used for the CSD measurement by a Malvern Mastersizer 2000 (Malvern Instruments Ltd., Worcestershire, U.K.). Crystallization Control System. The quality of crystallization experimental data is strongly dependent on the accuracy of the crystallization temperature control. Undesirable oscillation in the temperature of

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002 6185

Figure 2. Three-level cascade control system.

the crystallizer may make the solution highly supersaturated or shift that to the undersaturation region. The former may cause the solution to nucleate through primary homogeneous nucleation, and the latter provides the conditions for the dissolution of particles. Both of these phenomena are undesirable, and an appropriate control strategy is required to avoid their occurrence. In addition, fast control of the crystallizer temperature in tracking its setpoint is required. A simple proportional-integral-derivative (PID) controller usually does not provide the necessary performance to achieve these objectives. For kinetic experiments, a three-level cascade control system, shown in Figure 2, is employed. To accomplish this, the crystallizer temperature (Tc), the temperature of water flowing through the jacket (Tj), and the chiller temperature (Tch) are measured and transferred to the data acquisition system. The crystallizer temperature is compared with the desirable setpoint (Tc,sp), and the resulting error enters the first PID controller (TC1). The output from this controller is the setpoint for the jacket temperature (Tj,sp). The error between the jacket and its setpoint is used in a second PID controller (TC2) to compute the setpoint for the chiller temperature (Tch,sp). A third PID controller (TC3) in the chiller receives this setpoint from the computer and adjusts the cooling or heating power to bring the chiller temperature to its desired value. For this system the crystallizer and the jacket controllers act as master controllers for jacket and chiller controllers, respectively. Tuning of the controllers in the cascade control system was performed by trial and error starting from the chiller controller to the crystallizer controller. The sampling time for controllers was 10 s. Kinetic Experiments. The experimental setup used for the kinetic runs was basically the same as that shown in Figure 1 except that the 500-mL vessel was replaced with a 1.5-L double-jacketed glass crystallizer. The variable-temperature water bath was also replaced by a 5-kW, 30-L chiller (HX-150; Thermo Neslab Instruments, Inc., Portsmouth, NH), and the stirrer was driven by a motor; J-type thermocouples were used to measure the crystallizer and jacket temperatures. A typical run started with the preparation of 1.5 L of a saturated solution at 30 °C by the addition of appropriate amounts of deionized distilled water and ammonium sulfate particles. Vaporization of the solvent was reduced by covering the crystallizer top by Plexiglas. The density meter and thermocouples were scanned

every 10 s, which provided around 400-1000 sets of data during one batch. Before adding the seeds, the level of supersaturation was increased to avoid dissolution of small crystals due to possible undersirable control performance. However, this supersaturation was not sufficient to induce primary nucleation when seeds were added. Preweighed, presized seed crystals were added to the crystallizer at time zero, and the crystallizer temperature was set to follow the setpoint. Two to four samples were taken for size analysis in each run. To accomplish this, 70-120 mL of magma was withdrawn and weighed. The amount of magma taken for the CSD analysis was to satisfy the Malvern Mastersizer requirement. The experiments were continued until supersaturation in the crystallizer approached zero. The CSD samples were immediately filtered using a 20-µm filter paper. To avoid agglomeration upon drying, the crystals were then washed with 2-propanol saturated with ammonium sulfate. The wet crystals on the filter paper were collected and immediately used for size analysis. A 2-propanol solution saturated with ammonium sulfate was also used as the suspension liquid for the size measurement with a Malvern Mastersizer. Before each measurement, the 500-mL suspension chamber of the Mastersizer sensing zone and the connecting tubes were rinsed and then filled with the suspension liquid. Approximately 4-8 g of wet crystals was added to the chamber, vibrated for 20 s with ultrasound prior to size distribution measurement. Results and Discussion Developing a G, C, and T Correlation. Figure 3 shows the density of the solution at different concentrations vs temperature. As can be seen, the density decreases with temperature almost linearly. The lines are almost parallel, showing that changing the concentration over the range of experimental conditions does not considerably affect the slope of these lines. While it is possible to fit the data to a polynomial, fitting the data to the model

C)

k1F[1 + R1(T - 273)] + 1 k2F[1 + R2(T - 273)] - 1

(15)

gives a better accuracy. Using a nonlinear search

6186

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002

Figure 3. Density of ammonium sulfate at different concentrations measured at different temperatures.

Figure 4. Estimated vs experimental concentration.

optimization algorithm, the model constants k1, R1, k2, and R2 are estimated as -0.9599, 1.9343 × 10-4, 0.5858, 5.4613 × 10-4, respectively. The experimental and model concentrations are depicted in Figure 4. The R2 value is better than 99%. Developing a Solubility Correlation. Experimental data and eq 15 are used to estimate the saturation concentration of the solution at various temperatures. The results of this study and other work found in the literature are sketched in Figure 5. Mullin et al.4 and Tavare5 reported the solubility of ammonium sulfate. Tavare presented a correlation for the solubility without giving actual experimental measurements. The correlation is shown over its range of validity. The solubility measured in this work is slightly larger than those of

Mullin et al. and Tavare but is in the range of those of Seidel (stated in ref 4). The data generated in this work are fitted to a third-order polynomial equation:

C* ) -1 × 10-7(T - 273)3 + 1 × 10-5(T - 273)2 + 0.0022(T - 273) + 0.708 (16) The difference between the results may be attributed to the method by which supersaturation is measured. Mullin et al. measured the saturation temperature for known solution concentrations, while Tavare employed gravimetry for the concentration estimation. In the present work, the concentration is estimated from direct measurement of the density and temperature. This method did not include offline sample removal, which

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002 6187

Figure 5. Solubility concentration of ammonium sulfate at different temperatures. Table 1. Ranges of Experimental Variables concentration, kg/kg of solvent seed loading, kg/kg of solvent seed size, µm stirrer speed, rps supersaturation, kg/kg of solvent temperature, °C

0.70-0.80 0.013-0.025 20-525 13.25 0.002-0.011 0-30

usually contributes to the error in the analysis for the solution concentration, and therefore is expected to have better accuracy than others. Perhaps using each of these equations in a simulation program gives approximately similar results. However, if the simulated program is combined with experiments in which the supersaturation is calculated based on the measurement of concentration (such as the determination of the optimal temperature path based on a model), the solubility equation should be selected with great care. The reason is that supersaturation defined by the difference between the actual and saturation concentrations is very small compared to the actual concentration (∆C/C* ) 0.015), and a small error in the saturation concentration may lead to a large error in the supersaturation. Because supersaturation is the basic variable upon which other major crystallizer outputs are predicted, even a small error in that variable may lead to large errors in these outputs.

Figure 6. Typical temperature and supersaturation profiles.

Kinetic Experiments. Temperature profiles in batch crystallizers usually have an overall decreasing trend to provide supersaturation for crystal growth and nucleation. The temperature profile should be carefully designed to prevent dissolution of available crystals or formation of new crystals through primary nucleation. Table 1 lists the ranges of experimental variables in this study. Figure 6 displays a typical temperature and supersaturation profile used in the present work. The temperature is reduced from 29 to approximately 3 °C. This decreases the supersaturation from 0.0085 to 0.000 05 kg/kg of water. The temperature is carefully oscillated after t ) 4500 s to make variations in the supersaturation in order to provide the optimization program with more data. The highest supersaturation in all experiments is 0.011 kg/kg of water. Preliminary experiments showed that larger supersaturation would lead to a burst of nucleation which is an indication of primary nucleation. Similar observations were made by Youngquist and Randolph16 and Mullin et al. 4 for the upper limit of supersaturation for this system. Therefore, ammonium sulfate crystallizers may be categorized as class II crystallizers in which the supersaturation level for stable operation of crystallizers is very low. However, these findings do not agree with Tavare’s5,17 work, in which a value of 0.025 kg/kg of water is reported. On the basis of our observations, such a high supersaturation would lead to a shower of nuclei in the crystallizer which cannot be considered as secondary nucleation. Two simulation programs, one for size-independent growth rates and the other one for size-dependent growth rates, were employed in the optimization programs. The objective function in the optimization program minimized the difference between the measured and simulated variables: Ne Nm Ns

min Φ(θ) ) θ

∑ ∑ ∑ wj[yijk,m - yijk(θ)] i)1 j)1k)1

(17)

subject to eqs 7-10 where θ denotes the kinetic parameter vector and yijk,m

6188

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002

Figure 8. Comparison between experimental and simulated concentrations with a size-independent growth rate (eqs 18 and 19).

Figure 9. Comparison between experimental and simulated supersaturation with a size-independent growth rate (eqs 18 and 19).

Figure 7. Comparison between experimental and simulated CSDs with a size-independent growth rate (eqs 18 and 19): (a) t ) 1983 s; (b) t ) 4983 s; (c) t ) 8563 s.

and yijk(θ) represent the jth measurement and simulation prediction, respectively, taken at the kth sampling time in ith experiment. Ne is the number of experiments, Nm is the number of measured variables, and Ns denotes the number of samples taken for each variable in one run. The weight factor, wj, is assigned for each measured variable to recognize the relative importance of each variable. Two variables, supersaturaion and the percent volume for each size channel, are included in the optimization objective functions. Therefore, the objective is to find a set of kinetic parameters that minimizes the

difference between measured and simulated supersaturation and the volume distribution. The CSD measurements are performed for the size range between 30 and 1200 µm in 28 size channels. The reason for considering a lower size limit is that crystals smaller than this size have negligible volumes and removal of these small crystals from the filter paper and their subsequent accurate measurement is difficult. Including the volume percent of small size crystals is appropriate only when an in situ crystal size measuring device is used. The question may arise about how the nucleation rate can be obtained without measuring near nuclei size crystals. Experimental results in this study indicate that the new ammonium sulfate nuclei grow up to about 100 µm during an approximately 10 000-s batch experiment. In this work, the nucleation rate is inferred from the size distribution of crystals of 30 µm < L < 100 µm. Only the correct estimation of the nucleation and growth rates leads to the correct size distribution of these crystals. The optimization algorithm varies nucleation and growth parameters in the model until the error between simulated and experimental variables is minimized for crystals larger than 30 µm. A correct match between the experimental and simulated CSDs for small crystals is an indication of an accurate estimation of nucleation and growth rates.

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002 6189

Figure 11. Comparison between experimental and simulated concentrations with a size-dependent growth rate (eqs 20 and 21).

Figure 12. Comparison between experimental and simulated supersaturation with a size-dependent growth rate (eqs 20 and 21).

Figure 10. Comparison between experimental and simulated CSDs with a size-dependent growth rate (eqs 20 and 21): (a) t ) 1983 s; (b) t ) 4983 s; (c) t ) 8563 s.

An important issue in the formulation of the optimization objective function is the scaling of the variables. The errors for supersaturation are in the range of 0.003 kg/kg of solvent, while that for the percent volume in a size channel is around 1. The optimization on its own cannot recognize the importance of a unit change in a variable. The objective function without proper scaling or weighing factors may result in too much sensitivity to the CSD while ignoring the supersaturation. On the other hand, using large values for the supersaturation weighing factor ignores the kinetics and the population

balance and leaves the model with mass and energy balance equations. A second issue in this regard is the normalization of the errors for CSD in different size channels. Small crystals have very small volumes compared to the larger ones. As a result, the error for small crystals may be much smaller than that for large crystals; this may shift the sensitivity of the optimization to the larger crystals while neglecting the smaller crystals. If the CSD error is not properly normalized, the optimization algorithm treats the CSD error for all size channels of crystals equally and this may result in the large errors between simulated and experimental CSDs of small crystals. The CSD of fine crystals heavily influences the nucleation kinetics of crystals. In this study, the normalization was performed by dividing the error between the simulated and experimental CSDs for each size channel by the maximum error between these two size distributions. When experimental data were used with the sizeindependent growth rate equation, the optimization algorithm gave the following kinetic correlations:

G ) 1.0147 × 105∆C0.767 exp(-61751/RT) (18) B ) 1.0106 × 1013G0.221M3 exp(-15228/RT)

(19)

Parts a-c of Figure 7 display the measured and simulated CSDs at t ) 1983, 4983, and 8563 s, respec-

6190

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002

Figure 14. Comparison between experimental and simulated concentrations with a modified size-dependent growth rate.

Figure 15. Comparison between experimental and simulated supersaturation with a modified size-dependent growth rate.

for smaller crystals. From these plots, it was concluded that crystals of various size ranges do not grow at the same rate, and therefore a size-dependent growth algorithm was required. The new algorithm needed to prescribe larger growth rates for bigger crystals. Tavare’s5 simulation results using size-independent growth rates did not match his experimental data. Youngquist and Randolph16 also reported nonlinear log(n) vs L for the ammonium sulfate system, which indicates a sizedependent growth rate. In this study, the following sizedependent model was found to best fit the experimental data: Figure 13. Comparison between experimental and simulated CSDs with a modified size-dependent growth rate.

tively. Comparisons between the measured and simulated concentrations and supersaturation are presented in Figures 8 and 9. In general, the experimental results do not agree well with the simulated outputs. From parts b and c of Figure 7, it is noticed that sizeindependent growth kinetics overestimate the growth rate and the population density of small crystals while they underestimate those of large crystals. This problem cannot be fixed with changes in the available parameters in the size-independent growth kinetic correlations. Increasing the growth rate to decrease the difference between the measured and simulated plots for large crystals results in an overestimation of the growth

G ) 2108∆C2.42 exp(-32027/RT)(1 + 8.1 × 106L2) (20) B ) 3.5281 × 1014(G|L)0)0.55M0.58 exp(-25845/RT) (21) where both the coefficient and the order of L in eq 20 are obtained from the optimization algorithm. The growth rate was found to be second-order-dependent on the crystal size. Figures 10-12 depict the experimental and simulated results using the above kinetics. Both simulated supersaturation and CSD are in reasonable agreement with experimental data. The difference between associated plots could be due to the measurement error including experimental and random errors

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002 6191

Figure 16. Growth rate vs supersaturation level.

Figure 17. Growth rate vs nucleation rate.

and the simulation assumptions such as the absence of breakage and agglomeration. Accurate estimation of the supersaturation in crystallization processes is difficult and at the same time vital for the correct prediction of other process variables. A 0.6% error in the prediction of the concentration is almost equal to the maximum possible value of the supersaturation and therefore introduces a large error in the simulation results. The above model proved to be accurate enough for the reasonable prediction of the supersaturation. It is noted that the accuracy of the model and the instrument in the estimation and the measurement of the CSD is critical with class II crystallizers in which the supersaturation level is very low. From Figure 10-c it is noticed that the second-order size dependence in the model slightly overestimates the

volume percent of the large crystals. When the model was tested with a size-dependent order equal to 1.85 for crystals larger than 700 µm, better results were obtained. The CSD, concentration, and supersaturation plots for these kinetics are shown in Figures 13-15 for this modified kinetic. This suggests that a phenomenon such as breakage by the stirrer slightly decreases the growth of larger crystals. Figures 16 and 17 compare various kinetic correlations reported in the literature. In Figure 16, crystal growth vs supersaturation is sketched. Lower growth rates are given from the work of Tavare,5,17 who reports the supersaturation order as 0.91 in the growth kinetics. The highest growth rate is reported by Mullin et al. 4 for {001} and {100} faces. These investigators concluded that growth orders of 1 and 2 for {001} and {100} faces best fit the experimental

6192

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002

data. Growth rate predictions in the present work for crystal sizes of 0, 300, and 500 µm are somewhere between Tavare’s and Mullin et al.’s growth rate values. On the basis of these graphs, crystals with a size of 300 µm grow almost twice as fast as newly generated crystals when the solution is maintained at a supersaturation of 0.006 kg/kg of solvent. It is also noticed that the size-dependent growth kinetic varies almost linearly with supersaturation. Figure 17 shows a nucleation vs growth rate developed in this work and those found in the literature. Tavare5 reports the lowest nucleation rate, while the highest is due to Youngquist and Randolph16 and Jager et al.3 The nucleation rate from this work is higher than those of other investigators at low growth rates, but its gradual increase with the growth rate locates it somewhere between the other correlations at higher growth rates. It is noted that Youngquist and Randolph16 and Jager et al.3 report their results for continuous crystallizers, and this may somehow contribute to the similar trend observed in their work. In the development of the chart for the present work, the size-dependent growth rate for nuclei was used. The validity of the various kinetics shown in Figures 16 and 17 may be tested by employing the kinetic correlations in the simulation program and comparing them with the experimental results. Unfortunately, the work by Mullin et al.4 could not be used in our simulation program because the growth model predicts different growth rates for two major faces of the crystals, while the simulation used in this work assumes a spherical crystal shape with a shape factor to correct for the deviation from sphericity. Kinetic relations presented by Tavare,17 Youngquist and Randolph,16 and Jager et al.3 give only one of two correlations and therefore could not be tested in the simulation work. Tavare’s work5 with both nucleation and growth kinetic correlations was tested in the simulation program. The comparison of the results indicated that his model prediction for large crystals is low, and as a result, large errors for CSD and supersaturation are obtained. Conclusions Crystallization kinetics of ammonium sulfate were studied using a cooling batch seeded crystallizer. To accomplish this, the densities of near-saturated solutions of ammonium sulfate at different concentrations and temperatures were determined. A model was developed to fit the data. The solubility of the crystals was determined using this model and the experimental data. Crystallization kinetics were estimated from a nonlinear optimization program and a process simulation program. Supersaturation and CSD were employed as the major process outputs in the optimization objective function. The PDE in the simulation program was solved using a combination of the Lax-Wendroff and Crank-Nicholson methods. Both size-independent and -dependent growth models were developed for ammonium sulfate. The results indicated that a secondorder growth kinetic in the process model reasonably predicts the behavior of the process, while a sizeindependent growth model used by other workers is clearly inadequate. Fitting the experimental data with a size order of 1.85 for large crystals resulted in the lowest error between the measured and simulated outputs.

Acknowledgment The authors are grateful for the financial support provided by the Natural Sciences and Engineering Council (NSERC) of Canada and the Ontario Research and Development Challenge Fund (ORDCF). Nomenclature AN ) specific crystal surface area (m2/kg of solvent) B ) nucleation rate (#/kg of solvent/s) B0 ) nucleation rate at size zero (#/kg of solvent/s) C ) solution concentration (kg/kg of solvent) C* ) saturation concentration (kg/kg of solvent) ∆C ) supersaturation (kg/kg of solvent) Eb ) activation energy of birth (J/mol of crystal) Eg ) activation energy of growth (J/mol of crystal) G ) growth rate (m/s) g ) supersaturation exponent in eqs 4 and 10 G|L)0 ) growth rate evaluated at size zero (m/s) G ¨ l ) growth rate derivative with respect to size in eqs 1214 i ) growth rate exponent in eqs 5 and 11 j ) magma density exponent in eqs 5 and 11 k1 ) parameter in eq 15 k2 ) parameter in eq 15 ka ) area shape factor kb ) nucleation rate coefficient (#/kg of solvent/s/(m/s)i(kg/ kg of solvent)j) kg ) growth rate coefficient (m/s/(kg/kg of solvent)g) kv ) volume shape factor l ) index for size domain L ) size of crystals (m) ∆L ) computational size interval (m) LS ) seed mean size (m) LS0 ) seed mean size (m) at initial time M ) magma density (kg/kg of solvent) m0 ) zeroth moment of the size distribution mj ) jth moment of the size distribution defined in eq 1 n ) population density of crystals (#/kg of solvent/m) n0 ) nuclei population density (#/kg of solvent/m) Ne ) number of experiments Nm ) number of measured variables Ns ) number of samples taken for each measured variable in each run r ) number of dimensions in the definition of a simplex R ) gas constant (J/mol‚K) S ) mass of the solvent (kg) T ) solution temperature (K) t ) time (s) ∆t ) computational time interval (s) W ) weight of the crystals (kg) wj ) weight factor for the jth measured variable WS0 ) weight of the seed crystals at initial time (kg) y ) variable in the objective function Greek Letters R1 ) parameter in eq 15 R2 ) parameter in eq 15 ∆ ) difference θ ) optimization parameter vector λ1, λ2 ) variables in eq 12 λ3, λ4, λ5 ) variables in eq 13 λ6, λ7, λ8 ) variables in eq 14 F ) density of the solution (kg/m3) Fc ) density of the crystals (kg/m3) τ ) crystallizer time constant (s) Φ ) objective function Subscripts and Superscripts c ) crystals

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002 6193 i ) index for the number of experiments j ) index for the number of measured variables j ) index for moment j in eqs 1 and 3 k ) index for the number of samples m ) measured variable N ) nuclei, new crystals p ) index for the time domain q ) magma density exponent S ) seed crystals

Literature Cited (1) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes, 2nd ed.; Academic Press: New York, 1988. (2) Tavare, N. S.; Garside, J. Simultaneous Estimation of Crystal Nucleation and Growth Kinetics from Batch Experiments. Chem. Eng. Res. Des. 1986, 64, 109-118. (3) Jager, J.; De Wolf, S.; Kramer, H. J. M.; De Jong, E. J. Estimation of Nucleation Kinetics from Crystal Size Distribution Transients of a Continuous Crystallizer. Chem. Eng. Sci. 1991, 46 (3), 807-818. (4) Mullin, J. W.; Chakramborty, M.; Mehra, M. Nucleation and Growth of Ammonium Sulphate Crystals from Aqueous Solutions. J. Appl. Chem. 1970, 20, 367-371. (5) Tavare, N. S. Ammonium Sulfate Crystallization in a Cooling Batch Crystallizer. Sep. Sci. Technol. 1992, 27 (11), 14691487. (6) Garside, J.; Gibilaro, L. G.; Tavare, N. S. Evaluation of Crystal Growth Kinetics from Desupersaturation Curve Using Initial Derivatives. Chem. Eng. Sci. 1982, 37 (11), 1625-1628. (7) Jones, A. G.; Mullin, J. W. Programmed Cooling Crystallization of Potassium Sulphate Solutions. Chem. Eng. Sci. 1974, 29, 105-118.

(8) Roache, P. J. Computational Fluid Dynamics; Hermosa Publishers: Albuqurque, NM, 1982. (9) Ames, W. F. Numerical Methods for Partial Differential Equations; Academic Press: Boston, MA, 1992. (10) Bennett, M. K.; Rohani, S. Solution of Population Balance Equations with a New Combined Lax-Wendroff/Crank-Nicholson Method. Chem. Eng. Sci. 2001, 56, 6623-6633. (11) Miller, S. M.; Rawlings, J. B. Model Identification and Control Strategies for Batch Cooling Crystallizers. AIChE J. 1994, 40 (8), 1312-1327. (12) Tadayyon, A. Modelling, Control and Measurement in Continuous KCl Crystallizers. Ph.D. Thesis, University of Saskatchewan, Canada, 2000. (13) Bennett, M. K. Software for the Simulation of Industrial Crystallizers. M.E.Sc. Thesis, The University of Western Ontario, Canada, 2001. (14) Chapra, S. C.; Canale, R. P. Numerical Methods for Engineers; McGraw-Hill: New York, 1999. (15) Nelder, J. A.; Mead, R. A. Simplex Method for Function Minimization. Comput. J. 1965, 7, 308-313. (16) Youngquist, G. R.; Randolph, A. D. Secondary Nucleation in a Class II system, Ammonium Sulphate-Water. AIChE J. 1972, 18 (2), 421-429. (17) Tavare, N. S. Growth Kinetics of Ammonium Sulfate in a Batch Cooling Crystallizer Using Initial Derivatives. AIChE J. 1985, 31 (10), 1733-1735.

Received for review February 27, 2002 Revised manuscript received July 22, 2002 Accepted August 9, 2002 IE020163R