Ind. Eng. Chem. Res. 2003, 42, 6949-6955
6949
Nonlinear Oscillations in Ammonium Sulfate Crystallization: A Comparison of Different Model Predictions Prasanna Kumar Pathath† and Achim Kienle*,†,‡ Max-Plank-Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, D-39106 Magdeburg, Germany, and Lehrstuhl fu¨ r Automatisierungstechnik und Modelbildung, Otto-von-Guericke-Universita¨ t Magdeburg, D-39106 Magdeburg, Germany
In this paper, the focus is on nonlinear oscillations of ammonium sulfate crystallization. This system has been studied intensively in the past and can be viewed as a standard test system for crystallization processes. The focus is on crystallizers with fines dissolution and classified product removal. Using methods from numerical bifurcation analysis, we investigate how stability in such a system depends on the various operational parameters, like cut sizes and recycle rates of fines dissolution and product classification. Theoretical predictions of three different population models are compared with each other: (i) a model with simple kinetics and ideal fines classification, (ii) a model with complex kinetics involving breakage due to attrition and ideal fines classification, and (iii) a model with simple kinetics and nonideal fines classification. It is found that predictions of the instability regions with models i and ii are rather close. In contrast, the nonideal fines classification can change the stability rather drastically. 1. Introduction Crystallization is used in a wide range of industries to produce solids and provide high-purity separations. Nonlinear oscillations in continuous industrial crystallization processes are a well-known phenomenon leading to practical difficulties in design and operation. In our previous study,1 a simple population model with simple kinetics, including ideal fines classification and classified product removal, was considered. To gain a deeper insight in the origin and parameter dependence of the process behavior, numerical bifurcation analysis was done. For a given physical system, emphasis was made on the fines dissolution rate, classified product removal rate, fines cut size, and classified product cut size. Because of the simplicity of the model, the results that were presented were only of a qualitative nature. In the present study, predictions of the simple model (model 1) are compared to those of more detailed models: (i) a model with complex kinetics involving breakage due to attrition and ideal fines classification (model 2) and (ii) a model with simple kinetics and nonideal fines classification (model 3). As an application, ammonium sulfate crystallization is considered. The overview of the paper is as follows. In the first section, the models that will be used for the study are briefly explained. The next section consists of a short description of complex kinetics such as the size-dependent growth rate, nucleation, and attrition used for model 2. In the Results section, initially the models that will be used for the study are compared with each other by comparing the kinetics. By means of numerical bifurcation analysis, the instability region is predicted for model 1 in the space of adjustable operating parameters * To whom correspondence should be addressed. Tel.: +49-391-6110369. Fax: +49-391-6110515. E-mail: kienle@ mpi-magdeburg.mpg.de. † Max-Plank-Institute for Dynamics of Complex Technical Systems. ‡ Otto-von-Guericke-Universita¨t Magdeburg.
such as the fines dissolution rate and classified product removal rate. For a finite number of operating points, dynamic simulations are carried out with model 2. The results are discussed and compared to the predictions of model 1 in terms of stability, evolution of mass density function (MDF), and mass median crystal size L50. L50 is defined by a characteristic length up to which half of the overall crystal mass is reached. Finally, the influence of nonideal fines classification as proposed by Mitrovic´2 is studied. 2. Description of the Models 2.1. Model 1. Model 1 for the crystallization process has the following assumptions: (i) ideal mixing; (ii) isothermal operation; (iii) constant overall volume (liquid + solid); (iv) nucleation of crystals at negligible size; (v) no particle breakage, attrition, or agglomeration. For the modeling, the crystallizer is divided into continuous and dispersed phases. The dispersed phase is modeled by the population balance as introduced by Randolph and Larson.3 2.1.1. Dispersed-Phase Balance Model. The population balance equation is obtained as
∂F(L,t) ∂[GF(L,t)] q )- [hf(L) + hp(L)]F(L,t) ∂t ∂L V
(1)
with the following initial and boundary conditions:
F(L,t)0) ) 0, F(L)0,t) )
B(w) G(w)
(2)
The quantity F(L,t) is the number density function. The classification functions specifying the fines dissolution and product classification1 are given by
hf(L) ) R1[1 - h(L - Lf)]
(3)
hp(L) ) 1 + R2h(L - Lp)
(4)
10.1021/ie030200u CCC: $25.00 © 2003 American Chemical Society Published on Web 11/20/2003
6950 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
Figure 1. Continuous crystallizer with fines dissolution and classified product removal according to Gerla.18
where h(L) is the unit step function and R1 and R2 are the recycle ratios of the fines dissolution and classified product removal. These are defined as R1 ) qf/q and R2 ) (qq - q)/q. For the definition of the flow rates qf, qq, and q, we refer to Figure 1. The growth and nucleation rates that are dependent on the solute concentration are given by the following empirical power laws:4
G(w,t) ) kg[w(t) - wsat]g
(5)
B(w,t) ) kbmt[w(t) - wsat]b
(6)
where kg, kb and g, b are the coefficients and exponents for the growth and nucleation rate expressions. wsat is the weight fraction of the solute at saturation at the operating temperature, and mt is the suspension density based on the third moment defined by
mt ) kvFs
∫LL F(L,t) L3 dL ∞
0
(7)
In our previous study,1 the focus was on crystallization of KCl, and a nucleation rate that only depends on the concentration difference was used. In the present study, we focus on ammonium sulfate crystallization, and a nucleation rate that depends on the suspension density and concentration difference is used. 2.1.2. Continuous-Phase Balance Model. The mass balance of the solute in the liquid phase is given by
FL
dw q(FS - FLw) FS - FLw d qFLwin ) + + dt V dt V qFS L∞ (1 + kv L [hp(L) - 1)F(L,t) L3 dL] (8) 0 V
∫
where V is the volume of the crystallizer, FL is the density of the liquid phase, w is the weight fraction of the solute, FS is the density of the crystals, and is the void fraction, which is given by
∫LL FL3 dL
) 1 - kv
∞
0
(9)
with kv as the volume shape factor of the crystals and L0 and L∞ as the minimum and maximum lengths of the crystals. For the numerical solution, the given population model 1 is converted into a system of differential algebraic equations with some suitable method of lines. For that, a simple finite-volume approximation on a nonequidistant but fixed grid is used. The nonequidistant grid is used here, because it provides more grid points at smaller crystal lengths where a higher resolution of the solution is required. A fixed grid is used in contrast to an adaptive grid to avoid any influence of the grid adaption on the stability of the solution. For details, the reader is referred to work by Pathath and Kienle.1 2.2. Model 2. Model 2 is taken directly from the contributions of Mitrovic´2 and Gerstlauer et al.5 For detailed derivation of the mathematical models and technical aspects, the reader is referred to work by Mitrovic´2 and Gerstlauer et al.5 Compared to model 1, additional population phenomena are considered like the birth of nuclei whose sizes are different from zero in model 2. The nucleation rate is considered as a boundary condition for model 1 compared to model 2, where the nucleation is considered as a source term. Furthermore, the size-dependent growth rate and attrition caused by crystal stirrer collisions are considered in model 2. A short overview regarding the population phenomena is given in the following. According to classical nucleation theory by Volmer6 dating back to 1939, molecules form stable and unstable unions and clusters. Depending on the size and energy of these molecular clusters, the stable ones are able to serve as a basis for the formation of new particles within a pure supersaturated liquid. This type of nucleation is only affected by a high concentration of solute in the liquid phase. In general, a pure homogeneous solution is rarely encountered in crystallization from solution. The solution inside a crystallizer is usually contaminated with foreign particles. These particles cause a heterogeneous nucleation and occur at a reduced supersaturation in comparison to homogeneous nucleation. A common assumption is that the same kinetic equations are used for both types of nucleation. One important property to know in nucleation is the size of the newly formed nuclei, the so-called critical crystal length. The nucleation rate used in this paper as well as in refs 2 and 5 was proposed by Mersmann et al.7 and is mainly dependent on the supersaturation and temperature inside the crystallizer. It also depends on natural constants such as Avogadro’s number and Boltzmann’s constant. Furthermore, it is also influenced by properties of a specific chemical system like the diffusion coefficient (which is also dependent on temperature), the saturation concentration of the solute, the heterogeneous nucleation factor, the molecule diameter of the solute, and the surface tension. The growth of crystals is a quite complex chemical and physical process. Depending on the chemical system and the actual conditions inside a crystallizer, different growth mechanisms are prevalent. For model 2, a growth rate expression accounting for diffusion- and integration-limited growth derived by Mersmann et al.7 is used. This growth rate is a function of crystal length L and is dependent on a number of additional factors such as the mass-transfer coefficient, integration coefficient, supersaturation, and crystal strain energy.
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 6951
According to work by Gahn,8 a modified supersaturation definition accounting for the individual strain energy of crystals is used to describe the growth rate of crystals. As a result, the driving force may become negative at high-strain-energy values. For negative growth rates (dissolution of crystals), only the mass-transfer-limited part is considered. Between the positive and negative growth regimes of the population, a specific crystal length with zero growth rate is identified. The crystals of this specific length do not grow at all. For more details pertaining to this phenomenon, the reader is referred to work by Mersmann et al.7 and Gahn.8-10 The next phenomenon to analyze is the breakage of particles. Because we are considering a draft tube baffle (DTB) crystallizer, the dominating process is attrition of crystals with the stirrer. In crystallization literature, the breakage or attrition of crystals is also referred to as secondary nucleation.11 The collision of crystals with the stirrer leads to generation of an abraded original crystal with a length somewhat smaller than the removed crystals and a number of small attrition fragments resulting from the crystal stirrer collision. Therefore, source and sink terms due to attrition are considered in model 2. For further details concerning the equations of this attrition phenomenon, the reader is referred to work by Mitrovic´2 and Gerstlauer et al.5 So, the main difference between models 1 and 2 is that model 2 has a detailed kinetics. For the numerical solution of model 2, an equidistant finite-volume discretization method with a flux limiter is used. Because of the attrition of crystals in model 2, the length of the crystals is confined to a smaller domain and, hence, an equidistant grid can be used instead of the nonequidistant grid used in model 1 to obtain solutions with similar accuracy. 2.3. Model 3. Model 3 is similar to model 1 except the fines classification function hf is replaced by nonideal fines classification hfd, which is explained in the following. The nonideal fines classification was proposed in work by Mitrovic´.2 It accounts for the dissolution of large crystals in the settling zone of the crystallizer due to increased temperature of the recycle from the fines dissolution. Therefore, a dissolution factor12 is added to the ideal fines classification function, which is given by
Figure 2. Comparison of growth rates for models 1 and 2 at steady state.
Figure 3. Instability region in the R1/R2 plane for the ammonium sulfate system predicted by model 1.
hfd(L) ) R1[1 - h(L - Lf) + kfdh(L - Lfd)] (10) where Lfd is the cut size of the large crystals that are being dissolved at the rate of kfd times the fines dissolution rate. The discretization scheme for obtaining the numerical solution of model 3 is the same as that applied for model 1. 3. Results 3.1. Comparison between Models 1 and 2. For the study of nonlinear oscillations, crystallization of an ammonium sulfate-water [(NH4)2SO4-H2O] system in a 1100 L DTB crystallizer is considered. The physical, chemical, and operational parameters of this system for both models are taken from work by Mitrovic´2 and Kind and Lieb.4 To compare models 1 and 2, as a first step, it is good to compare the kinetics of the two models and see how far they fair with each other. For comparison, the growth rates at steady state for both models are considered. Figure 2 is a plot of growth rate versus crystal length. From Figure 2, it is shown that the growth rate of model 1, which is size-independent, has
Figure 4. Comparison of instability regions predicted by model 1 with predictions by model 2.
a steady-state value of around 0.032 µm/s for all crystal lengths. For model 2, the curve shown in Figure 2 indicates how the growth rate of the crystals changes with length. Figure 2 shows that the growth rate is small for smaller crystal lengths and increases steadily as the length of the crystals increases and thereafter remains almost constant after some specific crystal length. When the growth rates of the two models are compared, it is shown that there is no significant difference between the two. From this, it can be concluded that the two models are reasonably fair enough to compare. To study the origin and parameter-dependent behavior of the oscillations, a numerical bifurcation analysis of model 1 is done using the DIVA simulation environment.13 DIVA provides continuation methods for the direct calculation of parameter-dependent steady state and periodic solution branches. During the calculation,
6952 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
Figure 5. Periodic behavior of MDF for model 1 for R1 ) 20 and R2 ) 4.
the stability of the corresponding solution is monitored and bifurcation points are detected and located. In the present case, oscillatory behavior appears/disappears at Hopf bifurcation points. The continuation can be restarted from these Hopf bifurcation points to trace out the Hopf point curves in a two-dimensional parameter plane. For the details of the numerical methods and the
underlying basic concepts, the reader is referred to work by Mangold et al.,14 Kienle et al.,15 and Kubicek and Marek.16 Emphasis in this contribution is on the stability of systems with fines dissolution and classified product removal. Therefore, a bifurcation diagram is obtained using the two-parameter continuation method available
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 6953
Figure 6. Periodic behavior of MDF for model 2 for R1 ) 20 and R2 ) 4.
in the DIVA simulation environment1,14 between the operating parameters R1 and R2 (Figure 3). The classified product size and fines cut size are Lp ) 800 µm and Lf ) 300 µm. The curve in Figure 3 represents the locus of Hopf bifurcation points in the R1/R2 plane, where the system changes from stable to unstable. Oscillatory behavior occurs for the parameter combinations within
the shaded region of Figure 3. Without product classification, i.e., R2 ) 0, the instability occurs for a fairly high reflux ratio of R1 ) 73, which is not shown in Figure 3 and which is an unrealisticly high value because industrial crystallizers are cycling at R1 values of 20-30 and below.4,17 Further, the figure reveals that, as the classified product removal rate is increased, the
6954 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
crystallizer tends to oscillate at low fines dissolution rates. Because the results are produced with model 1, they are expected to be only of a qualitative nature. In the following, they will be compared with the results from model 2. For that, a detailed crystallizer model having attrition, nucleation, and size-dependent growth rate is considered. Pointwise dynamic simulation is carried out for model 2 in the instability and stability regions of Figure 3, which was obtained by numerical bifurcation analysis of model 1. Several fines dissolution rates R1 and classified product removal rates R2 are chosen; the 2 points in Figure 4 represent simulation results of the system that are stable, and the b points represent simulation results showing sustained oscillations. From Figure 4, it is evident that the behavior of model 2 is in good agreement with the predicted behavior from model 1. To gain further insight in the dynamic behavior of the two models, dynamic simulation is done for fines dissolution rate R1 ) 20 and classified removal rate R2 ) 4. Figure 5 shows the simulation result of the mass density distribution for model 1 at different characteristic times during the cycle of oscillations. In this figure, time advances from left top to the bottom right in rows. In the beginning, there are a large number of nuclei due to high supersaturation represented by the first peak. A second peak represents some big crystals that are present in the crystallizer that are in the process being removed from the crystallizer. As time progresses, the crystals grow toward larger crystal lengths, and that is why the MDF of the first peak increases toward larger crystal lengths. The second peak is slowly fading away as a result of product removal. As the crystals reach the product cut size of Lp ) 800 µm, they are removed and the peak starts decreasing. The supersaturation increases, and again new nuclei are formed after one cycle. The time period of this oscillation is 7500 s. A similar dynamic simulation is done for model 2 with the same initial conditions that were used for model 1. The mass density distribution is shown in Figure 6 for the same operational parameters that were taken for the simulation of model 1. The dynamic behavior can be explained in the following way. Initially, crystals are present at smaller crystal length and larger crystal length. The smaller and larger crystals are formed as a result of the attrition process, which is represented by this bimodal distribution. In addition to this, there is a very small peak at very small lengths. These crystals are called crystals of zero growth rate, as explained in the previous section. As the time progresses, the small crystals grow toward larger crystal lengths, as we can see that a small peak moves and grows larger. After some time, when the crystals grow, the attrition process becomes significant and large crystals break and change the shape of the distribution. Further, as time passes, the peak of the larger crystal decreases continuously as a result of attrition and classified product removal. Again supersaturation increases, and the small crystals start growing. The period of these oscillations is 8000 s. Compared to the period of the oscillations that is obtained with model 2, the time period of oscillations for model 1 is smaller. The physical reason behind this is that the growth rate for model 1 is higher and crystals reach the product cut size faster than the crystals in model 2. Further, attrition in model 2 has a large influence on the evolution of the MDF in Figure 6.
Figure 7. Comparison of L50 for particular operating parameters for models 1 and 2 for R1 ) 20 and R2 ) 4.
Figure 8. Dark-gray region: self-sustained oscillations for models 1 and 3. Light-gray region: self-sustained oscillations for model 3.
L50 values of the two models are also compared. Figure 7 is a plot between L50 and time. The solid line in the figure represents L50 of model 1, and the dashed line represents L50 of model 2. From the figure, it is evident that L50 for model 1 and the period of oscillation is smaller than the corresponding quantities of model 2. The physical explanation for this is that in model 1 the growth rate is a little larger than that in model 2, as was already discussed above. As a result, the crystals grow faster, and L50 for model 1 is also larger compared to that for model 2. Further, L50 is lowered for model 2 as a result of attrition. 3.2. Comparison between Models 1 and 3. In this section, finally the influence of nonideal fines classification as proposed in ref 2 is discussed. For this purpose, the instability region in the R1/R2 plane is determined for model 3 and compared to that for model 1. The dark gray region in Figure 8 represents the instability region for model 1 with ideal fines classification. The light gray region, which is overlapped by the dark gray region, represents the instability region for model 3 with nonideal fines classification. From Figure 8, we conclude that nonideal fines classification has a significant influence on the size and location of the instability region. Without product classification at R2 ) 0, model 3 is getting unstable at R1 ) 14, whereas model 1 is getting unstable at R1 ) 73. 4. Conclusions In this paper, the stability of ammonium sulfate crystallization with fines dissolution and classified
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 6955
product removal was studied by means of three different models. From the comparison of the different model predictions, the following conclusions can be drawn: (1) From the comparison of a model with simple kinetics and a corresponding model with detailed kinetics, we conclude that detailed kinetics is surprisingly not that important for the prediction of the instability region in the space of the adjustable operating parameters. (2) However, detailed kinetics makes a big difference if we look at the mass density distribution, which is predicted by the different models and which is of major interest in view of product quality. (3) From a comparison of two models with simple kinetics and ideal as well as nonideal fines classification, we find that nonideal fines classification has a large impact on the size and location of the instability regions in the space of the adjustable operating parameters. A similar strong influence can be expected for the product classification function. Hence, for a reasonable prediction of stability, detailed modeling of the product and the fines classification seems to be important. Furthermore, we conclude that the stability of crystallization processes can be directly manipulated by manipulating the fines and product classification. Future work will focus on experimental verification of these findings. Acknowledgment The authors thank Stefan Motz from Institute for System Dynamics and Control, Universita¨t Stuttgart, Stuttgart, Germany, and Ulrich Vollmer from MaxPlanck-Institute for Dynamics of Complex Technical Systems, Magdeburg, Germany, for providing the detailed model (model 2). Notation b ) nucleation rate exponent B ) nucleation rate [1/mm3‚s] F ) number density function [1/mm3‚mm] g ) growth rate exponent G ) growth rate [mm/s] hf ) ideal fines classification function hfd ) nonideal fines classification function hp ) product classification function kb ) nucleation rate constant [1/g‚s] kfd ) large crystal dissolution factor kg ) growth rate constant [mm/s] kv ) volume shape factor L ) characteristic particle length [mm] Lf ) fines removal cut size [mm] Lfd ) large crystal dissolution cut size [mm] Lp ) product classification cut size [mm] L0 ) minimum crystal length [mm] L50 ) mass median crystal length [µm] L∞ ) maximum crystal length [mm] mt ) suspension density [g/L] q ) feed rate [L/s] qf ) fines removal rate [L/s] qp ) product removal rate [L/s] R1 ) fines dissolution rate R2 ) classified product removal rate
t ) time [s] V ) total volume of the crystallizer [L] w ) weight fraction of the solute win ) weight fraction of the solute going in wsat ) weight fraction of the solute at saturation Greek Letters FL ) density of the liquid [g/L] FS ) crystal density [g/L] ) void fraction
Literature Cited (1) Pathath, P. K.; Kienle, A. A numerical bifurcation analysis of nonlinear oscillations in crystallization processes. Chem. Eng. Sci. 2002, 57, 4391. (2) Mitrovic´, A. Population Balance based Modeling, Simulation, Analysis and Control of Crystallization Process. Ph.D. Thesis, Universita¨t Stuttgart, Stuttgart, Germany, 2002. (3) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes; Academic Press: San Diego, CA, 1988. (4) Kind, M.; Lieb, A. Simulation study on the effect of classified product removal on the dynamic behaviour of continuous crystallizers. In Proceedings of the 14th International Symposium on Industrial Crystallization; IChemE: Warwickshire, U.K., 1999. (5) Gerstlauer, A.; Motz, S.; Mitrovic´, A.; Gilles, E. D. Development analysis and validation of population models for continuous and batch crystallizers. Chem. Eng. Sci. 2002, 57, 4311. (6) Volmer, M. Kinetik der Phasenbildung. Ph.D. Thesis, Theodor Steinkopff Verlag, Dresden, Germany, 1939. (7) Mersmann, A.; Angerho¨fer, A.; Gutwald, T.; Sangl, R.; Wang, S. General prediction of median crystal sizes. Sep. Technol. 1992, 2, 85. (8) Gahn, C. Die Festigkeit von Kristallen und ihr Einfluss auf die Kinetik in Suspensionskristallisatoren. Ph.D. Thesis, Lehrstuhl B fu¨r Verfahrenstechnik T. U., Mu¨nchen, Germany, 1997. (9) Gahn, C.; Mersmann, A. Brittle fracture in crystallization processes. Part A. Attrition and abrasion of brittle solids. Chem. Eng. Sci. 1999, 54, 1273. (10) Gahn, C.; Mersmann, A. Brittle fracture in crystallization processes. Part B. Growth of fragments and scale-up of suspension crystallizers. Chem. Eng. Sci. 1999, 54, 1283. (11) Mersmann, A. Crystallization Handbook; Marcel Dekker: New York, 1995. (12) Vollmer, U.; Raisch, J. Population balance modelling and H-∞ controller design for a crystallization process. Chem. Eng. Sci. 2002, 57, 4401. (13) Kro¨ner, A.; Holl, P.; Marquardt, W.; Gilles, E. D. An open architecture for dynamic simulation. Comput. Chem. Eng. 1990, 14, 1289. (14) Mangold, M.; Kienle, A.; Gilles, E. D.; Mohl, K. D. Nonlinear computation in DIVA -methods and applications. Chem. Eng. Sci. 1999, 55, 441. (15) Kienle, A.; Lauschke, G.; Gehrke, V.; Gilles, E. D. On the dynamics of the circulation loop reactor- Numerical methods and analysis. Chem. Eng. Sci. 1995, 50, 2361. (16) Kubicek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Springer: New York, 1983. (17) Kind, M.; Nieken, U. On the dynamic simulation of mass crystallization with fines removal. Chem. Eng. Process. 1995, 34, 323. (18) Gerla, J. H. Modelling, Measurement and Manipulation of Crystallizers. Ph.D. Thesis, T. U. Delft, Delft, The Netherlands, 1995.
Received for review February 28, 2003 Revised manuscript received September 26, 2003 Accepted October 1, 2003 IE030200U
