ARTICLE pubs.acs.org/crystal
Internal Fines Removal Using Population Balance Model Based Control of Crystal Size Distribution under Dissolution, Growth and Nucleation Mechanisms Z. K. Nagy,* E. Aamir, and C. D. Rielly Loughborough University, Loughborough, Leicestershire, LE11 3TU, United Kingdom ABSTRACT: The paper presents a methodology for the dynamic optimization of the crystal size distribution (CSD) considering growth, nucleation, and dissolution mechanisms for batch cooling crystallization processes. The optimal temperature trajectories are obtained by solving the population balance model (PBM) using the quadrature method of moments in conjunction with the method of characteristics. The coupled algorithm allows a computational efficient approach for the estimation of the shape of the CSD during crystallization processes with generic apparent size-dependent growth (caused by the actual mechanism or growth rate dispersion) and dissolution and secondary nucleation mechanisms. The approach is evaluated for the crystallization of potassium alum in water. The model parameters for growth and dissolution are identified based on pilot scale industrial and laboratory experimental data, respectively, and are used to design the optimal dynamic temperature trajectories to obtain the desired monomodal target shape of the CSD at the end of the batch. Simulation results illustrate that the incorporation of the dissolution mechanism in the model allows development of optimal in situ fine removal policies, which programmatically drive the process in the undersaturated region offering more flexibility in shaping the CSD. The online application of the controlled dissolution-based optimal control method can adapt the operating policy in the case of accidental seeding, which is a common disturbance in industrial crystallization processes and can eliminate the need for additional equipment used for external fines removal loop.
1. INTRODUCTION Crystallization is a very complex operation, involving a number of mass transfer or chemical reaction steps. In a supersaturated solution, the first step is the creation of a new surface by nucleation, followed by diffusion of the solute to the surface, and its adsorption and integration into the crystal structure. All these steps are governed by different physical laws. The reverse process of crystal growth is dissolution, which can be described as occurring in two steps: (i) surface reaction and detachment of the surface species in which crystals are disintegrated at the crystal surface and (ii) mass transfer of these species from the surface toward the bulk solution across a diffusion layer, which surrounds the crystals.1,2 The overall process is controlled by the slowest step. The dissolution kinetics of the readily soluble mineral salt is limited by mass transfer, while the dissolution of a sparingly soluble salt is found to be controlled by the kinetics of the detachment of the species from the surface of the crystals. During the dissolution process the solution is undersaturated. The dissolution kinetics of crystals has been studied extensively in the pharmaceutical domain, because the rate of dissolution affects the bioavailability of drug crystals. Dissolution tests are generally performed on the pharmaceutical dosages and are mainly used for quality control. In industrial processes, crystallization is an important operation, which affects the physical characteristics of the product, such as the crystal habit, size r 2011 American Chemical Society
distribution, and polymorphism, and in turn these determine the in vivo solubilization kinetics.18 Several mechanisms have been proposed in the literature to describe dissolution. Noyes and Whitney5 expressed the dissolution rate by assuming that the process is diffusion controlled and involves no chemical reaction, considering that dissolution is directly proportional to the difference between the solubility and the solution concentration. Dressman and Fleisher6 developed an expression based on NoyesWhitney equation, where the dissolution rate was expressed as the function of the remaining surface area of the crystal and the concentration gradient across the boundary layer. Hintz and Johnson7 and Lu et al.8 used the same dissolution rate expression, assuming a mass transfer controlled process, to simulate the dissolution of drugs in polydisperse powder form. Lu et al.8 modified the Noyes Whitney equation by including a nonspherical geometry for the crystals and a variation of the diffusion layer thickness with particle size. Recently, Shan et al.9 took into account both surface mechanism and mass transfer to simulate the dissolution of drug crystals in suspension, in a stirred vessel. The authors concluded that the mechanism changed from a diffusion-controlled process at high undersaturation to a process controlled by surface Received: November 23, 2010 Revised: March 9, 2011 Published: April 01, 2011 2205
dx.doi.org/10.1021/cg101555u | Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design mechanism at lower undersaturation.10 Because of the complexity of capturing accurately the dissolution mechanism, often generic semiempiric apparent size-dependent dissolution kinetic laws can be used, with parameters identified based on experimental data for a particular process. Although this approach only provides an empirical kinetic expression, it can describe the observed change in solute concentration and crystal size distribution (CSD) when used in conjunction with the population balance equation (PBE). Dissolution mechanism is also important in other industrial sectors, outside pharmaceutical manufacture. For example, the dissolution of fine crystals (known as ripening) is effective for the production of large crystals with a narrow crystal size distribution.10,11 Dissolution can be used to dissolve small particles and hence can be used to suppress the effects of secondary nucleation, which can be a dominant phenomenon in industrial crystallization, for example, due to attrition. The CSD is a key factor in the production of high quality products and determines the efficiency of several downstream processes. Many problems in downstream processing can be attributed to poor particle characteristics established in the crystallization step.12 The main difficulty in batch crystallization is to accomplish a uniform and reproducible CSD. Online control during the process allows for improved crystalline product quality, shorter process times, and reduction or elimination of compromised batches. Combination of a dissolution mechanism along with growth gives the flexibility to successfully shape the CSD and achieve improved crystal properties. The generic approach applied in industrial crystallization processes is to use external fine removal equipment.13 This generally consists of a draft pipe, pump, and external heat exchanger through which the slurry is recycled. The maximum fines size to be removed through the vertical draft pipe is given by the upward velocity through the tube determined by Stoke’s law of particle settling. This upward velocity is controlled by the recycle flow rate so that particles above a certain size sediment in an internal draft pipe and only fine particles are eliminated from the system. The slurry with the fine particles is then recycled through a heat exchanger where fines are dissolved and the supersaturated solution is reintroduced in the crystallizer. This approach requires additional equipment yielding increased capital and operating costs. Additionally, in the case of batch crystallization processes during the time when the slurry is eliminated and the draft pipe no longer reaches into the solution nucleation can occur, compromising the quality of the final product CSD. The recent development in model free (or direct design) approaches for crystallization processes have demonstrated the advantages of controlled dissolution to obtain larger crystals with narrower CSD.1416 These approaches use real time measurements provided by process analytical technologies to adapt the operating policies during the crystallization process to achieve controlled nucleation and in situ dissolution events. These approaches are very robust, easily adaptable, and can eliminate the need for extra equipment required for external fine removal. Although the advantages of these techniques are now widely accepted and simplified fine removal approaches, such as temperature cycling are commonly used in industry, generally these operating profiles are chosen based on trial-and-error. The only case study that demonstrates the application of a batch nonlinear model predictive control approach for online controlled dissolution was proposed recently.17 However, the approach does not provide a comprehensive evaluation of the design of open-loop control approaches that take dissolution into account. This paper
ARTICLE
proposes a methodology for the systematic model-based design of the open-loop operating policies that incorporate controlled dissolution to achieve a desired shape of the target CSD. Population balance models have been widely used for modeling crystallization processes and prediction of the CSD. However, the solution of the generic PBE generally requires complex numerical solution techniques characterized by significant computational burden, making only a small portion of the available solution approaches suitable for model based optimization and control.1822 The variety of solution approaches proposed for the PBE include the standard method of moments (SMM), the quadrature method of moments (QMOM),2327 the method of characteristics (MOCH),28 the method of classes (MOC),29,30 direct numerical solution (DNS) techniques, such as finite volume and finite difference schemes31,32 and dynamic Monte Carlo (DMC) simulation approaches.24,33 Both SMM and QMOM provide efficient solutions of the PBE and have been widely used in the literature for optimization and control purposes.34,35 However, these approaches only provide the moments of CSD and not the entire distribution. The method of characteristics (MOCH) in combination with the SMM has been used successfully for processes with size-independent growth and nucleation,28 including the modeling and control of polymorphic transformations.3638 Recently, a novel application of the method of characteristics (MOCH) in combination with the QMOM has been used successfully for processes with size-dependent growth and nucleation.39,40 The method has been extended for size-dependent dissolution mechanism,41,42 and according to the authors’ knowledge this is the first time when a comprehensive optimal distribution shaping control approach is proposed based on the solution of the PBE with combined size-dependent growth and dissolution together with nucleation. The experimental applications of both the model-free and model-based approaches rely on the application of different process analytical technologies (PATs), which can provide in situ real-time information about the crystallization process.4247 The paper illustrates a novel application of the method of characteristics (MOCH) in combination with the QMOM method for the prediction of the evolution of the entire CSD for a generic size-dependent dissolution mechanism in addition to size-dependent growth and nucleation. The method has been applied to identify the kinetic parameters of dissolution for the batch crystallization of potash alum in water using laboratory experimental data and to simulate the fine removal mechanism during dissolution in this system. The paper demonstrates that using model-based optimization optimal temperature profiles can be developed, which can lead to controlled dissolution (optimal fine removal) and can produce product distributions which would not be achievable operating within the saturated region of the phase diagram only.48,49 The proposed in situ fine removal approach based on an optimal design of a fine removal loop in the operating curve in the phase diagram allows achieving the desired monomodal target distributions without the use of external fine removal equipment.
2. EFFICIENT SOLUTION OF THE PBE FOR BATCH CRYSTALLIZATION PROCESSES CONSIDERING SIZE-DEPENDENT DISSOLUTION, SIZE-DEPENDENT GROWTH, AND NUCLEATION MECHANISMS For a one-dimensional dissolution, growth, and nucleation mechanisms characterized by a characteristic length, L, and a 2206
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
well-mixed crystallizer with dissolution, growth, and nucleation phenomena, the PBEs can be written in the following forms: Dissolution : ∂fn ðL, tÞ ∂ðDð S, L; θd Þfn ðL, tÞÞ þ ¼ 0, S < 0 ∂t ∂L
ð1Þ
function of the overall crystal mass (volume) expressed by the third moment of the distribution μ3(t). The PBEs are solved together with the overall mass balance, which is also calculated using the third moment μ3(t). Hence the third moment can be used to obtain the undersaturation S(t) and/or the supersaturation S(t) from the solute concentration: CðtÞ ¼ Cð0Þ kv Fc ðμ3 ðtÞ μ3 ð0ÞÞ
Growth and nucleation : ∂fn ðL, tÞ ∂ðGðS, L; θg Þfn ðL, tÞÞ þ ¼ BðS; θb ÞδðL r0 Þ, S g 0 ∂t ∂L
ð2Þ
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, D(S, L; θd) is the rate of crystal dissolution, G(S, L; θg) is the rate of crystal growth, B(S; θb) is the nucleation rate, S = (Csat C) is the absolute undersaturation, 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, δ(L r0) isR the Dirac delta 6 r0, with ¥ function (δ = ¥ if L = r0 and δ = 0 if L ¼ ¥ δ(x) dx = 1) and θd, θg, and θb are vectors of dissolution, growth, and nucleation kinetic parameters, respectively. The solution of eqs 1 and 2 is an initial value problem, with initial condition given by the size distribution of seed, fn(L,0) = fn,0(L0). Because of the undersaturated conditions crystals reaching a limit size (ro) disappear. Hence for the dissolution PBE the left boundary condition was left undefined. The right boundary condition is given by fn(¥, t) = 0 where the “infinite size” represents any limit size that is larger than the size of any crystal.50 For modeling the combined growth, nucleation, and dissolution events the solution of the two PBEs (1) and (2) is performed in a mutually exclusive way, depending on whether the system is in supersaturated or undersaturated region. However, since the operating trajectory can pass between the two regions, the states are propagated from one system to the other, when the solution algorithm switches between the two PBEs. Note that the discontinuity introduced by the switching between mechanisms may create numerical problems during the integration. However, during the simulations performed in this study no numerical difficulties have been encountered, using the generic variable order solver based on the numerical differentiation formulas (NDFs) for stiff, ordinary differential equation (function ode15s) in Matlab. In the case of numerical problems smooth, continuous switching functions (e.g., logsig) can be used between the mechanisms, instead of the discrete selection based approach. The generic relationship used in this study for apparent sizedependent dissolution, apparent size-dependent growth, and secondary nucleation are expressed as follows: Dð S, L; θd Þ ¼ kd ð SÞd ð1 þ ζLÞq GðS, L; θg Þ ¼ kg Sg ð1 þ γLÞp BðS; θb Þ ¼ kb Sb μ3
where Fc is the density of crystals and kv is the volumetric shape factor. Note that in the case of dissolution generally μ3(t) e μ3(0) and eq 4 shows an increase of the solute concentration. The moments (including μ3 required for the nucleation term and mass balance) can be calculated using the standard method of moments (SMOM) in the case of size-independent growth expression, or by using the quadrature method of moments (QMOM) for size-dependent growth rate expressions. Both the SMOM and the QMOM calculate the moments of the distribution defined by Z ¥ fn ðLÞLj dL, j ¼ 0, 1, 2, :::, ¥ ð5Þ μj ¼ 0
The quadrature method of moments (QMOM) is a generic solution approach for the PBE2325,51 which employs a quadrature approximation of the distribution function fn ðL, tÞ
Nq X
wi ðtÞδðLi ðtÞ, LÞ
ð6Þ
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) algorithm52 or via direct solution of a differential-algebraic (DAE) system,53 based on the idea of minimizing the error committed by replacing the integral from the moment definition with its quadrature approximation, Z ¥ Nq X j fn ðLÞLj dL wi Li ð7Þ μj ¼ 0
i¼1
Applying the moment transformation to eqs 1 and 2 with the quadrature approximation of eq 7 the resulting moment equations have the following form, with the right-hand side depending on whether the system is in the undersaturated (dissolution) or supersaturated (growth and nucleation) region of the phase diagram: ( 0, if S < 0 dμ0 ¼ ð8Þ BðS; θb Þ, if S g 0 dt
ð3Þ dμj dt
The size-dependent growth mechanism here is called apparent, to emphasize the fact that this is an empirical interpretation of the experimental observations. These observations may be caused by the actual growth kinetics, by diffusion controlled growth,48 or simply by growth rate dispersion, or a combination of these. It is also possible that in different stages of the crystallization processes one or the other mechanism is more dominant.42,48 The nucleation in the case of the generic seeded system a generic secondary nucleation mechanism is considered, which is a
ð4Þ
¼
8 Nq X > > j1 > > j wi Li DðS, Li ; θd Þ, j ¼ 1, 2, 3::: if S < 0 > > < i¼1 Nq > > X > j1 j > > j wi Li GðS, Li ; θg Þ þ BðS; θb Þr0 , j ¼ 1, 2, 3::: if S g 0 > : i¼1
ð9Þ The first equation in (8) indicates that according to the approach the number of particles in the dissolution part is considered constant (μ0 = constant). However, all particles eventually will reach a size of zero or r0 (i.e., below the detection limits) and hence will disappear. This formulation neglects the kinetics of 2207
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
“disappearance” of particles (the opposite mechanism to nucleation), which would describe how quickly each characteristic line fn decays to zero when the corresponding Li reaches r0. In this approach, it is considered that the particles disappear instantaneously, by having a size below the detection limit r0. The QMOM (or MOM) is only able to give the moments of the distribution, which are needed for the mass balance and nucleation term, but are not sufficient for the distribution shaping control. To obtain the shape of the CSD, the PBEs are solved using the method characteristics (MOCH), which can be used to reduce the generic PBEs given by eqs 1 and 2 to a system of ODEs. The characteristic equations are given by the following system of ODEs: dL ¼ dt
(
Dð S, L; θd Þ, GðS, L; θg Þ,
if S < 0 if S g 0
8 dDð S, L; θd Þ > > , f ðL, tÞ > < n dL
dfn ðL, t Þ ¼ > dt dGðS, L; θg Þ > > : fn ðL, tÞ þ BðS; θb ÞδðL r0 Þ, dL
if S < 0 if S g 0
ð10Þ with initial conditions L = L0 and fn(L, 0) = fn,0(L0), that is, the seed CSD at initial time t0 = 0. To obtain the dynamic evolution of the crystal size distribution fn(L, t), eqs 10 with a prescribed dissolution, growth, and nucleation rate expression can be integrated repeatedly for different initial values [L0, fn,0(L0)]. The initial conditions start with values calculated by choosing a discretization interval ΔL0 and using t0 = 0, L0 = max(0,L0,max kΔL0) and k = 0,1, ..., N, where N is the number of discretization points for the seed distribution and L0,max is chosen to be larger or equal to the maximum size range of the seed crystals. In the case of the dissolution mechanism, only the solution of the PBE (1) is an initial value problem only and hence the discretization interval ΔL0 will determine the number of integrations and thus the resolution of the dynamic evolution of the CSD. 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. The disappearance of fines is assumed to occur for L e r0, as shown in Figure 1. As shown on the figure the characteristic lines are integrated continuously even after the particles “disappear”. This approach is needed for consistency with the QMOM formulation where a constant total number of particles is assumed during dissolution (see eq 8). In this case, although the total number of particles virtually does not change during dissolution, some of them become of size zero (or below the detectable size, r0) and hence will not appear in the reconstruction of the CSD. This formulation allows the incorporation of the disappearance of fines in the model, and a smooth switching between growth and dissolution, as some of the particles which become smaller than the detectable size (r0) during dissolution may grow back to a detectable size when the system becomes again supersaturated, without requiring a nucleation event to occur. This scenario corresponds more closely to the real situation. 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. In the case of nucleation and growth, the initial conditions are the same as for the dissolution problem. However, in this case the boundary condition along the t axis of the L t plane also needs
Figure 1. Evolution of characteristic lines obtained from the method of characteristics in the case of dissolution mechanism.
to be considered, which determines the appearance of new characteristic lines due to nucleation, as shown in Figure 2. The solution approach for this case is described in detail by Aamir et al.39,40 for size-dependent growth and secondary nucleation mechanisms. The initial conditions for the integrations along the t axis of the L t plane, start 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. The iterations are stopped when tnext g tf where tf is the end time of the batch. 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 δ* = 1 if L ∈ [0, r0] and δ* = 0 if L ˇ [0, r0], with a switching defined over an interval. Thus, nucleation events are assumed to occur for L e r0. In this study, the same r0 = 1 μm is used for the size of detectable nuclei as the limit at which particles disappear during dissolution. When nucleation, growth, and dissolution are considered together the integrations are initialized according to the common initial conditions for both cases, and the boundary condition is used depending on whether the state of the system is under- or supersaturated. The right-hand sides of the characteristic differential equations change according to whether the system is supersaturated or undersaturated; however, the states are integrated continuously. The kinetic expressions described by eqs 3 give the following characteristic equations, which are used in the model: 8 d q > < kd ð SÞ ð1 þ ζLÞ ,
dL ¼ > dt : kg Sg ð1 þ γLÞp ,
if S < 0 if S g 0
dfn ðL, t Þ 8dt d q1 > , < fn ðL, t Þkd ð SÞ ζð1 þ ζLÞ ¼ > : fn ðL, t Þkg Sg γð1 þ γLÞp 1 þ BðS; θb ÞδðL r0 Þ,
if S < 0 if S g 0
ð11Þ 2208
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
Figure 2. Evolution of characteristic lines obtained from the method of characteristics in the case of nucleation and growth mechanisms.
Note that the third moment of the CSD required for the mass balance and evaluating the secondary nucleation in principle could be calculated from the discretized distribution using the characteristic lines computed by eqs 10 or 11. However, this method would require a high discretization of the distribution and an analysis of the discretization method for an accurate closure of the mass balance. The proposed combined QMOM-MOCH method first solves the population balance eqs 1 and 2 by transforming these into systems of ODEs by applying the standard method of moments (SMOM), in the case of size-independent growth and dissolution and nucleation, or the quadrature method of moments (QMOM), in more generic cases including size-dependent growth and dissolution, breakage, and aggregation. The advantage of using the combined QMOM/SMOM-MOCH is that the calculation of the moments is decoupled from the computation of the discretized CSD; hence, highly accurate moments can be obtained very efficiently and the resolution of the computed CSD can be chosen to fulfill the optimization requirements and to tailor computational requirement. Figure 3 shows a schematic representation of the combined QMOM-MOCH approach when dissolution, growth, and nucleation mechanisms are used together. In this case, both PBEs, one for the dissolution (1) and one for growth and nucleation (2), are included in the model, and are applied to the CSD, fn, depending on whether the operating curve is in the supersaturated (S > 0) or undersaturated region (S < 0). After the
Figure 3. Flowchart of the combined QMOM-MOCH approach for the solution of PBE using size-dependent growth and dissolution and nucleation mechanisms for supersaturated and undersaturated regions.
initialization of the method with initial conditions L = L0 and fn (L,0) = fn,0 (L0), the moments are calculated using eqs 8 and 9. Using the quadrature method of moments (QMOM) the dynamic evolution of the third moment μ3(t) and the supersaturation S(t) are computed. First the QMOM is applied with the dissolution or the growth and nucleation mechanisms based on whether the S < 0 or S > 0, respectively, and the dynamic evolution of the concentration, supersaturation, and third moment is computed by a single simulation for the entire batch. The dynamic evolution of the supersaturation suggests whether the system is in the supersaturated or undersaturated region. In the next part of the algorithm, the systems of pairs ODEs resulting from the MOCH is solved iteratively (characteristic line and corresponding fn). If S > 0, the system is in the supersaturated region and the MOCH is applied to the PBE with growth and nucleation kinetics, and if 2209
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
Figure 4. A schematic representation of the experimental setup.
S < 0 the system is in the undersaturated region and the MOCH is applied to the PBE with dissolution mechanism until the size of particles reaches L < r0, as shown in Figure 1. To obtain the dynamic evolution of the CSD, fn(L,t), the system of eqs 11, with the corresponding nucleation, growth, and dissolution expressions (given by eqs 3) are integrated repeatedly for different initial values [L0, fn,0(L0)]. In this way, the approach is able to consider all three mechanisms, that is, growth, nucleation, and dissolution, based on the supersaturated or undersaturated state of the system at a particular time of the batch.
3. EXPERIMENTAL SETUP 3.1. Materials. Potassium aluminum sulfate dodecahydarte (KAl2 (SO4)2 3 12H2O) (>99.95% purity, Fisher Bio Reagents) compound was used in the dissolution experiments. Deionized water was used as solvent. The solution was prepared using 8.4 g of potash alum in 100 g of water, corresponding to a saturation temperature of 30 C.48,54 3.2. Apparatus. The temperature in a 0.5 L jacketed glass vessel was controlled with a Pt100 thermocouple using a Huber VPC CC3 450 thermostat. An overhead stirrer with a four-blade pitch type impeller was used to agitate the system at 350 rpm. This agitation speed was chosen to be high enough to guarantee that particles were well suspended throughout the process, but low enough to avoid attrition or entrainment of bubbles due to vortex formation. A FBRM probe (model D600, Lasentec) was inserted into the solution to measure chord length distributions. The concentration was measured in situ using a conductivity probe. Conductivity (voltage) was measured using a CM 35 meter with WPA-35 conductivity probe. The CSD at the end of experimental runs, with different time durations, was measured using laser-diffraction equipment (Mastersizer 2000 with a hydro 2000 SM dispersion unit). Hexane was used to disperse potash-alum crystals to measure the CSD. The samples were analyzed off-line. Instead of performing one experiment during which samples would be withdrawn at different time intervals to evaluate the dynamic evolution of the CSD, experiments were performed by repeating the same experimental conditions (same seed, initial concentration, and temperature profile) but stopping the experiments at different stages of the batch and analyzing the sample obtained at the end. This approach was chosen to avoid unmodeled disturbances in the dynamic evolution of the process generated by taking the samples during a single experiment. The reproducibility of the experiments was confirmed by comparing the data from the conductivity and focused beam reflectance measurement (FBRM) probes for the overlapping parts for the experiments with progressively increasing duration. Images of crystals were also taken
ARTICLE
using a Leica DM LM microscope equipped with a Leica PFC 350 FX camera. A schematic representation of the experimental setup is shown in Figure 4. The operating conditions for the experiment are summarized in Table 1. 3.3. Seed Preparation. Seeds were prepared using sieve analysis. The sieve sizes used were 500, 355, 300, 300, 250, 200, 150, and 125 μm. The run time for the sieving operation was set to 120 min, and the rotation and shaking caused the crystals to distribute throughout the sieve stack. The product obtained on the sieve size of 300 μm was collected for seeding and was stored in a desiccator. Quantity of each sieving batch was 50 g and the shaking intensity was set to 7. 3.4. Method. A solution of potash alum and water was prepared to obtain the kinetic parameters for dissolution. A 0.5 L jacketed crystallization vessel equipped with thermocouple, conductivity, and FBRM probes was used. The saturation temperature used for these experiments was 30 C (8.4 g of potash alum dissolved in 100 g of water). Potash alum was dissolved in water by heating to 40 C at a rate of 0.8 C/min. The solution was equilibrated at 40 C for 30 min and then the temperature of the solution was reduced to 29 C at a rate of 0.5 C/ min. The temperature of the solution was maintained at 29 C prior to the start of experiment. After 10 min, 19 g of sieved seed in the size range between 300 and 355 μm (CSD determined using Malvern Mastersizer) was added to the solution and the temperature was maintained at 30 C for 10 min. During this period, the conductivity and FBRM readings were monitored to check if any amount of the seed had dissolved. The conductivity meter showed constant readings before and after the addition of the seed, which verified that the seed had not dissolved and the solution was saturated at 30 C. The FBRM number of counts increased suddenly corresponding to the seed addition and remained constant after that, confirming again that no seed dissolution occurred. In the next step, the process temperature was increased linearly to 35 C over a duration of 80 min. The same initial procedure was used for a sequence of nine runs using batch times of t = 0, 10, 20, ..., 80 min. At the end of each run, the temperature profile was stopped and the product was removed for CSD measurement.
3.5. Concentration Measurement Using Conductivity Meter. Conductivity (voltage) was measured using a CM 35 voltmeter with WPA-35 conductivity probe. For calibration, the conductivity was measured for several concentrations over a range of temperatures, as shown in Figure 5, to cover the metastable zone as well as a range of undersaturated conditions, where dissolution occurs. The temperature was decreased by 1 C at each concentration until the system nucleated. Hence, the nucleation points, as shown in Figure 5, were determined experimentally. The solubility curve was obtained from the literature.48,5456 The sensitivity of the conductivity measurement to the solid content in the slurry was evaluated by measuring the conductivity at constant temperature in equilibrated slurry, and changing the agitation speed to vary the solid fraction passing through the electrodes of the conductivity probe. Because of the relatively large surface area of the electrodes the sensitivity was practically negligible; hence the conductivity probe was considered to be suitable for concentration measurement in the slurry. Figure 6a,b shows the dependence of the conductivity on temperature at various concentrations, and the variation of the conductivity with concentrations at different constant temperatures. Both dependences appear to be linear with relatively constant slopes over the tested temperature and concentration ranges; hence, a simple bivariate linear regression based calibration model was used: Ck ¼ a0 Vk þ a1 Tk þ a2
ð12Þ
where T is the temperature in C, V is the measured voltage (output from the conductivity meter), and Ck is the experimental concentration values at the discrete measurement points k = 1, ..., K. The parameters for the 2210
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
Table 1. Operating Conditions for Dissolution Experiment for Seeded-Batch Cooling Crystallization operating conditions
units
experiment values
saturation temperature
C
C30
seed mass, (mseed)
kg
1.9 102 80
batch time, (tbatch)
min
initial solute concentration (Ci)
g of solid/g of water
0.084
initial temperature at seeding and start of profile, (T0)
C
29
final temperature, (Tf)
C
35
temperature profile followed
linear, Tlinear = T0 (T0 Tf) (t/tbatch)
points for smooth profile, (N) sampling time for FBRM measurement and conductivity measurement
s
60 10
sieve sizes for seed, (L)
μm
300355
agitation speed
rpm
350
density of crystals46 (Fc)
kg/m3
1750
volumetric shape factor, (kv)
0.62
mass of slurry, (mslurry)
kg
number of samples for CSD measurement
0.52 9
Figure 5. Measurement points for conductivity for the used concentrations and temperature ranges including solubility curve46 and the detected nucleation points. calibration model were estimated by a least-squares optimization approach solved using the fmincon function in MATLAB. The optimization problem for the parameter estimation is given by min ai
K X k¼1
exp
ðCk Ck Þ2
ð13Þ
where Ck and Cexp k are the estimated and experimental concentration values at the discrete measurement points k = 1, ..., K. The fitted parameters with their 95% confidence intervals for the calibration model are given in Table 2. The measured and the estimated concentrations using the calibration model showed very good agreement (Figure 7). The calibration eq 12 was validated against literature solubility data for anhydrous potash alum. An experiment was conducted in which the temperature of slurry containing excess solid was increased in several steps as shown in Figure 8a. Figure 8b shows a comparison between the estimated concentrations for the solubility curve using the calibration model and the literature data.48,5456 Good agreement between the measured results and the literature data was observed, with a sum square error of only 0.002. Since the
Figure 6. Experimentally observed relationship between conductivity and (a) temperature and (b) concentration. concentration measurements are expressed in weight fraction of anhydrous potassium alum, the solubility curve is calculated for the anhydrous potash alum and is given by the equation, Csol, anh ¼ Csol, hyd
Mw, anh Mw, hyd
! ð14Þ
where Mw,anh = 258.21 and Mw,hyd = 474.39 are the molecular weights for the anhydrous and hydrous forms, respectively. 2211
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
Table 2. Calibration Model Parameters for Concentration Measurement Using a Conductivity Probe for Potash AlumWater System parameters
values
error bound at 95% confidence interval
a0
0.2994
(0.0106
a1
0.0013
(0.0028
a2
0.0091
(0.0017
Figure 7. Comparison of measured and estimated concentrations using the calibration parameters shown in Table 2.
4. SIMULATION AND MODEL IDENTIFICATION OF THE SIZE-DEPENDENT DISSOLUTION FOR POTASH ALUMWATER SYSTEM The QMOM-MOCH approach has been validated for the potash alum (KAl(SO4)2) and water system using laboratoryscale experiments. The dissolution parameters were determined to capture the dynamic evolution of the shape of the size distribution, as well as the experimental concentration profile. Size-dependent growth and secondary nucleation parameters for potash alumwater system have already been identified and presented by Aamir et al.39 The optimization problem for the parameter estimation for the dissolution kinetics is given by ( ) Nd K X K X X exp exp 2 2 min wf ðfv, k ðLl Þ fv, k ðLl ÞÞ þ wc ðCk Ck Þ θ k¼1 l¼1
k¼1
ð15Þ subject to θmin e θ e θmax
ð16Þ
with θ = [kd, d, ζ, q] being the model parameter vector with the dissolution kinetic parameters, θmin and θmax are vectors with specified minimum and maximum bounds for each parameter, are the simulated and experimental respectively, Ck and Cexp k concentration values at the discrete time steps k = 1, ..., K, fv,k and exp are the values of the simulated and experimental volume fv,k distribution functions, corresponding to the discretized size Ll, l = 1, ..., Nd, with Nd being the number of experimental size bins, and wf and wC are weighting factors. The weighting factors generally can be used as the inverse of the standard deviations of the corresponding experimental data. However, in this study these were chosen arbitrarily to achieve similar magnitudes for the two terms in the objective function (wf = 500 and wC = 1). The bounds on the kinetic parameters (θmin and θmax) are often
Figure 8. (a) Determination of the solubility curve to validate the calibration parameters by increasing the temperature from 15 to 45 C in 5 C steps while containing excess solid in the slurry throughout the process. (b) Comparison between the experimental solubility curve using conductivity and literature data.
not easy to choose. In this study, first several optimizations were performed using an unconstrained optimization algorithm based on a derivative free simplex search algorithm (implemented in the Matlab function fminsearch), from various initial guesses. This optimization method provides very slow convergence but has the advantage of minimal information needed. On the basis of preliminary results from the unconstrained optimization, next the constrained optimization problem was solved using a sequential-quadratic-programming (SQP) method (fmincon function in Matlab), with the bounds selected to be within a broad range around the values obtained from the unconstrained optimization and adjusted if the optimal values were on the bounds (see Table 3). P The volume distribution function is 3 d calculated as fv,l = fn,lL3l / N l = 1(fn,lLl ΔLl). The optimization problem is solved using a sequential quadratic programming (SQP) based optimization approach implemented in the Matlab function fmincon. The resulting model parameters for the potash alum system are summarized in Table 3. To evaluate the robustness of the identified model, the confidence intervals of the estimated parameters were also calculated by the method described in detail by Nagy et al.57 Repeated dissolution experiments were conducted under identical experimental conditions, but stopping the batches at 2212
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
Table 3. Kinetic Parameters for Dissolution Nucleation and Growth, for Potash AlumWater System parameters
value
error bound at 95% confidence interval
bounds used in the optimization
Dissolution dissolution rate constant (kd), μm s1
1.2804
(0.072
[0.005, 18]
dissolution constant (ζ), μm1
0.0202
(0.063
[0.0001, 5]
dissolution constant (q), --
0.8604
(0.081
[0.0001, 10]
dissolution order constant (d), --
0.9801
(0.062
[0.0001, 10]
nucleation rate constant (kb), #/kg slurry1 s1
0.0380
(0.044
[0.0001, 5]
nucleation order constant (b), --
3.4174
(0.037
[0.005, 5]
growth rate constant (kg), μm s1 growth constant (γ), μm1
8.5708 0.0050
(0.036 (0.0035
[0.005, 18] [0.0001, 10]
Nucleation
Growth
growth constant (p), --
1.5778
(0.079
[0.0001, 5]
growth order constant (g), --
1.0000
(0.095
[0.0001, 10]
Figure 9. Dynamic evolution of the CSD throughout the batch for simulated and experimental CSD during dissolution mechanism.
different times (10, 20, ..., 80 min) for off-line CSD measurement using the Malvern Mastersizer unit. This approach was chosen to avoid disturbances in the mass balance due to samples taken during a single batch run. Details of the operating conditions of the experiments were already shown in Table 1. The same seed was introduced in each experiment shortly after the process temperature was stabilized at 29 C within 10 min. Figure 9 illustrates the dynamic evolution of the modeled and the experimental CSDs which are overall in good agreement during the entire batch. The proposed dissolution model with the identified parameters and solution algorithm based on the combined QMOM-MOCH method is able to describe the main features of the CSD throughout the entire batch. In the first part of the process (between 16 and 32 min) the deviation from the model prediction and experimental measurement is larger, which may be caused by the differences in the dissolution mechanisms at different stages of the batch, and by potential agglomeration during the experimental CSD measurement. Both explanations are possible and are facilitated by the presence of fine particles during this part of the process.
Figure 10. Evolution of characteristic lines (a) and number distribution function (b) for the simulated results for dissolution with the identified kinetic parameters.
The dissolution of the fine particles and the decrease of the size of the particles together with the narrowing of distribution are well captured by the model. Figure 10 illustrates the evolution of the characteristic lines and the number distribution function predicted by the simulation for dissolution, using the combined QMOM-MOCH. The evolution of the characteristic lines shows the narrowing of the distribution function due to the size-dependent dissolution kinetics. According to the approach, the particle numbers remain constant after the corresponding characteristic sizes reach the disappearance limit (r0 = 1 μm in this study). Since this size bin (e1 μm) acts as a sink for particles, they disappear and are not counted for the distribution. The results corresponding to the QMOM part of the combined QMOM-MOCH are shown in Figure 11. A quadrature approximation of only Nq = 2 was used, and Figure 11 represents the weights, abscissas, and the first four moments, obtained from the QMOM using eqs 8 and 9 with the product difference approach and the quadrature approximation given by eq 7. On the basis of the approach the weights (w) remain constant; 2213
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
also indicates that the seed was bimodal. As the temperature was increased from 30 to 35 C, the smaller particles dissolved, while the size of the larger particles decreased, hence shifting toward smaller sizes and narrowing the associated distributions. Figure 13 shows microscopic images of samples collected during the dissolution experiment. These images provide visual evidence of the disappearance of small particles and the reduction in size of the larger crystals compared to the seed. A total of eight batches were conducted, leading to nine CSD measurements (with the seed CSD at t = 0, and at t = 10, 20, ..., 80 min).
Figure 11. The abscissas (L), weights (w), and first four moments obtained from the quadrature method of moments (QMOM) part of the QMOM-MOCH method for size-dependent dissolution mechanism, corresponding to the simulation results presented in Figures 8 and 9.
however, the difference between the abscissas (L) decreases corresponding to the narrowing of the distribution also observed in Figures 9 and 10. The results indicate that a low order quadrature approximation (Nq = 2) can provide accurate calculation of the moments. Hence, in the MOCH part of the approach arbitrary discretization of the CSD can be used, based on the trade off between computational requirement and desired resolution of the CSD but with no effect on the accuracy of the moments and on the mass balance. 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. The seed used in the experiments (and simulation) was bimodal indicated by the two peaks in the distribution at the initial time in Figure 9. As the temperature increased according to the temperature profile in Figure 12b, the fine particles were dissolved, while the size of the larger particles decreased with the concomitant narrowing of the distribution, confirming the sizedependent mechanism of the dissolution. Figure 12 shows the concentration and FBRM results for the longest batch (80 min). The measurements from the shorter runs were consistent following the trend of the longest batch. The FBRM probe used in the experiments also detected a decrease in the square weighted mean chord length (SWMCL), and total number of counts per second, as shown in Figure 12a, caused by the disappearance of the smaller particles and the reduction in the size of the larger particles. The change in the concentration throughout the batch is shown in Figure 12b indicating a continuous increase with the dissolution of particles. The chord length distribution, as shown in Figure 12c,
5. MODEL-BASED DYNAMIC OPTIMIZATION OF TEMPERATURE TRAJECTORIES FOR THE DESIGN OF MONOMODAL CSD The combined QMOM-MOCH approach described in Section 2 was used to solve the population balance model (PBM) in a model-based dynamic optimization scheme. The aim was to determine the optimal temperature trajectories, which yield desired monomodal target CSD at the end of the batch. The final CSD is dependent on the supersaturation profile created over the batch time, and hence the heating/cooling trajectory is of critical importance. During optimization, both the temperature trajectory and the batch time were optimized. The batch time horizon [0, tf] was divided into Nb equally spaced time intervals of duration Δt (stages), with discrete time steps tk = jΔt, j = 0,1, ..., Nb for the solution of the dynamic optimization problem. The temperature trajectory is approximated by a piecewise linear function determined by the fixed initial temperature at t = 0, T(0) and the slopes RT(j) in each discretized period Δt. Since the batch time is also optimized, the duration of the time interval Δt is changing during the optimization, but the number of discretizations Nb is fixed. This formulation allows easy incorporation of the temperature rate constraints as bounds on the decision variables RT(j), which are important to obtain a practically implementable temperature trajectory. The optimization problem is formulated as follows: min RT ðjÞ, tf
Nd K X X target ðfv, k ðLl Þfv, k ðLl ÞÞ2
ð17Þ
k¼1 l¼1
subject to RT, min e RT ðjÞ e RT, min , j ¼ 0, 1, :::, Nb
ð18Þ
0 e tf e tf , max
ð19Þ
Cðtf Þ e Cf , max
ð20Þ
where RT(j) are the elements of the vector containing the slopes (dT/dt) for the temperature trajectories depending on the implementable heating and cooling capacity of the system, tf is the total batch time, C(tf) is the solute concentration at the end of the batch, Cf,max is the maximum acceptable concentration at the are the end of the batch to achieve the required yield, fv,k and ftarget v,k values of the 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 2214
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
Figure 12. (a) Measured total number of counts (#/s) and square weighted mean chord length detected by FBRM throughout the experiment; (b) temperature profile and measured concentration throughout the batch; (c) chord length distribution throughout the batch for the dissolution experiment conducted to determine the kinetic parameters.
Figure 13. Microscopic images for crystal size distribution (a) seed at t = 0 min, (b) t = 16 min, (c) t = 32 min, (d) t = 48 min, (e) t = 64 min, and (f) t = 80 min at the end of the experiment.
quadratic programming (SQP) approach implemented using the Matlab function fmincon. The kinetic parameters used are given in Table 3. Case 1: Designing Monomodal Distribution while Operating within Supersaturated Region. The temperature trajectory was optimized to achieve a fictitious monomodal distribution,
with the aim to try to eliminate completely the effect of secondary nucleation. The target monomodal distribution can be expressed as 2 1 2 target fn, monomodal ¼ pffiffiffiffiffiffi eðL 378Þ =ð2:50 Þ 2π2:50
2215
ð21Þ
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
Figure 14. (a) Comparison of simulated and target monomodal CSD (b) optimal temperature profile with 30 discretization points.
The monomodal seed was used, that is, a Gaussian distribution with a mean of 54 μm and a standard deviation of 15 μm. Figure 14a shows a comparison of the target and the simulated CSD at the end of the batch when the optimal temperature profile was implemented, shown in Figure 14b. The simulated product CSD is a bimodal distribution despite the monomodal target distribution. For these simulations, the number of discretizations was Nb = 30 and Cf,max = 0.4C(0) = 0.052 (kg/kg slurry). The bounds on the cooling rate were 0.5 C/min e R e 0, and the maximum batch time tf,max = 100 min. The value of Δt was calculated as Δt = (tf t0)/Nb. The optimal CSD captures the large peak very well; however, the optimal operating trajectory, which is constrained in this case to be within the metastable zone (R e 0) is unable to eliminate the development of the secondary peak due to secondary nucleation. The temperature trajectory was designed using growth and secondary nucleation. Both mechanisms occur in the supersaturated solution within the metastable zone. The secondary nucleation produces fines even when the supersaturation is not very large. This is in correlation with findings by other researchers58 who showed that in the case of particular seed loading in many cases it is practically impossible to avoid the formation of bimodal CSD, due to secondary nucleation. This problem generally is overcome by designing an additional fine removal loop, in which the slurry is recycled through an external heat exchanger in which the fines are dissolved. The size of fines eliminated from the system is controlled by the diameter of the draft tube and/or the flow
ARTICLE
rate. Without the external dissolution equipment, the only way to remove the fine crystals is to cross the solubility curve (Csat) into the undersaturated region, where eventually the fine crystals should preferentially dissolve leading to a monomodal CSD of large-sized crystals. Thus, to obtain the desired monomodal target distribution, an optimal temperature profile can be designed, while using both supersaturated and undersaturated regions. Case 2: Designing Monomodal Target Distribution via in Situ Fine Removal Using Controlled Dissolution. In this case, the temperature profile was optimized for the same monomodal target distribution described in case 1. The undersaturated region and supersaturated region both were used to design the optimal temperature profile, by allowing an increase in the temperature (0.5 C/min e R e 0.5) and using the model that incorporates growth, nucleation, and dissolution mechanisms. The model was solved using the combined QMOM-MOCH technique described in Section 2. The kinetic parameters used are presented in Table 3. Figure 15 shows the main results of the optimization considering growth, nucleation, and dissolution mechanisms. The target distribution is very well achieved with no secondary peak (see Figure 15a), which was impossible to eliminate while operating within the metastable zone (MSZ) only. It can be observed in Figure 15b that the temperature trajectory can be divided into three phases. During the first phase of 40 min of the batch, the temperature decreased from 40 to 28 C. In this period, seed crystals grow but the second peak of the distribution also develops due to secondary nucleation. In the second phase, the temperature increased from 28 to 32 C from 40 to 74 min, during which the fines are removed, but the larger particles also decrease in size due to dissolution. In the last phase, the temperature decreased again from 32 to 17 C for the remaining 26 min of the batch during which the remaining seed crystals grow and achieve the desired target CSD. The resulting batch time in this case was 100 min, which is longer than in case 1. The extra batch time is needed to give extra growth time to crystals to compensate for the decrease in size during the dissolution phase. Figure 15c shows the variation of the growth/dissolution rate during the three phases of the temperature profile. The supersaturation profile corresponding to the optimal temperature trajectory is shown in Figure 15e. The phase diagram showing the solubility curve along with the optimal temperature trajectory is presented in Figure 15f, and the complete dynamic evolution of the CSD is shown in Figure 16. In Figure 15f, a dissolution loop can be clearly seen, which indicates the cooling down, heating up, and cooling of the system again according to temperature profile in Figure 15b. When the temperature is decreased from 40 to 28 C the solution was supersaturated. During this time, the supersaturation reached its maximum value at around 22 min and then decreased to zero within the next 18 min. During this process, the seed has grown, and the apparition of a smaller peak is also observed due to secondary nucleation, as shown in Figure 16 at t = 32 and t = 39 min. In the next phase, the temperature was increased from 28 to 32 C, and the system entered the undersaturated region as shown in Figure 15e. When the temperature is increased, the smaller particles dissolve and the size of the larger particles decrease. The dissolution of particles during this heating phase leads to the formation of the dissolution loop, shown in Figure 15f. It can be seen in Figure 16 that between t = 59 and t = 74 min, the small peak in the distribution has disappeared and the distribution has 2216
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
Figure 15. (a) Comparison of simulated and target monomodal distribution, (b) optimal temperature profile, (c) growth and dissolution rates, (d) nucleation rate profile within the supersaturated and undersaturated regions for the monomodal target distribution, considering dissolution along with growth and nucleation, (e) supersaturation/undersaturation profile for monomodal target distribution, (f) phase diagram showing solubility curve and optimal operating profile with the dissolution loop for in situ fine removal to achieve the monomodal target distribution.
become narrower. During this period, the supersaturation started increasing and the system re-entered the supersaturated region from the undersaturated region, as shown in Figure 15e. During the last 26 min the temperature decreased from 32 to 17 C within the supersaturated region. The crystals have grown in size due to the available supersaturation and due to the size-dependent growth the CSD broadened, as shown in Figure 16 at t = 90 and t = 100 min. Some secondary nucleation was also observed in simulation during the last 26 min. However, secondary nucleation was not significant and the total volume of the generated particles is small; hence, no fines peak is shown in the final volume CSD in Figure 15a. Figure 17a shows the characteristic lines (L), which correspond to the combined phenomena of growth, nucleation, and dissolution. The characteristic lines correspond to the size-dependent growth in the first 40 min, which is indicated by the broadening of these lines. Following
this, the lines start narrowing down as the dissolution takes place. In the last phase, the lines broaden again to achieve the final shape of the distribution. The secondary nucleation can also be observed in Figure 17b, between 74 to 90 min, which is also indicated in Figure 15d. The described case study shows that by operating both within the metastable zone and in the undersaturated region, it is possible to obtain monomodal CSD. Using both phases to design an optimal profile provides the possibility to dissolve the fine particles produced by secondary nucleation. The simulation results show that in certain cases a monomodal target CSD is difficult to achieve, while just operating within the metastable zone and considering only growth and nucleation mechanisms. Operating above and below the metastable zone gives more flexibility to design an operating profile to achieve a target CSD using the kinetics of growth/dissolution and nucleation 2217
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design
ARTICLE
Figure 16. Dynamic evolution of CSD using size-dependent growth, secondary nucleation, and size-dependent dissolution kinetics for the potash alumwater system.
Figure 17. Evolution of characteristic lines (a) and number distribution function (b) for the simulated monomodal target distribution.
identified for the system. The above case study illustrates that the proposed method can be used to obtain optimal temperature trajectories by considering growth, dissolution, and nucleation mechanisms. 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 by including an in situ fine removal via controlled dissolution loops. Tailoring the final product CSD so that an optimal dissolution rate (heat release rate) is achieved for the case of exothermic reactive processes for which the dissolution is followed by reactions, can improve the safe operation of these systems. Similarly, the same strategy can also be applied for reactions that produce a significant amount of gases and induce excessive liquid rise, which can lead to the explosion of reactors (two-phase flow occurs in the connection pipe and pressure builds up).59,60
6. CONCLUSIONS The paper describes an approach for the population balance modeling of system with a generic size-dependent dissolution, size-dependent growth, and nucleation mechanism. The PBE is solved using a combined quadrature method of moments and method of characteristics approach, which provides a computationally efficient framework for reconstructing the dynamic variation of the whole CSD during the dissolution, growth, and nucleation processes. The approach is evaluated in the case of a seeded system of potash alum in water, for which the sizedependent dissolution kinetics was identified using laboratory experimental data consisting of off-line CSD measurement and in situ concentration measurement based on a conductivity probe. The process was also monitored using an FBRM probe. The technique was used to design an optimal temperature profile to achieve a monomodal target distribution using all three phenomena. Using both supersaturated and undersaturated regions gives an additional degree of freedom to design optimal temperature trajectories to dissolve small size particles and achieve monomodal target CSD, by designing controlled in situ dissolution loops for fines removal. The proposed approach can also be used for improved distribution shaping control, which can have significant safety benefits in the control of solidliquid reactions. ’ AUTHOR INFORMATION Corresponding Author
*E-mail: z.k.nagy@lboro.ac.uk; phone: þ44 1509 222516; fax: þ44 1509 223923.
’ 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 2218
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219
Crystal Growth & Design for providing experimental data for model identification of growth and nucleation kinetics.
’ REFERENCES (1) Sahin, O.; Ozdemie, M.; Kendirci, H.; Bulutcu, A. N. J. Cryst. Growth 2000, 219, 75–82. (2) Dokoumetzidis, A; Macheras, P. Int. J. Pharm. 2006, 321, 1–11. (3) Shan, G.; Igarashi, K.; Ooshima, H. Chem. Eng. J. 2002, 88, 53–58. (4) Myerson, A. S. Handbook of Industrial Crystallisation; ButterworthsHeinemann: Boston, USA, 1993. (5) Noyes, A. A.; Whitney, W. R. AIChE J. 1987, 19, 930–934. (6) Dressman, J. B.; Fleidher, D. J. Pharm. Sci. 1986, 75, 109–116. (7) Hintz, R. J.; Johnson, K. C. Int. J. Pharm. Sci. 1989, 5, 19–17. (8) Lu, A. T. K.; Frisella, M. E.; Johnson, K. C. Pharm. Res. 1993, 10, 1308–1314. (9) Shan, G.; Igarashi, K.; Ooshima, H. Chem. Eng. J. 2002, 88, 53–58. (10) Garcia, E.; Hoff, C.; Veesler, S. J. Cryst. Growth 2002, 237239, 2233–2239. (11) Mangin, D.; Garcia, E.; Hoff, C.; Klein, J. P.; Veesler, S. J. Cryst. Growth 2006, 286, 121–125. (12) Ward, J. D.; Mellichamp, D. A.; Doherty, M. F. AIChE J. 2006, 52 (6), 2046–2054. (13) Zipp, G. L.; Randolph, A. D. Ind. Eng. Chem. Res. 1989, 28 (9), 1446–1448. (14) Abu Bakar, M. R.; Nagy, Z. K.; Saleemi, A. N.; Rielly, C. D. Cryst. Growth Des. 2009, 9 (3), 1378–1384. (15) Woo, X. Y.; Nagy, Z. K.; Tan, R. B. H.; Braatz, R. D. Cryst. Growth Des. 2009, 9 (1), 182–190. (16) Abu Bakar, M. R.; Nagy, Z. K.; Rielly, C. D. Org. Process Res. Dev. 2009, 13, 1343-1356, 2009. (17) Nagy, Z. K. Comput. Chem. Eng. 2009, 33, 1685–1691. (18) Ramkrishna, D. Population Balances, Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, USA, 2000. (19) Gerstlauer, A.; Gahn, C.; Zhou, H.; Rauls, M. Chem. Eng. Sci. 2006, 61, 205–217. (20) Majumder, A.; Kariwala, V.; Ansumali, S.; Rajendran, A. Chem. Eng. Sci. 2010, 65 (13), 3928–3936. (21) Majumder, A.; Kariwala, V.; Ansumali, S.; Rajendran, A. Ind. Eng. Chem. Res. 2010, 49, 3862–3872. (22) Grosso, M.; Galan, O.; Barratti, R.; Romagnoli, J. A. Chem. Eng. Sci. 2009, 27, 291–296. (23) McGraw, R. Aerosol Sci. Technol. 1997, 27, 255–265. (24) Rosner, D. E.; McGraw, R. L.; Tandon, P. Ind. Eng. Chem. Res. 2003, 42, 2699–2711. (25) Marchisio, D. L.; Ptkturna, J. T.; Vigil, R. D.; Fox, R. O. AIChE J. 2003, 49, 1266–1276. (26) Qamar, S.; Noor, S.; Ul Ain, Q.; Seidel Morgenstern, A. Ind. Eng. Chem. Res. 2010, 49, 11633–11644. (27) Qamar, S.; Mukhtar, S.; Seidel-Morgenstern, A. J. Cryst. Growth 2010, 312, 2936–2945. (28) Hounslow, M. J.; Reynolds, G. K. AIChE J. 2006, 52, 2507– 2517. (29) Marchal, P.; David, R.; Klein, J. P.; Villermaux, J. Chem. Eng. Sci. 1988, 57, 1107–1119. (30) Puel, F.; Fevotte, G.; Klein, J. P. Chem. Eng. Sci. 2003, 58, 3715–3727. (31) Gunawan, R.; Fusman, I.; Braatz, R. D. AIChE J. 2008, 54, 1449–1458. (32) Zhang, T.; Tade, M. O.; Tian, Y. C.; Zang, H. Comput. Chem. Eng. 2008, 32, 2403–2408. (33) Haseltine, E. L.; Patience, D. B.; Rawlings, J. B. Chem. Eng. Sci. 2005, 60, 2627–2641. (34) Fujiwara, M.; Nagy, Z. K.; Chew, J. W.; Braatz, R. D. J. Process Control 2005, 15, 493–504.
ARTICLE
(35) Nagy, Z. K.; Braatz, R. D. AIChE J. 2003, 49, 1776–1786. (36) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallization; Academic Press: New York, USA, 1971. (37) Hermanto, M. W.; Chiu, M.-S.; Woo, X. Y.; Braatz, R. D. AIChE J. 2007, 53, 2643–2650. (38) Ono, T.; Kramer, H. J. M.; Ter Horst, J. H.; Jansens, P. J. Cryst. Growth Des. 2004, 4, 1161–1167. (39) Aamir, E.; Nagy, Z. K.; Rielly, C. D.; Kleiner, T.; Judat, B. Ind. Eng. Chem. Res. 2009, 48, 8575–8584. (40) Aamir, E.; Nagy, Z. K.; Rielly Chem. Eng. Sci. 2010, 65, 3602– 3614. (41) Aamir, E.; Nagy, Z. K.; Rielly, C. D. Proceedings of the 8th International Workshop of Industrial Crystallisation, 911th September, 2009, Lappeenranta, Finland; University of Lappeenranta: Finland, 2009; pp 6168. (42) Aamir, E.; Nagy, Z. K.; Rielly, C. D. Cryst. Growth Des. 2010, 10, 4728–4740. (43) Sarkar, D.; Doan, X.-T.; Ying, Z.; Srinivasan, R. Chem. Eng. Sci. 2009, 64, 9–19. (44) Yu, L. X.; Lionbergera, R. A.; Rawa, A. S.; D’Costa, R.; Wub, H; Hussain, A. S. Adv. Drug Delivery Rev. 2004, 56, 349–369. (45) Abu Bakar, M. R.; Nagy, Z. K; Rielly, C. D. Cryst. Growth Des. 2010, 10 (9), 3892–3900. (46) Simon, L. L.; Nagy, Z. K.; Hungerbuehler, K. Chem. Eng. Sci. 2009, 64, 3344–3351. (47) Kempkes, M.; Eggers, J.; Mazzotti, M. Chem. Eng. Sci. 2008, 63, 4656–4675. (48) Mullin, J. W. Crystallisation; Butterworth-Heinemann: London, UK, 2001. (49) Tadayyon, A; Rohani, S. Can. J. Chem. Eng. 2000, 78, 663–673. (50) Fevotte, G.; Alexander, C.; Nida, S. O AIChE J. 2007, 53, 2578–2589. (51) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. J. Colloid Interface Sci. 2003, 258, 322–334. (52) Gordon, R G. J. Math. Phys. 1968, 9, 655–663. (53) Gimbun, J.; Nagy, Z. K.; Rielly, C. D. Ind. Eng. Chem. Res. 2009, 48, 7798–7812. (54) Barrett, P.; Glennon, B. Trans IChemE, Part A. 2002, 80, 799– 805. (55) Xie, W.; Rohani, S.; Phoenix, A. Chem. Eng. Commun. 2001, 187, 229–249. (56) Zhang, G. P; Rohani, S. Chem. Eng. Sci. 2003, 58, 1887–1896. (57) Nagy, Z. K.; Fujiwara, M.; Woo, X. Y.; Braatz, R. D. Ind. Eng. Chem. Res. 2008, 47, 1245–1252. (58) Doki, N.; Seki, H.; Takano, K.; Asatani, H.; Yokota, M.; Kubota, N. Cryst. Growth Des. 2004, 4, 949–953. (59) Bhattacharya, A. Chem. Eng. J. 2008, 137, 347–360. (60) Simon, L. L.; Intrivigne, A.; Fischer, U.; Hungerbuhler, K. Chem. Eng. Sci. 2008, 63, 770–781.
2219
dx.doi.org/10.1021/cg101555u |Cryst. Growth Des. 2011, 11, 2205–2219