DOI: 10.1021/cg100122y
Modeling the Crystallization Kinetic Rates of Lactose via Artificial Neural Network
2010, Vol. 10 2620–2628
Shin Yee Wong,† Rajesh K. Bund,‡ Robin K. Connelly,†,‡,§ and Richard W. Hartel*,†,‡ †
Department of Biological Systems Engineering, University of Wisconsin, Madison, 460 Henry Mall, Madison, Wisconsin 53706, and ‡Department of Food Science, University of Wisconsin, Madison, 1605 Linden Drive Madison, Wisconsin 53706. §Current affiliation: Solae LLC, 4300 Duncan Avenue, St. Louis, Missouri 63110. Received January 26, 2010; Revised Manuscript Received April 22, 2010
ABSTRACT: A mixed-suspension, mixed-product removal (MSMPR) draft-tube-baffle crystallizer was used to estimate lactose crystallization kinetic rates under different processing conditions. The processing parameters investigated were the degree of supersaturation (S), agitation speed (rpm), the amount of seed (seed), and residence time (τ). Crystallization kinetic rates, including nucleation, growth, and aggregation rates, were calculated from the experimental crystal size distributions (CSD) via a moment transformation method. Two approaches were used to develop a correlation between the processing conditions and kinetic rates. The conventional power law approach estimated the model parameters of an assumed functional form via regression analysis. An artificial neural network (ANN) approach was also used. ANN does not require an assumed function, it develops a function by learning from the experimental processing conditions and kinetic pairs. The results showed that the ANN approach gave a much better prediction of the kinetic rates. Compared to the power law model, the R2 and rootmean-square error (RMSE) of the ANN approach improved by at least 49% and 33%, respectively. 1. Introduction The fundamental knowledge of lactose crystallization kinetics is beneficial in the analysis, design, and operation of industrial processes. Specifically, the crystal size of the final product depends upon the control of the rate of nucleation (creation of new crystal), growth, and aggregation (intergrowth of crystal). Numerous studies on the crystallization of alpha-lactose monohydrate have been reported with different levels of complexity. For example, one model incorporated the effect of growth rate dispersion (GRD) for a continuous cooling crystallization process,1 while another used one single kinetic model that included mutarotation, nucleation, and crystal growth for an unseeded batch crystallization process.2 Kinetic expressions for nucleation and/or growth have generally been reported; however, the development of kinetic expressions for aggregation is lacking. In addition, most kinetic rates expressions reported2-5 are obtained from experimental data using an assumed functional form. However, the mechanisms of all crystallization kinetics are extremely complex, and the high nonlinearity and parameter interactions might not be easily reflected by the assumed correlations.6 Modeling by artificial neural networks (ANN)7 provides a new alternative to the conventional approach. ANNs are self-organized,7 meaning they do not require any predefined function between the inputs and outputs. The relationship between inputs and outputs is determined by learning from the input and output pairs. Recently, there has been an increased application of ANN for reaction process systems.6-10 Among these applications, most of the applications in food systems focus on the growth kinetics of sucrose crystals in batch processes.8,9 Besides crystallization kinetic rates, ANN has also been used in other aspects of sucrose crystallization. For example, it was used to *Corresponding author. Phone: (608) 263-1965. Fax: (608) 262-6872. E-mail:
[email protected]. pubs.acs.org/crystal
Published on Web 05/10/2010
predict the mean crystal size from crystal images,10 to estimate the boiling point elevation of a sucrose-water system,7 and to function as a process model for model predictive control.7 In this study, the potential use of ANN to model a lactose crystallization process was investigated. First, a mixed-suspension, mixed-product removal (MSMPR) crystallizer was used to estimate the crystallization kinetic rates of a lactose cooling crystallization process experimentally. Then, the complex relationship between the process parameters and crystallization kinetic rates were modeled both by the conventional power law relationship and by the feed-forward ANN architecture. 2. Methods 2.1. Materials and Methods. A mixed-suspension, mixed-product removal (MSMPR)11 crystallizer was used for the continuous crystallization process. The schematic diagram of the crystallizer is shown in Figure 1. The crystallizer consisted of a tall-form 250 mL jacketed beaker, a polycarbonate draft tube-baffle set, and a 1” diameter three blade propeller attached to a stainless steel shaft. The draft tube baffle system was constructed so that both mother liquor and suspension are well mixed. This system allowed axial flow (of solution) inside the crystallizer: that is, flow is down inside the draft tube, and up in the annular region. In addition, the slower flow in the annulus reduces the tendency of large crystals to segregate along the outer wall. The temperature of the crystallizer was controlled by a circulating water bath (Thermo Electron Neslab RTE7) connected to the jacketed beaker. Supersaturated lactose solutions were prepared by dissolving pharmaceutical grade lactose (310, Foremost Farms USA) in deionized water. The solutions were heated to 80-85 C and vacuum filtered to remove any particles (6 μm and above). During the experiments, the supersaturated solution at the desired temperature (5 C higher than crystallizer’s temperature) was pumped into the crystallizer at a specific volumetric flow rate to a fixed location near the impeller. The solution was seeded (at varying amount) with the same pharmaceutical lactose (D[4,3] ≈ 83 μm) at time t = 0, which is at the beginning of the experiment. The lactose slurry, which developed in the crystallizer according to the crystallization kinetics at the experimental conditions, was r 2010 American Chemical Society
Article
Crystal Growth & Design, Vol. 10, No. 6, 2010
continually removed from the crystallizer at the same volumetric flow rate as the inlet. The refractive index (RI) was monitored by a Bausch & Lomb refractometer. After approximately 8-10 residence times, when there was no change in the RI, steady state was reached. The crystallizer was emptied and the lactose crystals were then collected by vacuum filtration. The crystal size distribution (CSD) was analyzed by Malvern Master Sizer 2000. The analysis method was modified from the method reported by Adi et al.12 First, approximately 0.08 g of lactose crystal was dispersed in 4 mL of 2-propanol (HPLC grade) with the aid of sonication (Fisher Scientific, FS30 ultrasonic bath) in a water bath for 3 min at room temperature. The sonicated sample was added dropwise into the sample cell (of Master Sizer 2000) containing 250 mL (approx) of 2-propanol until an obscuration of 4-7% was obtained. Particle size measurement of each sample was performed in the range of 6-1000 μm and analyzed with the reference refractive index of lactose (1.633) and 2-propanol (1.378). To establish the crystallization kinetic rates model, a series of experiments were conducted based on varying supersaturation (S), agitation speed (rpm), the amount of seed (seed), and residence time (τ). The residence time (τ) was defined as the ratio of the total volume (206 mL) and volumetric flow rate. The superaturation (S) of the R-lactose solution was calculated by a model developed by Visser,13 as shown in eqs 1-4. C ð1Þ S ¼ C s - FK m ðC - C s Þ C s ¼ 12:23 þ 0:3375T þ 0:001236T 2 þ 0:00007257T 3 þ 5:188 10 - 7 T 4
ð2Þ
Figure 1. Schematic diagram of the experimental setup: The crystallizer consisted of a baffle-draft-tube setup that allowed maximum mixing along the axial (z) direction.
2621
- 2374:6 þ 4:5683 F ¼ exp Tk
ð3Þ
K m ¼ - 0:002286T þ 2:6371
ð4Þ
The experimental conditions were selected based on minimum encrustation and sufficient yield for CSD measurement. Encrustation is defined as the coating of lactose crystals on the surface of the draft-tube/baffle or crystallizer wall. When it occurred, severe encrustation blocked the flow channel completely. To measure the level of encrustation, the emptied crystallizer (after vacuum filtration of the slurry) was dried in a 60 C oven overnight, and the weight of encrusted lactose crystals was measured. Encrustation levels less than 7.5 g were considered acceptable for crystallizer operation. In addition, a minimum of 0.08 g of lactose crystal was required for CSD analysis by Malvern Mastersizer. On the basis of these operating constraints, the range of operating conditions was S = [1.8 -1.97], seed (g) = [1.5 - 4.5], τ (min) = [21 - 28.6], and rpm = [837.2 - 2511.6] (tip speed (m/s) = [1.11 3.34]). Thirteen conditions were selected with each experiment repeated at least twice, and a total of 27 experiments were conducted. 2.2. Calculation of Crystallization Kinetic Rates from MSMPR Data. The CSD results obtained experimentally were used to calculate the crystallization kinetic rates via the approach14 summarized in Figure 2. First, the CSD data were transformed into the volume-based population density via eq 5 (see Figure 2). The first three moments (eq 6) were calculated through integration via the trapezoidal method. The growth (Gv), nucleation (Bo), and aggregation (βo) rates were calculated using eqs 7-9 in Figure 2, respectively. Equations 7-9 were derived from the population balance equation (PBE), with details on the derivation found in the Appendix. Unlike most reported studies on lactose crystallization kinetics,2-5 aggregation was also considered in this study. As shown in Figure 3, the microscopic photo of the experimental lactose crystals showed a very different crystal morphology between the samples with high and low aggregation rates. Figure 3a is an indication of the occurrence of aggregation during the crystallization process. Therefore, aggregation kinetic rates were accounted for when solving the PBE, resulting in eqs 7-9 for the crystallization kinetic rates. 2.3. Regression Analysis. Many empirical models describing the rate of secondary nucleation (Bo) and growth (G) have been proposed. In most of the models, nucleation rate was expressed as a power law function of supersaturation (S) and suspension density (MT, m3).15-17 Growth rate is commonly expressed as a function of supersaturation.15,17,18 The influence of temperature was expressed via Arrhenius relationship or direct function if the solubility curve is known. In this study, these models were modified or simplified to improve the goodness of fit of the assumed model. First, the temperature exponential term was removed, because all experiments were conducted at a fixed temperature (40 C). Second, the exponent of MT on nucleation rate was set to 1, which assumes the dominance of collisions between crystals and the vessel walls, or
Figure 2. Postprocessing MSMPR data: The experimental CSD data was transformed into crystallization kinetic rates via the moment transformation method.14
2622
Crystal Growth & Design, Vol. 10, No. 6, 2010
Wong et al.
more probably, between crystals and the stirrer over those between any two random crystals.16 In addition, to model the relation of the kinetic rates with respect to all four processing parameters (S, seed, rpm, τ), the kinetic rates were expressed as power law functions of all parameters (eq 10-12). Gv ¼ k 1 τk 2 S k 3 rpmk 4 seedk 5
ð10Þ
Bo ¼ k 6 τk 7 Sk 8 rpmk 9 seedk 10 MT
ð11Þ
βo ¼ k 11 τk 12 S k 13 rpmk 14 seedk 15
ð12Þ
Equations 10-12 were fitted to the experimental values by nonlinear least-squares curve fitting analysis using Matlab R2008a. Each experimental runs were considered as a single data point, and a total of 27*3 data points were used to estimate 15 model parameters in eqs 10-12. The model parameters, k1 to k15 were obtained when the sum of squares of the difference between predicted (pv) and experimental value (ev) was minimized (i.e., min Σi(pvi - evi)2). 2.4. Artificial Neural Network (ANN). An ANN is made up of three layers, including an input layer, a hidden layer, and an output layer. A neural network with three nodes in the hidden layer is shown in Figure 4. Each layer is connected by nodes. The architecture of a node is illustrated in Figure 5. In each node, the inputs are weighted and summed (eqs 13 and 15). The sum is then processed by a transfer function (eqs 14 and 16) before the output is passed to the next node or an output node. A tan-sigmoid transfer function (eq 14) was used in the hidden layer, and a linear transfer function (eq 16) was used in the output layer. A list of commonly used transfer functions can be found in Peacock.7 The choice of transfer function is dependent on the characteristics of the raw data to be fitted. If the transfer function closely approximates the data to be fitted, then fewer neurons will be required to achieve the desired accuracy. However, in this study, the raw data followed an arbitrary function. It might be best described by a nonlinear, continuous transfer function, but it is difficult to know what transfer function would be optimal. In this study, the default (tan-sigmoid) transfer function in the Neural Network Toolbox Fitting Tool in Matlab R2008a was used. This structure is useful for the function approximation (or regression) problem, and it was also used in other ANN kinetic studies.6 N h ¼ W h xn þ bh ð13Þ f ðN h Þ ¼
Figure 3. Microscopic photo for crystals collected from runs with (a) high aggregation rate (Run #5 in Table 1); (b) low aggregation rate (Run #7 in Table 1).
2 -1 1 þ e - 2N h
ð14Þ
N o ¼ W o f ðN h Þ þ bo
ð15Þ
f ðN o Þ ¼ N o
ð16Þ
Three ANNs, one each for growth, nucleation, and aggregation rates, were developed using Matlab R2008a. Rotational speed (rpm), amount of seed, supersaturation (S), and residence time (τ) were included as the input parameters. The structure of ANN can be summarized by eqs 13-16. The input vector (xn) was weighted and transformed to f(Nh) in the hidden layer, via eqs 13-14. This information was then transformed to the output layer, after which, the final output, f(No), was calculated by eqs 15-16. All weights (W) and biases (b) in the network were obtained by training the network with experimental data. The 27 experimental data points were split into three categories (training, validation, and testing) randomly by the ratio of 70%, 15%, and 15%. Training an ANN is accomplished by adjusting the weight and bias value
Table 1. Average Crystallization Kinetic Rates Obtained Experimentally Run #
rpm
seed (g)
S
τ (min)
Gv 1019 (m3 /s)
Bo 10-13 (#/m3-s)
βo 1012 (m3/#-s)
no. of experiments
1 2 3 4 5 6 7 8 9 10 11 12 13
1255.8 1255.8 2093 2093 2093 2093 837.2 1255.8 2511.6 1674.4 1674.4 1674.4 1674.4
2 3.5 3.5 2 2 3.5 2.5 2.5 2.5 1.5 4.5 2.5 2.5
1.817 1.812 1.812 1.809 1.952 1.941 1.919 1.913 1.873 1.943 1.915 1.913 1.904
21.21 25.19 21.75 26.55 21.49 23.45 22.93 21.93 22.56 21.60 22.65 25.19 26.67
8.74 ( 0.12 0.03 ( 0.0009 7.57 ( 0.7 0.03 ( 0.0018 2.35 ( 0.46 0.04 ( 0.003 5.88 ( 0.15 6.35 ( 0.21 0.92 ( 0.001 7.52 ( 1.36 0.03 ( 0.006 0.98 ( 0.051 6.56 ( 0.69
0.02 ( 0.005 432 ( 58 0.02 ( 0.0001 1705 ( 455 1.68 ( 0.07 3423.8 ( 342 0.01 ( 0.0006 0.03 ( 0.0018 1.82 ( 0.06 0.03 ( 0.01 1333.7 ( 429 1.44 ( 0.038 0.05 ( 0.0005
5.72 ( 0.36 4.13 ( 2.89 3.02 ( 0.56 1.98 ( 0.49 11.69 ( 4.02 10.21 ( 0.48 0.29 ( 0.014 1.34 ( 0.16 2.48 ( 0.10 4.83 ( 3.17 4.18 ( 1.86 4.87 ( 0.37 1.33 ( 0.25
2 2 2 2 2 2 2 2 2 2 3 2 2
Article
Crystal Growth & Design, Vol. 10, No. 6, 2010
2623
Figure 4. Schematic of the neural network (each node is represented as a circle).
Figure 5. Schematic of each neural network node. according to Levenberg-Marquardt optimization.9 At each iteration, the weights are modified based on the error in the training set. In addition, the error on the validation set is also calculated to crossvalidate the new set of weights and to stop training before overfitting. The testing data are used as a completely independent test of network performance. Finally, the performance of the network was measured by the root-mean-square error (RMSE) and correlation coefficient (R2) between all predicted (pv) and experimental (ev) data. sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 1X RMSE ¼ ðpv - evÞ2 ð17Þ n i¼1 P ðpv - pvÞðev - evÞ R2 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P P ðpv - pvÞ2 ðev - evÞ2
ð18Þ
3. Results and Discussion 3.1. Crystallization Rates. Two typical crystal size distribution (CSD) profiles obtained experimentally are shown in Figure 6. The distribution pattern of the CSD depended strongly on the aggregation rate. At high aggregation rate (Figures 3a and 6a), an equally distributed bimodal peak was observed, whereas, a single peak dominated in runs with low aggregation rates (Figures 3b and 6b). The PBE model (Figure 2) was used to generate nucleation, growth, and aggregation kinetic rates for each experimental conditions. The average kinetic rates (and their standard deviation) for each experimental condition are shown in Table 1. To determine which experimental operating parameter and their interactions are statistically important, the data (in Table 1) were fitted to a general linear model using Minitab 15. Then, the parameters were ranked according to the p-value of each term, with the five most important operating parameter or parameter interactions summarized in Table 2. For lactose crystallization, the growth rate is usually surface integration controlled.5 From Table 2, the amount of seed (seed) and its interactive effect with the three other processing parameters were identified as being important
parameters. This observation is similar to those reported in the literature. In a recent patented process for lactose crystallization, the authors stressed the importance of controlling the initial number of nuclei in addition to other processing parameters.19 Others20 argue that when large additional numbers of crystal units are introduced by seeding, the available S is divided equally between crystals resulting in less growth for individual crystals. When operating in the metastable zone, a highly supersaturated solution is more prone to secondary nucleation compared to a more dilute solution.21 Secondary nucleation is the dominant nucleation mechanism in a MSMPR crystallizer.22 Contact nucleation was reported as the most common mechanism of nucleation,17 which is induced by crystal contacts with other crystals or the crystallizer parts. Therefore, as shown in Table 2, the parameters τ, seed and rpm, which either increase the number of particles in suspension or the frequency of crystal contacts, were identified as significant parameters. During bulk crystallization, aggregation or intergrowth of crystals might occur via formation of nucleus-bridges between the crystals.23 This may be induced by Brownian motion, laminar shear flow, turbulent flow, or a combination of these mechanisms.24 As shown in Table 2, the rate of aggregation was affected mainly by rpm, τ, S, and their interactive effects, all of which might have promoted the motion of crystals in the crystallizer. 3.2. Crystallization Kinetic Rates. The crystallization kinetic rates (Gv, Bo, βo) were modeled using either the power law model (eqs 10-12) or ANN. Equations 19-21 are the kinetic equations obtained from the power law regression analysis. The relative value of the model parameters indicates the relative effect of each parameter on the crystallization kinetic rates. For example, according to this model, the growth rate (Gv) (eq 19) increases strongly with S, while changing the level of rpm, seed, and τ has minimum impact. For nucleation (eq 20), each of the operating parameters is predicted to have a large effect. From eq 21, the model shows that τ, rpm, and S have a large effect on aggregation, while the seed level does not. Gv ¼ 8:3 10 - 25 τ1:510
- 12
S 18:7 rpm2:310
- 14
seed2:210
- 14
ð19Þ
Bo ¼ 1:3 1010 τ9:6 S 2:8 rpm6:4 seed6 MT
ð20Þ
βo ¼ 1:2 10 - 11 τ - 4:2 S - 1 rpm1:4 seed - 0:09
ð21Þ
2624
Crystal Growth & Design, Vol. 10, No. 6, 2010
Wong et al.
Figure 6. Crystal size distribution (CSD) of runs with (a) high aggregation rate (Run #5 in Table 1); (b) low aggregation rate (Run #7 in Table 1). Table 2. Five Most Important Operating Parameters on the Kinetic Rates (*Denotes Interactive Effect) Gv (p-value)
Bo (p-value)
βo (p-value)
seed*S (0.068) seed (0.1) seed*τ (0.102) S (0.161) rpm*seed (0.171)
S (0.000) S*τ (0.000) seed (0.001) τ (0.001) rpm*seed (0.002)
rpm*S (0.077) rpm (0.106) S*τ (0.179) τ (0.223) S (0.238)
Table 3. Statistical Parameters for Model Parameters in eqs 19-21
eq 19
eq 20
The 95% confidence intervals of all model parameters are summarized in Table 3. The relatively large confidence intervals are indicative of the poor fit for the power law model. In addition, the range of the model parameters estimated were very different from values commonly reported in the literature. For example, the exponent of S for growth (eq 19) is very high compared to the typical value of 2-3 in the literature.5 From the large confidence intervals, it is clear that the power law regression fits do not match the data very well. This is probably because the crystallization kinetic rates do not follow the assumed functional form as suggested in eqs 10-12. Figure 7 shows the plot of predicted growth rate (Gv) versus the experimental value. Data estimated by the ANN lies close to the dotted line (y = x), which indicates a small difference between the predicted and experimental value. This is also illustrated by the relatively small RMSE (Table 4). Prediction from the power law model has significantly higher RMSE and lower R2 (Table 4). Thus, the power law model did not predict the experimental growth rate values very well, whereas the ANN gave a good fit. A similar plot for nucleation rate (Bo/MT) is shown in Figure 8. Secondary nucleation is highly complex, and the complete mechanisms are still not well understood. In the literature, modeling for secondary nucleation rate has met with only limited success.5,17 As illustrated in Figure 8 and Table 4, the estimation from power-law model did not match the experimental results, yielding an extremely high RMSE. However, model fitting by ANN was able to capture the relationship between input parameters and Bo, resulting in a low RMSE and high R2. Because of the complex nature of the phenomenon, the aggregation rate was difficult to model. However, similar to
eq 21
parameters
estimated value
standard error
95% confidence bounds
k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15
8.34 10-25 1.54 10-12 18.7 2.28 10-14 2.22 10-14 1.29 1010 9.6 2.8 6.4 6.0 1.20 10-11 -4.2 -1.0 1.4 -0.1
1.02 10-15 5.3 16.8 1.6 1.6 1.18 1024 8.5 27.4 2.6 2.7 9.31 10-11 2.1 6.7 0.6 0.6
(2.04 10-15 (10.7 (33.7 (3.2 (3.1 (2.35 1024 (17.0 (54.8 (5.2 (5.4 (1.86 10-10 (4.2 (13.4 (1.3 (1.3
Figure 7. Plot of growth rate (Gv) predicted from ANN and power law (eq 10) (for perfect prediction; all points should lie along the dotted line of y = x).
Gv and Bo, the prediction from ANN is significantly better than that obtained from power law (Figure 9 and Table 4). From the regression analysis, it was observed that all crystallization kinetic rates were weakly correlated to the processing parameters by the power law function and large errors were observed in eqs 10-12. Alternative functional
Article
Crystal Growth & Design, Vol. 10, No. 6, 2010
2625
Table 4. R2 and Root Mean Square Error (RMSE) for All Modelsa R2
growth nucleation aggregation a
RMSE
ANN
power law
ANN
power law
0.99 1.00 0.87
0.42 0.76 0.31
5.09 51.82 1.87
50.74 1108.34 3.72
ANN (artificial neural network).
Figure 9. Plot of aggregation rate (βo) predicted from ANN and power law (eq 12) (for perfect prediction, all points should lie along the dotted line of y = x).
Figure 8. Plot of nucleation rate (Bo) predicted from ANN and power law (eq 11) (for perfect prediction, all points should lie along the dotted line of y = x).
forms may improve the fit; however, due to the interactive nature of various processing parameters, the generalization of an expression correlating all the involved operating parameters is currently not possible.9 In contrast, ANNs are good at fitting the kinetic rates and recognizing patterns. Compared to regression analysis, the ANN approach can be easily implemented in Matlab R2008a, and it requires less formal statistical training than regression analysis. However, the network development is empirical in nature. It is highly dependent on the experimental data, and any extrapolation outside the training range should be avoided. 3.3. Application of ANN Kinetic Models. To demonstrate a potential application of the ANNs developed in this study, the ANNs were used to simulate the change in crystallization kinetic rates for different operating conditions. The predicted crystallization kinetic rates corresponding to the change in processing parameters are shown in Figure 10. These show a one-factor at a time simulation with the nominal conditions as rpm=1500 rpm, seed=2.5 g, S=1.85 and τ=22.5 min. In Figure 10, each curve shows the change in the crystallization kinetic corresponding to the changing level of one processing condition (x-axis), while holding the other three at nominal conditions. The level of the processing parameters investigated were equal to the range of operating conditions reported in Section 2.1; specifically, S=[1.8-1.97], rpm=[837.2 - 2511.6], seed (g) = [1.5 - 4.5], and τ (min) = [21 - 28.6]. For each parameters, 10 changing levels were tested. Like most kinetics studies, the ANNs predicted an exponential increase in growth rate with respect to the degree of supersaturation, as shown in Figure 10c. The effects of the other operating parameters, however, were not so clear, showing complex behavior depending on nucleation and aggregation. Figure 10a,b shows that the growth rate has a
quadratic dependence on rpm and seed, with Gv highest circa the mid range, and lowest in the high and low levels. The rate of convective mass transfer increases as the rpm increases, resulting in a higher Gv. However, above ca. 1500 rpm, the high rpm promoted aggregation, resulting in a decreased growth rate. At ca. 2300 rpm, an exponential increase in the nucleation rate was observed. As the number of secondary nuclei increased, increased growth rate dispersion would also occur, resulting in wider variation in the growth rates of very small crystals.11 Thus, the final product had a higher number of crystals in the small size range, resulting in a lower average growth rate. Similarly, the growth rate pattern observed in Figure 10b might also be due to the phenomenon of growth rate dispersion. In the case of residence time, Figure 10d, the growth rate decreased with increasing residence time until a point where the Gv was oscillating in a narrow range. The observed trend indicated that (at the nominal condition) the crystal growth was limited by the rate of mass transfer. When the volumetric flow rate in the feed and withdrawal tube decreased with an increasing level of residence time, it resulted in a decrease in the rate of mass transfer,21 thus lowering the Gv. The rate of secondary nucleation is very sensitive to the corresponding change in processing conditions. Secondary nuclei originated from the seed crystals.21 In general, a larger number (amount) of seeds is equivalent to a larger number of nucleation sites, resulting in an increased nucleation rate. However, this correlation is not always linear. For example, as shown in Figure 10b, the nucleation rate increases only after the amount of seed increases beyond 3.5 g. In addition, the nuclei are most likely generated by crystal contact with the crystallizer environment,11 for example, collisions with stirrers, impellers, and crystallizer walls or crystal-crystal contact. The frequency/probability of contact can be enhanced by agitation and τ, resulting in a higher Bo, as evidenced in Figure 10a,d). Supersaturation had the highest impact on nucleation rate with nucleation increasing exponentially when S was greater than about 1.85 (Figure 10c). At high supersaturation, the size of critical nucleus decreases, and thus the probability of the nuclei surviving to form crystal is also larger.21 According to Figure 10, the ANN model predicted an increasing aggregation rate with higher seed (b) and lower τ (d), while very high rpm (a) caused a decrease in aggregation and S had no effect (c). In the literature, there were limited reports that detail the mechanism of aggregation during
2626
Crystal Growth & Design, Vol. 10, No. 6, 2010
Wong et al.
Figure 11. The normalized average % change in crystallization kinetic rates with respect to all processing conditions.
the time during which the crystals are in contact and the probability that the nucleus bridges formed between crystals are not destroyed at subsequent collisions. From the results shown in Figure 10, at a fixed flow condition, aggregation increased slightly with an increasing amount of seed (increase in the total number of crystal in suspension). At high rpm, the high velocity of the crystals might promote the frequency (probability) of collisions, but at the same time, the contact time at collision might be too short for the formation of nucleus bridge, which results in a lower aggregation rate (Figure 10a). Similar to the impact of residence time on Gv, the flow rate in the crystallizer decreased with an increasing level of residence time, resulting in the decrease in convective flow, thus lowering the aggregation rate. To further clarify the impact of each parameter, the information in Figure 10 can be transformed to a chart where the % change in crystallization kinetic rates is normalized with respect to a 1% change in the processing conditions, as shown in Figure 11. The % change in processing parameters was calculated by ðlevel - lower limitspecific condition Þ %change ¼ lower limitspecific condition 100%
Figure 10. Changes in the crystallization kinetic rates corresponding to changes in the (a) agitation speed (rpm), (b) amount of seed (g), (c) supersaturation (S), and (d) residence time (τ) at nominal condition of rpm = 1500 rpm; seed = 2.5 g; S = 1.85; and τ = 22.5 min.
lactose crystallization. An explanation to the simulated trend can be speculated based on the aggregation studies for other crystals. Linnikov23 developed an equation to describe the aggregation of crystals by modifying the Smoluchowski equation. In the equation, aggregation was correlated to the total number of crystals in a solution, the hydrodynamics, the probability of collision of suspended crystals, the proper orientation of crystals at the moment of collisions,
ð22Þ
The total % change in the processing conditions were approximately 9.5%, 200%, 200%, and 37% for S, rpm, seed, and τ, respectively. At the nominal conditions investigated, the nucleation rate is most sensitive to the changing level of operating conditions. In addition, among all other operating conditions, Figure 11 indicated that nucleation rate can be best controlled by the degree of supersaturation. Both Figures 10 and 11 can be use to estimate a new set of processing conditions that improve the quality of lactose obtained from the nominal operating condition. For example, to produce lactose crystals with a larger average size, one would want to promote growth, and minimize nucleation and aggregation rate. From Figures 10 and 11, this can be achieved by operating at the optimum level of seed or rpm, which are 2.8 g or 1600 rpm (tip speed = 2.13 m/s), respectively. 4. Conclusions In a cooling lactose crystallization process, the quality of the lactose crystals is governed by the crystallization kinetic rates (nucleation, growth, and aggregation rates). In the context of process modeling and optimization, a good model that correlates the processing parameters to crystallization
Article
Crystal Growth & Design, Vol. 10, No. 6, 2010
kinetic rates is essential to the analysis, design, and operation of the crystallization process. However, the mechanisms of crystallization kinetics are complex and highly nonlinear. As indicated in the experimental results, the crystallization kinetic rates are influenced by the highly interactive effects of most processing parameters. Therefore, they cannot be easily modeled by an assumed function form (e.g., power law kinetics). The ability of ANN to predict the highly complex and interactive nature of crystallization kinetic rates results in a significantly better match with experimental data, as shown in this study. Acknowledgment. This project was supported by National Research Initiative Grant #2007-55503-18448 from the USDA Cooperative State Research Education and Extension Service. Derivation of Equations 7-9 To post-process the MSMPR data obtained from Mastersizer analysis, equations 7-9 were used to calculate the crystallization kinetic rates by taking the first three volume based moments. These equations were derived from the population balance equation;14 however, the details on derivation were not available in the reference. To justify the use of equations 7-9, their derivations are presented in this section. The population balance equation is given as11 Dn þ r 3 ðv ðA1Þ BnÞ - birth þ death ¼ 0 Dt where n is the population density, vB is the set of internal and external coordinates, where vB = vBe þ vBi birth ¼ birthagg þ birthbr
As shown in Randolph & Larson, eq A1 for a macroscopic external coordinate region (with volume V) having an arbitrary number of inputs and outputs of flow rate Qk and population density nk can be written as
ðA2Þ
For a continuous, constant volume, isothermal, wellmixed, MSMPR crystallizer, assuming (1) The feed stream is free of suspended solids; (2) the crystallization process can be described satisfactorily with one internal coordinate (v, particle volume) in particle phase space, eq A2 can be written as Dn D n þ ðGv nÞ - birth þ death þ ¼ 0 ðA3Þ Dt Dv τ The ith volume-based moment of the distribution was defined as Z ¥ μi ¼ vi nðvÞ dv ð6Þ 0
i
0
dμi μi þ þ dt τ
0
Z
Z ¥ D vi ðGv nÞ dv vi birthdv Dv 0 0 Z ¥ i þ v deathdv ¼ 0 ¥
Multiplying eq A3 by v dv and integrating from zero to infinity gives Z ¥ D n i Dn v þ ðGv nÞ - birth þ death þ dv ¼ 0 ðA4Þ Dt Dv τ 0
ðA5Þ
0
The third term can be integrated by parts to give Z ¥ Z ¥ D vi ðGv nÞ dv ¼ ½vi Gv n¥0 - i vi - 1 Gv n dv Dv 0 0
ðA6Þ
As v f ¥, n(v) f 0, and at v = 0, n(0)Gv = Bo (which is the nucleation rate). So eq A6 can be written as, Z ¥ D vi ðGv nÞ dv ¼ - 0i Bo - jGv μi - 1 ðA7Þ Dv 0 Substituting eqs A7 into A5, eq A5 becomes, dμi μi þ - 0i Bo - jGv μi - 1 - birthi þ deathi ¼ 0 dt τ where
Z birthi ¼
¥
ðA8Þ
vi birthðvÞ dv
ðA9Þ
vi deathðvÞ dv
ðA10Þ
0
deathi ¼
11
Dn dðlog vÞ þ r 3 ðv Bi nÞ þ n dt - birth þ death Dt X Qk nk þ ¼ 0 v k
Reversing the order of integration and differentiation for the first terms, the equation becomes, Z Z Z ¥ D ¥ i n ¥ i D v n dv þ v n dv þ vi ðGv nÞ dv Dt 0 τ 0 Dv 0 Z ¥ Z ¥ vi birthdv þ vi deathdv ¼ 0
Z
death ¼ deathagg þ deathbr
2627
¥
0
Assuming the particle generated from the breakage of particles is negligible, that is, in eq A1, birthbr ≈ 0 and deathbr ≈ 0. To account for aggregation, the birth and death function of the aggregation of two particles of volume u and v-u into a particle volume v can be represented by25 Z 1 ¥ birthagg ðvÞ ¼ Kðu, v - uÞnðu, tÞnðv - u, tÞ du ðA11Þ 2 0 Z deathagg ðvÞ ¼ nðv, tÞ
¥
Kðu, vÞnðu, tÞ du
ðA12Þ
0
The aggregation kernel, K(u,v), is a measure of the frequency with which a particle of volume u aggregates with one of volume v. The formulation of the aggregation kernel can be chosen to correspond to the particular mechanism of aggregation.25,26 Assuming the collisions rates (leading to aggregation) do not depend on the sizes of the colliding particles but only on the product of their number densities, n(v) and the rate (βo), the aggregation kernel K(u,v) can be expressed as βo. Then, eqs A9 and A10 can be written as, Z ¥ Z v β birthi ¼ o vi nðu, tÞnðv - u, tÞ du dv ðA13Þ 2 0 0 Z deathi ¼ βo
¥
Z i
0
¥
nðu, tÞ du dv
v nðv, tÞ 0
ðA14Þ
2628
Crystal Growth & Design, Vol. 10, No. 6, 2010
Wong et al.
Then applying eq 6 to eqs A13 and A14, eq A8 becomes27 dμi μi þ - 0i Bo - jGv μi - 1 dt τ
8 0
