On the Design of Optimally Informative Experiments for Dynamic

Ind. Eng. Chem. Res. 2004, 43, 4889-4902

4889

PROCESS DESIGN AND CONTROL On the Design of Optimally Informative Experiments for Dynamic Crystallization Process Modeling Bing H. Chen,† Sean Bermingham,‡ Andreas H. Neumann,‡ Herman J. M. Kramer,‡ and Steven P. Asprey*,†,§ Centre for Process Systems Engineering, Imperial College London, London SW7 2AZ, United Kingdom, Laboratory for Process Equipment, Faculty of Design, Engineering and Production, Delft University of Technology, Leeghwaterstraat 44, 2628 CA, Delft, The Netherlands, and Orbis Investment Advisory Ltd., Orbis House, 5 Mansfield Street, London W1G 9NG, United Kingdom

In this paper, we present the challenging application of now well-established general and systematic procedures for model development, statistical discrimination, and validation to a published large-scale dynamic crystallization process model. Because of the model’s size, this represents, to our knowledge, the first application of such statistical methods to such a largescale dynamic model. For completeness, a brief review of both the model development procedures and the dynamic model are included in the paper. In reviewing the model development procedures, we cover such methods as parametric identifiability testing (to determine whether the parameters, as they appear in the model, can in fact be identified), as well as optimal design of dynamic experiments for both model discrimination among three crystallization models (differing in their kinetics only) and parameter precision improvement within the single “best” dynamic model. Because of the relatively large scale of the model, an optimization-based approach is used for testing of model parameter identifiability that involves semi-infinite programming (SIP) to ensure that the entire control (or input) space has been explored. The problem of designing dynamic experiments is cast as an optimal control problem that enables the calculation of optimal sampling points, experiment durations, fixed and variable external control profiles, and initial conditions of a dynamic experiment subject to general constraints on inputs and outputs. Within this framework, methods are presented to provide experiment design robustness, accounting for parametric uncertainty and subsequently model prediction uncertainty. The paper details the progression of the three crystallization models through the model development procedures and shows the Gahn and Mersmann model (Chem. Eng. Sci. 1999, 54, 1273) to be superior to its competitors. 1. Introduction Crystallization is a well-known process that is widely used in the chemical, and now biochemical, industries to obtain very high purity chemicals that are soluble in a particular solvent. Crystallization is second only to distillation in importance as a separation and purification technique. However, unlike distillation, the mechanisms describing the various physical and kinetic processes occurring in the crystallizer are very poorly understood, and for this reason, until recently, crystallization has received little attention from other fields of research such as equipment design, control, and process modeling. Four types of solution crystallization processes are found in industry and can be differentiated by the way in which the supersaturation is generated. This can be * To whom correspondence should be addressed. E-mail: [email protected]. Tel: +44 (20) 7462-0838. † Imperial College London. ‡ Delft University of Technology. § Orbis Investment Advisory Ltd.

done via evaporation, cooling, precipation/reaction, and salting/drowning out for these cases. By making use of one of these crystallization techniques, a solid crystalline product exhibiting a characteristic crystal size distribution (CSD) is produced from a crystallizer. A crystallizer is designed to create an environment conducive to the formation of crystalline material of a desired quality and quantity. The product quality is determined by the efficiency of the solid/liquid separation in downstream processes such as filtering and washing. The most important quality constraint is the CSD, and the most important considerations in determining the CSD in a crystallizer are those associated with the crystallization kinetics. To tackle the difficult problem of modeling such processes, in this paper, we briefly review and follow the general, systematic procedure to support the development and statistical verification of dynamic models described by Asprey and Macchietto1 and Asprey and Mantalaris.2 Subsequently, we show the application of these procedures to a published dynamic crystallization model13 with three competing kinetics models embedded

10.1021/ie030649n CCC: $27.50 © 2004 American Chemical Society Published on Web 06/29/2004

4890 Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004

Figure 1. Overall model-building scheme.

within a general dynamic crystallizer model framework. As this paper presents the union of two well-established sidessthe model-building procedure on one side and crystallization modeling on the othersthe interested reader is referred to the more extensive reviews of the model-building procedures found in Asprey and Macchietto1 and very extensive review of crystallization processes and their dynamic modeling in Bermingham.13 In general, we consider models of the form

f(x3 (t),x(t),u(t),w,θ,t) ) 0 y(t) ) g(x(t))

(1)

where x(t) is an ns-dimensional vector of time-dependent state variables; u(t) is an nu-dimensional vector of timevarying controls or inputs to the process; w is an nidimensional vector of constant controls; θ is a P-dimensional vector of model parameters to be determined within a continuous, realizable set Θ; and y(t) is an Mrdimensional vector of measured response variables that are functions of the state variables, x(t). In most cases, g(x(t)) will simply be a “selector” matrix, selecting those state variables that are in fact measured. In this general model formulation, we can envision a model taking the form of strictly algebraic equations (AEs), strictly differential equations (ODEs), or a mixed set of differential/ algebraic equations (DAEs). For models such as the crystallization model presented here, partial differential and algebraic equation (PDAE) systems can be entertained through proper discretization of the spatial or crystal number domain. As such, in the case of our crystallization model, the state variables include such variables as solvent and solute mass, enthalpy, and a discretized version of the crystal size distribution. Timevarying controls to the system include such variables

as impeller speed and residence time (in turn affected through a change in flow rate through the crystallizer). Finally, parameters to be estimated from experimental data include such variables as the rate constant for surface integration and surface-related energy increase within the embedded crystallization kinetics model. The measured variables include three quantiles of the crystal size distrubtion (CSD). The overall model-building strategy devised by Asprey and Macchietto1 is depicted in Figure 1. As they point out, the first step of this strategy involves the user specifying one or more models to describe the process at hand. In stage I, the strategy then uses these models to perform preliminary parametric identifiability and distinguishability tests before any data are collected. The question of identifiability is to determine, on the basis of the mathematical behavior of the model itself and a proposed set of variables to be measuredsbefore supporting experimental data have been gathered and analyzedswhether all of the parameters within the model can, in principle, be uniquely identified. The question of distinguishability is to determine whether, given two or more model structures, each can be distinguished from one another. In stage II, the next step in the strategy involves designing experiments for model discrimination, followed by parameter estimation and model adequacy checking to determine whether any of the proposed models can be rejected, with the final purpose of selecting the “best” model from the given set. Once we are left with the best model, stage III is then applied to design experiments for improving the precision of the parameters within that model, eventually arriving at a final model formulation and its most precise parameter values. It must be stressed that, while verifying the model(s) in stages II and III, we always check for

Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 4891

statistical lack-of-fit. Should we find that the model taken forward to stage III has significant lack-of-fit at any point, stage II is revisited with the newly gathered experimental data, and the other models are once again entertained (including any newly devised models). In this paper, we show the application of these statistical techniques to the development of a published dynamic predictive crystallization process model comprising several partial and ordinary differential and algebraic equations.13 To obtain good predictive behavior, a high degree of segregation between (local) kinetics and (overall) hydrodynamics is essential. This separation is achieved by using a compartmental modeling approach. Embedded within a first-principles-based model, three crystallization kinetic models are proposed to model the process, for which we use model discrimination techniques to select the best model from the candidate set. As such, the crystallization application is used to demonstrate the overall procedure of the model development techniques to a relatively large-scale dynamic process model. In the remainder of the paper, in section 2, we first briefly review the mathematical problems faced when implementing the model-building strategy. Next, in section 3, we review the modeling approach used for crystallization processes, from which three competing models arise. In section 4, we cover the complete application of the aforementioned model-building strategy to the dynamic crystallization models devised. Finally, we provide some concluding remarks.

2.1. Preliminary Analysis: Parametric Identifiability and Model Distinguishability. Generally speaking, to test for parametric identifiability, we wish to ensure that, when the same input u(t) is applied to the same model, but with different parameter vectors θ and θ* (such that θ * θ*), the two predicted model output trajectories do not overlay one another. This has been mathematically stated by Asprey and Macchietto1 as follows. Definition: Parametric Identifiability.1 A model giving output trajectory y(θ,u(t),t) is globally identifiable if, for any two parameter sets, θ ∈ Θ and θ* ∈ Θ*, a time horizon of interest, t ∈ [0,tf], for all system inputs u(t) ∈ U, and the same initial conditions y(t)0), the global maximum

max (θ - θ*)T(θ - θ*)

(2)

θ∈Θ,θ*∈Θ*

subject to

∫0T(y(u(t),θ) - y(u(t),θ*))TW(y(u(t),θ) y(u(t),θ*)) dt < y

∀u(t) ∈ U

f(y3 (t),y(t),x(t),u(t),w,θ,t) ) 0 (the model equations)

(3)

l θLi e θi e θU i / /U θ/L i e θi e θi

T (y(1)(t) ∫ 0 u(t)∈ Uθ∈Θ,θ*∈Θ

ΦD ) max min

y(2)(t))TW(y(1)(t) - y(2)(t)) dt (4) subject to

f(1)(y3 (t),y(t),x(t),u(t),w,θ,t) ) 0

2. Review of the Model-Building Scheme

ΦI )

positive values and W is a weighting matrix so as to avoid scaling problems. In words, ΦI is the largest distance between two parameter vectors that still give (essentially) the same response trajectories. If such a distance is arbitrarily small, the model can be deemed globally identifiable. Fortunately, there is a systematic way in which the tolerances ΦI and y can be chosen. In particular, y should be chosen as the integral of the expected value of the noise (or error) that would be experienced during an actual experiment, while ΦI can be chosen to be commensurate with the numerical accuracy to which the model can be solved. To test for distinguishability, we wish to ensure that, when the same input u(t) is applied to two different models with parameter vectors θ and θ*, the two predicted model output trajectories do not overlay one another. Again, this has been mathematically stated by Asprey and Macchietto1 as follows. Definition: Distinguishability.1 Two models giving output trajectories y(θ,u(t),t) and y*(θ*,u(t),t) are distinguishable if, for any θ ∈ Θ, θ* ∈ Θ, t ∈ [0,tf] for all u(t) ∈ U, and the same initial conditions y(t)0), the global minimum

i ) 1, ..., P i ) 1, ..., P

gives ΦI e ΦI, where ΦI and y are arbitrarily small

f(2)(y3 (t),y(t),x(t),u(t),w,θ*,t) ) 0 l θLi

(5) e

θi eθU i

/ /U θ/L i e θi e θi

i ) 1, ..., P i ) 1, ..., P

gives ΦD g ΦD for some ΦD g 0. In words, ΦD is the largest value of the minimum distance between the responses of the two models over the time horizon of interest. If such a value is sufficiently large, the two models can be deemed distinguishable. Again, a value for ΦD g 0 can be chosen to reflect the expected value of the integral of the noise (or error) that would be experienced during an actual experiment. As an advantage, the optimization-based approach is not limited to any particular model form or size. On the other hand, we are now faced with solving semi-infinite or max-min dynamic optimization problems, and we must choose a priori a tolerance for each test (ΦI and ΦD, best chosen to be commensurate with the integral of the expected variance of experimental measurements). However, it is possible to simplify the situation for special model forms. In particular, for nonlinear control-affine models (and models that can be reformulated to this form), we can use the input-output representation by using local functional expansions. By doing so, the input space U is mapped to a domain, X, of the state variables, and we are left with a semiinfinite programming or max-min problem without “dynamics”. A more complete derivation of the optimization-based approach to testing identifiability and distinguishability, as well as solution methods, was developed by Asprey and Macchietto.1


2.2. Designing Dynamic Experiments for Model Discrimination. The problem of model discrimination occurs when NM rival models have been proposed to describe a system and it is not certain which of the NM models is best. For any experiment, it is expected that the best model will provide the most accurate prediction of the observed measurements. Experimental designs optimized for model discrimination are termed Toptimal. Although there have been several approaches to designing dynamic experiments for model discrimination (including such works as refs 3 and 4), we take the approach outlined in Chen and Asprey,5 who used a dynamic extension of the Buzzi-Ferraris criterion6 for the design of dynamic experiments. In the case of dynamic design, we sum the model divergence over all sampling times

× Mr). To design a dynamic experiment for model discrimination, the objective function in eq 7 is embedded into the dynamic optimization framework presented by Asprey and Macchietto1 and included in Appendix A for completeness. 2.3. Designing Dynamic Experiments for Improving Parameter Precision. The aim of experimental design for improved parameter precision is to decrease the size of the inference regions of each of the parameters in any given model. This is equivalent to making the elements of the parameter variance-covariance “small”. In general, all previous methods for designing dynamic experiments for improving parameter precision start with the basis of the highest posterior density (HPD) region, where, in particular, the approximate marginal posterior density covariance is given by (see, e.g., Bard7)

nsp

max φ

∑ Tl,l′,sp(φ,θˆ l,θˆ l′) sp)1

(6)

in which

Tl,l′,sp(φ,θˆ l,θˆ l′) ) (yˆ l(x,θ,tspk) yˆ l′(x,θ,tspk))TSl,l′-1(yˆ l(x,θ,tspk) - yˆ l′(x,θ,tspk)) (7) where Sl,l′ ) S + Wl + Wl′ and S is the covariance matrix of experimental error. Subsequently, we can define

where

[

Wl ) Vl,spΣθ,l-1Vl,spT

∂yˆ l,1(tsp1,θˆ l) ∂yˆ l,1(tsp1,θˆ l)

∂ˆ θl,1 ∂θˆ l,2 ∂yˆ l,2(tsp1,θˆ l) ∂yˆ l,2(tsp1,θˆ l)

Vl,sp )

∂θˆ l,1

∂θˆ l,2

... ...

∂yˆ l,1(tsp1,θˆ l) ∂θˆ l,Pl ∂yˆ l,2(tsp1,θˆ l) ∂θˆ l,Pl

l l l l ∂yˆ l,M(tsp1,θˆ l) ∂yˆ l,M(tsp1,θˆ l) ∂yˆ l,M(tsp1,θˆ l) ... ∂θˆ l,1 ∂θˆ l,2 ∂θˆ l,Pl

]

(8)

T -1 σ˜ rs ∑ ∑ 1 Sr Ss] r)1s)1

V(θˆ ,φ) ) [

(12)

where Sr is the matrix of partial derivatives of the rth response variable in the model with respect to the parameters θ calculated at all experimental points

Sr )

|

∂y(φ,θ,t) ∂θ

θ)θˆ

(13)

σ˜ rs ˜ 1 is the rsth element of the inverse of the estimate Σ of the variance-covariance matrix of the residuals Σ ) cov(yr, ys), estimated by N

σ˜ rs )

(9)

Vl,sp is defined as the sensitivity matrix for model l at sampling time tsp1, and Σθ,l is covariance matrix of the parameters for model l. The parametric sensitivity information for models comprising a set of differential and algebraic equations can be obtained by the solution of the first-order derivative of the model equations

∂y ∂g(x) ∂x ) ∂θ ∂x ∂θ

Mr M r

(10)

where, by partial differentiation of eq 1 with respect to the parameter vector θ

[(yri - fri(θˆ ))(ysi - fsi(θˆ ))] ∑ i)1 N-1

(14)

To compare the magnitudes of different matrices, various real-valued functions have been suggested as metrics. Three common criteria are D-optimality (maximizing the determinant of the information matrix MI), E-optimality (maximizing the smallest eigenvalue of the information matrix MI), and A-optimality (maximizing the trace of the information matrix MI). Various authors have extended these concepts to the dynamic case, with the first attempt made by Espie and Macchietto,8 followed by the works of Zullo9 and the more recent work of Ko¨rkel et al.10 Here, we use the unique criterion put forth by Zullo,9 who defined a dynamic information matrix for experiment design for the improvement of parameter precision in nonlinear, multiresponse, dynamic situations and who developed the optimization framework described elsewhere (Asprey et al.11), also included in Appendix A for completeness. The design criterion requires maximization of a metric of the information matrix, defined in the following way for dynamic models Mr Mr

∂x3 ∂x ∂f fy3 + fy + )0 ∂θ ∂θ ∂θ

(11)

As in the case of the original Buzzi-Ferraris6 method, “in order to discriminate among more than two rival models, the criterion requires the choice for the new experimental design vector, φ, that maximizes eq 7 over any pair (l, l′) of models (l * l′).” The procedure is stopped when no setting of φ maximizes eq 7 above (nsp

MI(θˆ ,φ) ≡

T σ˜ rs ∑ ∑ 1 Qr Qs r)1s)1

(15)

where θˆ is the vector of the best available estimates of the model parameters and φ is the vector of experiment design variables (defined in detail below). The matrix Qr is the matrix of sensitivity coefficients of the rth response variable in the model computed at each of the nsp sampling points and is defined as


Again, to design a dynamic experiment, this objective function is embedded into the dynamic optimization framework of Asprey and Macchietto,12 as detailed in Appendix A. 3. Crystallisation Process Model 3.1. Review of the Crystallizer Model. We use a compartmental modeling approach to model the crysCODE 22-L DTB crystallizer, an evaporative crystallizer with a fines destruction loop (Figure 2). The compartmental modeling approach was selected to account for the spatial distribution of process conditions, such as the supersaturation, energy dissipation, and CSD, throughout the crystallizer, thus allowing separation of the local intrinsic kinetics and the overall hydrodynamics. Essential in this approach is that all compartments within the compartmental model are described with the same compartment model, i.e., the same equations of conservation, physical and thermodynamic property relations, kinetic rate expressions, and parameters. Differences between compartments with respect to nucleation, growth, dissolution, attrition, breakage, and aggregation rates are therefore purely a result of spatial variations in process conditions. One

Figure 3. Compartment model: variables and structure.

of the key features of this approach is the choice of the kinetic model and the estimation of the kinetic parameters in these models to be used as the basic building block for the construction of compartmental models (for a further discussion on why this is so and the compartmentation procedure, refer to Bermingham13). In this work, we focus on the choice of the kinetic model and the estimation of the kinetic parameters for the wellmixed single-compartment model. For this, experimental data were used from a 22-L DT evaporative crystallizer (see Neuman et al.18). 3.1.1. Model of a Single Compartment. A compartment is defined as a well-mixed zone/volume with respect to process conditions such as supersaturation, energy dissipation rate, CSD, temperature, and pressure. Consequently, a compartmental model for a crystallizer must be set up such that negligible gradients in process conditions exist within the compartments. Considerable gradients might, of course, exist between the various compartments. A single-compartment model is schematically depicted in Figure 3, showing the states within the compartment.

Figure 2. Comprehensive compartmental model for the DTB crystallizer: (a) location of compartments and (b) connectivity diagram.


An arbitrary number of inlet streams and outlet streams, respectively, enter and exit the compartment. Each flow is characterized by an overall flow rate, component mass flow rates, crystal number densities, and an enthalpy flow rate. Here, we present the equations of conservation for a single-solute single-solvent system, such as the crystallization of ammonium sulfate from water. This system is used for all experimental and modeling work in the remainder of this paper. Hence, we have two liquidphase components (NCL ) 2), one solid phase (NPS ) 1), and no liquid-phase reactions (NRL ) 0). We also have liquid-phase component 1 as the solute (A), liquidphase component 2 as the solvent (B), and solid phase as the crystals (C). Using these definitions, the crystallization reaction in the above-mentioned system can be written as follows

υS,1A + υS,2B a C

and particle length, respectively, can be written as follows

(17)

For the crystallization of ammonium sulfate (A) from water (B) to crystal (C), stoichiometric coefficients υS,1 and υS,2 are equal to -1 and 0, respectively. 3.1.2. Mass Balance for Solvent. With these newly introduced definitions, we can write the liquid-phase component mass balance for the solvent as

with initial condition n(L,t)0) ) n0(L). The initial condition can be used to express a trivial situation such as a clear liquid (no crystals) or a size distribution of seed crystals (often used in batch processes). The classical boundary condition with respect to crystal size for the population balance equation is19

n(L)0,t) )

B0(t) G(L)0,t)

G(L) g 0, ∀L

(21)

This boundary condition is applicable and sufficient when the crystal growth rate is positive for all crystal sizes. B0 denotes the birth rate of crystals with size zero. Depending on the kinetic model employed, nucleation is modeled either using a birth rate at the boundary (e.g., Ploss et al.14) or as a distributed process with respect to crystal size (e.g., O Ä Meadhra et al.,15 Gahn and Mersmann16). If the growth rate is negative for all sizes, a boundary condition at infinite length is required and sufficient with initial condition mL,2(t)0) ) mL,2,0. 3.1.3. Mass Balance for Solute. By assuming that the vapor flow consists of only solvent, we obtain the following equation for the mass accumulation rate of solute in the liquid phase

with initial condition mL,1(t)0) ) mL,1,0. 3.1.4. Population Balance for Crystals. The time evolution of a crystal size distribution (CSD) is given by the population balance equation (PBE). The PBE for a uniformly mixed volume, with the amount and the size of particles expressed in terms of number density

n(L)∞,t) ) 0

G(L) < 0, ∀L

(22)

If the growth rate is negative for certain crystal size intervals and positive for others, the PBE requires more than one boundary condition with respect to crystal size. 3.1.5. Enthalpy Balance. Finally, the dynamics of the temperature are given by the energy balance

with initial condition H(t)0) ) H0.


Figure 4. Model structure of the framework of Gahn and Mersmann.29

3.1.6. Crystallization Kinetic Models. The notion of modeling the CSD in a crystallizer using the model assumptions associated with the mixed-suspensionmixed-product-removal (MSMPR) concept has been used to simulate the CSD present within a steady state in a crystallizer.21,22 According to MSMPR, the population density function n(x) is the function of the kinetic growth rate of the crystal (Gkin) and residence time (τ) of suspension in the crystallizer

ln[n(x)] ) ln(n0) -

x Gkinτ

(24)

where n0 is the population density in the first size class formed by a secondary nucleation process. However, MSMPR cannot follow the so-called downward and upward curvatures in the large and small size ranges of real distribution measurements. Because of the lack of accurate knowledge of the birth function and initial growth rate of secondary nuclei, the problem was overcome by assuming that all particles grow at the same rate, and an “effective” nucleation rate was used whereby crystals are born at the same size close to zero. Many correlations have been proposed for the effective rate of secondary nucleation (B0,eff), which is dependent on a number of factors expressed as a power-law equality, namely, a term describing the state of the fluid phase in the crystallizer, such as the supersaturation (σ), which relates to the kinetic growth rate of the crystals (Gkin(σ)); the total mass of solid in slurry (Mp); and a term for the impeller frequency (Nimp) or power input via impeller ()

B0,eff ∼ σbNimp hMpj ∼ GkinikMpj

(25)

Over the years, quite a number of kinetic models have been developed for secondary nucleation and growth dominated systems. O Ä Meadhra20 discussed the models of Ottens et al,23 Ploss et al,14 Jager,24 and van der Heijden et al.25 He found that only models with nucleation rates depending on the number and size of the parent crystals present were able to fit observed dynamics in the CSD. The two kinetic models developed in the previous phase of the crysCODE project at TU Delft contained similar relations between the nucleation rates and the CSD.26,15 Both also include the effect of size-dependent attrition on the effective crystal growth

rate. Thus

B0,eff ) n0Geff(σ,x)

(26)

In the large size range (downward curvature)

Geff(σ,x) ) Gkin(σ) - Gattr(x)

(27)

Gattr(x) ) K′attr(1 - {1/[1 + (x/xg)n]})

(28)

whereas in the small size range (upward curvature)

Geff(σ,x) ) Gkin(σ,x) )

{[ ( ) ]}

p1σp2(1 - (1 - p3) exp -

x p4

p5

(29)

where x is the length of crystal and K′attr, xg, and n) are unknown parameters evaluated from experimental CSD curves obtained under growing conditions. However, all of these models are largely of an empirical nature: they do not contain the particle mechanics of importance for attrition, and they miss the vital link between the attrition process and the hydrodynamic conditions inside the crystallizer. As a result, many of the parameters need to be redetermined for different crystallizer and/or pump geometries and sometimes even for different operating conditions, thus making these models unreliable for scale-up. More recently, Gahn and Mersmann16 developed a model framework for attrition, secondary nucleation, and growth. This framework contains a more physical basis than any previous work: the secondary nucleation rate caused by crystal-impeller collisions is calculated largely from first principles. Consequently, the number of kinetic parameters is significantly lower. More importantly, these parameters are not of an empirical nature, i.e., they have physical meaning. One can therefore expect this model to have the predictive properties required for design purposes. As this model is expected to exhibit the best predictive properties at this moment in time, it is selected for implementation in the general crystallization process modeling framework. The model framework of Gahn and Mersmann,16 shown in Figure 4, can be seen as a synthesis of three submodels:27 (1) a procedure to determine the total number of crystals colliding per second with the faces

4896 Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 Table 1. Implemented Kinetic Models kinetic model

expected predictive value

nature

number of kinetic parameters

power law purely empirical limited O Ä Meadhra15 partly first principles reasonable Gahn16 largely first principles good

5 9 2

and edges of the impeller blades and the corresponding impact energy per collision;14,28 (2) a relation between the impact energy and the attrition volume produced due to a single collision of a crystal corner with a hard, flat surface and the resulting number and size distribution of resulting attrition fragments;29 and (3) a relation to describe the growth rate of the fragments formed by the attrition.30-32 The calculation of the growth/dissolution within the PBE is based on the following set of equations

G(L) ∆c(L) kd(L) + ) cs 2krcs 2kd(L)

x( )

G(L) ∆c(L) ) cs 2kd(L) kd(L) )

[

kd(L) 2krcs

2

div(G, O)

div(O, P)

div(G, P)

1 (640 rpm) 2 (775 rpm) 3 (910 rpm) optimal

19 038 1766 13 663 45 134

6209 5943 5660 51 861

4077 5486 25 614 72 560

a All three experimental data sets were collected under a heat input of 120 kW/m3 and at the impeller speeds indicated.

kd(L) ∆c(L) krcs cs ∆c(L) > 0 (30) (31)

( )( )] 1/5

data seta

+

∆c(L) e 0

DAB jL4 2 + 0.8 3 L υL

Table 2. Objective Function Values for Model Discrimination Using the Extended Buzzi-Ferarris Criterion

υL DAB

1/3

(32)

∆c(L) ) c - c/real(L)

Figure 5. Model predictions for optimal model discrimination between the Gahn (G) and O Ä Meadhra (O) models. Table 3. Optimally Designed Inputs for Model Discrimination between the Gahn and O Ä Meadhra Models

(33)

tswi (s) (i ) 1, ..., 5)

residence time (s)

impeller speed (rpm)

0 600 1200 53 400 54 000

2000 2000 6500 6500 6500

600 600 600 600 1080

where

j )

NeNimp3Dimp5 V

and

( )

c/real(L) ) c* exp

ΓS RTL

For this study, the crystallization kinetics are described using the three aforementioned competing models (see Table 1 for a summary). The reader is referred to the original papers for further details; in particular, the Gahn and Mersmann16 model is too lengthy to include here. 3.2. Variable Definitions. Here, we provide the reader with the mapping between the variable names used in describing the model-building techniques and the variable names used in describing the crystallization process model. In particular

{mL,1,mL,2,H,V,n(L,t)} ∈ x {n(L,t)} ≈ {X10,X50,X90} ∈ y {kr,ΓS} ∈ θ

(34)

{Q,Nimp(t),φV(t)} ) {u(t),w} ∈ φ Note here that the continuous crystal size distribution is approximated by the 10%, 50%, and 90% quantiles of the distribution. 4. Application of the Model-Building Procedures to the Dynamic Crystallization Model 4.1. Identifiability of the Crystallizer Model. Using the dynamic optimization approach for testing parametric identifiability, as outlined in eqs 2 and 3,

Table 4. Lack-of-Fit Analyses of the Three Models under the Optimally Designed Experiment Conditionsa,b quantile

Gahn

O Ä Meadhra

power law

X10 X50 X90

737.1 (4.03 × 10-11) 688.3 (1.78 × 10-15) 656.3 (0)

27 517.5 (1) 10 224.9 (1) 7025.2 (1)

25 410.4 (1) 11 215.2 (1) 6413.6 (1)

a

Reference χ2 ) 925.7. b Test’s p values are in parentheses.

we tested the identifiability of the dynamic crystallization model with the Gahn kinetics model using ΦI ) 1.0 × 10-5 and y ) 1.0 × 10-4 and W set to the identity matrix (the scaling of the three responsessthe quantiles of the crystal size distributionsare already commensurate with one another). Thus, the model involved two kinetic parameters, kr and ΓS. The solution to this semiinfinite programming problem was found to be

max (θ - θ*)T(θ - θ*) ) 2.528 × 10-12 (eΦI)

θ∈Θ,θ*∈Θ*

showing that the model is indeed identifiable due to the fact that, to ensure that the same model output is obtained, the two instances of the parameter vector, θ and θ*, must be close to one another. 4.2. Designing Dynamic Experiments for Model Discrimination. As previously mentioned, we study the problem of the evaporative crystallization of ammonium sulfate from water. In this case, an accurate measurement of the supersaturation is extremely difficult; however, detailed CSD measurements are available in 2-min intervals for each experiment by use of


Figure 6. Overlay plots of the predictions from the Gahn model with experimental data.

Figure 7. Overlay plots of the predictions from the O Ä Meadhra model with experimental data.

light-scattering measurements. By using more characteristics of the CSD than the median size and/or by using the observed dynamics in the CSD, the absence of a supersaturation measurement can (in principle) be countervailed (see Bermingham13 for a detailed expla-

nation). For parameter estimation purposes, to characterize the CSD, the median size (X50) and the 10% and 90% quantiles (X10 and X90, respectively) were chosen. Although more quantiles of the CSD could have been used, this would also increase the dimensionality of the


Figure 8. Overlay plots of the predictions from the power-law model with experimental data.

experiment design and parameter estimation problems; it is not clear a priori whether including more quantiles would warrant this increase. Three existing data sets were analyzed using the extended Buzzi-Ferraris divergence measure (eq 6), results of which are shown in Table 2. The three existing data sets were generated on a 22-L DT crystallizer to ensure that the relevant process variables are uniformly distributed over the crystallizer. For details of the data sets, the reader is referred to Neumann18 and Bermingham.13 As can be seen in Table 2, experiment 1 is relatively informative for discriminating between the Gahn (G) and O Ä Meadhra (O) models, whereas experiment 3 is relatively informative for discriminating between the Gahn (G) and power-law (P) models. All existing experiments are noninformative for discriminating between the O Ä Meadhra (O) and power-law (P) models. The final row of Table 2 shows results of optimal experiment designs carried out for model discrimination using the extended Buzzi-Ferraris criterion (eq 6) and piecewise-constant control inputs, with a maximum experiment duration of 60 000 s. An a priori selection of the number of sampling points (based on the fixed 2-min intervals of the light-scattering measurements), along with the number of switching times/constant levels of the input profiles, is required. We arbitrarily chose five time intervals for the control vector parametrization (nswi ) 5, i ) 1, ..., nu, where nu ) 2, corresponding to impeller speed and residence time). We impose the constraints

tswi,j - tswi,j-1 g 10.0 h

i ) 1, ..., nu; j ) 1, ..., nsw (35)

As can be seen in Table 2, these optimal designs are at least three times more informative than the existing

Table 5. Objective Function Values for Parameter Precision Using A-, E-, and D-Optimal Criteria for the Gahn Model data seta impeller residence speed (rpm) time (s) A-optimal E-optimal 1 2 3

640.2 4500 755.2 4500 910.2 4500 optimal designs

86 667 92 077 89 519 101 634

4.36 3.22 1.07 42.27

D-optimal 5.016 × 108 1.918 × 108 4.003 × 108 8.028 × 108

a All three experimental data sets were collected under a heat input of 120 kW/m3 and at the impeller speeds indicated.

data sets for the purpose of discrimination. For illustrative purposes, the model predictions for the optimal (G, O) model discrimination design are shown in Figure 5, and the designed inputs are listed in Table 3. Under the optimally designed conditions for discrimination between the Gahn (G) and O Ä Meadhra (O) models, we subsequently simulated additional data using a “true” Gahn model under these experimental conditions, with a homoskedastic normal pure error variance of 6.0% for all three responses. Subsequent parameter estimation was carried out, giving the results in Table 4, in which we report a univariate χ2 lack-offit (LOF) value for each model, along with the “p value” of the test. A small p value indicates statistically insignificant LOF, whereas a high p value (>0.20) indicates significant LOF. These LOF results are supported by the inspection of overlay plots of the model predictions with the experimental data, as presented in Figures 6-8, showing that the Gahn model indeed fits the collected data much better than the other two models. As can be seen in Table 4, as expected, the Gahn model is the only model to pass the lack-of-fit test

Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 4899 Table 6. Optimally Designed Inputs for Parameter Precision Using the A-, E-, D-Optimal Criteria with the Gahn Modela A-optimal

E-optimal

D-optimal

tswi (s) (i ) 1, ..., 5)

residence time (s)


tswi (s) (i ) 1, ..., 5)

residence time (s)


tswi (s) (i ) 1, ..., 5)

residence time (s)


0 12 006 46 870 47 470 48 070

2000 6500 2000 2000 2000

600 600 600 1080 1080

0 12 997 13 597 32 129 53 546

4187 4273 4873 4310 4492

600 600 600 1080 1080

0 53 200 53 800 54 400 55 000

3766 6500 2000 2000 6012

1080 600 1080 600 600

a

Time duration ) 60 000 s.

Table 7. Improvement of Confidence Intervals via Parameter Precision Design Using A-, E-, and D-Optimal Criteria for the Gahn Modela estimate (precision) E

parameter

estimate (initial)

confidence interval

A

kr ΓS

1.51 × 10-4 1.53 × 10-5

(1.037 × 10-6 (3.94 × 10-6

1.52 × 10-4 1.51 × 10-5

a

1.51 × 10-4 1.53 × 10-5

confidence interval (precision) E D

D

A

1.52 × 10-4 1.50 × 10-5

(8.89 × 10-7 (2.36 × 10-7

(8.01 × 10-7 (2.03 × 10-7

(1.59 × 10-6 (3.97 × 10-7

Value of Student’s t-test ) 1.645. Table 8. Correlation Matrices of Parameter Re-estimation via Parameter Precision Design Using A-, E-, and D-Optimal Criteria for the Gahn Model kr

ΓS

kr ΓS

A-optimal 1.000 -0.9967

-0.9967 1.000

kr ΓS

E-optimal 1.000 -0.9968

-0.9968 1.000

kr ΓS

D-optimal 1.000 -0.9991

-0.9991 1.000

Table 9. Lack-of-Fit Analyses of Parameter Re-estimation via Parameter Precision Design Using A-, E-, and D-Optimal Criteria for the Gahn Modela,b Figure 9. Model predictions for the parameter precision experiment design (D-optimal) for the Gahn model.

(calculated LOF < reference χ2) among the three candidates, even after just one optimally designed dynamic experiment. 4.3. Designing Dynamic Experiments for Parameter Precision. Because the Gahn model was determined as the best model in the model discrimination study, we turn to improving the precision in the parameter estimates for this model. We use the A-, E-, and D-optimal criteria outlined in section 2.3 to design experiments for this activity. Three existing data sets are again used to analyze their parametric information content according to the three different criteria; the results are reported in Table 5. It can be seen in Table 5 that experiment 1 gives the largest amount of parametric information according to the E- and D-optimal criteria, whereas the A-optimal criterion shows that experiment 1 is less informative than, but close to, experiment 3. The optimal experiment designs, listed in the final row, give much higher parametric information, as expected. The designed piecewise-constant control inputs are reported in Table 6, in which the parameter precision experiment design starts from the parameter point estimate θ0 ) [1.512 × 10-4, 1.534 × 10-5], as obtained by parameter estimation based on an existing experimental data set. As can be seen in Table 6, each different criterion gives different optimally designed inputs for parameter precision, although the A- and E-optimal designs

quantile X10 X50 X90 a

A-optimal

E-optimal

D-optimal

742.6 (1.09 × 10-10) 762.5 (5.53 × 10-12) 703.6 (5.31 × 10-14) 731.7 (1.48 × 10-11) 679.2 (0) 677.0 (0) 729.6 (1.00 × 10-11) 662.4 (0) 764.8 (4.53 × 10-9)

Reference χ2 ) 925.7. b Test’s p values are in parentheses.

are closer. The model predictions for the precision design using the D-optimal criterion are illustrated in Figure 9. Table 7 reports values of the estimated model parameters and confidence intervals obtained by using data collected at the designed experiment conditions for improving parametric precision, designed by each of the three criteria. As can be seen, the confidence intervals decrease for the A- and E-optimal criteria, whereas the confidence interval of kr increases for the D-optimal design. With the D-optimal criterion, in trying to shrink the confidence interval of one parameter, it is possible that the remaining parameters are punished. Because the D-optimal criterion is trying to minimize the area of the confidence ellipsoid, it sometimes tries to “squeeze” the ellipse along the minor axis, which, in turn, can cause stretching in the major axis. This also means a possible increase in the correlation of the parameter estimates, as can be found in Table 8. Finally, Table 9 shows the lack-of-fit analyses of the parameter reestimation via parameter precision design using the A-, E-, and D-optimal criteria for the Gahn model. All cases show no significant LOF, as indicated by the small p values.


5. Concluding Remarks In this paper, we have briefly reviewed the general, systematic procedure put forth by Asprey and Macchietto1 to support the development and statistical verification of dynamic process models. Within this procedure, methods were presented to address several key aspects, such as parametric identifiability and model distinguishability testing, as well as optimal design of dynamic experiments for both model discrimination and parameter precision improvement. We showed the application of these statistical techniques to the development of an already-published dynamic predictive crystallization process model.13 Embedded within a firstprinciples-based model, three crystallization kinetic models were proposed to model the process, for which we used model discrimination techniques to select the best model from the candidate set. We subsequently used dynamic experiment design for parameter precision with the Gahn kinetics model to show how the precision of the parameter estimates could be improved through use of our methods. As such, the crystallization application was used to demonstrate the overall procedure, applied to a relatively large-scale and challenging dynamic process modeling exercise. Acknowledgment The authors acknowledge funding of this project by the EPSRC. Nomenclature A ) ammonium sulfate B ) water B0 ) birth rate of crystals at L0 (s-1) C ) crystal c ) molar concentration (mol‚m-3) c* ) molar saturation concentration of stress-free crystals (mol‚m-3) / creal(L) ) molar saturation concentration of crystals of size L (mol‚m-3) cs ) molar concentration in the solid phase (mol‚m-3) DAB ) binary diffusion coefficient (m2‚s-1) Dimp ) impeller diameter (m) f, g ) analytic vector fields (ns-dimensional) G ) linear crystal growth rate (m‚s-1) Gkin ) kinetic growth rate of the crystal (m‚s-1) H ) enthalpy (J) kd ) mass-transfer coefficient (m‚s-1) kr ) rate constant for surface integration (m‚s-1) L ) crystal length (m) m ) mass (kg) M ) molar mass (kg‚mol-1) MI ) information matrix Mp ) total mass of solid in slurry during crystallization Mr ) the number of measured response variables N ) number of experimental data sets Ne ) Newton number or power number NM ) number of rival models n(L,t) ) crystal number density (m-3‚m-1) n(x) ) population density function NC ) number of components NI ) number of solid/liquid inlet streams Nimp ) impeller speed (rpm) NM ) number of rival models in model discrimination NP ) number of solid phases nsp ) number of sampling points (or times) nsw ) number of switching points (or times) nu ) number of time-varying process inputs or controls

P ) number of parameters Qr ) sensitivity matrix R ) ideal gas constant (8.314 J‚mol-1‚K-1) S ) covariance matrix of experimental error t ) time (s) tspi ) ith sampling point of the response variables tswi,j ) jth switching time of the ith time-varying process input T ) temperature (K) T(φ) ) model discrimination criterion u(t) ) nu-dimensional vector of time-varying process inputs V ) size of control volume (m3) Vl,sp ) dynamic sensitivity matrix for model l at sampling time tsp w ) ni-dimensional vector of constant controls W ) weighting matrix WC ) critical work of indentation for crack formation (J) x(t) ) ns-dimensional vector of state variables yˆ ) estimated or predicted response variables y(t) ) M-dimensional vector of system outputs or response variables y0 ) M-dimensional vector of initial conditions Greek Symbols φ ) flow rate in crystallization Σ ) M × M covariance matrix of the response variables θ ) P-dimensional vector of model parameters θˆ ) present best point estimate of the parameter vector θ φ ) vector of experiment design variables Σθ ) M × M covariance matrix of the parameters χ2 ) χ2 statistic ΦD ) defined in eq 4 ΦI ) defined in eq 2 ) specific power input (W‚kg-1) y ) arbitrarily small positive value ΦD ) small positive value best chosen to be commensurate with the integral of the expected variance of experimental measure ΦI ) small positive value best chosen to be commensurate with the integral of the expected variance of experimental measure ΓS ) surface-related energy increase or condition of deformation (J‚m‚mol-1) φH ) enthalpy flow rate or enthalpy flux (J‚s-1) φm ) mass flow rate or mass flux (kg‚s-1) φmol ) molar flow rate or molar flux (mol‚s-1) φn ) number density production rate (m-3‚m-1‚s-1) φn,coll ) crystal collision rate at a certain radial position (m-1‚m-1‚s-1) φV ) volumetric flow rate or volumetric flux (m3‚s-1) φV,dt ) volumetric flow rate through draft tube (m3‚s-1) ηcoll ) collision efficiency (m-1) ηgeom ) geometric target efficiency (m-1) ηL ) dynamic viscosity of the liquid (N‚m-2‚s) ηL/S ) apparent dynamic viscosity of the slurry (N‚m-2‚s) ηtarget ) target efficiency µ ) shear modulus (N‚m-2) νL ) kinematic viscosity of the liquid (m2‚s-1) ν ) stoichiometric coefficient matrix FL ) material density of the liquid phase (kg‚m-3) FS ) material density of the solid phase (kg‚m-3) σ ) relative supersaturation Superscripts L ) lower bound U ) upper bound Subscripts aggl ) agglomeration attr ) attrition break ) breakage

Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 4901 dis ) dissolution edge ) edge of the impeller blade face ) face of the impeller blade grow ) growth L ) liquid phase, i.e., the continuous phase nucl ) primary nucleation S ) solid phase V ) vapor phase

x0Li e x0i e x0Ui

g(x0) g 0 wLi e wi e wU i

AE ) algebraic equation CFD ) computational fluid dynamics CSD ) crystal size distribution DAE ) differential/algebraic equation DT ) draft tube (crystallizer) DTB ) draft tube baffle (crystallizer) MSMPR ) mixed-suspension-mixed-product removal ODE ) ordinary differential equation PBE ) population balance equation PDAEs ) partial differential and algebraic equations

Appendix A. Dynamic Optimization Framework1 When using either model discrimination criteria or parameter precision criteria to design a dynamic experiment, we are faced with solving a nonlinear dynamic optimization problem. Within this optimization framework, a design variable, tsp,l, is associated with each discrete sampling point l ) 1, ..., nsp, of the measured responses. The time-varying inputs to the process, u(t), are assumed to be piecewise-constant, piecewise-linear, or piecewise-quadratic functions of time defined over a number of control intervals. Zeroth- or higher-order continuity can be enforced at the interval boundaries. For the piecewise-constant case, each control variable, ui(t), is defined in terms of its value (zi,j) over each interval j, as well as the “switching” times delineating these intervals (tswi,j)

tswi,j-1 e t e tswi,j

i ) 1, ..., nu; j ) 1, ..., nswi

(A1)

i ) 1, ..., nu; j ) 1, ..., nswi (A2)

where nu is the number of time-varying controls and nswi is the number of switching points (defined a priori) associated with the ith control. All zi,j and tswi,j are design variables, and so are the initial conditions of the experiment, x0, and any timeinvariant controls w. Overall, then, the vector of design variables, φ, is defined as

φ)

[

tspi, i ) 1, ..., nsp; (tswi,j, zi,j, i ) 1.., nu; j ) 1, ..., nswi); x0T; w

]

(A4) (A5)

i ) 1, ..., nw

min max e tspl - tspl-1 e ∆tsp ∆tsp l l

Abbreviations

ui(t) ) zi,j

i ) 1, ..., M

l ) 1, ..., nsp

(A6) (A7)

min min and ∆tsp are, respectively, the minimum where ∆tsp l l and maximum allowable intervals between consecutive sampling points. To allow maximum versatility, independent limits are allowed for each interval between consecutive sampling points. As it turns out, these constraints are active almost all of the time, as the experimental design criteria attempt to place all sampling points as closely as possible to the maxima of the sensitivity coefficients. To avoid mathematical singularities caused by the collapse of one or more control intervals, the following linear constraints on the time-varying control switching points are introduced

min max e tswi,j - tswi,j-1 e ∆tsw ∆tsw i,j i,j

i ) 1, ..., nu; j ) 1, ..., nswi (A8) Simple bounds are also imposed on the levels of the piecewise-constant controls within each control interval U z,Lij e zi,j e zi,j

i ) 1, ..., nu; j ) 1,..,nswi (A9)

Finally, the constraint

tswi,nsw e tspnsp

i ) 1, ..., nu

(A10)

indicates that there is no point in effecting any changes to the control variables past the last sampling point. The above formulation defines, in a very comprehensive and versatile manner, a large range of dynamic experimental conditions. Although this is an attractive feature, very large optimization problems might arise, even for small dynamic systems. For instance, the experiment design problem for a model with only 2 response variables, 10 sampling times, 2 variable initial conditions, 2 time-varying controls, 5 switching times for each time-varying control, and the corresponding piecewise levels involves ca. 30 design variables. These dynamic optimization problems are solved using a control vector parametrization approach.17 Literature Cited

(A3)

The above formulation allows for a different number of intervals and control levels of each time-varying control variable, which can be switched at different times. Here, it is assumed that the sampling points are the same for all response variables; however, the formulation could easily be altered to permit different sampling times for each response variable. A variety of constraints can be imposed on the design variables. There are simple bounds and optional nonlinear equality/inequality constraints on the initial conditions, simple bounds on the constant inputs (w), and linear constraints on the sampling point distribution

(1) Asprey, S. P.; Macchietto, S. Statistical Tools for Optimal Dynamic Model Building. Comput. Chem. Eng. 2000, 24, 1261. (2) Asprey, S. P.; Mantalaris, A. Global Parametric Identifiability of a Dynamic Unstructured Model of Hybridoma Cell Culture. In Proceedings of the 8th IFAC International Conference on Computer Applications in Biotechnology (CAB-8); Dochain, D., Perrier, M., Eds.; Elsevier Science & Technology: New York, 2001. (3) Espie, D. M. The Use of Nonlinear Parameter Estimation for Dynamic Chemical Reactor Modeling. Ph.D. Thesis, University of London, London, U.K., 1986. (4) Cooney, M. J.; McDonald, K. A. Optimal Dynamic Experiments for Bioreactor Model Discrimination. Appl. Microbiol. Biotechnol. 1995, 43, 826. (5) Chen, B. H.; Asprey, S. P. On the Design of Optimally Informative Dynamic Experiments for Model Discrimination in Multiresponse Nonlinear Situations. Ind. Eng. Chem. Res. 2003, 42, 1379.

4902 Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 (6) Buzzi-Ferraris, G.; Forzatti, P. Sequential Experimental Design for Model Discrimination in the Case of Multiple Responses. Chem. Eng. Sci. 1984, 39, 81. (7) Bard, Y. Nonlinear Parameter Estimation; Academic Press: London, 1974. (8) Espie, D. M.; Macchietto, S. The Optimal Design of Dynamic Experiments. AIChE J. 1989 35, 223. (9) Zullo, L. Computer Aided Design of Experiments. An Engineering Approach. Ph.D. Thesis, University of London, London, U.K., 1991. (10) Ko¨rkel, S.; Bauer, I.; Bock, H. G..; Schlo¨der, J. P. A Sequential Approach for Nonlinear Optimum Experimental Design in DAE Systems. In Scientific Computing in Chemical Engineering II: Simulation, Image Processing, Optimization, and Control; Keil, F., Mackens, W., Voss, H., Werther, J., Eds.; Springer: Berlin, 1999; Vol. 2, pp 338-345. (11) Asprey, S. P.; Macchietto, S.; Pantelides, C. C. Robust Optimal Designs for Dynamic Experiments. In Proceedings of the IFAC Advanced Control of Chemical Processes (ADCHEM) 2000; Biegler, L. J., Bramilla, A., Scali, C., Marchetti, G., Eds.; Elsevier Science Ltd.: New York, 2001; Vol. 3, pp 845-850. (12) Asprey, S. P.; S. Macchietto. Designing robust optimal dynamic experiments. J. Process Control 2002, 12, 545. (13) Bermingham, S. K. A design procedure and predictive models for solution crystallization processes: Development and application. Ph.D. Thesis, Delft University of Technology, Delft, The Netherlands, 2003. (14) Ploss, R.; Mersmann, A. A New Model of the Effect of Stirring Intensity on the Rate of Secondary Nucleation. Chem. Eng. Technol. 1989, 12, 137. (15) O Ä Meadhra, R. S.; Kramer, H. J. M.; Rosmalen, G. M. Model for Secondary Nucleation in a Suspension Crystallizer. AIChE J. 1996, 42, 973. (16) Gahn, C.; Mersmann, A. Brittle Fracture in Crystallization Processes Part A. Attrition and Abrasion of Brittle Solids; Part B. Growth of Fragments and Scale-up of Suspension Crystallizers. Chem. Eng. Sci. 1999, 54, 1273. (17) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a Class of Multistage Dynamic Optimization Problems 1: Problems without Path Constraints. Ind. Eng. Chem. Res. 1994, 33, 2111. (18) Neumann, A. M. Characterisation of Industrial Crystallisers of Different Scale and Type. Ph.D. Thesis, Delft University of Technology, Delft, The Netherlands, 2001. (19) Bermingham, S. K.; Verheijen, P. J. T.; Kramer, H. J. M. Optimal design of solution crystallization processes with rigorous models. Trans. Inst. Chem. Eng. A 2003, 81, 893.

(20) O Ä Meadhra, R. Modeling of the kinetics of suspension crystallizerssA new model for secondary nucleation. Ph.D. Thesis, Delft University of Technology, Delft, The Netherlands, 1995. (21) O’Dell, F. P.; Rousmalen, R. W. Magma density and dominant size for size dependent growth. AIChE J 1987, 24, 738. (22) Rojkowski, Z. H. Crystal growth rate models and similarity of population balances for size dependent growth rate and for constant growth rate dispersion. Chem. Eng. Sci. 1993, 48, 1478. (23) Ottens, E. P. K.; Janse, A. H.; DeJong, E. J. Secondary nucleation in a stirred vessel cooling crystallizer. J. Crystal Growth 1972, 13/14, 500. (24) Jager, J. Control of Crystallizers: The Physical Aspects. Ph.D. Thesis, Delft University of Technology, Delft, The Netherlands, 1990. (25) Van der Heijden, A. E. D. M.; Elwenspoek, M. Contact nucleation: In-situ and ex-situ observations of surface breeding. J. Crystal Growth 1990, 99, 1087. (26) Eek, R. A. Control and dynamic modeling of industrial suspension crystallizers. Ph.D. Thesis, Delft University of Technology, Delft, The Netherlands, 1995. (27) Neumann, A. M.; Bermingham, S. K.; Kramer, H. J. M.; Rosemalen, G. M. van. The effect of the impeller speed on the product crystal size distribution (CSD) in a 22-liter draft tube (DT) crystallizer. J. Crystal Growth 1999, 198/199, 732. (28) Mersmann, A.; Sangl, R.; Kind, M.; Pohlisch, J. Attrition and secondary nucleation in crystallizers. Chem. Eng. Technol. 1988, 11, 80. (29) Gahn, C.; Mersmann, A. Theoretical prediction and experimental determination of attrition rates. Trans. Inst. Chem. Eng. A 1997, 75, 125. (30) van der Heijden, A. E. D. M. Secondary nucleation and crystallization kinetics. Ph.D. Thesis, University of Nijmegen, Nijmegen, The Netherlands, 1992. (31) Gahn, C. Die Festigkeit von Kristallen und ihr Einfluss auf die Kinetik in Suspension-kristallisatoren. Ph.D. Thesis, Mu¨nchen University of Technology, Herbert Utz Verlag Wissenschaft, Mu¨nchen, Germany, 1997. (32) Zacher, U. Die Kristallwachstumsdispersion in einem kontinuierlichen Suspensions-kristallisator. Ph.D. Thesis, Mu¨nchen University of Technology, Mu¨nchen, Germany, 1995.

Received for review August 8, 2003 Revised manuscript received May 12, 2004 Accepted May 17, 2004 IE030649N

On the Design of Optimally Informative Experiments for Dynamic

Recommend Documents