Ind. Eng. Chem. Res. 1993,32, 1275-1296
1276
Model Identification and Control of Solution Crystallization Processes: A Review James B. Rawlings,’ Stephen M. Miller, and Walter R. Witkowskit Department of Chemical Engineering, The University of Texas at Austin, Austin, Texas 78712 The modeling of crystal quality (size, shape, and purity) is discussed. Numerical methods for the model solution are reviewed including collocation and Galerkin’s method on finite elements. Measurement technologies for on-line crystal size distribution are presented. The difficulties with laser light scattering methods are discussed in detail. Parameter estimation is reviewed and new methods based on nonlinear optimization are presented. Sample calculations of parameters and their uncertainty from experimental data are included. The research literature covering control of continuous and batch crystallization is reviewed, and sample calculations of optimal cooling profiles for batch crystallizers are presented. 1. Introduction
The objective of this article is to review the area of model identification and control of crystallization processes. The discussion is focused on the following topics: model development, model solution, measurement techniques, parameter estimation, and control. The final section draws conclusions and suggests possible future directions. In this review we attempt to assess the current state of understanding of this field. We identify the issues that are well understood and those that are not. Although effort has been expended to be reasonably thorough in the references, omissions undoubtedly have occurred. These are not intentional, and the authors would appreciate having them brought to their attention. 2. Model Development 2.1. Population Balance. The objective of crystallization modeling is to provide a description of the population of the crystals dispersed in the continuous phase. The mathematical framework that has been developed to describe general dispersed-phase systems is known as population balance modeling. The concept of a population balance arose as a generalization of the residence time distribution. In early work, Rudd 11611 demonstrated that one cannot always formulate a model based on particle age, and it is necessary to consider a generalized describing variable (catalyst particle activity in his example). Other contributions during this time were made byRandolph [136,137,141] Curl [38],andBehnken et al. [81. These early efforts were followed by the classic paper by Hulburt and Katz 1821 and book by Randolph and Larson 11421. These authors extended the concept of particle size distribution to include an arbitrary number of “internal” coordinates necessary to fully describe the state of a particle. Tsuchiya et al. [1811,Eakman et al. 1481, and Ramkrishna [132-1351 apply this methodology to microbial populations. Hidy and Brock [77l, Friedlander 1561,and Hidy [ 761 review population balance methods in aerosol science applications. Rawlings and Ray 1146-1491 review the contributions in emulsion polymerization that employ population balance models, and describe the advantages of this description over other models. Ramkrishna’s 11351
* To whomcorrespondenceshould be addressed. Phone: (512) 471-4417; email: jbraw @ che.utexas.edu. + Current address: Sandia National Laboratories, Division 1434, Albuquerque, NM 87185-5800. 08SS-58S5/93/2632-1275$04.00/0
review article provides an extensive discussion of population balance modeling, applications, and solution methods. 2.1.1. Stochastic versus Deterministic Formulation. One of the first decisions that has to be made in modeling a particulate system is whether to use astochastic or deterministic framework. This decision is usually based on whether one considers the particle population fluctuations to be important. If fluctuations are large and need to be modeled for accurate predictions of observable5 of interest, they can be described in a stochastic framework. Ramkrishna 11341 and Fan and coworkers [53, 80, 811 present stochastic population balances. Ramkrishna is concerned with modeling small populations of cells, for example in sterilization processes. At small population levels, the fluctuations about the mean number of cells become of the same order of magnitude as the mean itself, necessitating the more complex stochastic machinery. Fan and co-workers,however, recommend using the stochastic approach even for large populations of crystals. They attempt to model crystal size distribution (CSD) data displaying deviations from the decaying exponential size distribution expected of well-mixed crystallizers. It is not clear why a stochastic approach is required since the data can be fit by other deterministic models that include sizedependent growth or growth rate dispersion. It also seems plausible that, given the difficulty of measuring CSD, fluctuations in observed CSD are likely due to measurement noise rather than true crystal population fluctuations. In this case it would be reasonable to model the crystallizer with a deterministic population balance and simply append additive stochastic variables to account for the measurement noise. One could then design filters and regulators for CSD control that account for the expected CSD measurement fluctuations. Therefore in this article we focus exclusively on deterministic population balances. The key property of the particulate system is some general number distribution function. Let f ( z , t ) denote the crystal distribution function, which depends on the time, t , and other suitably chosen coordinates, z. The z coordinates include the usual three spatial (external) coordinates as well as the additional internal coordinates which are required to completely specify the state of the crystal. Let dV, denote an infinitesimal volume in the crystal phase space. The total number of crystals per volume of crystallizer in some 0 1993 American Chemical Society
1276 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993
subregion of crystal phase space, a, is
One derives a population balance by accounting for changes in the crystal number distribution due to changes in the coordinates z and any crystal formation and removal mechanisms. The mathematical arguments are analogous to those used in the development of conservation equations in transport phenomena 11061. The resulting population balance is
where B and D are the density functions of formation of new crystals and removal of existing crystals, respectively. If one considers the ideally stirred reactor, which has no spatial dependence, it is convenient to split the z coordinates into the external, x, and internal, y, parts,
Since neither f, B, nor D have any x dependence, the microscopic population balance can be integrated to give
in which fk(y) is the crystal distribution function of the kth stream, Qk is the volumetric flow rate of the kth stream (positive sign denotes flow out), and V is the crystallizer slurry volume. With suitable choices of the Q k , eq 3 describes batch, semibatch, and continuous crystallizer operation. The crystallizer is clearly a distributed parameter system regardless of whether it is spatially distributed as well. To complete the description, one has to choose the internal coordinates, y, specify the functional forms of ay/at, fk, B, and D, and supply the necessary initial and boundary conditions. Specifying these functions and conditions requires a detailed understanding of the transport and kinetic events taking place, and is the major hurdle to overcome in modeling any particulate system. 2.1.2. Crystal Quality: Size, Shape, and Purity. Equation 3 is a general starting point, and has been used to address issues such as growth rate dispersion, coalescence and breakage of particles, and a variety of other mechanisms. The choice of internal coordinates in vector y is dictated by the intended use of the model. In industrial manufacture of crystals, the crystal quality is usually determined by the crystal size, shape (habit or morphology), and purity (composition). The vast majority of the engineering research literature is concerned with modeling the crystal size and size distribution. In this article we too shall review and formulate models intended to enable improved control of crystal size distribution. However, a brief review of the progress being made in crystal habit and purity modeling and control is worthwhile and points out several new opportunities for improved crystal quality control. The fact that crystal size is important in applications is rather obvious: crystals that are too small are not easily separated from the mother liquor and may cause handling and environmental problems due to dust formation. A poorly shaped CSD may promote agglomerationand caking and pose storage and processing problems. However, crystal morphology may play an equally important role in many applications. For example, precise control of the twinning grain morphology and size distribution in high-
speed films is vital for achieving their superior photographic properties. The aspect ratio of the particles in powdered reinforcement influences the strength and toughness of composite materials [1521. If the final product is a suspension of crystals in a liquid, the fluid rheology is strongly affected by the particle morphology [1111.
The relative growth rates of the different faces of a crystal determine its shape. The crystal habit is determined by the slowest growing faces. This phenomenon serves as the basis for attempts to control habit by developingtailor made additives to serve as habit modifiers [I, 21. A tailor-made additive is similar in chemical structure to the crystallizing species, so it is incorporated into the crystal structure. The difference between the additive and the crystallizing species, however, inhibits the growth of the crystal face containing the additive. This can produce a dramatic effect on the crystal morphology. A comprehensive recent review of this area is provided by ref 40. Another recent development that may have a significant effect on industrial practice in the near future is the development of public domain software for morphology prediction [36,44,451. The use of this software for quantitative prediction of crystal morphology for certain ideal organic crystals is well established, and progress is being made for the more difficult ionic crystals as well [1111. Attempts have been made to calculate the effect of impurities and different solvent mixtures on final crystal habit [1211. This area has been reviewed by Myerson [1191. His conclusions are that although there are not as yet reliable methods to predict habit or select solvents and habit modifiers without extensive trial-and-error experimentation, recent progress may make habit and purity control a realistic goal for the near future. 2.1.3. CSD Modeling. Therefore, although states of a crystal such as purity and shape can be included as elements of y, we restrict ourselves to modeling the crystal size distribution (CSD),in which case thevector of internal coordinates is simply L, a characteristic crystal length. If one additionally makes the following assumptions: no crystal breakage or coalescence,and crystals are nucleated at negligibly small size, LO,then eq 3 reduces to
in which G is the crystal growth rate, Bo is the nucleation rate density, and 6(L- LO)is the Dirac delta function acting at L = LO.The functions fk(L) can be defined to include washout from a CSTR (mixed product removal) and sizedependent removal (classified product removal and fines destruction). If LO 0, it can be shown that eq 4 is equivalent to
-
(5)
with the boundary condition BO (6) f(0,t) = GI,=, 2.1.4. Mass and Energy Balances. As illustrated by the constitutive relations for G and Bodiscussed in section 2.2, growth and nucleation rates are functions of supersaturation (a measure of the difference between concentration and saturation concentration) and temperature. Therefore, mass and energy balances are required to complete the description of the crystallization process.
Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 1277 The solute mass balance in the continuous phase must account for flow of solute into and out of the system and mass transfer to the solid phase via formation and growth of crystals. The mass balance can be written as
where C is the solute concentration in the continuous phase, and pc and k , are the crystal density and volumetric shape factor, respectively [I 141. The fraction of V that is solidsfree (continuous phase) is represented by E , which is equal to 1 - k , p3 in which p3 is the third moment of the CSD. The energy balance given below accounts for enthalpy differences between the flow streams, the heat of crystallization, and heat removal by cooling systems, evaporation, etc. pc,v-
dT
d(PV) = -dt dt
Tis the crystallizer temperature, cp is the specific heat of the slurry, p is the slurry density, A Hc is the heat of crystallization, Hk is the enthalpy of the kth stream, and Hextis the net heat removal, incorporating effects of the cooling system, evaporation, loss to surroundings, and physical mixing. P is the system pressure; the pressure effect term is usually neglected for cooling crystallization. In reactive crystallization, the supersaturation is generated by chemical reactions. To model such cases, one appends the chemical reaction rates multiplied by appropriate stoichiometric coefficientsto the right-hand sides of the species’ mass balances. Expressions for the reaction rates in terms of species concentrations are also necessary. One then requires a vector of mass (or mole) balances to describe the system. In order to simplify the discussion, the case of reactive crystallization is not presented explicitly. However,the modifications required to handle this case are, in principle, straightforward. The major difficulty is coming up with appropriate reaction rate expressions. Equations 5-8 can be used to describe evaporative or cooling crystallizers in a variety of configurations,including batch and continuous modes with size-dependent removal. The methods of model identification and control discussed in this article are illustrated in terms of cooling crystallization but are easily extended to apply to evaporative crystallization. 2.2. Constitutive Relations. In order to complete the modeling of the crystallization process, constitutive relations are required for growth and nucleation, the two competing phenomena in the phase change from solution to solid phase. I t is well-known that neither growth nor nucleation can occur unless the solute concentration exceeds its saturation value (i.e., the solution is supersaturated). However, there are many factors that affect formation and growth of crystals, including crystallizer hydrodynamics, temperature, and presence of suspended or dissolved impurities. A brief discussion is given below of some of the theories and models used to describe growth and nucleation. It should be noted that, for use in modelbased control systems, the validity of the chosen models should be examined over a range of operating conditions to ensure predictive capability. 2.2.1. Growth. At the microscopic level, solute molecules moving from the bulk solution adsorb on the crystal solid surface and are incorporated into the crystal lattice.
Figure 1. Crystalincorporation sites. Flat faces (A), step sites (B), and kink sites (C).
A well-defined,smooth crystal face is planar and new solute molecules must migrate across the surface to find energetically favorable incorporation sites. Three such potential sites are shown in Figure l. Site A is thermodynamically unfavorable compared with B, a step site, or C, a kink site. Surface adsorption and diffusion determine whether a solute molecule is incorporated into the crystal or returns to the bulk phase. The observed growth rate is then caused by the flow of steps across the surface. The origin of the steps on the surface was originally postulated as the result of collisions between solute molecules that had adsorbed onto the surface. These collisions allow the formation of a critical-sized twodimensional nucleus on the surface. However, the predicted rate of step formation with this mechanism is smaller than the experimentally observed growth rates, especially at low supersaturation. Burton, Cabrera, and Frank (BCF) [25] overcame this deficiency by considering the step to be formed by a screw dislocation, a commonly observed crystal imperfection. The crystal continues to grow without the formation of new steps as solute is fed to the steps by surface diffusion flux governed by Fick‘s law and is incorporated along the spiral. The simplest form of the BCF model predicts the macroscopic growth rate to have the form G = k ’ - S2 t a n h ( T )~ B C F ~BCF
S=- c - cent cent
in which S is the relative supersaturation. If the surface diffusion and incorporation rates are fast enough, the mass transfer of solute molecules from the bulk to the surface also affects the macroscopic growth rate [101,1891.Most of the chemical engineering literature treats this problem with an empirical mass-transfer coefficient that must be correlated with the fluid dynamics and fluid properties. Other approaches include treating the diffusion through a stagnant film, or applying some form of boundary layer theory [60].Garside [591provides a thorough survey and discussion of growth rate theories. For engineering purposes, the semiempirical power law has become the standard representation of the growth rate, G = k , exp(-EglRT)Sg (10) Note that g = 1 and g = 2 are the high and low
supersaturation asymptotes of the BCF model, respectively. Also notice that the temperature dependence has been incorporated using an Arrhenius-type expression. Neither of the growth rate models given above indicate dependence of G on L. This size-independent growth rate assumption is usually referred to as McCabe’s AL law. There are several examples of systems that violate
1278 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 Nucleation
Secondary
Primary
Homogeneoua
Heterogeneous
with Externd surfsce
Crystal-Crystal Initial Contact Breeding
fracture
Dendritic Separation
Figure 2. Nucleation mechanisms.
McCabe’s AL. law; Canning and Randolph 1261 give a discussion of such systems, and Budz et al. 1221 give an expression for the size-dependent growth rate of potassium sulfate. Size-dependent growth rate is usually attributed to either bulk diffusion effects or the Gibbs-Thomson effect, which suggests an increasing growth rate with increasing size because of an inverse relationship between solubility and size. Garside et al. 1641 present a theory of size-dependent surface integration kinetics. Modeling the growth rate is further complicated by a phenomenon known as growth rate dispersion. Growth rate dispersion describes the situation in which not $1 of the crystals grow at identical or constant rates even though the crystallizer conditions remain constant. As Girolami and Rousseau 1683 demonstrate, many cases that have been thought to be occurrences of size-dependent growth are actually manifestations of growth rate dispersion. Considerable attention recently has been focused on growth rate dispersion, and Zumstein and Rousseau 11991 provide a review of these results. Two mechanisms have been put forward to explain deviations in the CSD due to growth rate dispersion. Randolph and White 11441assume that the crystals have the same time-averaged growth rate, but individual crystals experience stochastic fluctuations in the growth rate causing the growth rate dispersion. Janse and de Jong 1891and Berglund and Larson 1121postulate that each crystal has a constant, deterministic growth rate but that these growth rates differ from crystal to crystal; there has been some work that supports this hypothesis by relating mechanical stresses and strains of a crystal to its growth rate [13,661 . Modifications to the population balance, eq 5, required to described either of these cases are available. As demonstrated by White and Wright [1921 the spread of the seed population CSD during batch experiments is a good experimental test to determine whether growth rate dispersion is an important factor. 2.2.2. Nucleation. The nucleation of crystals from solution has been reviewed extensively [61,107,142,1741. The complexity of mathematically modeling the formation of new crystals is largely due to the variety of mechanisms that are collectively termed “nucleation”. Nucleation mechanisms are commonly lumped into one of two categories-primary and secondary nucleation-and can be further classified as shown in Figure 2. Mechanisms of formation of crystals that are independent of the presence of other suspended crystals are classified as primary or spontaneous nucleation. Primary nucleation is associated with high levels of supersaturation and is usually partitioned as homogeneous nucleation, which occurs in the pure bulk solution, and heterogeneous nucleation, which is induced by foreign surfaces such as impurities. Classical thermodynamic free energy minimization is used to derive the rate of homogeneous nucleation 11901.
This theory postulates the production of embryos from the combination of solute molecules in a series of bimolecular reactions. The free energy of the embryos achieves a maximum at a critical size particular to the chemical system. Once an embryo exceeds this critical size, the free energy decreases with further growth, leading to spontaneous nucleation. The nucleation rate is then determined by the rate at which embryos cross the free energy maximum. Except in the case of precipitation crystallization, the very high supersaturation required to overcome interfacial tension between embryo and solution makes homogeneous nucleation an unlikely mechanism for crystal formation for most chemical systems under industrial conditions 135, 1901. Foreign surfaces and particles promote nucleation as a result of an ordering process caused by interactions across the interface [1901. The result of this catalysis is that primary, heterogeneous nucleation can occur at supersaturation levels significantly lower than required for homogeneous nucleation and is, therefore, the dominant mechanism of primary nucleation when impurities are present. Nevertheless, the supersaturation levels of heterogeneous nucleation are often still too high for good crystal growth and production of crystals of desirable morphology 1351, Also, classical theory suggests that primary, heterogeneous nucleation is characterized by a process that is either starved for nuclei or overwhelmed by a burst of new crystals [351, making CSD control difficult. In addition to the issues of poor crystal quality and control difficulties associated with primary nucleation, attempts to predict and scale-up the primary nucleation rate have been largely unsuccessful. This is a result of the difficulty of accurately characterizing the low levels and types of impurities, and describing all the locations of high local supersaturation (i.e., near cooling coils or boiling surfaces). Fortunately, secondary nucleation, which is more easily controlled than primary nucleation and occurs at supersaturation levels conducive to good crystal quality, always accompanies primary nucleation and is the dominant mechanism in most industrial crystallizations [17,69,141I. Secondary nucleation describes the nucleation that takes place due to the presence of other solute crystals. There are a variety of proposed mechanisms whereby the crystals promote formation of new crystals. Botsaris 1171 organizes the promotion around two questions: what are the sources of the potential nuclei, and how are the potential nuclei extracted from the source and displaced into the bulk solution to initiate new crystals. Proposed sources of potential nuclei include protrusions from growing crystal surfaces, solute clusters on or near the crystal surface, and embryos in the supersaturated solution. Several mechanisms for the conversion of “potential nuclei” into nuclei
Ind. Eng. Chem. Res., Vol. 32, No. 7,1993 1279 have been proposed; these include spontaneous removal of dendrites due to free energy driving force, fluid shear, and contact nucleation (contact of crystals with an external surface or other crystals) [ I 7, 35, 92, 1741. The variety and complexity of the mechanisms of secondary nucleation have frustrated attempts to use theory for quantitative prediction. Therefore the modeling of nucleation is usually simplified by assuming an empirical functional form. The following expression is commonly used to described secondary nucleation:
Bo = kb eXp(-Eb/T)Sbp{
(11)
in which kb, b, E b , and j are considered to be empirical constants, Pk is the kth moment of the CSD, and S is the relative supersaturation. The mechanistic effect of temperature on Bo is complicated and, as Garside and Shah [65] show with a summary of several kinetic studies, nucleation can actually decrease with temperature, which corresponds to a negative activation energy,E,. A common hypothesis is that higher temperatures lead to increased growth rates, involving greater efficiency of molecular diffusion on and integration into crystal surfaces; thus, fewer potential nuclei are available for secondary nucleation [75]. Nevertheless, evidence suggests that b is independent of temperature [651 and an Arrhenius-type expression is probably adequate for characterizing temperature dependence. While secondary nucleation rate is usually correlated with p3, the nature of the phenomenon suggests that the surface area, k&z, would be a reasonable choice and has been used by some researchers 11891. In many cases, contact of crystals with external surfaces is the predominant initiator of secondary nucleation and a first-order dependence on p3 is expected; several studies validate this assumption [651. On the other hand, Sherwin et al. [I671 claim that, in the presence of a large amount of crystals, the dependence of nucleation rate on supersaturation is very strong, while the dependence on crystal surface (or mass) is a small effect and can be ignored 0’ = 0). As shown in the summary of kinetic studies by Garside and Shah [651, j = 0 is commonly used. Some researchers even have found the correlation between Bo and p3 to be such that j is negative; Winter [I941 surmises that this corresponds to occurrence of substantial agglomeration. For scale-up purposes either the stirrer speed, power input, or tip speed are also included as factors in eq 11165, 1653. 2.3. Model Parameters. When using power law or BCF constitutive relation for the nucleation and growth rates, the parameters b, g, kb, k, (or k’, and ~ B C Ffor the BCF relation), E b , and E, need to be estimated on the basis of experimental data. As discussed above, k and j are often chosen using mechanistic considerations. All other parameters such as the crystal‘s density, area, and volume shape factors, etc. are considered known or available from independent experiments. As discussed further in section 5.2, another goal of experimental model validation is to assess the limitations of the nucleation and growth constitutive relations and determine their predictive capability outaide of the narrow operating region in which they are identified. 3. Model Solution
The crystallizer model described in the previous section consists of a coupled set of nonlinear integro-differential equations. The population balance itself is a hyperbolic partial differential equation. The integrals appear because
the transfer of mass from the solution to the solid phase must be represented as a sum over the rates for each individual crystal. Thus the overall mass and energy balances, eqs 7 and 8, involve integrals over the distribution. The model needs to be solved rapidly and accurately in order to achieve the goals of the parameter estimation and feedback control sections discussed subsequently. Analytical solutions of the nonlinear models of interest do not exist, and although some stability issues can be treated with linearized models, the numerical solution of the full model is also necessary. Unfortunately for our purposes, the numerical methods to treat this class of equations are much less developed than the methods available for more standard engineering problems such as ordinary differential equations and optimization or solution of linear and nonlinear algebraic equations. Standard software does not exist to solve the models of interest, although this will hopefully change as more numerical analysts turn their attention to this class of problems [1021. Recent efforts, such as the book by Delves and Mohamed 1431,indicate that progress is being made on this front. In the meantime, applied scientists and engineers have employed a variety of numerical methods to solve specific models on a case-by-case basis. This literature is reviewed subsequently. Attention is focused on efficient methods, since model-based control algorithms must be solved in real time. Certain classes of methods, such as finite difference, are not discussed for this reason. 3.1. Global Methods. The method of weighted residuals has been applied most often to solve integrodifferential equations. Ramkrishna 11351 has recently reviewed many of these efforts. In the method of weighted residuals the population density function is approximated as a linear combinations of chosen basis functions, \ki(L),
f a t ) = C U i ( t ) *i(L)
(12)
i=l
The \ki(L)are known functions on the interval IO,-). The methods are called global because the basis functions are used to describe f over the entire L domain. The unknown ai are determined by substituting eq 12 into the population balance to define a residual, r(L,t). If the residual is identically zero, then one has the true solution. The sense in which the residual is made small determines the type of weighted residual method being used. The idea of weighted residuals is to find the ai that force the residual to be orthogonal to a chosen set of weighting functions, W i ( L ) . If the basis functions can represent f well with a small number of terms, then it makes sense to force the residual to be orthogonal to these same functions,
c * i ( L ) r(L,t) dL = 0
(13)
This choice of weighting functions is known as Galerkin’s method. Witkowski and Rawlings [I981solve a crystallizer model using this approach with Laguerre polynomials as the basis functions. They calculate accurate solution of both stable and oscillating cases with as few as four basis functions. The modest computation time required to solve the model allows one to consider using it as the basis for an on-line control scheme. Chang and Wang [30,311 use Galerkin’s method with shifted Legendre polynomials. Because CSDs emanating from a continuous crystallizer have an exponential hi1 at large L, Legendre polynomials do not represent f well unless one uses many basis functions. Laguerre polynomials have the exponential weighting function and are therefore a more natural and
1280 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993
efficient choice. Hwang and Shih [83] use the Galerkin method with block pulse functions as basis functions, and again a very large number of basis functions are required for accurate solution. Wang et al. [I911 expand this approach by using generalized orthogonal polynomials and state that the flexibility of this technique leads to accurate results for a wide array of crystallization problems. If the integrals in eq 13 must be computed numerically, then one should consider the collocation method. The collocation method makes the residual small by forcing it to vanish at selected Li locations, called the collocation points. The choice of these locations is vital. Using the is known as orthogonal collocation. Villadsen zeros of and Stewart [I871 show that orthogonal collocation is equivalent to Galerkin’s method if one uses an optimal Gaussian quadrature rule to evaluate the integrals in eq 13. Ramkrishna [1331 discusses problem-specific polynomials and suggests their application in the method of weighted residuals in solving integro-differentialequations. Sampson and Ramkrishna [I631 and Singh and Ramkrishna [169,170] use collocation with problem-specific polynomials for the basis functions to solve population balances describing coagulation. Finally Lakatos et al. 11051 solve a crystallizer model using collocation with Laguerre polynomials. The method of moments, discussed thoroughly by Hulburt and Katz [82] and Randolph and Larson [1421, is probably the most popular method for solvingpopulation balances. Ramkrishna points out that the method of moments is equivalent to the method of weighted residuals if one chooses the weighting functions to be polynomials,
w,(L) = Li-l
(14)
Since polynomial weights are usually a poor choice for population balances on semi-infinite intervals, the method of moments is not recommended. Randolph and Larson [I421 also demonstrate that reconstructing f from its moments is numerically unstable. 3.2. Finite Element Methods. In finite element methods, the solution is approximated over subdomains of the independent variables. Specifically for problems described by eq 5, the [O,.D) L domain is subdivided into elements, and f is approximated by basis functions on each element. Usually the order of the approximation in the elements is lower than in the global methods. Thus, if a suitable basis can be found for a global method, it usually converges faster to the true solution than a finite element method 1431. A limitation of global methods is that one cannot always find an efficient set of basis functions, and the convergence is quite slow and oscillatory if there are discontinuities in the solution. Discontinuities in df/dL arise in crystallizer systems with product classification or fines destruction. Finite element methods can be tailored to handle discontinuities and sharply changing CSDs. Peterson et al. [12Aand Gelbard and Seinfeld 1671use collocation on finite elements to solve population balance equations with coalescence, condensation, nucleation, and removal terms. The semi-infinite L domain is truncated at some large value. The finite domain is then divided into elements. Low order shifted Legendre polynomials are used to approximatef on each element. The collocation method is used to minimize the residual. In these studies, the polynomial approximations off and dflcU are forced to be continuous at element boundaries. One can obviously change the continuity conditions for problems with discontinuities in f o r dfldL. The choice of the element locations is important for efficient solution. Gelbard and
Seinfeld use a logarithmic scaling to more effectively place elements at small sizes where large numbers of particles occur. The best choice of element location is also obviously problem dependent. Villadsen and Michelsen [I861 provide a general discussion of collocation on finite elements for boundary value problems. Sastry and Gaschignard [I641 also solve a coalescence population balance by polynomial approximation on finite elements. They also truncate the semi-infinite domain and use problem-dependent heuristics to choose the element locations. Rawlings and Ray [I#, 1491 use collocation on finite elements to solve population balances arising in emulsion polymerization. Discontinuities in df/ dL are handled efficiently by allowing the element boundaries to move with the discontinuity. Steemson and White [I 731 solve the steady-state continuous crystallizer model using cubic spline approximation. The flexibility of the method allows them to solve a complex system with several crystallizers, classifiers,and recycle streams. They also demonstrate improved efficiency over the method of Chang and Wang [301 described previously. Witkowski [I963applies collocationand Galerkin’smethod with finite elements to solve dynamic crystallizer models. He demonstrates that accurate solutions to batch models can be generated with modest amounts of computation time. None of the presently available methods produce an efficient and accurate solution for a broad class of models. Most solution methods have been tailored to specific problem types. The difficulty in generalizingthe approach arises because the best element locations and polynomial orders within elements may change with time. One might naturally consider adapting or refining the grid as the solution changes. The adaptive finite element method is an active area of current research. Russell and Christiansen 11621 and Delves and Mohamed 1431 provide reviews of adaptive methods in boundary value problems and integral equations, respectively. Most of these approaches have not yet been teated on population balance models, but could be applied if the global and fixed finite element methods are insufficient. 4. Measurement Techniques for Crystallization Processes
The discussion of this section focuses on the possible measured variables in crystallization processes, and is important background for understanding the parameter estimation and control problems discussed in subsequent sections. 4.1. Concentration Measurement. A measurement of the concentration of the solute in the continuous phase is needed to characterize the state of a crystallizer. There are several methods for the rapid measurement of continuous-phase concentration. It has been shown that quick and accurate determination of concentration from the measurement of refractive index is possible[75,117,1681. Garside and Mullin [63] suggest the use of an on-line densitometer to determine solution concentration. Qiu and Rasmuson [I311 report successful use of density measurements in assessing concentration of samples extracted from a crystallizer, and Witkowski 11961,Riebel et al. [1511, and Miller 11141 have demonstrated the use of densitometry for accurate, on-line concentration measurements. 4.2. Dynamic Crystal Size Distribution Measurement. Obtaining on-line information about the CSD is necessary for the development and implementation of crystallizer control algorithms, and as demonstrated in section 5, accurate measurement of solid-phase properties
Ind. Eng. Chem. Res., Vol. 32, No. 7,1993 1281 Photodiode Detector
Sample Cell
L
F
4
Figure 3. Schematic diagram of particle aizer based on light scattering theory.
is required for uniquely determiningthe model’s nucleation and growth parameters. This section provides an overview of the available methods for particle sizing with emphasis on the laser light scattering method. Sieve analysis is the oldest particle sizing technique. While it provides a goodcheck of the quality of a crystallizer product, it is time consuming and not suitable for control purposes. Bunville [24] and Lieberman [IIO] point out that there are now several commercially available instruments, based on a wide variety of physical phenomena, that provide on-line particle size information. Coulter counters are based on electrical resistance differences between the particle and solution. Sedimentation and hydrodynamic chromatography are based on transport properties. Optical properties give rise to techniquesbased on diffraction, light blockage, and direct imaging. Rohani and Paine [I571suggest a device for measuringsuspension density based on temperature change required to dissolve crystals in a sample cell. Unfortunately, no method developed to date is without restrictions and problems. Coulter counters suffer from orifice-blockingproblems that prohibit industrial application, or at least restrict their application to streams that have been classified by size. Coulter counters and optical-property techniques are sensitive to the presence of foreign particles and entrained bubbles. The temperature-difference technique of Rohani and Paine is not affected by foreign particles and bubbles, but requires substantial calibration to allow absolute suspension density measurement and is better suited for detecting changes in suspension density. Particle sizing with light scattering is an attractive candidate in crystallization applications since it generally avoids the pore-plugging problems of Coulter counters, requires little calibration, and is easily automated for possibleuseasafeedbacksignalfor acontroller. Randolph et al. [I451 use an on-line laser particle sizer to monitor the nucleation and growth rate of a continuous KC1 crystallizer. Brown et al. 1201 and Felton and Brown 1551 use light scattering to monitor sugar and KI crystal growth rates, respectively. Boxman’s recent thesis [I81 provides a thorough discussion of the issues that must be resolved in order to use light scattering instruments as reliable sensors for crystallizer control. 4.2.1. Particle Sizing Using Light Scattering Theory. Particle sizing with laser light scattering is based on the interaction between an incident light source and solid particles. This interaction is discussed thoroughly by Borhen and Huffman [161. The amounts of reflection, refraction, and absorption of the incident wave depend on the refractive index of the particle, its geometry, and the wavelength of the incident light. The full Maxwell equations for the electromagnetic field at all spatial points
due to interaction with an optically smooth sphere can be solved analytically [1121. This solution, known as Mie scattering, allows one to predict the intensity of the scattered light as a function of scattering angle; it involves complex, oscillation functions that can be evaluated efficiently using numerical methods given by Dave [39] and Wiscombe [1951. In most crystallization applications, the crystals of interest are large compared to the wavelength of the incident light (typically a He-Ne laser, X = 0.632 pm). For these large particles, the phenomenon primarily responsible for scattering light in the “near-forward direction” (small scattering angles) is Fraunhofer diffraction. In studies by Jones 1983,Fymat and Mease 1581,and de Boer et al. [421, the full Mie scattering solution is compared with Fraunhofer diffraction results, and it is shown that, for many cases of interest, a valid assumption is that the near-forward scattering pattern is solely due to Fraunhofer diffraction. An advantage to the Fraunhofer approximation is that the scattering depends on the projected area of the particle and is independent of refractive index. Figure 3 is a schematic of an instrument developed by Swithenbank et al. [I 751, which is produced commercially by Malvern Instruments. The scattering intensity pattern is difficult to measure, so the instrument measures the scattering energy pattern. The scattering intensity decreases markedly with increasing s (increasing scattering angle),so areas of the annular photodiodes, which are used as detectors, increase by about 3 orders of magnitude from the center to the periphery to compensate and reduce problems related to the limited dynamic range of the photodiodes. The Fraunhofer approximation allows the derivation of the following simple expression for the energy that is scattered by a particle of radius L to the ith annular detector of inner and outer radii SI and 8’1, respectively: e(L,si) = KTL’{[ J(,
T)
2TLSi
+ J:(
T 2ULSi )-]
where F is the focal length of the lens, and K is an optical constant that depends on the power of the laser and the detector sensitivity. If the suspensionof particles passing through the sample cell is dilute enough that multiple scattering can be ignored, the energy of scattered light at the detector is a simple summation of the energies of the individual particles times the number density distribution, e,(si) = K e ( L , s i )f ( ~ d~ ) (16) Since it has been assumed that the particles are spherical,
1282 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993
20
Tc,
-
25
-
original distribution recovered distribution ---.--
2o
=8
10 15t
,
i
. . . . . . . . . . . . . . . ,. original distribution recovered distribution ------
I
r
i
L.I
1
10 100 Particle Size (microns)
1000
Figure 5. Recovery of the CSD from light scattering data with 5 7% noise using penalized and constrained least squares with y = lo-‘.
an expression for et in terms of amass density distribution instead of number density distribution is et(si)= K e ’ ( L , s i )m ( L ) d~
(17)
where J Lm ( ~L ) dL ~ is ~the mass of particles of density p with L1 5 L 5 Lz, and e’(L,si) is given by e’(L,si)= ,e(L,si) 4lrpL Light scattering theory provides eq 17, which calculates the scattering pattern for a given CSD. However, particle sizing requires the solution of the inverse problem-determining the CSD based on a finite number of measurements of the scattering pattern. The conversion from number to mass distribution is made because the common strategy is to obtain weight fraction distribution instead of the number CSD; the normalization removes the optical constant K. One approach to the inverse problem is to employ a quadrature rule to evaluate the integral in eq 17. One then derives a matrix equation for the values of m at the p quadrature points, P
et(si)= x w j e ’ ( L j , s i )m(Lj)
i = 1, .... n
(18)
1=1
in which Lj are the quadrature points and W j are the quadrature weights. As is the case for many inverse problems, the inverse operator for the integral equation and, thus, eq 18 is illconditioned. This ill-conditioning is due to the fact that significantlydifferent CSDs can give very similar scattering patterns; therefore, small errors in the measurement of et can lead to large errors in the CSD inference. Wing [I931 reviews this class of problems and discusses the methods that have been developed to handle the ill-conditioning. Since the matrix equation is in general overdetermined (p < n),one can solve the least squares problem for the m(Lj). This fails completely in practice because of the ill-conditioning of the E matrix (Eij = wj e’(Lj,si));CSDs with wild variations, including negative components, are common. For example, Figure 4 gives the results of the least squares solution for a scattered light pattern predicted from theory for a known particle size distribution with no measurement noise. This is an illustration of the instability due to ill-conditioning, because the substantial error in the recovered distribution is for a case in which the only error in the scattered light pattern is on the order of the machine precision of the computer. To avoid these instabilities, statisticians often use a biased estimator called ridge regression [791,which simply uses a least
squares objective function that is modified to include a penalty function on the variations in the solution. This is also sometimes referred to as constrained linear inversion (an unfortunate choice since constraints are not actually imposed). Phillips [I281 and Twomey 11831 suggest using a penalty on the second derivative of the recovered distribution. The modified or penalized least squares solution is given by
m = (ETE+ yH)-lETe,
(19) in which H is the smoothing matrix (e.g., coefficients for the finite difference approximation to the second derivative). A further modification to ensure a physicallymeaningful solution is to incorporate the constraint that the solution must be everywhere positive; the result is the following quadratic programming problem min cTm+ l/z mTBm m subject to
(20)
.....p
m j 2 0 , j=1,2
where B = ETE + yH, and c = -2ETet. A comparison of Figures 4 and 5 demonstrates the benefit derived from the constrained, penalized inversion algorithm given by eq 20. The distribution is inferred with greater success using this inversion algorithm then when using simple least squares, even though the data contain 5% random noise while the data used in the least squares inversion contains no added noise. While Figure 5 illustrates that the constrained, penalized least squares approach can be effective, the choice of y is important. The tradeoff is that y must be large enough to avoid instabilities, but not so large that the solution is strongly biased. The approach used by many researchers is atrialand-error procedure for determining y 1341, but other methods are available. Crump and Seinfeld [37J,suggest a cross-validation method for determining an appropriate value for y. Shaw [I661demonstrates that y can be treated as a parameter of a low-pass filter that removes the oscillatory modes of the solution; the eigenvalues of ErE can then be used as a guide for the selection of y. Boxman [I81demonstrates the use of a weighted least squares technique to solve the inverse problem. The weighting matrix is the inverse of the covariance matrix of the light scattering measurements; it is obtained from observations of fluctuations of light intensity over time. This weighting and a nonnegativity constraint help to stabilize the solution. An iterative approach is used to minimize a measure of the goodness of fit to the data. Boxman also demonstrates the use of a statistical analysis
Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 1283 According to the Beer-Lambert law, 2.520
.
15
-
10
-
where T is the slurry turbidity, 1 is the flow cell width, A,(L) is the projected area for a crystal of characteristic size L, and Q(L) is the extinction efficiency factor. Assuming Fraunhofer diffraction is the dominant scattering phenomena, Q = 2 for all L. An extension to cases of nonspherical particles can be made by applying the theorem discussed above relating geometrical cross section to surface area, transforming eq 22 to the following:
5 -
to assess the reliability of a solution and to guide the way in which the size domain is discretized. A further complication of inferring the CSD from scattered light intensity measurements is due to the assumption that the particles are spherical. Brown and Felton [19]suggest making particle shape corrections using a theorem proved by van de Hulst [1851, that states that the average geometrical cross section of a convex particle with random orientation is one-fourth its surface area. The light scattering results are taken to be for equivalent spheres, and the theorem is used to shift the equivalentsphere diameter to a characteristic length of the actual particle. Jones et al. [951 demonstrate that agreement with sieve analysis is obtained when this method is used to correct particle size to correspond to the characteristic length of potassium sulfate crystals. Many crystals have needle or platelet morphological structures. Methods based on a shift of equivalent-sphere diameters could lead to large errors for such large aspect ratio particles. To illustrate, Figure 6 presents the results of a Malvern particle sizer for a monodispersion of cellulose acetate cylinders (45-pm diameter, 195-pm length). Note that instead of the description of a monodispersion, the distribution is bimodal with modes located near the primary dimensions of the cylinders. Obviously, no method that simply shifts the CSD would be adequate for crystals of large aspect ratios. Another assumption made in developing eq 17 is that there is no multiple scattering. For high-density slurries that are common in industrial crystallizers, this poses a problem. Some studies have been made of the effect of multiple scattering and on methods to extend the light scattering theory to account for this effect 146, 781, but most of these methods are empirical and require calibration. Jager et al. [851 present a different approach-the use of an automatic dilution unit. The challenge of this approach is maintaining a diluent of proper quality to assure crystal stability and a representative sample. 4.2.2. Transmittance. As discussed by Witkowski et al.[197], information about the CSD can be obtained from a measurement of the transmittance, the fraction of light that is transmitted by the crystal slurry: transmittance = I / I o where I is the intensity of undiffracted light that passes through the suspension of crystals and IO is the intensity of the incident light. Most particle sizers based on light scattering allow measurement of transmittance, but a relatively inexpensive spectrophotometer could also be used to make this measurement.
where k,is the surface area shape factor (theratio of surface area of a crystal of size L to a cube with sides of length L). Therefore, there is a one-to-one relationship between transmittance and the second moment of the CSD. The expressions given above assume that there is no multiple scattering. The rule-of-thumb given by van de Hulst 11851is that multiple scattering effects on scattering patterns can usually be ignored for T 1 < 0.3. The validity of the Beer-Lambert law for predicting transmittance for T 1 < 0.3 is confirmed by Miller [1141. 4.2.3. Digital Image Processing. Kane et al. [99] used a strobe and a camera to obtain images that, with considerable tedium, provided a CSD measurement. The recent advances in image processing and analysis algorithms and computer hardware may make such an approach feasible for on-line application. The advantages of direct observation techniques are that they do not suffer from the ill-conditioning of inference techniques, and assumption about particle shape are not necessary. The primary remaining challengesof imagingtechniques are obtaining quality images that are representative of the CSD and compensating for flaws in the images. Fantini et al. [541 give an algorithm that automatically handles images that are out of focus, a problem in analysis of microscopic images. They compare their results for characterizing droplets in a spray with those of light scattering techniques, Using digital image processing for particle sizing is still in the developmental stage, but there have been several demonstrations of the application of digital image analysis to determining particle size distributions 172, 1761. 5. Parameter Estimation
The most common method for simultaneously determining nucleation and growth kinetics involves operating a continuous mixed suspension, mixed product removal (MSMPR) crystallizer. The removal function for an MSMPR crystallizer is given by h(L) = 1. With this removal function and the assumption that G # G(L), the steady-state solution of eq 5 is f,(L) = exp(-L/G,T)(B,'/G,) (24) The plot of the logarithm of the CSD versus L should give a straight line, (25) log(f,) = log(B,"/G,) - L/G,T in which T = V/Q is the crystallizer residence time, and f,,B,' and G, are the steady-state CSD, nucleation rate, and growth rate, respectively. The growth rate can be
1284 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993
obtained from the slope of the semi-log plot, and the nucleation rate can then be calculated from the intercept. One can determine the orders of the nucleation and growth rates with respect to supersaturation by measuring the steady-state CSD at several different operating conditions, and making log-log plots of G,and B," versus S. This method of parameter estimation predominates the crystallization literature; Chen and Larson 1321 provide a representative sample of papers describing application of the method. Timm and Larson 11801 also use this method, but use the additional measurement of the slurry density to fix the intercept and determine Bso. Hartel et al. [731 use the measurement of the slurry density and the cumulative number density, the semilog plot of which should provide a straight-line relationship similar to that of eq 25, to determine G, and BEo. N*lt and hEek 11243 claim that the two most common sources of errors of the graphical, steady-state MSMPR method are uncertainties in the determination of the CSD from sieve data and the extrapolation of the CSD to L = 0. In other words, the data are always noisy and the intercept and, thus, the nucleation rate estimate are sensitive to this noise. They suggest that the use of cumulative weight distribution data will abate both of these difficulties. They show that an expression for the cumulative weight distribution can be obtained from eq 24; the term Gsr shows up as the only unknown in the cumulative weight distribution and can be estimated using nonlinear regression on sieve data. After obtaining Gar, Bso can be determined using a measurement of the slurry density in the following expression,
where M Tis the slurry density. Zumstein and Rousseau [200] also use this method. There are several disadvantages to determining kinetics with MSMPR experiments. First of all, the system must achieve steady state before the CSD is measured. This can be time consuming, especially if one obtains sufficient data to ensure statistical confidence in the parameters [I 711, and it is often not feasible for producers of fine and specialty chemicals due to lack of sufficient material to run continuous experiments. Second, it is sometimes difficult to achieve a true steady state. Crystallizers are often quite sensitive to small upsets and can actually be unstable and oscillate continuously. The stability issue is discussed further in section 6. Further problems are encountered with steady-state MSMPRexperiments when the semilog plot of the steadystate CSD is not a straight line. Significant upward curvature is often reported at the small crystal sizes. The upward curvature is a signal that the model structure is incorrect-the mixed product withdrawal assumption is not met, growth rate dispersion is occurring, and/or growth rate is size dependent. JanEiE and Garside [881 and Hartel et al. [731 provide attempts to extend the MSMPR technique to cases of size-dependent growth rate. As discussed in section 5.3, the information content of the data of an experiment is an important consideration in parameter estimation. Another disadvantage of the steady-state MSMPR method is that it throws away any dynamic data. Timm and Larson [I801 use transient responses to confirm kinetic characterization performed on steady-state MSMPR data of three different chemical systems and suggest the use of transient data for parameter estimation. Tavare [I781 outlines attempts to use the dynamic CSD measurement from a MSMPR crystallizer
in three simplified parameter estimation procedures, one of which is based on a Laplace transformation of the population balance. While the methods seemed successful when using simulated experimental observations, there was substantial scatter and disagreement between the different algorithms when actual data were used. After noting the sensitivity of the results to noise, Tavare suggests that MSMPR transient analysis may be used to confirm kinetics derived from steady-state analysis, but that additional work on experimental design would be needed for use in parameter estimation. Jager et al. 1863 investigate the Laplace transformation method and two other methods that have been suggested for parameter estimation from continuous crystallizer transient data. They show the ineffectiveness of these methods with the currently available measurement technology and demonstrate some success with an alternative approach-one based on nonlinear parameter estimation using light scattering measurements. Parameter estimation based on transient batch crystallization is an attractive alternative to the MSMPR method. Tavare 11791 summarizes the advantages of the batch crystallizer: a large number of operating variables can be studied in a relatively short time and a batch crystallizer can be successfully operated with a minimum development time and investment. Garside et al. 1621 and Palwe et al. [I251 propose using transient solute concentration data from a batch experiment in an "initial derivatives" technique to determine growth parameters. After assuming a simple power law expression for G, the kinetic parameters are expressed in terms of the first and second derivatives of the supersaturation-time curve at zero time. The derivatives are obtained from a polynomial fit to the data. Tavare [I 77J extends this approach to non-isothermal crystallizers. As Biegler et al. [I41 demonstrate in a case study of a chemical reactor, nonlinear parameter estimation is a challenging problem, but optimization schemes provide an intuitive approach that is general and that does not require model simplification. As discussed subsequently, nonlinear optimization can be applied successfully to determine crystallization kinetic parameters, but it requires careful experimental design and selection of objective function. Palwe et al. [I251 report that convergence is difficult to achieve for straightforward parameter estimation with a rigorous optimization of the s u m of the squared errors between predicted and measured supersaturation. It is not clear from their discussion why this method fails. Winter 11941 uses an optimization scheme with dynamic slurry density measurements from a batch crystallizer (obtained by analyzing the final product from runs that are replicates except for the batch duration time), but does not address the uncertainty in his optimal parameter values. While Qiu and Rasmuson [I311 do not give a great deal of detail about their estimation technique, they successfully use an optimization routine with measurements of concentration and final CSD to estimate growth and nucleation parameters and confidence intervals. Qiu and Rasmuson show that results simulated using the estimated parameters allowpredictions of weight mean sizes that agree with experimental data to within 25 % . Uncertainty bounds are not reported in most of the crystallization parameter estimation studies, which casts doubt on the validity of the kinetic parameter values. Tavare [I781 demonstrates that the scatter in parameter estimates using different estimation algorithms on the same set of data is significant, so it is not surprising that several different investigators have published quite dif-
Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 1286 ferent kinetic parameters for identical systems. Causes for this disagreement have been discussed by Jones [951, Qian et al. [l30], N*lt and Broul[1231, Budz [23], JanEit and Garside 1881, and Timm and Larson [1801. Subsequent sections demonstrate how uncertainties for typical crystallization kinetic studies can be large enough to account for much of the disagreement in previous studies. As shown in section 5.3, this problem can be overcome by carefully designing the experiments, including the choice of measurements. 5.1. Parameter Estimation Posed as an Optimization Problem. For general nonlinear parameter estimation problems, analysis techniques for assessingthe validity of the model structure and the parameter uncertainty are not global (i.e., values of the parameters are required). This section describes how the crystallization kinetics of a system can be efficiently determined by casting the estimation problem as a nonlinear optimization problem. Rippin [I531 demonstrates that the choice of the experimental configuration can strongly influence the information content of the data used to characterize a system, and he points out that it should be the first step in setting up a parameter estimation problem. Because of the dynamic nature of batch crystallization, much information about a system can be obtained from a single run and, as discussed previously, batch crystallization offers advantages of operating convenience over continuous crystallization. In addition, as will be demonstrated in section 5.3, solid-phase information is very important in parameter estimation; batch crystallization has the advantage that, assuming the seed CSD is known, a measure of the solute concentration gives a measure of the third moment of the CSD. Therefore, a batch crystallizer is chosen to illustrate the parameter estimation techniques. To make the discussion explicit, we list the modeling equations describing the batch crystallizer,
d&?
= -3p,k,GhKfL2 dL dt pc,,V-d T = -3AH~ck,VG$omfL2dL - UAc(T- T,) (29) dt where is the solute concentration on a per-mass-solvent basis and h is a conversion factor equal to the volume of slurry per mass of solvent. These equations can be derived from eqs 5 , 7 , and 8 by assuming that growth rate is size independent and that V and p are constant. The conversion of the basis for solute concentration is done to be consistent with convention for batch crystallizers. The boundary condition and initial conditions are as follows: f ( L , t )= B o / G , L = 0 f ( L , t )= f,(L), t = 0 C ( t ) ,=
c,,
(30)
Operating variables include the transient cooling profile and the initial conditions for the experiment-the initial concentration and the seed distribution. Possible measurements include &), the final CSD, and the dynamic CSD (or some function of the CSD, such as transmittance). Section 5.3 discusses the consequences on the estimation results of the assignment of operating variables and the choice of measurements, but assume for the moment that these decisions have been made. One then chooses an objective function that measures how well the model solution with a given set of parameter values predicts the experimental data. The least squares objective function is a common choice because it does not require consideration of the experimental data's error structure. As discussed by Witkowski et al. E 197J,weighted least squares allows the scaling of measurements that have different dimensions and weighting so that parameter estimates will be most influenced by the measurements felt to be most reliable; in this case, the objective function is
where 9ij and yij denote the ith measurement and model prediction, respectively, taken at the j t h sampling instant during the experiment, ai is the weighting assigned to the ith measured variable, N , is the number of measured variables (Le., concentration, transmittance, etc.), N d is the number of samples of each measured variable, and BT= [k,, kb, b,gl represents the unknown parameter vector. While the weights of eq 32 can be chosen to incorporate the data error structure (e.g., letting oi equal the inverse of the variance of the ith measurement), it would require replicate experiments to assess this statistical information. Therefore, other than accounting for scaling of the measurements, the weights are usually set somewhat arbitrarily. The maximum likelihoodmethod, on the other hand, enables the estimation of weights for the data along with the estimation of the parameters. If one assumes that the model structure is correct, then the deviations between the data and the model predictions are due to random (and possibly biased) errors in the measurements. The maximum likelihood method can be thought of as determining the parameters that maximize the probability that a given set of data could be obtained, considering the error structure of the measurements. Consider the following additional assumptions about the error structure of the measurements: 1. Errors at the j t h sampling instant are normally distributed with zero mean (no bias) and are homoscedastic (variance of errors is independent of the magnitude of the measured variable). 2. Errors at different sampling instants are uncorrelated. 3. Measured variables are independent (measurement covariance matrix is diagonal). With these assumptions, it can be shown 161 that the maximum likelihood method is equivalent to the minimization of
t =0
From section 2.2, reasonable choices for constitutive relations for operation in which the temperature variation is small are G =kPg Bo = kbsb;, (31) where P 3 is the third moment of the distribution in units of cmVg solvent.
where eij 9ij - yij. The method also yields the following estimate for the variance of the ith measured variable:
.
N.i
(34)
where 0' is the set of optimal parameters.
1286 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993
Whether the objective function is chosen for weight least squares or for maximum likelihood, the optimal parameter values are found by solving the nonlinear optimization problem, min W.W,W);@
e
subject to eqs 27-31. 5.2. Interpretation of Estimates. The estimation process can be thought of as a nonlinear operator that transforms a data set into an estimate of the parameters. Due to the random errors of the measurements, the parameters are also random variables with a probability distribution that depends on the nature of the estimator and on the distribution of the measurement errors; Bard [6]refers to the probability distribution of the parameters as the sampling distribution. If many replicates of the experiment used to estimate the parameters were available, the sampling distribution could be rigorously characterized. Since replicate data from crystallization experiments are difficult to obtain, an alternative is to perform a Monte Carlo study in which simulated replicate data are derived by adding random errors to the original data set in accordance with the assumed measurement probability distribution. However, Monte Carlo studies are too time consuming for many purposes as well, because the estimation process can be computationally intensive. An approximation of the confidence region can be obtained by assuming that the model can be represented by linear functions in the vicinity of the estimate, yj(e) = yj(e*) + Bj(e - e*)
(35)
where yy(f3) = [yu(e),...,YNA(e)lT, and Bj is the N , X Np matrix (whereNpis the number of parameters) given below Bj -
=P
These partial derivatives can be obtained by finite difference; however, since the model involves differential equations, they can be efficiently and accurately computed by integrating the sensitivity equations along with the model equations 16,271. With the linearized model and assumption that the measurement errors are normally distributed, the estimates are normally distributed and it can be shown that the quantity Nd
(0 - 8 * ) T ( ~ B ~ V - 1 B j-) B*) (B j=l
is distributed as x2 with Np degrees of freedom 161, in which V is the diagonal measurement covariance matrix with elements Vii = ui2 = Qi.2 Therefore, letting V8-I = BTV-lBj (V8 is the parameter covariance matrix for the linearized problem), and letting Pr(E) denote the probability of some event E,
CFl
Pr{(e- e*)Tv;l(e- e*) IxNP2(a)} = 1- a
(36)
In other words, the approximate l O O ( 1 - a)% confidence region from the linearization method is the ellipsoidal region defined by
(e - e*)Tv;l(e- e*) = xNP2(a)
(37)
This region can be characterized by the eigenvectors and eigenvalues of V0-l. The eigenvectors give the direction of the ellipsoid's axes, which have lengths that depend on
0.5-
__-_-_----__-_-, Figure 7. Confidence region demonstrating significant correlation between parameters. Conservativeconfidence interval for Bi is defmed by 8ei.
the eigenvalues-the ith semiaxis lies along the ith eigenvector of V8-I and has length [ x n $ 2 ( c ~ ) / X ~ l ~ where /~, hi is the ith eigenvalue. The shape and orientation of the ellipsoid indicate the degree of correlation between the parameters. As can be seen in Figure 7 from the elongated ellipse for a hypothetical two-parameter problem, 81 and 02 may be varied individually by only f0.5; however, if they are adjusted simultaneously, they can be varied significantly and still be in the confidence region. Since it is difficult to report succinctly a confidence region for cases in whichNp > 2, confidence intervals-PrI8 E[e' f Se(a)]j= (1- a)-are often desired. Conservative estimates of these can be obtained by locating the box in the parameter space that circumscribes the ellipsoid (cf. Figure 7). Bard [6]gives several alternative approximate techniques for obtaining approximate confidence regions and intervals, and Donaldson and Schnabel 1471 give an extensive discussion of confidence regions for nonlinear least squares problems. As illustrated by several Monte Carlo studies [4,47J, approximate techniques can overestimate confidence regions, but there is no guarantee that the true confidence region will not be grossly underestimated. Nevertheless, the approximate technique outlined above does give a convenient method for assessing parameter uncertainty that at least gives qualitative information about how well-determined the parameters are. 5.3. Experimental Design. The goal of the experimental design can now be stated precisely. One wishes to choose operating conditions and measurements such that the correlation between parameters and the volume of the uncertainty region are minimized. As discussed, the.fact that the crystallizer model is nonlinear in the parameters requires that the parameter estimates be obtained before the uncertainty can be assessed. A related complication exists for the experimental design problem-a model-based procedure requires the values of the parameters to minimize the parameter uncertainty. Whether the approach is trial-and-error or quantitative and model-based, experimental design is sequential in nature: collect data, estimate parameters and uncertainty region, design new experiments. 5.3.1. Measured Variables. Physical arguments provide some insight into the choice of measured variables. For example, intuition suggests that since the third
Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 1287 32.5
0.3 I 0.25
8
I
'I
0.2
,d
0.15
*
E
0.1 0.05 28.5
1
0
10
20
30 40 50 Time (minutes)
60
70
80
Figure 8. Temperature profile for batch crystallization of KNOs from water. 0.495 0.49
8 8
h
0.485
-
0.475
-
0.47
0.46
Ip
0.455
8
0.45
20
30
40 50 60 Time (minutes)
70
80
90
Figure 10. Transmittance measurements from experiment with temperature profile given in Figure 8.
I
0.48
0.465
'a
f-\ -
10
0
-
E
\
0.485 0.48 0.475 0.47
-
0.465
.
0.46
-
0.445 0.44 .
0.455
-
-
0.45
-
0.445
-
0.44 0
10
20
30
40
60
50
70
80
90
Time (minutes)
Figure 11. Fit to the data of the model with nominal parameters given in Table I. Table I. Parameter Estimation Results Using Concentration Measurements estimate interval
8.014 k0.423
1.10 k0.12
14.710 16.268
Table 11. Parameter Estimation Results Using Concentration and Transmittance Measurements ~ ~
estimate interval
W g )
8
h(kb)
8.167 k0.205
1.09 kO.06
15.779 11.475
0.85 12.48
~
_
_
b 1.39 i0.406
The concentration measurement does not contain sufficient information to determine the parameters uniquely. This result is the central point that is overlooked in much of the crystallization literature and is a possible explanation for the disagreement between various investigators' reported parameters. The practical significance is that crystallizer control based on the parameters identified using concentration measurements alone is not likely to achieve the product CSD that is desired. As discussed in section 4.2.2, there is a one-to-one relationship between transmittance and the second moment of the CSD. It is therefore reasonable to expect that the measurement of transmittance would provide independent information for the estimation process. A comparison of Tables I and I1 demonstrates that the transmittance measurement is in fact effective in reducing parameter uncertainty. The reduction in parameter uncertainty has important implications for controller design. Given the identified model in Table I, the controller design problem is difficult since one does not even know the sign of the nucleation order. In other words, the controller must take control action without knowing whether increasing supersatura-
_
_
1288 Ind. Eng. Chem. Res., Vol. 32,No. 7, 1993 0.5 32
0.49
31.5
0.48
31
0.47
30.5 30
0.46
29.5 0.45
29 28.5
1
0
10
20
40 50 Time (minutes)
30
60
70
I
0.44 0
80
Figure 12. Temperature profile for batch crystallization of KNOs from water.
estimate interval
8.918 10.127
1.32 10.03
ln(kd 17.207 *0.356
b 1.78 10.09
tion will increase or decrease the nucleation rate. This obviously extreme amount of model uncertainty is a consequence of the poor experimental design used for identification. The controller design problem using Table I1 is much simplified since the uncertainty is so small. It would be a much simpler matter to detune the nominal feedback controller to be robust to this small amount of uncertainty. 6.3.2. Operating Variables. The inspection of the kinetic expressions, eq 31,provides considerable insight for the choice of operating variables. For example, it is obvious that if one wishes to identify k b and b independently, then the variation in the supersaturation should be made as large as possible during the course of the experiment. The justification is simply that if the supersaturation is constant during the experiment, one can only estimate the product, k b Sb. Similarly, if temperature dependence is being characterized, the temperature variation should be made as large as possible. The parameter estimation results obtained in the previous section were based on a single batch experiment using a temperature profile with a rapid initial drop and a slow approach to the final temperature. It is possible that the choice of a different temperature profile would lead to improved parameter estimation (a smaller parameter uncertainty region). As discussed previously, the experimental design procedure for nonlinear problem is sequential. With the results obtained above, simulations could be used to determine a temperature profile that would be expected to reduce the parameter uncertainty. The process is actually often trial-and-error-analysis is performed using data from experiments with new operating conditions that were selected using intuition. For example, intuition suggests that additional information could be obtained by performing an experiment that is a replicate of the one described above, with the exception of the use of gradual cooling like that shown in Figure 12 instead of the rapid cooling shown in Figure 8. The results from parameter estimation using concentration and transmittance data from both experiments are given in Table 111. As expected, the additional information contained in the data of the new experiment leads to a decrease in the parameter uncertainty. It should be noted that the nominal parameters of Table I11 do not fall within the previous confidence region (given in Table 11);again, that
20
30
40 SO 60 Time (minutes)
70
80
90
Figure 13. 13: Fit to concentration data of model with nominal parameters given in Table 111.
Table 111. Parameter Estimation Results Using Concentration and Transmittance Measurements from Two Batch Runs R
10
0.3
I
I
0.25
P 'I c
I
0.2 0.1s 0.1
0.05
0
10
20
30
40
SO
60
70
80
90
Time (minutes)
Figure 14. Fit to transmittance data of model with nominal parameters given in Table 111. Table IV. Comparison of Linear Statistics with Monte Carlo Analysis of 600 Replicate Data Sets linear statistics confidence value, %
50 70 80 90 95 98
freq that param eat from replicate data are within elliptical confidence region, % 49 67 76 87 93 97
is the nature of experimental design for nonlinear problems. The fit of the model to the data is given in Figures 13 and 14. The seed load used in the experiments discussed above was quite light. This initial condition is another operating variable that could be used to affect the quality of parameter estimation. As pointed out in section 5.2, the validity of the approximate parameter confidence region obtained by linearizing the model and using linear statistics is problem specific but can be tested using a Monte Carlo study. The validity of the approximate confidence intervals for this model is confirmed by a Monte Carlo study in which replicates of the data used in the discussion above were produced. Table IV illustrates that the elliptical confidence regions from linear statistics adequately characterize the parameter uncertainty over a range of confidencevalues [115].
6. Control Issues
The control of batch and continuous industrial crystalkers recently has been reviewed by Rawlings et al. [1501. They discuss the following issues in detail: the
Ind. Eng. Chem. Res., Vol. 32, No. 7,1993 1289 *
254microns
I
v)
for cases in which the moments do not close and the model cannot be converted into a set of ODES 157,1981. As summarized by Randolph 11371,the two explanations that have been developed for CSD instability have been termed high-order cycling and low-order cycling. Highorder cycling corresponds to a high-order kinetic relationship between the nucleation rate and the supersaturation. Low-order cycling is due to nonrepresentative product removal (Le., f k not equal to fi. For an MSMPR crystallizer, the value of b required to cause cycling is much higher than has been estimated for most chemical systems. Therefore, low-order cycling is the more likely explanation for the commonly observed CSD instabilities. All of these stability analyses are local, and no one has determined Lyapunov functions to ascertain the global stability. There have not been any simulations of the nonlinear equations, however, that indicate the global stability is markedly different from the local stability. Numerous investigators have addressed the problem of stabilizingan oscillatingcrystallizerusing feedback control. Many of the early studies proposed measuring some moment of the CSD on line and using a proportional feedback controller to adjust the flowrate of a fines destruction system [70, 71, 1081. After reviewing these studies, Rousseau and Howell 11591 proposed using supersaturation as the measured variable, since accurate on-line measurements of properties of the CSD are difficult to make. They demonstrated through simulation that a proportional feedback controller could stabilize a cycling crystallizer, although the sensitivity to measurement error was more severe. Randolphandco-workers [7,11,137,140,143,145,1601 proposed inferring the nuclei density from a measurement of the CSD in the fines stream in conjunction with a proportional feedback controller. Significantly improved disturbance rejection is demonstrated experimentally for the proposed control scheme. Recently, Rohani [I551 reported preliminary tests on using a turbidimeter for an on-line measurement of slurry density in the fines loop. Again, stabilization and improved disturbance rejection were demonstrated by simulation using a proportional feedback controller to manipulate the fines removal rate. One of the primary motivations for feedback is to overcome model uncertainty. For the controller to work well in practice, it should be tuned to be robustly stabilizing. That is, the controller shouldnot only stabilize the nominal model, but also all models within some uncertainty region that reflects how well the system has been identified. The areas of robust stability and robust performance are current topics of control research. To date, essentially none of this research has been applied to crystallization. There have been relatively few studies and no experimental implementations of multivariable controllers for continuous crystallizers. The multiple-input, multipleoutput control systems that have been suggested are based on linear state-space models. A variety of methods to obtain these models have been put forth. Tsuruoka and Randolph [I821 use the method of lines to cast the population balance in state-space form, and de Wolf et al. [41] and Jager et al. 1871 demonstrate that time-series analysis can be used to obtain a state-space model. Hashemi and Epstein 1741 use the method of moments to transform the population balance into a set of ODES and then linearize them. They apply singular value decomposition analysis to develop measures of the system observability and controllability. More recently Eek et al. 1511 have studied the observability and controllability
::K, A,. , .
0
0
5
10
15
20
25
30
35
40
45
50
Time (hours)
Figure 15. Oscillationsin an industrialcrystallizer,after Randolph t1361.
interaction between crystallizer design and control, available measurements and manipulated variables, and guidelines for choice of sensors and final control elements. The state-of-the-art in practiced industrial crystallizer control is summarized for numerous types of common crystallization process designs. The focus of this review, therefore, is on the open research problems and areas in which research advances may lead to further improvement in industrial practice. 6.1. Continuous Crystallization. The need for largetonnage fertilizers created the demand for continuous crystallizers that could produce large and uniform crystals, and operate stably for long periods of time [9,113,154]. Although general estimates for the minimum production rate justifying continuous crystallization vary from 1ton/ day [I223 to 50 tons/day 1101, the scale of the equipment unavoidably creates nonhomogeneous regions of temperature, product concentration, and slurry density within the crystallizer at locations such as the heat exchanger or boiling surface, feed entrance, and in areas of lower mixing intensity or within baffled regions. The design intent is achieved if the vessel operation is controlled in a manner such that the vessel locations of these regions are time invariant and the process variables within the regions do not violate desired bounds. Operation outside these bounds generally results in fouling, deviations in crystal size, and lowered product purity. Additionally, the rate of change of the variables within acceptable bounds must be moderated to the extent that a higher supersaturation, or driving force, is not created than can be relieved by crystal growth and secondary nucleation. Such rapid changes would trigger spontaneous nucleation, producing a cyclic swing of smaller crystals for the downstream product recovery equipment to handle 11501. Sustained oscillations in continuous crystallizers are an important industrial problem and can lead to off-specification products, overload of dewatering equipment, and increased equipment fouling. Figure 15 shows this oscillatory behavior for an industrial KC1 crystallizer. Similar results have been observed in academic investigations with bench-scale crystallizers 1138, 1721. 6.1.1. Stability Issues. The stability of the open-loop, continuous crystallizer has received considerable attention 15, 7, 52, 84, 90, 109, 138, 139, 1421. The analytical techniques employed are either Laplace transform of the linearized population balance, conversion to a set of ordinary differential equations (ODES)by taking moments, or spectral analysis of the linear operator obtained from the population balance 11671. Recently analytical techniques have been developed to enable the calculation of the linearized integro-differential operator’sspectrum even
1290 Ind. Eng. Chem. Res., Vol. 32, No.7, 1993
issues using forward light scattering to measure the CSD and withdrawal rate of both fine and large crystals as manipulated variables. Myerson et al. 11201 outline a scheme for multivariable, optimal control of an MSMPR crystallizer. The suggested strategy uses a quadratic optimal controller to manipulate heat exchanger duty and feed flow rate. The controller is based on the linearized moment equations, and mass and energy balance equations, which are relinearized at each sampling period, and uses state estimates determined by an extended Kalman filter. On-line densitometry was proposed as a means of measuring the third moment of the CSD. 6.2. Batch Crystallization. Batch crystallizers are used extensively in the chemical industry, often for small volume, high value-added specialty chemicals. Batch operation offers a flexible and simple processing step for plants with frequently changing recipes and product lines. The control objectives for batch operation are often quite different from those of continuous crystallization. In the production of specialty organic materials, such as pharmaceuticals and photomaterials, product purity and CSD are of prime importance. If acceptable CSD and purity standards are not met, a batch of crystals has to undergo further processing steps, such as milling or recrystallization. Most of the previous process control studies of batch crystallizers have dealt with finding the open-loop, supersaturation versus time trajectory that optimizes some performance criterion derived from the CSD. The “natural” cooling curve displays the greatest cooling in the initial period when the temperature difference between the hot crystallizer fluid and the surroundings is the greatest. This causes a large degree of supersaturation at early times that leads to excessive nucleation, smaller final crystals, and aggravated operability problems such as fouling. Mullin and Nfvlt [I181 show that adjusting the supersaturation to maintain a constant level of nucleation would increase crystal size and decrease fouling of heat exchanger surfaces. Jones and Mullin 1971 calculate a cooling profile that maintains a constant supersaturation. They use simulations to conclude that this operating policy leads to a product CSD that is superior to that from natural or linear cooling policies. It was experimentally verified that the constant-supersaturation cooling policy improves both the mean size and the variance of the CSD compared with natural cooling. Ulrich [I841 tests various coolingprofiles on a system of an organic product in an acetonelwater solution and finds that a temperature curve qualitatively similar to that suggested by Jones and Mullin gives products with larger crystals and less variation than that corresponding to natural cooling or linear cooling. Ajinkya and Ray [31suggest that it would be desirable to determine the optimal operating policy that maximizes mean particle size and minimizes CSD variance in minimum time. Nevertheless, with the simplified model examined (assumptions included monodispersion of seeds and negligible nucleation), the only problem that could be addressed was that of manipulating cooling rate to maximizethe size of the crystals. Morari [I161determines the analytical solution of this problem and tests the algorithm using a simulation study, but the model used is unrealistically simple and no experimental verification was performed. Jones 1931 applies optimal control theory and control vector iteration to find the cooling curve that maximizes the final size of the seed crystals. The temperature of the size-optimal control profile decreases slowly at early times and then rapidly at the end of the batch time, in direct
32 .k
I
31.5 31
30.5 30
29.5 29 28.5
28 0
10
20
30
40
50
60
70
80
Time (minutes)
Figure 16. Temperature profiles for batch crystallizer. Cooling policies illustrated include natural cooling, constant-rate or linear cooling, and cooling controlled to maximize final-time seed size.
contrast to the natural cooling curve. Chang and Epstein [281 also compute optimal temperature profiles for a variety of objectivefunctions based on average size,volume, or variance of the final CSD. They use a gradient method to solve the two point boundary value problem for the optimal profile. Surprisingly, their results show an initial increase in supersaturation followed by a slow decrease, which is much closer to natural cooling than Jones’ result. Although these works employ slightly different objective functions and crystallizer models, the source of the disagreement is not immediately clear. This disagreement who provides a general is discussed further by Jones [94], review of attempts to control batch crystallizers. It should be noted that the open-loopoptimal controllers of Jones and of Chang and Epstein are based on the ODES that result from applying the method of momenta to the population balance, so the objective function for the optimization must be a function of the momenta of the CSD. The moment equations are a closed set only for the cases of no classified removal and for growth rates that are size-independent or linearly dependent on size; therefore, there are limitations of the applicability of these methods. As explained by Rawlings et al. [1501,the open-loop optimal control problem can be posed as a nonlinear programming problem. This formulation is flexible with respect to the crystallizer configuration, the form of the kinetic expression, and the choices of objective function and manipulated variables. I t also allowsthe incorporation of input, output, and final-time constraints (e.g., constraints on cooling rate, supersaturation level, and yield, respectively). To illustrate the method, the parameters given in Table I11and the corresponding model were used to compute the temperature profile for maximizing the terminal seed size. The temperature profile was parameterized as piecewise linear, and the constraints imposed included a minimum temperature bound and a constraint that the yield must be greater than or equal to that resulting from a linear cooling profile. The optimal cooling profile is shown in Figure 16, and the corresponding supersaturation profile is given in Figure 17. Optimal cooling leads to a predicted 9% increase in the terminal seed size over that of linear cooling, and a 18% increase over that of natural cooling. Jones et al. 1961 investigate the feasibility of using fines destruction in a batch crystallizer. While the energy costs of using fines destruction would have to be taken into account, this may provide an additional manipulated variable. Their experimental results indicate that fines destruction can remove the crystals formed by secondary
Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 1291 0.06 0.05
I
I-
0
-
,
optimal linear ---.-natural ........
'
10
20
30
40
50
60
70
80
Time (minutes)
Figure 17. Supersaturation profile corresponding to the cooling profiles in Figure 16.
nucleation, increase the average crystal size, and decrease the coefficient of variation of the CSD. Again, although we are not covering this area explicitly, one should note that the species feed flow rates can be used as manipulated variables in reactive crystallization, providing other handles for influencing supersaturation. Much insight is provided by computing the optimal open-loop cooling profiles for batch crystallizers. One can determine from such calculations general operating guidelines for producing CSDs with selected characteristics (cf. Figure 16). The optimal cooling profile by itself, however, is not sufficient to control the crystallizer. As demonstrated by Bohlin and Rasmuson 1151 via simulation and experimentally by Chianese et al. 1331,even the qualitative benefits expected from optimal cooling profiles may not be realized, at least not reproducibly. Poorly characterized kinetics and inaccurate process control are the probable causes of the shortcomings. Feedback is necessary to compensate for deviation from the nominal trajectory, to reject disturbances, and to account for the errors in the model on which the open-loop policy is based. Chang and Epstein [291 suggest a strategy for incorporating feedback into the optimal control scheme discussed above that is based on the assumption that the system state trajectory is in the vicinity of the nominal optimal trajectory. A one-step iteration in a gradient search technique is used to find a new trajectory with an improved performance index. Simulation studies were used to demonstrate the application of this control scheme, but it has not been experimentally verified, and it assumes that the moments of the CSD can be measured or estimated. A natural way to incorporate feedback into the controller is to periodically compute an optimal control trajectory after using available measurements to estimate the unmeasured states and possibly update the model's parameters. Eaton et al. 149,501 and Patwardhan et al. [I261 have applied on-line nonlinear estimation and control strategies to several chemical engineering processes. A more conventional approach that could be used as a basis of comparison is to use state estimation with a linear quadratic (LQ) optimal controller based the linearized ODES obtained after using some lumping technique on the crystallizer model (method of moments, method of lines, etc.). Kwakernaak and Sivan 11041 and Bryson and Ho [213 are standard references for LQ design. A feedback control scheme that is not rooted in optimal control theory has been suggested by Rohani et al. [1581. The strategy is based on the hypothesis that the final CSD can be improved by maintaining the fines slurry density at some constant value over a batch run. They experi-
mentally verify that by maintaining the fiies slurry density at a set point using a conventional PI controller to manipulate the fines destruction rate, the final mean crystal size is increased and the coefficient of variation is reduced. A simulation study by Rohani and Bourne 11561 shows a self-tuning regulator to be more effective in setpoint tracking and disturbance rejection than the PI controller. Another variable that can be used to affect the CSD is seeding. As shown by Bohlin and Rasmuson 1151, and Chianese et al. 1331, seeding increases average crystal size and enables more reproducible operation than unseeded crystallization. However, because of the inherent bimodal nature of seeded crystallizerproducts, seedingmay actually increase the coefficient of variation, depending on the cooling policy, the seed loading, and the crystallization kinetics. Crystallizers are usually seeded a t or very near the saturation temperature, but Karpinski et al. [I001 demonstrates experimentally that seeding can take place at significantly higher supersaturation levels without leading to lower product quality. The timing of the seeding of a batch sucrose crystallizer is one of the steps that is automated in a control scheme discussed by Virtanen 11883. The goals of automating the existing batch equipment were to shorten the crystallization time, produce a more uniform CSD, and reduce the labor and energy costs. In the control strategy, a microcomputer adjusts the steam flow rate to maintain constant supersaturation based on on-line density and refractive index measurements. The computer also regulates the timing of the feed addition, concentration, seeding, and growing steps for each batch. The control scheme was implemented at the Finnish Sugar Co. Ltd. and the biggest gains were achieved in decreasing batch to batch variations, decreasing the labor costs, and improving power plant economy. 7. Conclusions and Suggested F u t u r e Directions 7.1. Model Development a n d Solution. The population balance methodology is well established and successful for macroscopic modelingof the dynamics of crystal size distribution. Quantitative microscopic prediction of crystal shape and purity is not as well developed, but as discussed in section 2.1.2, rapid and exciting progress is being made. A feasible and worthwhile goal for future research is to bridge these two levels of models so that the macroscopic design and control decisions can be made with models that contain the best available knowledge of the microscopic chemistry. For engineering design and control purposes, the crystal nucleation and growth expressions are empirical and will likely remain so in the near future due to the complexity of these phenomena. Perhaps what is required in the near term for developing models with an intended control purpose is a more flexible empiricism. The use of flexible nonlinear functions for nucleation and growth, such as neural networks, radial basis functions, or multivariate splines, within the fundamental population balance framework to describe the CSD may prove to be a useful approach. The use of fundamental models in conjunction with flexible empirical models is a developing area of research 191, 103, 1291. The solution of these models for prediction and control in real time is feasible with today's computing capabilities. Although efficient numerical methods have been developed for particular crystallizer configurations, it would be reasonable to expect progress on the development of a unified approach to handle the entire class of problems,
1292 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993
much as a large class of ODES are efficiently handled by implicit predictor/corrector methods with variable step size. The benefit would be that effort currently being expended on case-by-casesolution methods could be spent on other issues, and the organization and dissemination of the software to researchers and industrial practitioners would be simplified and streamlined. 7.2. Measurement and Parameter Estimation. The single area in which improvements are in hand but not completely disseminated and used is parameter estimation. An effort should be made to use the latest and most powerful methods to perform parameter estimation because population balance models pose, in general, difficult parameter estimation problems. It is not obvious whether a set of data are sufficient to identify a crystallization model without carefully and systematically analyzing a nonlinear optimization problem. Investigators should report parameter uncertainty as well as best fit parameters. Without some measure of parameter uncertainty, it is essentially impossible for an investigator to deduce how informative someone else’s experimental data are. It is now possible to implement efficient experimental design. Models can be identified in some cases without time-consuming and expensive steady-state experiments. The choice of measurements is important in experimental design. It seems that laser light scattering does not as yet provide a reliable CSD measurement for the purposes of parameter estimation. The inherent assumptions of single scattering by optically smooth spheres seems to be the biggest problem with the technology, in spite of recent research efforts. The use of less informative but more reliable measurements, such as transmittance and turbidity, seems to be adequate for accurate parameter estimation. Direct digital imaging may prove to be the most effective technology for the future. Although the method is still problematic, progress is being made and the instrument cost is decreasing. 7.3. Control. Although considerable research effort has been devoted to CSD control, we are only now seeing the advances in measurement and computing technologies necessary for successful industrial implementation of the ideas. It is reasonable to expect closed-loop CSD control to become part of accepted industrial practice in the near future. In most process control applications, one tries to develop a model with minimal effort that is suitable for control purposes. The main reason for feedback control is, after all, the mitigation of model error and unmeasured disturbances in the closed loop. Hence the popularity of linear, blackbox models in process control. Crystallizers, however, have complex dynamics such as open-loop instability and long time delays caused by the interaction between growth of the CSD and new crystal nucleation. Common control objectives include optimal regulation of a batch crystallizer over a large operating region. It seems likely that low-order, linear models will prove inadequate in some of these applications and the nonlinear population balance models will be necessary. Research is needed to address in which applications these models are necessary and how they best can be used in the on-line control of crystallizers. Acknowledgment Financial support from Eastman Kodak Company and the National Science Foundation under Grant CTS8957123is gratefully acknowledged. The cylinder samples for Figure 6 were made available by Chester W. Sink of Tennessee Eastman Co.
Notation ai = coefficients in basis function expansion off A, = heat-transfer area for cooling crystallizer B = rate of crystal formation mechanisms Bj = matrix for linearized model at jth sampling instant Bo = rate of crystal nucleation at size LO C = solute concentration (mass of solute/volume of slurry) C = solute concentration (mass of solute/mass of solvent) CWt= solute saturation concentration cp = heat capacity D = rate of crystal removal mechanisms eij = error between measurement and prediction, 4 i j - yij e&) = energy of scattered light measured by detector of inner radius si E b = nucleation rate activation energy E, = growth rate activation energy f = crystal distribution function F = focal length of lens in particle sizer g = crystal growth rate order G = crystal growth rate h = conversion factor, volume of slurry per mass of solvent H = smoothing matrix AH = heat of crystallization I = intensity of scattered light IO = intensity of incident light j = slurry density power law coefficient for nucleation rate expression JO = Bessel function of the zeroth kind J1 = Bessel function of the first kind K = optical constant of particle sizer k, = crystal area shape factor kb = nucleation rate coefficient kBCF = Burton, Cabrera, and Frank growth rate parameter k’, = growth rate coefficient in the BCF model k, = crystal volume shape factor L = characteristic particle or crystal size m(L) = particle or crystal mass density distribution m = vector of mass density distribution values at quadrature locations MT = slurry density, pCkvpa Nd = number of measurementsamples during an experiment N , = number of quantities measured Np = number of parameters estimated P = system pressure Q = volumetric flow rate r = residual of the weighted residual method s = radius of annular scattered energy detector S = supersaturation, (C - CWt)/Cmt R = gas constant t = time T = temperature T, = coolant temperature U = heat-transfer coefficient V = volume of the crystallizer’s contents V = measurement covariance matrix VB= parameter covariancematrix for the linearized problem wi = ith residual weighting function or quadrature weight at the ith quadrature point x = vector of external coordinates y = vector of internal coordinates 9ij = ith measurement taken at the jth sampling instant yij = model prediction for the ith measurement taken at the jth sampling instant z = vector of particle’s coordinates Greek Letters a = joint confidence limit rn; = chi-square statistic with Npdegrees of freedom 6 = Dirac delta function y = penalty weight in constrained linear inversion e = volume fraction of continuous phase
Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 1293 0 = vector of parameters X = wavelength of light Xi = ith eigenvalue pk = kth moment of the CSD p = density of the crystallizer contents pc = density of the crystal solid ui2 = variance of ith measured variable 7 = turbidity or crystallizer residence time wi = weight for the ith measurement in the parameter
estimation objective function CP = parameter estimation objective function qi = ith basis function for the weighted residual method Subscripts k = kth flow stream 0 = initial conditions Literature Cited (1) Addadi, L.; Berkovitch-Yellin, Z.; Domb, N.; Gati, E.; Lahav, M.; Leiserowitz, L. Resolution of Conglomerates by Stereoselective Habit Modifications. Nature 1982,296,21-34. (2) Addadi, L.; Berkovitch-Yellin, Z.; Weissbuch, I.; Lahav, M.; Leiserowitz, L. Morphology Engineering of Organic Crystals with the Assistance of Tailor-Made Growth Inhibitors. Mol. Cryst. Liq. Cryst. 1983,96,1-17. (3) Ajinkya, M. B.; Ray, W. H. On the Optimal Operation of Crystallization Processes. Chem. Eng. Commun. 1974,l, 181-186. (4) Alper, J. S.;Gelb, R. I. Standard Errors and Confidence Intervals in Nonlinear Regression: Comparison of Monte Carlo and Parametric Statistics. J. Phys. Chem. 1990,94, 4747-4751. (5) Anahus,B. E.; Ruckenstein, E. On the Stability of a Well Mixed Isothermal Crystallizer. Chem. Eng. Sci. 1973,28,501-513. (6) Bard, Y. Nonlinear Parameter Estimation; Academic Press: New York, 1974. (7) Beckman, J. R.; Randolph, A. D. Crystal Size Distribution Dynamics in a Classified Crystallizer: Part 11. Simulated Control of Crystal Size Distribution. AIChE J. 1977,23 (4),51&520. (8) Behnken, D.; Horowitz, J.; Katz, S. Particle growth processes. Ind. Eng. Chem. Fundam. 1963,2 (3),212-216. (9) Bennett, R. C. Advances in Industrial Crystallization Techniques. Chem. Eng. Prog. 1984,89-95. (10) Bennett, R.C. Matching Crystallizer to Material. Chem. Eng. 1988,95 (8),118-127. (11) Bennett, R. C.; Randolph, A. D. Application of On-Line CSD Control to Industrial KC1 Crystallizers. In Potash 83; Mckercher, R. M., Ed.; Pergamon: Willowdale, Ontario, 1983. (12) Berglund, K. A.; Larson, M. A. Modeling of Growth Rate Dispersion of Citric Acid Monohydrate in Continuous Crystallizers. AIChE J. 1984,30,280-287. (13) Bhat, H. L.; Sherwood, J. N.; Shripathi, T. The Influence of Stress, Strain and Fracture of Crystals on the Crystal Growth Process. Chem. Eng. Sci. 1987,42 (4),609-618. (14)Biegler,L. T.; Damiano,J.J.; Blau, G. E. Nonlinear Parameter Estimation: a Case Study Comparison. AIChE J. 1986,32 (l), 2945. (15) Bohlin, M.; Rasmuson,A. C. Applicationof ControlledCooling and Seeding in Batch Crystallization. Can. J. Chem. Eng. 1992,70, 120-126. (16) Bohren, C. F.; Huffman, D. R. Absorption and Scattering of Light by Small Particles; Wiley: New York, 1983. (17) Botsaris, G. D. Secondary Nucleation-A Review. In Industrial Crystallization; Mullin, J. W., Ed.; Plenum Press: New York, 1976. (18)Boxman, A. Particle Size Measurement for Control of Industrial Crystallizers. Ph.D. Dissertation, Delft University of Technology, 1992. (19) Brown, D. J.; Felton, P. G. Direct Measurement of Concentration and Size for Particles of Different Shapes Using Laser Light Diffraction. Chem. Eng. Res. Des. 1985,63, 125-132. (20)Brown,D. J.; Gregory, T.; Weatherby,E. J. The Development of On-line Measurement of Crystal Growth by Laser Diffraction. In Control and Instrumentation-The Changing Scene, Proceedings of The Institute of Measurement and Control Conference; 1985. (21) Bryson, A. E.; Ho, Y. Applied Optimal Control; Hemisphere Publishing: New York, 1975.
(22) Budz, J.; Jones, A. G.; Mullin, J. W. On the Shape-Size Dependence of Potassium Sulfate Crystals. Ind. Eng. Chem. Res. 1987,26,82&824. (23) Budz, J. A. Kinetic Measurements in Crystallizing Systems. Seminar presented at Tenneasee Eastman Co., Kingsport, TN, 1987. (24) Bundle, L. G. CommercialInstrumentation for Particle Size Analysis. In Modern Methods of Particle Size Analysis; Barth, H. G., Ed.; Wiley: New York, 19W,pp 1-42. (25) Burton, W.; Cabrera, N.; Frank, F. The Growth of Crystals and the Equilibrium Structure of Their Surfaces. Philos. Trans. R. SOC.1961,243,299-358. (26) Canning, T. F.; Randolph, A. D. Some Aspects of Crystallization Theory: Systems that Violate McCabe’sDelta LLaw. AIChE J. 1967,13 (l), 5-10, (27) Caracotaios,M.; Stewart, W. E. Sensitivity Analysisof Initial Value Problems with Mixed ODEsand AlgebraicEquations. Comput. Chem. Eng. 1985,9 (4),359-365. (28) Chang, C.; Epstein, M. A. Identification of Batch Crystallization Control Strategies Using CharacterieticCurves.AIChE Symp. Ser. 1982,78 (215),68-75. (29)Chang, C.; Epstein, M. A. F. Simulation Studiesof a Feedback Control Strategy for Batch Crystallizers. AIChE Symp. Ser. 1987, 83 (253),110-119. (30) Chang, R.; Wang, M. Modeling the Batch Crystallization Process via Shifted Legendre Polynomials. Ind. Eng. Chem. Process Des. Dev. 1984,23,463-468. (31) Chang, R.;Wang, M. Shifted Legendre Function Approximation of Differential Equations: Application to Crystallization Processes. Comput. Chem. Eng. 1984,8(2),117-125. (32) Chen, M.; Larson, M. A. Crystallization Kinetics of Calcium Nitrate Tetrahydrate from MSMPR Crystallizer. Chem. Eng. Sci. 1985,40 (7),1287-1294. (33) Chianese, A.; Cave, S. D.; Mazzarotta, B. Investigation on Some Operating Factors Influencing Batch Cooling Crystallization. In Industrial Crystallization 84; JanEib, S . J., de Jong, E. J., Eds.; Elsevier: Amsterdam, 1984;pp 443-446. (34) Chow, L.; Tien, C. Inversion Techniques for Determining the Droplet Size Distribution in Clouds: Numerical Examination. Applied Opt. 1976,15 (2),378-383. (35) Clontz, N. A.; McCabe, W. L. Contact Nucleation of Magnesium Sulfate Heptahydrate. AIChE Symp. Ser. 1971,67 (110), 6-17. (36) Clydesdale, G.; Docherty, R.; Roberta, K. J. HABIT-A Program for Predicting the Morphology of Molecular Crystals. Comput. Phys. Commun. 1991,64,311-328. (37) Crump, J. G.; Seinfeld, J. H. A New Algorithm for Inversion of Aerosol Size Distribution Data. Aerosol Sci. and Technol. 1982, 1, 15-34. (38) Curl, R. Dispersed Phase Mixing: I. Theory and Effect in Simple Reactors. AIChE J. 1963,9,175. (39) Dave, J. V. Coefficients of the Legendre and Fourier Series for the Scattering Functions of Spherical Particles. Appl. Opt. 1970, 9 (8),1888-1896. (40) Davey, R.J.; Polywka, L. A,; Maginn, S. J. The Control of Crystal Morphologyby Additives: Molecular Recognition, Kinetics, and Technology. In Advances in Industrid Crystallization; Butterworth-Heinmann; 1991;pp 150-166. (41) de Wolf, S.;Jager, J.; Kramer, H. J. M.; Eek, R.; Bosgra, 0. H. Derivation of State-Space Models of Continuous Crystallizera. Presented at the IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, Maastricht, The Netherlands, 1989. (42) deBoer, G. B. J.; deWeerd, C.; Thoenes, D.; Goossens, H. W. J. Laser Diffraction Spectrometry: Fraunhofer Diffraction Versus Mie Scattering. Part. Charact, 1987,4,14-19. (43) Delves, L. M.; Mohamed, J. L. Computational Methods for Integral Equations. Cambridge University Press: Cambridge, UK, 1985. (44) Docherty, R.; Clydesdale, G.; Roberta, K.; Bennema, P. Application of Bravais-Friedel-Donnay-Harker Attachment Energy and Ising Models to Predicting and Understanding the Morphology of Molecular Crystals. J. Phys. D: Appl. Phys. 1991,24,89-99. (45) Docherty, R.; Roberta, K. J.; Dowty, E. MORANG-A Computer Program Designed to Aid in the Determinationa of Crystal Morphology. Comput. Phys. Commun. 1988,51,423-430. (46) Dodge, L. G. Change of Calibration of Diffraction-Based Particle Sizers in Dense Sprays. Opt. Eng. 1984,23 (5),626-630. (47) Donaldson, J. R.; Schnabel, R. B. Computational Experience with Confidence Regions and Confidence Intervals for Nonlinear Least Squares. Technometrics 1987,29 (l), 67-82.
1294 Ind. Eng. Chem. Res., Vol. 32, No.7, 1993 (48) Eakman, J. M.; Fredrickson, A. G.; Tsuchiya, H. M. Statistics andDynamicsof MicrobialCell Populations. Chem.Eng.Prog. Symp. Ser. 1966,62, 37-49. (49) Eaton, J. W.; Rawlings, J. B. Feedback Control of Chemical Processes Using On-line Optimization Techniques. Comput. Chem. Eng. 1990, 14 (4/5), 469-479. (50) Eaton, J. W.; Rawlings, J. B.; Edgar, T. F. Model-Predictive Controland SensitivityAnalysis for Constrained Nonlinear Processes. In Proceedings of the IFAC Workshop on Model Based Process Control; 1988. (51) Eek, R. A.; Boxman, A.; Dijkstra, S. Observability and Controllability Aspects of Continuous Industrial Crystallisers. Automatica 1992, in press. (52) Epstein, M. A. F.; Sowul, L. Phase Space Analysis of Limit Cycle Development in CMSMPR Crystallizers Using Three-Dimensional Computer Graphics. AIChE Symp. Ser. 1980,76 (193),6-17. (53) Fan, L. T.; Chou, S. T.; Hsu, J. P. Transient Analysis of Crystallization: Effect of the Size-Dependent Residence Time of ClassifiedProduct Removal. AIChE Symp. Ser. 1987,83 (253), 120129. (54) Fantini, E.; Tognotti, L.; Tonazzini,A. Drop Size Distribution in Sprays by Image Processing. Cornput. Chem. Eng. 1990,14 (ll), 1201-1211. (55) Felton, P. G.; Brown, D. J. Measurement of Crystal Growth by Laser Light Diffraction. Inst. Chem. Eng. Symp. Ser. 1980,59, 7:1/ 1-7:1/ 9. (56) Friedlander, S. K. Smoke, Dust, and Haze: Fundamentals of Aerosol Behavior; Wiley: New York, 1977. (57) Fuentes, Y. 0.;Rawlings, J. B. Stability and Bifurcation of Population Balance Models. Submitted for publication in Chem. Eng. Sci. 1993. (58) Fymat, A. L.; M e w , K. D. Mie Forward Scattering: Improved Semiempirical Approximation with Application to Particle Size Distribution Inversion. Appl. Opt. 1981, 20, 194-198. (59) Garside, J. Advances in the Characterisation of Crystal Growth. AIChE Symp. Ser. 1984, 80 (240), 23-38. (60) Garside, J. Industrial Crystallization From Solution. Chem. Eng. Sci. 1985, 40 (l),3-26. (61) Garside, J.; Davey, R. Secondary Contact Nucleation: Kinetics, Growth and Scale-up. Chem. Eng. Commun. 1980,4393424. (62) Garside, J.;Gibilaro, L. G.;Tavare, N. S. Evaluation of Crystal Growth Kinetics from Desupersaturation Curve Using Initial Time Derivatives. Chem. Eng. Sci. 1982, 37, 1625-1628. (63) Garside, J.; Mullin, J. W. Continuous Measurement of Solution Concentration in a Crystalliser. Chem. Ind. 1966, 20072008. (64) Garside, J.; Phillips, V. R.; Shah, M. B. On Size-Dependent Crystal Growth. Znd. Eng. Chem. Fundam. 1976,15 (3), 230-233. (65) Garside, J.; Shah, M. B. Crystallization Kinetics from MSMPR Crystallizers. Ind. Eng. Chem. Process Des. Deu. 1980,19, 509-514. (66) Garside, J.; Webster, G.; Davey, R. J.; Ruddick, A. J. The Influence of Stress, Strain and Fracture of Crystals on the Crystal Growth Process. In Industrial Crystallization 84; JanEiE, S. J., de Jong, E. J., Eds.; Elsevier: Amsterdam, 1984; pp 459-462. (67) Gelbard, F.; Seinfeld,J. H. NumericalSolution of the Dynamic Equation for Particulate Systems. J. Comput. Phys. 1978,28,357375. (68) Girolami, M. W.; Rousseau, R. W. Size-Dependent Crystal Growth-A Manifestation of Growth Rate Dispersion in the Potassium Alum-Water System. AIChE J. 1985, 31 (ll),1821-1828. (69) Grootscholten, P. A. M.; de Leer, B. G.M.; de Jong, E. J.; Asselbergs, C. Factors Affecting Secondary Nucleation Rate of Sodium Choride in an Evaporative Crystallizer. AIChE J. 1982,28, 728-737. (70) Gupta, G.; Timm, D. C. PredictiveKorrective Control for Continuous Crystallization. Chem. Eng. Prog. Symp. Ser. 1971,67 (110), 121-128. (71) Han,C.D.AControlStudyonIsothermalMixedCrystallizers. Ind. Eng. Chem. Process Des. Dev. 1969,8 (2), 150-158. (72) Hanzevack, E. L.; McNeill, S. R.; Bowers, C. B.; Ju, C. H. An Inexpensive Experimental System for Study of Two Phase Flow by Laser Image Processing. Chem. Eng. Commun. 1988,65, 161-167. (73) Hartel, R. W.; Berglund, K. A,; Gwynn, S. M.; Schierholz, P. M.; Murphy, V. G. Crystallization Kinetics for the Sucrose-Water System. AIChE Symp. Ser. 1980, 76 (193), 65-72. (74) Hashemi, R.; Epstein, M. A. Obaervabilityand Controllability Considerations in Crystallization Process Design. AICh E Symp. Ser. 1982, 78 (215), 81-90.
(75) Helt, J. E.; Larson, M. A. Effect of Temperature on the crystallization of Potassium Nitrate by Direct Measurement of Supersaturation. AIChE J. 1977,23 (6), 822-830. (76) Hidy, G. Aerosols: A n Industrial and Environmental Science; Pergamon: Orlando, 1984. (77) Hidy, G. N.; Brock, J. R. The Dynamics of Aerocolloidal Systems; Pergamon: New York, 1970. (78) Hirleman, E. D. Modeling of Multiple Scattering Effects in Fraunhofer Diffraction Particle Size Analysis. Part. Charact. 1988, 5, 57-65. (79) Hoerl, A. E.; Kennard, R. W. Ridge Regression-1980: Advances, Algorithms, and Applications. Am. J. Math. Man.Sci. 1981, 1, 5-83. (80) Hsu, J.; Fan, L. Transient Analysis of Crystallization: Effect of the Size Dependent Growth Rate. Chem. Eng. Commun. 1987,56, 19-40. (81) Hsu, J. P.; Fan, L. T.; Chou, S. T. Transient Analysis of Crystallization: Effect of the Initial Size Distribution. Chem. Eng. Commun. 1988,69, 95-114. (82) Hulburt, H.; Katz, S. Some Problems in Particle Technology. Chem. Eng. Sci. 1964,19,555-574. (83) Hwang,C.; Shih,Y. Solutionof PopulationBalance Equations via Block Pulse Functions. Chem. Eng. J. 1982,2.5,39-45. (84) Ishii, T.; Randolph, A. D. Stability of the High Yield MSMPR Crystallizer with Size-Dependent Growth Rate. AIChE J. 1980,26, 507-510. (85) Jager, J.; de Wolf, S.; Klapwijk, W.; de Jong, E. J. A New Design for On-Line Product-Slurry Measurements. In Industrial Crystallization 87; N$vlt, J., ZdEek, S., Eds.; Elsevier: Bechyn6, Czechoslovakia, 1987; pp 415-418. (86) Jager, J.; de Wolf, S.; Kramer, H. J. M.; de Jong, E. J. Estimation of Nucleation Kinetics from Crystal Size Distribution Transienta of a Continuous Crystallizer. Chern. Eng. Sci. 1991,46 (3), 807-818. (87) Jager,J.;Kramer,H.J.M.;deJong,E.J.;deWolf,S.;Bosgra, 0. H.; Boxman, A.; Merkus, H. G.; Scarlett, B. Control of Industrial Crystallizers. Powder Technol. 1992,69 (l),11-20. (88) JanEiE, S.; Garside, J. A New Technique for Accurate Crystal Size Distribution Analysis in an MSMPR Crystallizer. In Industrial Crystallization; Mullin, J. W., Ed.; Plenum Press: New York, 1976. (89) Janse, A. H.; de Jong, E. J. The Occurrence of Growth Dispersion and Its Consequences. In Industrial Crystallization; Mullin, J. W., Ed.; Plenum Press: New York, 1976. (90) Jerauld, G. R.; Vasatis, Y.; Doherty, M. F. Simple Conditions for the Appearance of Sustained Oscillations in Continuous Crystallizers. Chem. Eng. Sci. 1983,38, 1673-1681. (91) Johansen, T. A.; Foss, B. A. Representing and Learning Unmodelled Dynamics with Neural Network Memories. In Proceedings of the 1992 American Control Conference; 1992; pp 30373043. (92) Johnson, R. T.; Rousseau, R. W.; McCabe, W. L. Factors affecting contact nucleation. AIChE Symp. Ser. 1972,68 (121), 3141. (93) Jones, A. G. Optimal Operation of a Batch Cooling Crystallizer. Chem. Eng. Sci. 1974, 29, 1075-1087. (94) Jones, A. G. Design and Performance of Batch Crystallizers. In Institutionof ChemicalEngineers: North WesternBranchPapers NO.2; 1982; pp 4.1-4.11. (95) Jones, A. G.; Budz, J.; Mullin, J. W. Crystallization Kinetics of Potassium Sulfate in an MSMPR Agitated Vessel. AIChE J. 1986, 32 (12), 2002-2008. (96) Jones, A. G.; Chianese, A.; Mullin, J. W. Effect of Fines Destruction on Batch Cooling Crystallization of Potassium Sulphate Solutions. In Industrial Crystallization 84; JanEiE, S. J., de Jong, E. J., Eda.; Elsevier : Amsterdam, 1984; pp 191-195. (97) Jones, A. G.; Mullin, J. W. Programmed Cooling Crystallization of Potassium Sulphate Solutions. Chem. Eng. Sci. 1974,29, 105-118. (98) Jones, A. R. Error Contour Charts Relevant toparticle Sizing by Forward-scattered Lobe Methods. J.Phys. D: Appl. Phys. 1977, 10, L163-L165. (99) Kane, S. G.; Evans, T. W.; Brian, P. L. T.; Sarofim, A. F. Determination of the Kinetics of Secondary Nucleation in Batch Crystallizers. AIChE J. 1974,20 (5), 855-862. (100) Karpinski, P.; Budz, J.; Naruc, Z. Effect of Seeding Moment on Quality of CSD from Batch CoolingCrystallizer.In Science Papers of the Institute of Chemical Engineeringand Heat Systems; Belina, D., Ed.; Technical Universityof Wroclaw,Wroclaw: Wroclaw, Poland, 1980; Book No. 38, Symposium No. 5, pp 172-179.
Ind. Eng. Chem. Res., Voi. 32,No. 7, 1993 1295 (101) Karpineki, P. H. Crystallization as a Mass Transfer Phenomenon. Chem. Eng. Sci. 1980,35,2321-2324. (102) Kostova, T. V. Numerical Solutions of a Hyperbolic Differential-Integral Equation. Comput. Math. Appl. 1988,15 (68),427-436. (103) Kramer, M. A.;Thompson, M. L.; Bhagat, P. M. Embedding Theoretical Models in Neural Networks. In Proceedings of the 1992 American Control Conference; 1992;pp 475-479. (104)Kwakernaak,H.; Sivan,R.Linear Optimal Control Systems; Wiley: New York, 1972. (105) Lakatos, B.; Varga, E.; Hal& S.; Blickle, T. Simulation of Batch Crystallizers. In Industrial Crystallization 84;JanEiC, S . J., de Jong, E. J., Eds.; Elsevier: Amsterdam, 1984; pp 185-190. (106) Landau, L.D.; Lifshitz, E. M. Fluid Mechanics; Pergamon: New York, 1959. (107) Larson, M. A. Advances in the Characterization of Crystal Nucleation. AZChE Symp. Ser. 1984,80 (240),39-44. (108) Lei, S. J.; Shinnar, R.; Katz, S. The Regulation of a Continuous Crystallizer with Fines Trap. Chem. Eng. Prog. Symp. Ser. 1971,67(110),129-144. (109) Lei, S.-J.; Shinnar, R.; Katz, S. The Stability and Dynamic Behavior of a Continuous Crystallizer with a Fines Trap. AZChE J. 1971,17 (6),1459-1470. (110) Lieberman, A. On-line Particle Counting for Crystallizer Operation Control. AZChE Symp. Ser. 1982, 78 (215),76-80. (111) Meenan, P.; Roberts, K. J.; Sherwood, J. N.; Yuregir, K. R. Understanding and Controlling the Crystal Morphology of Some Ionic Crystals. Powder Technol. 1991,65,219-225. (112) Mie, G. Beitrage zur Optik Ober Medien speziell kolloidaler Metalltisungen. Ann. Phys. 1908,25,377-445. (113) Miller, P.; Saeman, W. Continuous Vacuum Crystallization of Ammonium Nitrate. Chem. Eng. Prog. 1947,43 (12),667-690. (114) Miller, S.M. Modelling and Quality Control Strategies for Batch Cooling Crystallizers. Ph.D. Dissertation, The University of Texas at Austin, 1993. (115)Miller, S.M.; Rawlings, J. B. Control Strategies for Batch Cooling Crystallizers. Presented at the AIChE National Meeting, Miami, FL, 1992. (116) Morari, M. Some Comments on the Optimal Operation of Batch Crystallizers. Chem. Eng. Commun. 1980,4,167-171. (117) Mullin, J. W.; Leci, C. J. Desupersaturation of Seeded Citric Acid Solutions in a Stirred Vessel. AZChE Symp. Ser. 1972,68(121), 8-20. (118) Mullin, J. W.; N h l t , J. Programmed Cooling of Batch Crystallizers. Chem. Eng. Sci. 1971,26,369-377. (119) Myerson, A. S.Crystallization Research in the 1990s: An Overview. In Crystallization as a Separations Process; Myerson, A. S., Toyokura, K., Eds.; American Chemical Society: Washington, DC, 1990,pp 1-15. (120) Myerson,A. S.;Rush, S.; Schork, F. J.; Johnson, J. L. Control of Crystallization Processes. In Zndustrial Crystallization 87; Nfvlt, J., fXEek, S., Eds.; Elsevier: Amsterdam, 1987;pp 407-410. (121) Myerson, A. S.;Saska, M. Calculation of Crystal Habit and Solvent-Accessible Areas of Sucrose and Adipic Acid Crystals. In Crystallization as a Separations Process;Myerson, A. S., Toyokura, K., Eds.; American Chemical Society: Washington, DC, 1990; pp. 55-71. (122) Nfvlt, J. Industrial Crystallization: The Present State of the Art; Verlag Chemie: Weinheim, 1978. (123) Nfvlt, J.; Broul, M. Kinetic Exponents in Crystallization and the Accuracy of Their Determination. Znt. Chem. Eng. 1982,22, 543-548. (124) Njvlt, J.; & e k , S. Possible Inaccuracies Involved in the Crystal Population Balance Method and Their Elimination. Cryst. Res. Technol. 1981,16 (7),807-814. (125)Palwe, B.; Chivate, M.; Tavare, N. Growth Kinetics of Ammonium Nitrate Crystals in a Draft Tube Baffled Batch Crystallizer. Znd. Eng. Chem. Process Des. Dev.1985,24,914-919. (126) Patwardhan, A. A.;Rawlings, J. B.; Edgar, T. F. Nonlinear Model Predictive Control. Chem. Eng. Commun. 1990,87,123-141. (127) Peterson, T. W.; Gelbard, F.; Seinfeld, J. H. Dynamics of Source-Reinforced,Coagulating,and CondensingAerosols.J.Colloid Interface Sci. 1978,63 (3),426-445. (128) Phillips, D. L. A Technique for the Numerical Solution of Certain Integral Equations of the First Kind. J. Assoc. Comput. Mach. 1962,9,84-97. (129)Psichogios, D. C.; Ungar, L. H. A Hybrid Neural NetworkFirst Principles Approach to Process Modeling. AZChE J. 1992,38 (lo), 1499-1511.
(130) Qian, R.-Y.; Chen, Z.; Ni, H.; Fan, Z.; Cai, F. Crystallization Kinetics of Potassium Chloride from Brine and Scale-up Criterion. AZChE J. 1987,33(lo),1690-1697. (131) Qiu,Y.;Rasmuson,A. C.NucleationandGrowthofSuccinic Acid in a Batch Cooling Crystallizer. AIChE J. 1991,37(9),12931304. (132) Ramkrishna, D. Solution of Population Balance Equations. Chem. Eng. Sci. 1971,26,1136-1139. (133)Ramkrishna, D. On Problem-Specific Polynomials. Chem. Eng. Sci. 1973,28,1362-1365. (134) Ramkrishna, D. Statistical Models of Cell Populations. In Advances in Biochemical Engineering; Ghose, T. K., Fiechter, A., Blakebrough, N., Eds.; 1979; p 1. (135) Ramkrishna, D. The Status of Population Balances. Rev. Chem. Eng. 1985,3 (l), 49-95. (136)Randolph, A. D. Comments on Recent Advances in the Analysis of Crystallization Processes. Chem. Eng. Prog. Symp. Ser. 1971,67(110),1-5. (137) Randolph, A. D. CSD Dynamics, Stability, and Control (A Review Paper). AZChE Symp. Ser. 1980,76 (193),1-5. (138) Randolph, A. D.; Beckman, J. R.; Kraljevich, Z. I. Crystal Size Distribution Dynamics in a Classified Crystallizer: Part I. Experimental and Theoretical Study of Cycling in a Potassium Chloride Crystallizer. AZChE J. 1977,23,500-509. (139)Randolph, A. D.; Beer, G. L.; Keener, J. P. Stability of the Class I1 Classified Product Crystallizer with Fines Removal. AZChE J. 1973,19 (6),1140-1148. (140) Randolph, A. D.; Chen, L.; Tavana, A. Feedback Control of CSD in a KC1 Crystallizer with a Fines Dissolver. AZChE J. 1987,33, 583-591. (141)Randolph, A. D.; Larson, M. A. Transient and Steady State Size Distributions in Continuous Mixed Suspension Crystallizers. AZChE J. 1962,8 (5),639-645. (142) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes, 2nd ed.; Academic Press: San Diego, 1988. (143) Randolph, A. D.; Low, C. C. D. Some Attempts at CSD Control Utilizing On-line Measurement of Nucleation Rate. In Industrial Crystallization 81;JanEiC, S . J., de Jong, E. J., Eds.; North-Holland: Amsterdam, 1981;pp 29-34. (144) Randolph, A. D.; White, E. T. Modeling Size Dispersion in the Prediction of Crystal-Size Distribution. Chem. Eng. Sci. 1977, 32,1067-1076. (145) Randolph, A. D.; White, E. T.; Low, C.-C. D. On-line Measurement of Fine-Crystal Response to Crystallizer Disturbances. Znd. Eng. Chem. Process Des. Dev. 1981,20,496-503. (146)Rawlings,J. B.;Ray, W. H. EmulsionPolymerizationReactor Stability: Simplified Model Analysis. AIChE J. 1987,33(lo), 16631677. (147) Rawlings,J. B.; Ray, W. H. Stability of Continuous Emulsion Polymerization Reactors: a Detailed Model Analysis. Chem. Eng. Sci. 1987,42(ll), 2767-2777. (148) Rawlings, J. B.; Ray, W. H. The Modelling of Batch and Continuous Emulsion Polymerization Reactors. Part I Model Formulation and Sensitivity to Parameters. Polym. Eng. Sci. 1988, 28,237-256. (149) Rawlings, J. B.; Ray, W. H. The Modelling of Batch and Continuous Emulsion Polymerization Reactors. Part I1 Comparison With Experimental Data From Continuous Stirred Tank Reactors. Polym. Eng. Sci. 1988,28,257-273. (150) Rawlings, J. B.; Sink, C. W.; Miller, S. M. Control of Crystallization Processes. In Industrial Crystallization-Theory and Practice; Myerson, A. S., Ed.; Butterworth Boston, 1992;pp 179-207. (151) Riebel, U.; Kofler, V.; Loffler, F. Control of Supersaturation in Instationary Suspension Crystallization. In Industrial Crystallization 90,Mersmann, A., Ed.; 1990; pp 595-599. (152) Ring, T. A. Kinetic Effects on Particle Morphology and Size Distribution During Batch Precipitation. Powder Technol. 1991, 65,195-206. (153) Rippin, D. W. T. Statistical Methods for Experimental Planning in Chemical Engineering. Comput. Chem. Eng. 1988,12, 109-116. (154) Robinson, J. N.; Roberts, J. E. A Mathematical Study of Crystal Growth in a Cascade of Agitators. Can. J. Chem. Eng. 1957, 35, 105-112. (155)Rohani, S. Dynamic Study and Control of Crystal Size Distribution (CSD) in a KC1 Crystallizer. Can. J . Chem. Eng. 1986, 64, 112-116.
1296 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 (156) Rohani, S.;Bourne, J. R. Self-tuning Control of Crystal Size Distribution in a Cooling Batch Crystallizer. Chem. Eng. Sci. 1990, 45 (12),3457-3466. (157) Rohani, S.;Paine, K. Measurement of Solids Concentration of a Soluble Compound in a Saturated Slurry. Can. J. Chem. Eng. 1987,65,163-165. (158) Rohani, S.;Tavare, N. S.; Garside, J. Control of Crystal Size Distribution in a Batch Cooling Crystallizer. Can. J. Chem. Eng. 1990,68,260-267. (159) Rousseau, R. W.; Howell, T. R. Comparison of Simulated Crystal Size Distribution Control Systems Based on Nuclei Density and Supersaturation. Znd. Eng. Chem. Process Des. Dev. 1982,21, 606-610. (160) Rovang, R. D.; Randolph, A. D. On-Line Particle Size Analysis in the Fines LOOPof a KCl Crystallizer. AZChE Symp. . Ser. 1980; 76 (193),798-807. (161)Rudd.D.F.AGeneralizationoftheReaidenceTimeConceDt. Can. J: Chem.’Eng. 1962,40,197-202. (162) Russell, R. D.; Christiansen, J. Adaptive Mesh Selection Strategies for Solving Boundary Value Problems. SZAM J. Numer. Anal. 1978,15 (l), 59-80. (163) Sampson, K. J.; Ramkrishna, D. A New Solution to the Brownian Coagulation Equation through the Use of Root-Shifted Problem-Specific Polynomials. J. Colloid Interface sci. 1985,103 (l),245-254. (164) Sastrv. K. V. S.: Gaschienard. P. Discretization Procedure for’the’ Coalescence Equation of Particulate Processes. Znd. Eng. Chem. Fundam. 1981,20,355-361. (165) Schulz, V. J.; Baumann, K. H.; Weiss, S. Kristallwachstum und Keimbildung bei der Kristallisation von Harnstoff und Kaliumsulfat. Wiss. 2. Tech. Hochsch. “Carl Schorlemmer” LeunaMerseburg 1982,24,500-515. (166) Shaw, G. E. Inversion of Optical Scattering and Spectral Extinction Measurements to Recover Aerosol Size Spectra. Appl. Opt. 1979,18 (7),988-993. (167) Sherwin, M. B.; Shinnar, R.; Katz, S. Dynamic Behavior of the Well-Mixed Isothermal Crystallizer. AZChE J. 1967,13,11411153. (168) Sikdar, S. K.; Randolph, A. D. Secondary Nucleation of Two Fast Growth Systems in a Mixed Suspension Crystallizer: Magnesium Sulfate and Citric Acid Water Systems. AZChE J. 1976, 22 (l),110-117. (169) Singh, P. N.; Ramkrishna, D. Transient Solution of the Brownian Coagulation Equation by Problem-Specific Polynomials. J. Colloid Interface Sci. 1975,53(2),214-223. (170) Singh, P.N.;Ramkrishna, D. Solution of Population Balance Equations by MWR. Comput. Chem. Eng. 1977,I, 23-31. (171) Sink, C. W. Eastman Chemicals Division, Tennessee Eastman Co. Personal communication, 1987. (172) Song, Y. H.;Douglas, J. M. Self-Generated Oscillations in Continuous Crystallizers: Part 11. An Experimental Study of an Isothermal System. AZChE J. 1975,21,924-930. (173) Steemson, M.; White, E. Numerical Modelling of Steady State Continuous Crystallization Processes Using Piecewise Cubic Spline Functions. Comput. Chem. Eng. 1988,12 (l),81-89. (174) Strickland-Constable, R. F. The Breeding of Crystal Nuclei-A Review of the Subject. AZChE Symp. Ser. 1972,68(121), 1-7. (175)Swithenbank, J.; Beer, J. M.; Taylor, D. S.; Abbot, D.; McCreath, G. C. A Laser Diagnostic Technique for the Measurement of Droplet and Particle Size Distribution. In Erperimental Diagnostics in Gas Phase Combustion Systems; Zinn, B. T., Ed.;Progress in Astronautics and Aeronautics; A M Washington, DC, 1976;pp 421-447. (176) Tanaka, H.; Hayashi, T.; Nishi, T. Application of Digital Image Analysis to Pattern Formation in Polymer Systems. J.Appl. Phys. 1986,59 (111,3627-3643. (177) Tavare, N. S. Growth Kinetics of Ammonium Sulfate in a Batch CoolingCrystallizer Using Initial Derivatives. AZChE J. 1985, 31 (lo), 1733-1735.
(178)Tavare, N. S.Crystallization Kinetica From Transients of an MSMPR Crystallizer. Can. J. Chem. Eng. 1986,64,752-758. (179) Tavare, N. S.Batch Crystallizers: A Review. Chem. Eng. Commun. 1987,61,259-318. (180) Timm,D. C.; Larson, M. A. Effect of Nucleation Kinetics on the Dynamic Behavior of a Continuous Crystallizer. AIChE J. 1968,14,462-457. (181) Tauchiya, H.M.; Fredrickson, A. G.; Ark, R. Dynamica of Microbial Cell Populations. Adu. Chem. Eng. 1966,6,125-206. (182)Tsuruoka, S.;Randolph, A. D. State Space Representation of the Dynamic crystallizer Population Balance Application to CSD Controller Design. AIChE Symp. Ser. 1987,83 (253), 104-109. (183) Twomey,S.On the NumericalSolutionof FredholmIntegral Equations of the First Kind by the Inversion of the Linear Syetem Produced by Quadrature. J. Assoc. Comput. Mach. 1963, 10,97101. (184) Ulrich, M. Optimization of Batch Solution Crystallization. Ger. Chem. Eng. 1979,4,195-200. (185) vande Hulst, H. Light Scattering by SmallParticles;Dover Publications: New York, 1981. (186)Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Modeb by Polynomial Approximation; Prentice-Hall: Englewood Cliffs, NJ, 1978. (187)Villadsen, J. V.; Stewart, W. E. Solutionof Boundary-Value Problems by Orthogonal Collocation. Chem. Eng. Sci. 1967,22,14831501. (188) Virtanen, J. MicrocomputerControl of SugarCrystallization. Sugar J. 1986,8-9. (189) Wakao, H.;Hiraguchi, H.; Ishii, T. A Simulation of Crystallization from Aqueous Supersaturated Solutions in a Batch Isothermal Stirred Tank. Chem. Eng. J. 1987,s (7),169-178. (190) Walton, A. G. Nucleation in Liquids and Solutions. In Nucleation; Zettlemoyer,A. C.,Ed.;MarcelDekker: New York, 1969. (191) Wang, M.-L.; Horng, H.-N.;Chang, R.-Y. Modeling of Crystallization Systems Via Generalized Orthogonal Polynomials Method. Chem. Eng. Commun. 1987,57,197-213. (192) White, E. T.; Wright, P. G. Magnitude of Size Dispereion Effects in Crystallization. Chem. Eng. Prog. Symp. Ser. 1971,67 (110),81-87. (193) Wing, G. M. ‘Integral Equations of the First Kind”; Technical Report LA-UR-84-1234;Los Alamos NationalLaboratory, 1984. (194) Winter, B. Modell zur diskontinuierlichen Massenkristallisation und dessen Anwendungbei der Masstabstibertragung. Chem. Tech. 1984,36(9),375-378. (195) Wiscombe,W. J. Improved Mie Scattering Algorithms. Appl. Opt. 1980,19 (9),1505-1509. (196) Witkowski, W. R. Model Identification and Parameter Estimation of Crystallization Processes. Ph.D. Dissertation. The University of Texas at Austin, 1990. (197) Witkowski, W. R.; Miller, S. M.; Rawlings, J. B. Light Scattering Measurements to Estimate Kinetic Parameters of Crystallization. In Crystallization as a Separations Process; Myerson, A. S., Toyokura, K., Eds.; American Chemical Society: Washington, DC, 1990;pp 102-114. (198)Witkowski, W. R.; Rawlings, J. B. Modelling and Control of Crystallizers. In Proceedings of the 1987 American Control Conference;American Control Conference: Minneapolis, MN, 1987; pp 1400-1405. (199) Zumstein, R. C.; Rousseau, R. W. Growth Rate Dispersion by Inital Growth Rate Distributions and Growth Rate Fluctuations. AZChE J. 1987,33 (l),121-129. (200) Zumstein, R. C.; Rousseau, R. W. Utilization of Industrial Data in the Development of a Model for Crystallizer Simulation. AZChE Symp. Ser. 1987,83 (253),130-139. Received for review November 20, 1992 Accepted April 12, 1993
