Combined Quadrature Method of Moments and Method of

Ind. Eng. Chem. Res. 2009, 48, 8575–8584

8575

Combined Quadrature Method of Moments and Method of Characteristics Approach for Efficient Solution of Population Balance Models for Dynamic Modeling and Crystal Size Distribution Control of Crystallization Processes E. Aamir,† Z. K. Nagy,*,† C. D. Rielly,† T. Kleinert,‡ and B. Judat‡ Loughborough UniVersity, Loughborough, LE11 3TU, United Kingdom, and BASF SE, Ludwigshafen, D-67056, Germany

The paper presents a novel methodology for the estimation of the shape of the crystal size distribution (CSD) during a crystallization process. The approach, based on a combination of the quadrature method of moments (QMOM) and the method of characteristics (MOCH), provides a computationally efficient solution of the population balance equation (PBE) and hence a fast prediction of the dynamic evolution of the CSD for an entire batch. Furthermore, under the assumption that for supersaturation-controlled crystallization the main phenomenon is growth, an analytical CSD estimator is derived for generic size-dependent growth kinetics. These approaches are evaluated for the crystallization of potassium alum in water. The model parameters are identified on the basis of industrial experimental data, obtained using an efficient implementation of supersaturation control. The proposed methods are able to predict and reconstruct the dynamic evolution of the CSD during the batch. The QMOM-MOCH solution approach is evaluated in a model-based dynamic optimization study, which aims to obtain the optimal temperature profiles required to achieve desired target CSDs. The technique can serve as a soft sensor for predicting the CSD, or as a computationally efficient algorithm for off-line design or online adaptation of operating policies based on knowledge of the full CSD data. 1. Introduction Crystallization from solution is an industrially important unit operation due to its ability to provide high purity separation. It is widely used in industries such as petrochemicals, pharmaceuticals, microelectronics, and food processing. Batch cooling crystallization provides the advantages of being simple, flexible, and generally requiring less process development and investment than many other separation/purification techniques. However, despite the long history and widespread application of batch crystallization, there are a disproportionate number of problems associated with its control. The crystal size distribution (CSD) is a key factor in the production of a high quality product and determines the efficiency of the downstream processes.1,2 Therefore many problems in downstream processes can be attributed to poor particle characteristics established in the crystallization step. The control objectives for batch crystallization processes can be defined in terms of product purity, crystal habit, morphology, average particle size, crystal size distribution, bulk density, product filterability, and dry solid flow properties, which all depend on the CSD. The main difficulty in batch crystallization is to accomplish a uniform and reproducible CSD.3 Online control during the process allows for improved crystalline product quality, shorter process times, and reduction or elimination of compromised batches. Even if some of these objectives can be expressed in terms of the moments of the distribution, knowing and predicting the entire shape of the distribution allows the design and adaptation of operating policies to achieve improved product quality. A potential way to enhance the control of CSD is to use supersaturation control,4-10 or direct nucleation control,11,12 which drives the process within the metastable zone to avoid nucleation or to generate controlled nucleation/dissolution * To whom correspondence should be addressed. E-mail: z.k.nagy@ lboro.ac.uk. Tel.: +44 1509 222516. Fax: +44 1509 223923. † Loughborough University. ‡ BASF SE.

events. Although these approaches can provide improved consistency of the CSD, the actual prediction and estimation of the shape of the distribution at the end of the batch can provide useful information for monitoring or designing the operating curve for the supersaturation controller. Model-based approaches can be used for better control13-21 but also for product design by reverse engineering the process to achieve the desired CSD and shape.22-25 Population balance models have been widely used for modeling crystallization processes.26,27 The solution of the generic population equation (PBE) usually requires computationally expensive, complex numerical solution techniques.28,29 The variety of solution approaches proposed for the PBE include the standard method of moments (SMM), the quadrature method of moments (QMOM),30-32 the method of characteristics (MOCH),22 the method of classes (MOC),33,34 direct numerical solution (DNS) techniques, such as finite volume and finite difference schemes,35,36 and dynamic Monte Carlo (DMC) simulation approaches.31-37 The first two methods only provide average characteristics of the CSD expressed in terms of low order moments (e.g., the mean and standard deviation), whereas the last two methods are usually computationally too expensive for online real-time applications. Both SMM and QMOM provide efficient solutions of the PBE and have been widely used in the literature for optimization and control purposes.4,14 However, these approaches only provide the moments of the CSD and not the entire distribution. Several techniques are available to reconstruct the distribution from its moments, for example, using linear or nonlinear inversion approaches. These techniques have the disadvantage that they require a larger number of moments and generally suffer from solution multiplicity and ill-conditioning problems. Approximate distribution functions (e.g., polynomial, normal, gamma, or log-normal), or a weighted sum of distributions, using for example orthogonal polynomials as weighting functions, can also be used to approximate the shape of the distribution based on the

10.1021/ie900430t CCC: $40.75  2009 American Chemical Society Published on Web 07/31/2009

8576

Ind. Eng. Chem. Res., Vol. 48, No. 18, 2009

moments.27,38 However, the solution of such inverse problems is usually not unique. There is also a lack of systematic methodologies for the choice of the suitable type and number of base functions and distributions. Hence the approximate distribution functions resulting from both categories of reconstruction methods are subject to spurious oscillations and the correctness of the resultant shape of distribution is difficult to evaluate in most practical cases. The method of characteristics (MOCH) in combination with the SMM has been used successfully for processes with size-independent growth and nucleation,22 including the modeling and control of polymorphic transformations.10,39 However this approach does not apply in the more generic case of a PBE with size-dependent growth and dissolution, or when breakage and agglomeration mechanisms need to be considered. The approach presented in the paper combines the advantages of QMOM and the MOCH to provide a computationally efficient technique for the prediction of the entire CSD. The algorithm can be applied for the solution of population balance equations with generic size-dependent growth and nucleation kinetics and can provide a more general framework for the efficient solution of PBEs even in the case of breakage and agglomeration. According to the authors’ knowledge, this is the first time that these two approaches have been combined to provide an efficient solution approach for a model-based control of distribution shape, with the original idea first presented by the authors in Aamir et al.45 The method has also been used to identify the kinetic parameters for the batch cooling crystallization of potash alum from water using industrial experimental data. An analytical CSD estimator is also presented for the first time, which can be used in the case of supersaturation controlled, growthdominated processes. It is shown that the proposed approach provides a computationally efficient CSD estimation technique, which can be used for off-line parameter estimation, crystallization design, or for online estimation and control. The approach is applied in a model-based optimal control approach to determine optimal temperature trajectories to achieve desired target CSDs.

of moments (in more generic cases including size-dependent growth, breakage, and aggregation. Both methods calculate the moments of the distribution defined by, µj )

∫

∞

0

fn(L) Lj dL,

j ) 1, 2, ..., ∞

(2)

The QMOM is a generic solution approach for the PBE.30,32,40 It employs a quadrature approximation of the distribution function Nq

fn(L, t) ≈

∑ w (t) δ(L (t), L) i

i

(3)

i)1

where Nq is the number of quadrature points and the corresponding weights, wi, and abscissas, Li, can be determined through the product-difference (PD) algorithm41 or via direct solution of a differential-algebraic (DAE) system,42 based on the idea of minimizing the error committed by replacing the integral from the moment definition with its quadrature approximation, µj )

∫

∞

0

Nq

fn(L) Lj dL ≈

∑wL

j i i

(4)

i)1

Applying the moment transformation to eq 1 with the quadrature approximation of eq 4 the resulting moment equations have the form32,40 dµ0 ) B(S;θb) dt N

q dµj )j wiLj-1 i G(S, Li ;θg) + dt i)1

∑

B(S;θb)rj0,

j ) 1, 2, 3,...

(5)

Considering a single growth direction with one characteristic length L and a well-mixed crystallizer with growth and nucleation as the only dominating phenomena, the population balance equation (PBE) has the form

The generic PBE eq 1 can be reduced to a system of ODEs by applying the MOCH. The aim of the MOCH is to solve the PBE by finding characteristic curves in the L - t plane that reduce the partial differential equation to a system of ODEs. The L - t plane is expressed in a parametric form by L ) L(Z i) and t ) t(Z i ), where the parameter Z i gives a measure of the distance along the characteristic curve. Therefore, fn(L,t) ) i ),t(Z i )), and applying the chain rule gives fn(L(Z

∂fn(L, t) ∂(G(S, L;θg)fn(L, t)) + ) B(S;θb)δ(r0, L) ∂t ∂L

dfn dL ∂fn dt ∂fn ) + dZ i dZ i ∂L dZ i ∂t

2. Combined QMOM-MOCH Approach for the Efficient Solution of PBEs

(1)

where fn(L,t) is the crystal size distribution expressed as a number density function (number of crystal per unit mass of slurry), t is time, G(S,L;θg) is the rate of crystal growth, B(S;θb) is the nucleation rate, S ) (C - Csat) is the absolute supersaturation, C is the solute concentration, Csat ) Csat(T) is the saturation concentration with T being the temperature, r0 is the size of nuclei, δ(r0,L) is the Kronecker delta (δ ) 1 if L ) r0 and δ ) 0 if L * r0) and θg and θb are vectors of growth and nucleation kinetic parameters, respectively. The solution of eq 1 is an initial value problem, with initial condition given by the size distribution of seed, fn(L,0) ) fn,0(L0). Equation 1 can be transformed into a system of ODEs by applying the standard method of moments (in the case of sizeindependent growth and nucleation) or the quadrature method

(6)

In the case of generic growth kinetics, eq 1 can be rewritten in the form of ∂fn(L, t) ∂fn(L, t) dG(S, L;θg) + G(S, L;θg) ) -fn(L, t) + ∂t ∂L dL B(S;θb) δ(r0, L)

(7)

Comparing eqs 6 and 7 it can be shown that Z i ) t and the characteristic equations are given by the following system of ODEs: dL ) G(S, L;θg) dt

(8)


dG(S, L;θg) dfn(L, t) ) -fn(L, t) + B(S;θb) δ(r0, L) dt dL

(9)

with initial conditions L ) L0 and fn(L,0) ) fn,0(L0), that is, the seed CSD. To obtain the dynamic evolution of the crystal size distribution fn(L,t), eqs 8 and 9 with prescribed nucleation and growth expressions can be integrated repeatedly for different initial values [L0,fn,0(L0)]. The initialization of the integrations in the L-t plane are illustrated in Figure 1, showing typical evolutions of the characteristic lines during the integration. To simulate the growth of the seed, the initial conditions start from along the L axis of the L-t plane, with values calculated by choosing a discretization interval ∆L0 and using t0 ) 0 and L0 ) max (0,L0,max - k∆L0), k ) 0,1, ..., N, where N is the number of discretization points for the seed distribution and L0.max is chosen to be greater or equal to the maximum size range of the seed crystals. The discretization interval ∆L0 will determine the number of integrations (the number of characteristic lines) and hence the resolution of the dynamic evolution of the seed CSD. For this part of the integration the initial values for the probability distribution function are calculated from the seed distribution fn,0(L0) ) fseed(L0), and all integrations start from an initial time t0 ) 0. Figure 1 also shows the method employed to represent the contribution to the overall distribution function from nucleation events, which may occur during the batch. In this case, the characteristic lines for nucleation and growth of new particles, start from initial conditions along the t axis of the L-t plane, using L0 ) 0, fn,0(L0) ) 0, and t0 ) tnext. The initial time for the next integration, tnext, is calculated by interpolating the characteristic line for L ) r0, as shown in Figure 1. The number of integrations within this part of the algorithm is not predetermined and will depend on the evolution of the characteristic lines governed by the growth kinetics. This is an adaptive feature of the algorithm, which allows the high resolution prediction of the part of the CSD that results from nucleation events. The iterations are stopped when tnext g tf where tf is the end time of the batch. For the solution of eqs 8-9, it is considered that at the moment of nucleation, nuclei can have any size between 0 and r0. This is described by the modified delta function defined as δ(r0, L) )

{

1 if L ∈ [0, r0] 0 if L ∉ [0, r0]

8577

(10)

Thus nucleation events are assumed to occur for L e r0. For seeding crystallization, secondary nucleation is considered as

Figure 2. Flowchart of the combined QMOM-MOCH approach for the solution of PBEs.

the dominating nucleation phenomenon, which is generally expressed as a function of the supersaturation and the volume of the existing crystals, given by the third-order moment of the size distribution. Hence in the model, B ) kbSbµ3. The formulation given by eq 10 allows the direct consideration of apparent nucleation kinetics in the model, where r0 is the size of the particles when they are first detected with a particular measurement approach. In this study r0 ) 1 µm is used since it represents approximately the size of particles which can be detected by typical in situ process analytical tools based on image analysis or focused beam reflectance measurements. Both the growth and nucleation rates are functions of the supersaturation, S, which can be calculated from the material balance. The solute concentration is given by C(t) ) C(0) - kvFc(µ3(t) - µ3(0))

Figure 1. Evolution of characteristic lines with the generic approach of calculating the initial conditions for the method of characteristics in the case of growth and nucleation mechanisms.

(11)

where Fc is the density of crystals and kv the volumetric shape factor. The solution of eqs 8, 9, and 11 requires a priori knowledge of the dynamic evolution of the supersaturation, S(t) and/or the third moment µ3(t), which can be obtained by using the moment transformation of eq 1 via the SMM or QMOM. The main steps of the proposed algorithm are shown in Figure 2. In the case of secondary nucleation and size-dependent growth, the ODEs from the QMOM have to be integrated together with eq 11 once for the duration of the batch to predict the evolution of µ3 and the variation of supersaturation with time, then S(t) and µ3(t) are used in the nucleation and growth kinetic expressions during repeated integrations of eqs 8,9 with different initial conditions, to map out the complete evolution of the full CSD via the MOCH. When nucleation is included

8578


(

fn(L) ) fn,0(L0) 1 +

Figure 3. Operating curves represented in the phase diagram for seeded and unseeded batch cooling crystallizers.

in the model, an iterative integration of eqs 8-9 is needed, since the number of initial conditions along the t axis, for nucleating and growing particles, is not known a priori. Note that the method also applies for dissolution problems and can be extended for certain breakage and agglomeration mechanisms, which would appear as extra terms on the right-hand side of eq 9 and would be computed using the evolution of the weights and abscissas from the solution of the QMOM part of the algorithm. 3. Analytical CSD Estimator for Growth-Dominated Crystallization Systems The traditional way of controlling cooling crystallization processes is to follow a predetermined temperature profile in time. Recent developments in the direct design of crystallization systems have led to a widespread application of supersaturation control.4,5,7 The direct design approach is based on the idea of operating the system within the metastable zone bounded by the nucleation and solubility curves (see Figure 3). In this technique a supersaturation set point profile is chosen experimentally and it is followed in the phase diagram using a supersaturation controller based on concentration measurement. The supersaturation is usually chosen to be constant and with the application of properly designed control algorithms, in the case of seeded crystallization the process is maintained at the desired constant supersaturation throughout the entire batch. For constant supersaturation controlled seeded crystallization, seed is added to suppress nucleation, and the operation will be dominated by growth. Considering the generic case of sizedependent growth given by G ) kgSg(1 + γL)p

kgSgγ(1 - p)t (1 + γL0)1-p

)

p/(p-1)

(16)

Discretizing the initial seed size distribution fn,0(L0) for different values of L0, eqs 15, 16 can be used to compute the dynamic evolution of the CSD for a generic growth-dominated process (the solution is valid for p * 1 and γ * 0; for p ) 1 and/or γ ) 0 the analytical expressions can also be derived). Equations 15, 16 are of significant practical importance, because they show that if the initial seed distribution and growth kinetics of a supersaturation controlled crystallization process are known, the evolution of the CSD can be estimated/predicted by using concentration (supersaturation) measurements. 4. Model Identification and Validation The QMOM-MOCH approach has been validated for the batch cooling crystallization of the inorganic compound, potash alum (KAl(SO4)2) in water. The experimental data were obtained using an industrial pilot-plant crystallization system located at BASF (Ludwigshafen, Germany). A supersaturation controller was implemented in the process and seeded experiments were carried out at constant supersaturation set-points. The CSD was measured at different time intervals using laser-diffraction, whereas concentration measurements were obtained based on density measurements of the slurry. The results from two experimental runs (experiments A and B) are shown in Figure 4. Experiment A was conducted at a supersaturation set-point Ssp ) 0.6% (weight percent in kg solute/kg slurry), whereas experiment B was at Ssp ) 0.3% (weight percent). Seed was introduced in both cases shortly after the supersaturated state had been reached. In the case of experiment A it can be seen that the supersaturation controller exhibits an overshoot during the initial part of the operating curve, which leads to secondary nucleation. Since experiment A captures both the growth and the nucleation mechanisms, it was used for model parameter

(12)

where θg )[kg,g,γ,p] is the growth parameter vector, for growthdominated systems (B ) 0) the model eqs 8-9, reduce to the form: dL ) kgSg(1 + γL)p dt

(13)

dfn(L, t) ) -kgSgγp(1 + γL)p-1fn(L, t) dt

(14)

In the case of well-controlled constant supersaturation, which follows the desired set-point value, Ssp, eqs 13, 14, can be solved analytically, with the solution: L ) (((1 + γL0)1-p + kgSgγ(1 - p)t)1/

(1-p)

- 1)/γ (15)

Figure 4. Experimental results in the case of supersaturation controlled experiments. (a) Experiment A: Ssp ) 0.6 wt %, used for parameter identification. (b) Experiment B: Ssp ) 0.3 wt %, used for validation.


8579

Table 1. Size-Dependent Growth and Nucleation Parameters for the Crystallisation of Potash Alum in Water (Units for S ) kg/kg slurry) parameter

value

growth rate constant (kg), µm s-1 growth constant (γ), µm-1 growth constant (p), -growth order constant (g), -nucleation rate constant (kb), µm-3 s-1 nucleation order constant (b), --

8.5708 0.0050 1.5777 1.0000 0.0380 3.4174

identification using the QMOM-MOCH approach described in section 2, whereas experiment B was used for validation of the analytical estimator. For the potash alum system, size-dependent growth has been reported in the literature,44 and was observed experimentally. Hence a generic size-dependent growth expression given by eq 12 was used in the model identification. The nucleation and growth parameters were determined to capture the dynamic evolution of the shape of the size distribution, as well as the experimental concentration profile. The optimization problem for the parameter estimation is given by K

min{wf θ

Nd

∑ ∑(f

Figure 5. Dynamic evolution of the modeled and experimental CSD for experiment A.

K

exp 2 V,k(Ll)-fV,k (Ll)) + wC

k)1 l)1

∑ (C

k

2 - Cexp k ) }

k)1

(17) subject to (18)

θmin e θ e θmax

with θ ) [kg,g,γ,p,b,kb] being the model parameter vector with the growth and nucleation kinetic parameters, θmin and θmax are vectors with specified minimum and maximum bounds for each parameter, respectively, Ck and Ckexp are the simulated and experimental concentration values at the discrete time steps k exp are the values of the simulated and ) 1, ..., K, fv,k and fV,k experimental volume probability distribution functions, corresponding to the discretized size Ll, l ) 1, ..., Nd, with Nd being the number of experimental size bins, and wf, wC are objective function weighting factors. The volume distribution function is calculated as Nd

∑ (f

(19)

Figure 6. Experimental and simulated results: (a) concentration and (b) de Broucker mean diameter during the entire batch of experiment A.

The optimization problem is solved using a sequential quadratic programming (SQP)-based optimization approach implemented in the Matlab function fmincon. Note that finding the best kinetic parameters is generally a difficult optimization problem due to the strong correlation between the parameters and the generally nonconvex optimization problem 17, 18. Supersaturation controlled experiments can be used to design experiments which allow to decouple the identification of the kinetic parameters or quick metastable zone determination experiments can be used for providing experimental data based initial guesses for the parameter identification process.43 The resultant model parameters for the potash alum system are presented in Table 1.45 The dynamic evolution of the modeled and experimental CSDs are in very good agreement during the entire batch, as shown in Figure 5. It can be seen that due to the particular size-dependent growth kinetics of this system, the CSD broadens with decreasing height during the batch. The formation of a secondary CSD peak can also be observed, which is the result of secondary nucleation. The QMOM-MOCH

approach with the model using the identified growth and nucleation parameters is able to describe the main features of the CSD throughout the entire batch. Figure 6 shows the comparison between the experimental and modeled concentration and weight mean size throughout the batch, which are also in relatively good agreement. Figure 7 illustrates the evolution of the characteristic lines and the discretized number distribution function predicted by the simulation of experiment A, using the combined QMOMMOCH. The evolution of the characteristic lines show the broadening of the distribution function due to the size-dependent growth kinetics, as can be observed from Figure 5. The distribution function is initialized at t ) 0 with values obtained from the seed distribution, after which the values of fn decrease as the distribution broadens. At different time steps new nuclei and new characteristic lines appear using the methods described in section 2. The discretization intervals along the time axis depend on the growth kinetics according to the approach described in section 2 and illustrated in Figures 1 and 2. It can be observed that during 10-30 min into the batch, when the

fV,l ) fn,lLl / 3

3

n,lLl

∆Ll)

l)1

8580


Figure 9. Performance of the analytical estimator initialized with CSD at t ) 30 min (experiment B).

Figure 7. Evolution of characteristic lines (a) and number distribution function (b) for the simulated experiment A.

growth is faster (and the nucleation rate is also more significant), the discretization is finer, compared to the later stages of the batch characterized by slower growth. The combined QMOM-MOCH can be used not only for the model identification, but also for CSD prediction. The simulation time for the reconstruction of the entire evolution of the CSD during the batch, only takes a few seconds on a standard PC running Matlab. For even faster computational performance, the analytical estimator given by eqs 15 and 16 can be used in the case of growth-dominated processes. Figure 8 shows the CSD for the seed, selected times during the batch, and the final CSD for experiment B, compared to predictions from the analytical estimator with initial concentration, temperature profile, and seed distribution set to match the experiment, along with kinetic parameters deduced from experiment A (see Table 1). In this case the supersaturation was well maintained at its constant setpoint. The simulated and experimental CSDs are in good agreement overall. However, the experimental CSD shows evidence of secondary nucleation during the batch, indicated in Figure 8 by the secondary CSD peak developed during the

crystallization. The analytical estimator is derived based on the assumption of constant supersaturation and no nucleation. Therefore, initializing it with the seed CSD and applying it in open-loop, the analytical estimator is not able to predict the development of the secondary peak. However, in the case of many practical applications, an online measurement of the CSD is available typically with a sampling time in the range of 1-15 min (e.g., by using focused beam reflectance measurement coupled with inverse geometric modeling to transform chord length distributions into size distributions).46,47 In these cases the analytical estimator can be used in closed-loop, initializing it with the new CSD measurement every time it becomes available. Figure 9 illustrates the results when the estimator was initialized using the measured CSD after 30 min. The effect of the secondary nucleation, which occurred in the first 30 min of the batch, on the final CSD is partially predicted. The proposed method can be used as an efficient estimator for monitoring and predicting the CSD at the end of the batch, or in off-line or on-line optimization approaches for designing crystallization systems to produce consistently the desired final CSD. These results show that in practical applications the analytical estimator can be used for prediction of the CSD even if during the initial phase of the batch the supersaturation is not constant. As new CSD measurements become available and the supersaturation

Figure 8. Comparison between measured and simulated CSD using the analytical CSD estimator (experiment B).


8581

Figure 10. Simulated dynamic evolution of CSD with optimal temperature profile throughout the batch.

reaches its constant set point value, the predicted CSD will converge to the correct value. 5. Model-Based Dynamic Optimization of Temperature Trajectories for Shaping the CSD The combined QMOM-MOCH approach described in section 2 was used to solve the PBM in a model-based dynamic optimization scheme, with the aim to determine optimal temperature trajectories which yield desired target CSDs at the end of the batch. The final CSD is dependent on the supersaturation profile created over the batch time, hence the cooling trajectory is of critical importance for the final CSD. In the optimization both the temperature trajectory and the batch time were optimized. For the solution of the dynamic optimization problem the batch time horizon [0, tf] is divided into Nb equally spaced time intervals of ∆t (stages), with discrete time steps tk ) j∆t, j ) 0, 1, ..., Nb, and the temperature trajectory is approximated by a piece-wise linear function determined by the fixed initial temperature at t ) 0, T(0), and the slopes RT(j) in each discretized period ∆t. Note, that since the batch time is also optimized the duration of the time interval ∆t changes during the optimization, but the number of disretizations Nb is fixed. This formulation allows an easy incorporation of the temperature rate constraints (as bounds on the decision variables RT(j)), which are very important to obtain a practically implementable temperature trajectory. The optimization problem is formulated as follows: K

min

RT(j),tf

Nd

∑ ∑(f

target 2 v,k(Ll)-fV,k (Ll))

k)1 l)1

(20)

subject to RT,min e RT( j) e RT,max,

j ) 0, 1, ..., Nb

(21)

0 e tf e tf,max

(22)

C(tf) e Cf,max

(23)

with RT(j) the elements of the vector containing the slopes, (dT/ dt) for the temperature trajectories and tf is the total batch time, C(tf) is the solute concentration at the end of the batch and Cf,max is the maximum acceptable concentration at the end of the batch target are the values of the to achieve the required yield, fV,k and fV,k simulated and the target volume probability distribution functions at the discrete time steps k ) 1, ..., K, where measurement data were available, corresponding to the discretized sizes Ll, l ) 1, ..., Nd, with Nd being the number of experimental size bins. The optimization problem is solved using a sequential quadratic programming (SQP) approach implemented in the Matlab function fmincon. The kinetic parameters used are described in Table 1. Figures 10 and 11 show the results of the distribution shaping optimal control. The seed and target distributions are shown on the first plot (at t ) 0) in Figure 10. The target CSD is a fictitious bimodal distribution. For simulation purposes a monomodal seed was selected, represented by a Gaussian distribution with mean of 54 µm and standard deviation of 15 µm, chosen based on CSD measurement of the seed used in the experiments. Figure 10 shows the dynamic evolution of the CSD toward the target bimodal distribution, when the optimal temperature trajectory shown in Figure 11a is implemented in the simulation. The final CSD is in good agreement with the target distribution.

8582


Figure 11. Optimal control results of the simulations for bimodal target distribution with pronounced secondary peak: (a) optimized temperature profile with 50 discretization points, (b) phase diagram showing solubility and optimal operating curve, (c) concentration profile during the batch (d) supersaturation profile (kg/kg slurry) during the batch, (e) nucleation rate profile during the batch, (f) growth rate profile during the batch.

Figure 12. Optimal control results for bimodal target distribution with minimized fines due to secondary nucleation: (a) comparison of model and target distribution, (b) optimized temperature profile with 30 discretization intervals, (c) concentration profile during the batch, (d) supersaturation profile during the batch.

The first 10-15 min of the batch are devoted mainly to the growth of the seed crystals, which approach rapidly the larger size mode of the target distribution, due to the high driving force in this periodslarge supersaturation as can be seen on Figure 11b,d. The supersaturation achieves its peak value at about 10 min (see Figure 11d) yielding the appearance of a second mode due to secondary nucleation, which develops clearly by t ) 17 min. The optimal temperature profile, shown in Figure 11a, was

obtained using 50 discretization points in the optimization. The selection of the number of discretizations is a trade-off between computational time and accuracy. No significant improvement in the objective function was achieved when the discretization was increased from 30 to 50 and hence a finer discretization was not considered necessary. The target CSD can be achieved within the temperature range of 40 to 5 °C during a batch period of approximately 1 h with the desired yield of 60%, as shown


in Figure 11c. The temperature profile shows several distinguishing slope changes, which correspond to secondary nucleation events generated to achieve the desired shape of the target distribution. After the initial peak at around 10 min, the supersaturation raises again between 25-40 min and subsequently toward the end of the batch between 50-55 min, as can be seen in Figure 11d. These increases in the supersaturation lead to additional nucleation events and growth, required to achieve the target distribution, as shown in Figure 11e,f. Figure 10 shows that the first mode of the bimodal target distribution is skewed requiring the subsequent nucleation events to improve the shape of the obtained CSD. The third region of increased supersaturation has significantly smaller values and is mainly to facilitate growth of the larger particles to give a better fit to the second mode of the target CSD; Figure 10 shows no significant change in the first mode of the CSD after t ) 41 min, indicating mainly the growth of the second mode. The previous target CSD showed a pronounced secondary peak to illustrate the ability of the approach to achieve a bimodal distribution. In industrial scenarios, the major emphasis is to suppress the nucleation and minimize the formation of small particles during the crystallization process. Therefore the temperature profile was optimized for another fictitious target distribution, with a significantly smaller fraction of fine particles. Figure 12a shows that the CSD predicted by the model is very close to the target distribution. The slight deviation in the fitting of the first peak of the target distribution indicates that it is not possible to achieve any arbitrary target distribution by designing the cooling profile only. There are limitations on the potentially attainable CSD shapes given by the particular nucleation and growth kinetics, as well as the shape of the seed distribution. The optimized temperature profile ranges from 40 to 7 °C in Figure 12b and shows a sequence of gradually increasing cooling rate until t ) 60 min, when the nucleation of the particles required for the second small peak occurs. Figure 12d shows that this nucleation event leads to a decrease in the supersaturation, to a level which is maintained for the next 20 min of the batch, initially by reducing the cooling rate. The cooling rate then is gradually increased again to produce a slightly higher supersaturation level to complete the growth of the particles to the required size and to achieve the required yield of 60%; see Figure 11c. The batch time for this case has increased significantly compared to the previous case to approximately 95 min. The number of discretizations used is 30 since the previous case showed that no significant improvement was obtained by increasing the discretization to 50. Although the secondary peak is significantly reduced in the second case, note that with the desired shape of the primary mode of the target distribution and required yield it is not possible to eliminate completely the secondary peak from the resulting CSD. However the two case-studies illustrate that the proposed method can be used to obtain optimal temperature trajectories for designing different target distributions. The combined QMOM-MOCH, model can predict the CSD efficiently and can be used with off-line or online optimization approaches for designing crystallization systems to produce consistent target CSDs. 6. Conclusions The paper describes a new methodology of solving population balance equations. The approach combines the quadrature method of moments (QMOM) with the method of characteristics (MOCH), and provides a computationally efficient method of reconstructing the crystal size distribution (CSD). The approach can also be used to design optimal temperature profiles for

8583

crystallization product design with a desired shape of the final CSD. By applying the approach in conjunction with control at constant supersaturation, analytical expressions can be derived for the estimation of the CSD for growth-dominated processes. The QMOM-MOCH approach and the analytical CSD estimator are evaluated in the case of the seeded crystallization of potash alum in water with size dependent growth. Experimental and simulation results demonstrate the efficiency of the proposed approach. Optimal temperature trajectories have also been designed for various bimodal target CSDs, suggesting that the approach based on the combined QMOM-MOCH model has the ability to be used for off-line or online optimization of batch crystallization processes. Acknowledgment Financial support provided by the Engineering and Physical Sciences Research Council (EPSRC), U.K., (grant EP/E022294/ 1) is gratefully acknowledged. BASF, Ludwigshafen, Germany, is also acknowledged for partial financial support and for providing experimental data for model identification. Literature Cited (1) Chung, S. H.; Ma, D. L.; Braatz, R. D. Optimal model-based experimental design in batch crystallization. Chemom. Intell. Lab. Syst. 2000, 50, 83. (2) Mullin, J. W. Crystallisation, 4th ed.; Butterworth Heinemann: Oxford, U.K., 2001. (3) Braatz, R. D. Advanced control of crystallization processes. Annu. ReV. Control 2002, 26, 87. (4) Fujiwara, M.; Nagy, Z. K.; Chew, J. W.; Braatz, R. D. First-principles and direct design approaches for the control of pharmaceutical crystallization. J. Process Control. 2005, 15, 493. (5) Nagy, Z. K.; Chew, J. W.; Fujiwara, M.; Braatz, R. D. Comparative performance of concentration and temperature controlled crystallizations. J. Process Control. 2008, 18, 399. (6) Liotta, V.; Sabesan, V. Monitoring and feedback control of supersaturation using ATR-FTIR to produce an active pharmaceutical ingredient of a desired crystal size. Org. Process Res. DeV. 2004, 8, 488. (7) Zhou, G. X.; Fujiwara, M.; Woo, X. Y.; Rusli, E.; Tung, H. H.; Starbuck, C.; Davidson, O.; Ge, Z. H.; Braatz, R. D. Direct design of pharmaceutical antisolvent crystallization through concentration control. Cryst. Growth Des. 2006, 6, 892. (8) Doki, N.; Seki, H.; Takano, K.; Asatani, H.; Yokota, M.; Kubota, N. Process control of seeded batch cooling crystallization of the metastable alpha-form glycine using an in-situ ATR-FTIR spectrometer and an insitu FBRM particle counter. Cryst. Growth Des. 2004, 4, 949. (9) Gron, H.; Borissova, A.; Roberts, K. J. In-Process ATR-FTIR spectroscopy for closed-loop supersaturation control of a batch crystallizer producing monosodium glutamate crystals of defined size. Ind. Eng. Chem. Res. 2003, 42, 198. (10) Hermanto, M. W.; Chiu, M.-S.; Woo, X. Y.; Braatz, R. D. Robust optimal control of polymorphic transformation in batch crystallization. AIChE J. 2007, 53, 2643. (11) Abu Bakar, M. R.; Nagy, Z. K.; Saleemi, A. N.; Rielly, C. D. The impact of direct nucleation control on crystal size distribution in pharmaceutical crystallization processes. Cryst. Growth Des. 2009, 9 (3), 1378. (12) Woo, X. Y.; Nagy, Z. K.; Tan, R. B. H.; Braatz, R. D. Adaptive concentration control of cooling and antisolvent crystallization with laser backscattering measurement. Cryst. Growth Des. 2009, 9, 182. (13) Rawlings, J. B.; Miller, S. M.; Witkowski, W. R. Model identification and control of solution crystallization processes: A review. Ind. Eng. Chem. Res. 1993, 32, 1275. (14) Nagy, Z. K.; Braatz, R. D. Robust nonlinear model predictive control of batch processes. AIChE J. 2003, 49, 1776. (15) Larsen, P. A.; Patience, D. B.; Rawlings, J. B. Industrial crystallization process control. IEEE Control Syst. Mag. 2006, 26, 70. (16) Ward, J. D.; Mellichamp, D. A.; Doherty, M. F. Choosing an operating policy for seeded batch crystallization. AIChE J. 2006, 52, 2046. (17) Nagy, Z. K.; Fujiwara, M.; Braatz, R. D. Modelling and control of combined cooling and anti-solvent crystallization processes. J. Process Control. 2008, 18, 856.

8584


(18) Worlitschek, J.; Mazzotti, M. Model-based optimization of particle size distribution in batch-cooling crystallization of paracetamol. Crys. Growth Des. 2004, 4, 891. (19) Sheikhzadeh, M.; Trifkovic, M.; Rohani, R. Fuzzy logic and rigid control of a seeded semi-batch, anti solvent, isothermal crystallizer. Chem. Eng. Sci. 2008, 63, 991. (20) Sheikhzadeh, M.; Trifkovic, M.; Rohani, R. Adaptive MIMO neurofuzzy logic control of a seeded and an unseeded anti-solvent semi-batch crystallizer. Chem. Eng. Sci. 2008, 63, 1261. (21) Zhang, G. P.; Rohani, S. On-line optimal control of a seeded batch cooling crystallizer. Chem. Eng. Sci. 2003, 58, 1887. (22) Hounslow, M. J.; Reynolds, G. K. Product engineering for crystal size distribution. AIChE J. 2006, 52, 2507. (23) Lee, K., Lee, J. H., Fujiwara, M., Ma, D. L., Braatz, R. D. Runto-run control of multidimensional crystal size distribution in a batch crystallizer, Proceedings of the American Control Conference; IEEE Press: Piscataway, NJ, 2002; p 1013. (24) Rusli, E., Lee, J. H., Braatz, R. D. Optimal distributional control of crystal size and shape. Proceedings of the Fifth World Congress on Particle Technology; Orlando, FL, 2006; Paper 240f. (25) Nagy, Z. K. Model based robust control approach for batch crystallization product design. Comput. Chem. Eng., doi:10.1016/j.compchemeng.2009.04.012. (26) Hulburt, H. M.; Katz, S. Some problems in particle technologysA statistical mechanical formulation. Chem. Eng. Sci. 1964, 19, 555. (27) Randolph, A. D., Larson, M. A. Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallization; Academic Press: San Diego, CA, 1971. (28) Ramkrishna, D. Population Balances, Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, CA, 2000. (29) Gerstlauer, A.; Gahn, C.; Zhou, H.; Rauls, M.; Schreiber, M. Application of population balances in the chemical industryscurrent status and future needs. Chem. Eng. Sci. 2006, 61, 205. (30) McGraw, R. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci. Technol. 1997, 27, 255. (31) Rosner, D. E.; McGraw, R. L.; Tandon, P. Multi-variate population balances via moment- and Monte Carlo simulation methods. Ind. Eng. Chem. Res. 2003, 42, 2699. (32) Marchisio, D. L.; Ptkturna, J. T.; Vigil, R. D.; Fox, R. O. Quadrature method of moments for population balance equations. AIChE J. 2003, 49, 1266. (33) Marchal, P.; David, R.; Klein, J. P.; Villermaux, J. Crystallization and precipitation engineeringsI. An efficient method for solving population balance in crystallization with agglomeration. Chem. Eng. Sci. 1988, 57, 1107. (34) Puel, F.; Fevotte, G.; Klein, J. P. Simulation and analysis of industrial crystallization processes through multidimensional population

balance equations. Part 1: A resolution algorithm based on the method of classes. Chem. Eng. Sci. 2003, 58, 3715. (35) Gunawan, R.; Fusman, I.; Braatz, R. D. Parallel high resolution simulation of particulate processes with nucleation, growth, and aggregation. AIChE J. 2008, 54, 1449. (36) Zhang, T.; Tade´, M. O.; Tian, Y. C.; Zang, H. High-resolution method for numerically solving PDEs in process engineering. Comput. Chem. Eng. 2008, 2, 2403. (37) Haseltine, E. L.; Patience, D. B.; Rawlings, J. B. On the stochastic simulation of particulate systems. Chem. Eng. Sci. 2005, 60, 2627. (38) Flood, A. E. Thoughts on recovering particle size distributions from the moment form of the population balance. DeV. Chem. Eng. Mineral Process. 2002, 10, 501. (39) Ono, T.; Kramer, H. J. M.; Ter Horst, J. H.; Jansens, P. J. Process modeling of the polymorphic transformation of L-glutamic acid. Cryst. Growth Des. 2004, 4, 1161. (40) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 2003, 258, 322. (41) Gordon, R. G. Error bounds in equilibrium statistical mechanics. J. Math. Phys. 1968, 9, 655. (42) Gimbun, J., Nagy, Z. K., Rielly, C.D. Simultaneous quadrature method of moments for the solution of population balance equations, using a differential algebraic equation framework. Ind. Eng. Chem. Res., 10.1016/ j.compchemeng.2009.04.012. (43) Nagy, Z. K.; Fujiwara, M.; Woo, X. Y.; Braatz, R. D. Determination of the kinetic parameters for the crystallization of paracetamol from water using metastable zone width experiments. Ind. Eng. Chem. Res. 2008, 47, 1245. (44) Brecevic, L.; Garside, J. On the measurement of crystal size distributions in the micrometer size range. Chem. Eng. Sci. 1980, 36, 867. (45) Aamir, E.; Nagy, Z. K.; Rielly, C. D.; Kleiner, T.; Judat B. Efficient crystal size distribution estimation approach for growth dominated crystallisation processes. Proceedings of 17th International Symposium of Industrial Crystallisation, Maastricht, The Netherlands, 14-17th September, 2008; Vol 3 p 1733. (46) Ruf, A.; Worlitschek, J.; Mazzotti, M. Modeling and experimental analysis of PSD measurements through FBRM. Part. Part. Syst. Charact. 2000, 17, 167. (47) Hukkanen, E. J.; Braatz, R. D. Measurement of particle size distribution in suspension polymerization using in situ laser backscattering. Sens. Actuators, B 2003, 96, 451.

ReceiVed for reView March 16, 2009 ReVised manuscript receiVed June 13, 2009 Accepted July 10, 2009 IE900430T

Combined Quadrature Method of Moments and Method of

Recommend Documents