16th European Symposium on Computer Aided Process Engineering and 9th International Symposium on Process Systems Engineering W. Marquardt, C. Pantelides (Editors) © 2006 Published by Elsevier B.V.
Model-based optimization for operational policies in seeded cooling crystallization Ali Abbas^, S. Mostafa Nowee^, Jose A Romagnoli'^
'^ Nanyang Technological University, Singapore ^University of Sydney, Australia ^Louisiana State University, USA In this solution crystallization study, a population balance model that predicts the crystal size distribution (CSD) is used for targeting, by optimization, the product mean size and the coefficient of variation of the final CSD. The model is robust over a wide range of conditions and predicts the effects of heating during a batch since the kinetics of dissolution have been identified and incorporated into the model. The dynamic temperature profile and the initial seed size distribution are optimized while the initial seed induction time is conferred through knowledge of saturation conditions. This is quite a highly developed approach in contrast to existing anecdotal and/or ruleof-thumb experience currently dictating industrial operational methods. Results from optimizations under different objective functions are presented and validated experimentally. Keywords: Crystallization; Optimization; Seeding 1. Introduction Crystallization has been referred to as an art rather than a science. The main objective of crystallization systems engineering is to develop advanced understanding of crystallization phenomena enabling the efficient bulk processing of high quality crystalline products. Crystallization, in many industries, is the most common way of production of high value chemicals with high purity and desired size and shape. Crystal size is typically the most important factor to be controlled because of its close relation to the characteristics of the crystalline product, such as flowability and bulk density, and also different applications require different particle size. Crude methods are still being used in industrial operations to operate/control their crystallization processes and those industries that implement any control rarely do this at optimal conditions. Nevertheless, there have been a number of successfiil research investigations as evidenced in the literature on the control of the crystallization process with most of these studies aiming to control the product crystal size or size distribution. Crystallization processes are complex in nature and are known to be highly non-linear and thus pose difficulties in their modeling and control. They also carry liquid as well as solid (particulate) phases adding to the process's complexity. Their mathematical description is based on complex multi-dimensional population balance equations, coupled with integral DAEs. In addition, the kinetic history of the process largely influences the end product quality. Any unanticipated disturbances during process operation could lead to irreversible changes to product particle characteristics. It is thus evident and based on research observations that the use of the model is crucial for optimizing the crystallization operations allowing for end-product crystal size reproducibility. Researchers have recently recognised the significant role a crystallization model can play in optimizing operations. Manipulated variables that would be candidates for optimization decision variables are traditionally the working temperature and seeding. Considerable effort has been made on particle size control in batch cooling crystallization (Mullin, 2001). Most of these works consider supersaturation as the key variable to be optimized. This strategy which is based on the determination of the optimal temperature profile was first discussed by Mullin and Nyvlt (1971). It was though based on a key understanding namely that constant nucleation rate is required in order to keep the supersaturation low and within the metastable zone. Other works dealing with an
1347
1348
A.Abbasetal
optimal cooling profile generation intended to keep supersaturation low in the beginning stages of the process inhibiting nucleation but promoting crystal growth. Although seeding was believed to be an important factor for controlling crystal size (and allowing for the production of mono-dispersed CSDs), the amount of seed to be added was largely based on empirical knowledge. Only in recent years some empirical research appeared considering seeding as a factor to be optimized (Doki et al., 1999, Doki et a l , 2001, Doki et al., 2002, Jagadesh et al., 1996, Jagadesh et al., 1999, Kubota et al., 2001). In these investigations the seed chart was used to determine optimal seeding policies. These works demonstrated the impact seeding has on the crystallization showing even under un-controlled profile of cooling (eg. natural cooling), the rate of nucleation could be kept low. It could also be shown via simple mass balance calculations and the assumption that only seeds consume the generated supersaturation during the batch, the predicted final crystal mean size correlates well with the experimental ones. Other finds included the ability to produce an end-product with narrow monomodal CSD and further due to secondary nucleation, the final product could have a bimodal distribution. A model-based dynamic optimization of seeded cooling crystallization was discussed by Chung et al. (1999). Seeding parameters in addition to supersaturation were optimized and the objective function used was maximisation of final mean size with minimum coefficient of variance. They concluded that the effect of seeding variables on the final product characteristics are much more pronounced than that of supersaturation. Al-Zoubi and Malamataris (2003) and De Faveri et al. (2004) incorporated factorial experimental design using a surface response model in determining initial seed concentration and seed introduction time as well as the temperature profile. Their target was a specific crystal yield, size, shape and purity. Another investigation that dealt with model-based optimization of seeded cooling crystallizations was done by Choong and Smith (2004), who determined the optimal seeding characteristics, mean size and seed mass, by maximisation of product mean size and minimisation of CSD coefficient of variance as objective functions, respectively. Recent research trends, as evident in the literature, show an increase in the use of the dynamic crystaUization model for achieving various optimal objectives. The simultaneous use of the two variables, temperature and seeding, is evident in more recent works. The most comprehensive study appears to be that of Patience et al. (2004) who solve the population balance equation (PBE) for a growth-dependent dispersive crystal system and use temperature profile as well as two seeding variables as decision variables. They also validate their optimization resuhs experimentally. In this study, we extend on a previously developed dynamic crystallization model (Abbas, 2003) to include a mathematical description of (1) both super- and under-saturation crystallization regimes and (2) seeding. In (1) dissolution and disappearance kinetic rates are introduced to the model. We report the results of a model-based parameter estimation study undertaken to determine the kinetics parameters and thus identify the crystallization model. Finally an offline model-based optimization analysis was conducted to determine the optimal cooling profile in conjunction with the optimal initial seed size distribution. The latter being a unique determination that is novel to this work and extends on model-based optimization work by previous workers. We further present results of experimental validation of such dual decision variable offline optimization. 2. The crystallization model Since the crystallization process is a particulate one, the population balance equation (PBE) as proposed by Hulburt and Katz (1964) is used accounting for the evolution of crystal particles across temporal and size domains. For a batch crystallization system such as the ammonium sulfate one we are dealing with in this study, crystal growth is assumed to be non-dispersed and independent of crystal size and agglomeration and attrition are considered negligible, the population balance simplifies to
dt
dL
Where n(L,t) is the number density of crystals, t is time, L is the characteristic crystal size, G is the growth rate of the crystals while B is the nucleation rate. In the case where crystal dissolution by heating exists such as when the process is operating under temperature control giving rise to
Model-Based
Optimization for Operational Policies
1349
temperature increases, the PBE will change and terms denoting dissolution (Z)) and disappearance (Da) of particles will need to be incorporated (Equation 2). The structure of PBE remains the same as in Equation 1, but the dissolution and disappearance terms are introduced to replace the growth and nucleation ones respectively and with opposite signs.
M^=^M^_z)
(2)
dt dL "" The switch from supersaturation state to undersaturation state results in elimination of growth and nucleation processes, and in the occurrence of dissolution and disappearance processes. The reverse occurs when the process switches from under- to supersaturation (eg. switching from heating to cooling). In other words, where there are high supersaturation levels, one would expect to observe nucleation and growth of crystals and Equation 1 would hold true. Where the process is going in the direction of under-saturation, i.e. heating, one would observe dissolution and disappearance of crystals rendering Equation 2 to be true. Clearly, this is only true for those solute systems that have a solubility that increases with temperature. In crystallization processes some phenomena that may occur are nucleation, growth, agglomeration and breakage. These phenomena govern the process and ultimately dictate the product characteristics and particularly the shape of the product PSD. In this study, two additional kinetic rates are considered namely dissolution and disappearance. As the objective of this research is to identify optimal operational policies, the kinetic model structure is kept simple and confined to a power law structure.
R = k^C^
(3)
Where R is the vector of kinetic rates [B G D Da\ representing the rates of nucleation, growth, dissolution and disappearance, k is the vector of kinetic rate parameters [kb kg kd ka], AC is the supersaturation which is equal to the difference between the concentration of the solute C and the saturated concentration C . Finally 0 represent the vector of rate orders [b g d a] corresponding to each of the kinetic rates respectively. To complete the model, the mass and energy balances, after Rawlings et al. (1993), for the solute and crystallizer temperature, are incorporated. The model was solved by finite difference discretization method. For the supersaturated and undersaturated states, the backward and forward finite techniques were used respectively, provided stability. The initial condition for the system is the size distribution of seeds. The solution was implemented and simulated using the process modeling package gPROMS (Process Systems Enterprise Ltd, UK) which executed the solution within a timeframe of a few seconds using a computer running a Pentium 4 processor.
3. Model identification and optimization Using optimization based (maximum likelihood criteria) parameter estimation, it was possible to identify the parameter sets k and 0, that provide the highest probability of the model predicting the real data. The gEST function of gPROMS was used for parameter estimation here. The problem was split in two sections, cooling and heating. The parameters in nucleation and growth kinetics were estimated using the experimental data from the cooling phase while the dissolution and disappearance kinetics were estimated from the heating phase. The estimated parameters are log(A:) = [3.25 -6.55 -10.04 -6.83] and 0= [10.8 6.0 2.0 1.0]. Experiments for the purpose of parameter estimation were conducted in a glass jacketed crystallizer unit with a volume of 1 litre. Cooling/heating was implemented by means of a circulator (Julabo, USA) connected to the Honeywell industrial distributed control systems (DCS). The temperature profile was preprogrammed and realized using a standard PID controller using set-point tracking. The feed to the crystallizer was prepared from nano-pure water and Laboratory reagent grade ammonium sulfate.
1350
A.Abbasetal
The ammonium sulfate solution was prepared so as to have a concentration corresponding to the saturation temperature of 50°C and a volume of 0.85 litres. We use a conductivity (Orion conductivity probe model 018020A, USA) and temperature to track the solution concentration while CSD was determined offline using laser diffraction (Malvern Mastersizer-S, UK). The experiments designed for parameter estimation consisted of two phases of linear temperature change; one cooling and one heating. During the cooling phase supersaturation develops affecting the crystallization. This is then followed by a heating phase affecting the dissolution of the crystals formed during the cooling phase. Solution concentration and CSD data under these cooling/heating regimes are recorded providing information on the kinetic rates and allows the quantification of the parameter sets. For experiments leading to offline optimization validation, the temperature profiles and seed determinations from optimization results were implemented. The subsequent step of optimization is discussed next. The optimization of the model was performed using the gOPT function in gPROMS. Two objective functions were posed; one is to maximize the mean crystal size, L , in the shortest time while keeping the distribution variance {Var.) at a minimal value (Figure la) and the other is to minimize the distribution variance while constraining L to a specific end-point value again in minimum time (Figure lb). The control variables are the temperature of the crystallization and seed loading profile. The formulation of these objective optimization problems is summarized in Figures la and lb. Ten control intervals were implemented in this dynamic optimization. Temperature being the time-variant control variable was specified as piecewise-linear while the seed distribution was specified as a piecewise-constant profile for each class of crystal size. Finding the optimal seed distribution was possible by applying anecdotal knowledge, bounding the seed distribution to a size range known to provide adequate seeding conditions. For instance the seed size must be above a critical size determined from experimental investigations otherwise the seeds will disappear before realising their effect. A number of constraints were implemented and consisted of the path rate of cooling constraint and end-point constraints reflecting the end-product characteristics.
4. Results and Discussion The solution of this non-linear optimization problem represented the optimal combination of //, T(t) and n(L,t) initial that gives the maximum or minimum value of the objective function while satisfying the given constraints. The significance of the dynamic optimization problem is that a path to optimality is drawn. In this case the optimal path is described by a temperature trajectory and a corresponding mean size trajectory (Figures 2a & 2b). The optimal temperature profiles are very similar and comprise of three distinct regions. In region one, we observe an initial decrease in temperature that can be explained by the drive to cross over the saturation temperature to generate supersaturation. Region 2 sees a phase of slow cooling in which the supersaturation is being stabilized and growth of crystals is occurring under controlled nucleation towards the maximum size. In the third and final region, a sudden decrease in temperature is observed and is the result of the optimization attempting to minimize tf and so forcing the temperature down sharply. The corresponding mean size profiles simulated under the optimal temperature profiles are shown in Figure 2b. Again similarities between the two exist thought the results of the maximization of the mean size are slightly better. No size reduction of the seeds is because the induction time of the seeds is fixed at a point very close to saturation conditions rendering the rate of dissolution to be slow so that the reduction in seed size is not appreciable. The crystals continue to grow and stabilise at about the same size for both optimizations leading to suggest that the two optimization problems provide similar end-size solutions and thus are similar in the way they are posed. As for the time of the batch, no significant difference is observed between the two optimizations. The optimal seed distributions for the two optimization problems are shown in Figure 3 a. The two optimizations observed very similar seeding optima. The end product size distribution was found to vary for the two optimization problems which are shown in Figure 3b. The height and shape of the distributions are distinctly different suggesting the possibility of shaping the CSD to a specific shape desired. These outputs from the optimizations where implemented to the bench scale crystallizer and the results show a good agreement. The temperature profiles consist of four linear profiles reasonably close to the profile given by optimizations (Figures 4a & 4b). Mean size measurements were close to what it was expected
Model-Based Optimization for Operational Policies
1351
from simulation studies, although there is a very minor fluctuation from the mean size growing profiles presented from simulations (Figures 5a & 5b), which could be mainly due to the oscillatory behavior in the controlled measured temperature inside the crystallizer . However, as could be seen in Figure 6 the final size distribution has a very narrow shape compared with a process without seeding and with linear cooling profile, which proves the claim of controlling the size distribution by the addition of the seed. Objective function 1:
mcLX,L\tf) t,T.(^\o.t,\
Objective function 2:
min
Subject to:
Subject to:
Crystalliser model: Equations 1-3
CrystalHser model: Equations 1-3
Initial conditions: n(L,t) = 0
Initial conditions: n(L,t) = 0
C(t) = Co Batch time:
t=0
C(t) = Co
t=0
t"^'" < tj. < t'J'"
Control variables: T'"'" < T{t) < T"'"'
Batch time:
J""'"
End-point constraints:
T """ < T\ty j
cr
Figure la. Optimization formulation for maximum mean size.
t=0
t=0
t'J'" ^tf ^ tf"''
Control variables: T'"'" < T{t) < T'"'^
0
Var.ytj^j
^^^
Path constraint:
0
End-point constraints:
T'"'" < T\tj ^
cr
Figure 2(a) Optimal temperature profiles for the two optimization objectives (left) and 2(b) corresponding optimal mean size profiles (right).
Figure 3(a) Optimal seed distributions for the two optimization objectives (left) and 3(b) corresponding optimal end-product size distributions (right).
1352
A. Abbas et al.
Figure 4(a) Optimal temperature profile for objectiveftinction:maximum mean size (left) and 4(b) corresponding optimal mean size profile versus experimental (right).
Figure 5(a) Optimal temperature profile for objective function: minimum variance (left) and 5(b) corresponding optimal mean size profile versus experimental (right).
Size [micron]
Figure 6. Comparison between fmal product distribution of optimal and linear cooling profiles. References: A. Abbas, (2003), A model-based framework for advanced operation of crystallization process - A pilot plant study, PhD dissertation, The Department of Chemical Engineering, University of Sydney. N. Al-Zoubi and S. Malamataris, (2003) International Journal of Pharmaceutics, 260, 123-135. K.L. Choong and R. Smith, (2004) Chemical Engineering Science, 59, 313-327. S.H. Chung, D.L. Ma and R.D. Braatz, (1999) Canadian Journal of Chemical Engineering, 77, 590-596. D. De Faveri, P. Torre, P. Perego and A. Converti, (2004) Journal of Food Engineering, 61,407-412. N. Doki, N. Kubota, A. Sato and M. Yokota, (2001) Chemical Engineering Journal (Lausanne), 81, 313-316. N. Doki, N. Kubota, A. Sato, M. Yokota, O. Hamada and F. Masumi, (1999) AIChE Journal, 45, 2527-2533. N. Doki, N. Kubota, M. Yokota and A. Chianese, (2002) Journal of Chemical Engineering of Japan, 35, 670-676. H.M. Hulburt and S. Katz, (1964) Chemical Engineering Science, 19, 555-74. D. Jagadesh, N. Kubota, M. Yokota, N. Doki and A. Sato, (1999) Journal of Chemical Engineering of Japan, 32, 514-520. D. Jagadesh, N. Kubota, M. Yokota, A. Sato and N.S. Tavare, (1996) Journal of Chemical Engineering of Japan, 29, 865873. N. Kubota, N. Doki, M. Yokota and A. Sato, (2001) Powder Technology, 121, 31-38. J.W. Mullin, (2001) Crystallization, Butterworth-Heinemann, Oxford. J.W. Mullin and J. Nyvlt, (1971) Chemical Engineering Science, 26, 369-377. D.B. Patience, P.C. Dell'Orco and J.B. Rawlings, (2004) Organic Process Research & Development, 8, 609-615. J.B. Rawlings, S.M. Miller and W.R. Witkowski, (1993) Industrial & Engineering Chemistry Research, 32, 1275-96.