Modeling of Nucleation and Growth Kinetics for Unseeded Batch

Subscriber access provided by University of Newcastle, Australia

Article

Modeling of nucleation and growth kinetics for unseeded batch cooling crystallization Huayu Li, Yoshiaki Kawajiri, Martha A. Grover, and Ronald W. Rousseau Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b04914 • Publication Date (Web): 15 Mar 2017 Downloaded from http://pubs.acs.org on March 27, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Modeling of nucleation and growth kinetics for unseeded batch cooling crystallization Huayu Li, Yoshiaki Kawajiri,∗ Martha A. Grover, and Ronald W. Rousseau School of Chemical and Biomolecular, Georgia Institute of Technology, 311 Ferst Drive, Atlanta, GA 30332, USA E-mail: [email protected], Phone +1-404-894-2856, Fax +1-404-894-2866 Abstract Both primary and secondary nucleation rates are considered in a model developed for unseeded batch crystallization. A carefully designed strategy was employed to minimize the effects of the stochastic nature of induction time; nucleation was induced at designed supersaturations on known temperature plateaus. Crystallization kinetics of paracetamol from ethanolic solutions were extracted from measurements of in-situ solute concentrations and combined with sieve (ex-situ) data on the final product. Parameters in models for primary and secondary nucleation and for crystal growth rate were estimated by fitting a full population balance model to the measurements, and the evolution of the crystal size distribution was compared against in-situ estimation from focused-beam reflectance measurements using the technique that we previously developed. The resulting models suggest that primary nucleation produces fewer surviving crystals that had been 1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 49

expected, and that most of the product crystals from the process involving a temperature plateau result from secondary nucleation.

1.

Introduction

Numerous chemical species are separated and purified by crystallization from solution, and crystal quality, as defined by size, morphology, and purity, often distinguishes a successful process from one that is unsuccessful. The crystallization process can be operated either continuously or batch-wise, but high value-added chemicals are often manufactured in batch processes, which are simple, flexible, and less expensive to develop and integrate into a complete process 1. Even with current efforts to convert pharmaceutical manufacturing to continuous operations, the use of batch operations remains commonplace 2, 3. In batch crystallization processes, several methods can be used to generate a supersaturated solution, including cooling, evaporation of solvent, addition of anti-solvent, and change of pH. Nuclei are created in a supersaturated solution by either primary or secondary nucleation. Primary nucleation produces nuclei directly from dissolved solute molecules, which is often described by the classical nucleation theory (CNT) 1 or most recently by the two-step theory 4, 5. Primary nucleation usually requires high supersaturation to overcome the resistance to formation of a new phase. On the other hand, secondary nucleation can occur at low supersaturation; it is postulated that the resulting nuclei are detached from existing crystals through mechanisms such as fluid shear or crystal collision with other crystals or crystallizer internals 1.

2


Page 3 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Estimation of kinetic parameters for primary nucleation is usually carried out using data on either metastable zone width (MSZW) or induction time. Using linearized power-law expressions or CNT equations, the average MSZW or induction time is related to supersaturation and temperature 6. Another property of primary nucleation is that its onset is stochastic; i.e. it has a certain level of randomness. By investigating the distribution of induction times, it is found that the onset of primary nucleation is consistent with a Poisson process 7-9. The mean value and the standard deviation of the induction time depend on the nucleation rate and the volume of solution 8. To overcome the uncertainty associated with the stochastic nature of primary nucleation, the addition of seed crystals to supersaturated solutions is often used to facilitate controlled secondary nucleation. The product in such operations is expected to be of the same crystal form as the seeds and the reproducibility of the process is often improved 10, 11. The kinetics of seeded processes are often studied using a population balance equation together with a measurement of the processes. One common measurement is the solution concentration, obtained from in-situ spectra detection 1216

, or from another property, such as the density of the solution 17. However, concentration alone is

insufficient to decouple nucleation and growth 18, and thus additional information about the crystal phase must be provided. The traditional approach is offline measurement of the crystal size distribution (e.g. using sieving, image analysis, or laser diffraction 12, 17), which is often applied only to the final product. In-situ measurement of CSD is a challenging task but substantial progress has been reported recently. Online image analysis 19-22 can give both shape and size information and requires

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 49

dedicated and powerful algorithms to process the images. This method was used in several studies for parameter estimation 18, 23. Another technique is focused beam reflectance measurement (FBRM), which measures the chord-length distribution (CLD) of the crystal population. One use of the CLD is correlating the moments of the CLD with the moments of the CSD 16, 24, although the precise relationship between these two sets of distributions remains unclear. The relationship can be calibrated with some known CSD, such as the CSD of seed crystals and mass of crystals in some seeded processes 14. First-principles models that directly transform the CLD to CSD were also developed 25, 26, but their use in estimating kinetic parameters has been limited 15, 27. Fitting a population balance to process measurements is often applied to seeded runs for the estimation of secondary nucleation and growth rates 12, 14, 26-28. However, unseeded experiments are rarely employed with this approach 18. In unseeded runs, the nuclei that originate from primary nucleation may act as seed crystals for secondary nucleation 29, 30. To understand the crystallization kinetics, both primary and secondary nucleation mechanisms should be considered in modeling crystallization processes. The present work aims to estimate kinetic parameters of primary nucleation, secondary nucleation, and crystal growth by fitting the full population balance equation (PBE) to experimental data obtained from unseeded cooling crystallizations. The experimental procedure, which is described later, attempted to induce primary nucleation at a fixed constant temperature and, therefore at a fixed supersaturation. FBRM and attenuated-total-reflectance Fourier transform infrared (ATR-FTIR) measurements are collected during cooling and translated to CSD and

4


Page 5 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


supersaturation 31, 32. CSDs of the final products are also analyzed by a sieving test. The PBE is solved by the conservation element/solution element (CE/SE) scheme 33-35, and the fitting error to the experimental data is minimized to estimate the kinetic parameters. With the estimated parameters, the nucleation and growth rates are calculated to interpret the crystallization process in the unseeded batch process.

2.

Methods and Materials

2.1 Population balance equation and crystallization kinetics The batch cooling crystallization process can be described by a population balance equation in which the size of crystals is represented in one dimension and invariant crystal growth (no growth-rate dispersion and size-independent growth) is assumed: ∂n ∂n +G =0 ∂t ∂x

(1)

where n is the number density [no./(µm · kg of solvent)], G is the growth rate [µm/min], x is the one-dimensional size of a crystal [µm], and t is the time [min]. Boundary conditions, initial conditions and mass balance are n (t , x = 0) = B / G

(2)

n(t = 0, x) = n0

(3)

c(t ) = c0 − kv ρµ3 (t )

(4) 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 49

where B is the nucleation rate [no./min/kg of solvent], n0 is the initial CSD, c0 is the initial solute concentration in the solution, c(t) is the solute concentration [g solute/kg solvent], kv is the volume shape factor of the crystals, ρ is the density of crystals, and µ3 is the third moment of the crystal population density function. For unseeded experiments, the initial CSD n0 is zero for all size x. In this work, three kinetic phenomena are considered: primary nucleation, secondary nucleation, and crystal growth. The classical nucleation theory (CNT) was used to describe primary nucleation, where the difference between volume free energy and surface free energy is assumed to determine the nucleation rate 1. The nucleation model can be described by  16πν 2  σ3 B1 = k b1 exp  − 3 3 2  3k T (ln S )  

(5)

in which kb1 is a constant [no./(min·kg of solvent)], ν is the volume of one solute molecule [m3], k is the Boltzmann constant [m2 kg/(s2 K)], σ is the crystal-solution interfacial energy [J/m2], T is the temperature [K], and S is a dimensionless supersaturation ratio, c/cs. In this study, kb1 and σ were the parameters to be estimated, T and S were measured during experiments, and ν is approximated by molecular weight and density of the solute crystals. Observations of dense tiny liquid droplets or spherical particles, formed prior to occurrence of crystalline nuclei, rich in solute molecules, have led to postulates involving a two-step mechanism 35-37. According to this theory, primary nucleation is divided into two steps: (1) the formation of the dense liquid phase and (2) the transformation of the dense liquid phase to crystals.

6


Page 7 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Two-step nucleation has been observed in both inorganic 5 and organic systems 38, but has not been reported for the system in the present study. Since little is known about the transient liquid phase, we used CNT to model primary nucleation in this study. It also may be questioned if the primary nucleation is homogeneous or heterogeneous. It is very difficult to guarantee that foreign particles are completely excluded from the crystallizer, even if care is taken to provide a clean solution each time the experiment is run. In this study, the mechanism was assumed to remain the same across all the experiments since materials and experimental procedures were consistent. Under this assumption, the estimated surface tension in Equation (5) is an effective value representing the combination of homogeneous and heterogeneous nucleation. An empirical expression was used to model secondary nucleation: α

B2 = kb 2 ( S − 1) msβ

(6)

where kb2 is the nucleation constant [no./(min·kg of solvent)], and ms is the mass of crystals [g crystal/kg solvent]. In this model, the secondary nucleation rate is empirically related with mass of crystals, supersaturation, and collisions of crystals. The effect of crystal collisions is incorporated in the parameter kb2, which was assumed constant in our experiments since stirring speed (mixing intensity) was fixed. Growth rate is related to temperature and supersaturation by the expression

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

 E  G = kg exp  − a  ∆cγ  RT 

Page 8 of 49

(7)

where kg is the growth constant [µm/min], Ea is the activation energy [J/mol], and ∆c = c −cs [g solute/kg solvent] is the supersaturation that drives crystal growth. While S − 1 = (c − cs )/cs could alternatively be used as the driving force, ∆c was chosen herein to facilitate direct comparison of the present results with previous studies 15, 39.

2.2 Numerical method Solving the PBE with the kinetics models of Equations (5)–(7) requires an efficient numerical scheme, especially when multiple parameters are to be estimated by fitting the model to experimental data. The parameters influence the difficulty of the simulation, and hence, determine CPU time and accuracy. For example, when nucleation and growth occur slowly, the simulation can be accomplished with rough discretization and CPU time is short. On the other hand, it is possible that another parameter set causes fast nucleation and growth, which makes the PBE very difficult to solve and CPU time becomes excessive. Improper discretization could make the numerical scheme unstable, and thus a fine discretization of the spatial domain is required. Fine discretization generally increases the computational burden. An ideal numerical method should maintain an acceptable accuracy for different parameter sets. In this work, we used the conservation element/solution element (CE/SE) scheme to solve the PBE. Originally designed for aerodynamic problems, this scheme adopts a staggered way of discretizing spatial domains, and the conservation of mass is enforced locally and globally 35. The 8


Page 9 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


method has been applied to solve partial differential equations in the chemical engineering field, such as simulated moving bed 34, 40 and crystallization 41. In particular, Qamar et al. 42 compared CE/SE with the finite-volume method (FVM) and the finite-element method (FEM) for onedimensional population balance modeling, and concluded that CE/SE has “much better performance as compared to the other schemes.” The scheme was coded in MATLAB, in which the size range in our simulation was 0–2000 µm and the length of spatial finite element was set to 5 µm.

2.3 Material and apparatus Batch cooling crystallization of paracetamol (SigmaAldrich, 99%) dissolved in ethanol (KOPTEC, 200 proof anhydrous) was prepared in a crystallizer shown schematically in Figure 1. There were three main components in the experiment: FBRM, ATR-FTIR, and a crystallization system (OptiMax from Mettler Toledo, 1000 mL glass vessel), which was connected to a computer for data collection and equipment control. The vessel was equipped with bar-type baffles and agitated by four pitched-blade stirrers. The FBRM (D600L, Mettler Toledo) was set in the fine mode with the laser focused at 0 µm and set at a scanning speed of 2 m/s. The CLD histogram has 100 bins between 1 µm and 1000 µm, and data were recorded every ten seconds. The ATR-FTIR (ReactIR iC10, Mettler Toledo) measured the IR spectrum from 653 to 2998 cm-1 every minute. The React IR was purged by compressed air and cooled by liquid nitrogen. The OptiMax included a temperature probe, stirring system, 1-L glass crystallizer, and heating/cooling metal jacket. The temperature range in this work was from 0 to 70 ◦C and the stirring speed was fixed at 400 rpm,

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 49

ensuring good mixing in the crystallizer. These instruments communicated with the computer via the iC software provided by Mettler Toledo (iC FBRM 4.2.234, IC IR 4.3.27, iControl 5.1.29).

Figure 1: Experimental setup 2.4 Experimental procedure The conventional linear cooling strategy and proposed temperature plateau strategy are shown in Figure 2 (a), and (b), respectively. In the conventional linear cooling strategy, the temperature decreases linearly, and thus if the induction times in two replicate runs, tind,1 and tind, 2 are not identical due to stochasticity, there exist corresponding two nucleation temperatures, T1 and T2, and thus the nucleation temperature cannot be determined uniquely. On the other hand, in the proposed temperature plateau strategy43, the temperature profile can be designed in such a way that solution temperature is held at Tplateau for a sufficiently long time. This temperature is 10


Page 11 of 49

determined using the solubility data so that the solution becomes saturated at this temperature. In this strategy, even if there is some uncertainty in the induction time, nucleation occurs only at a single temperature Tplateau, which allow us to determine the nucleation temperature uniquely, and the cooling in the second stage begins only after nucleation subsides. Thus we anticipate this approach reduces the influence of stochasticity of induction time on the crystal size distribution of the product.

(a)

(b)

Temperature

Temperature

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


T1

Tplateau

T2

tind,1

tind,2

tind,1

time

tind,2

time

Figure 2 Illustration of temperature manipulation strategy: (a) conventional linear cooling; (b) proposed temperature plateau approach Employing the proposed temperature plateau strategy, three runs were carried out at conditions shown in Table 1. Each run began with the crystallizer containing 0.5 L of undersaturated paracetamol in ethanol. Prior to starting a run, the solution was maintained at 70 ºC for one hour to insure complete dissolution. The solution then was cooled (0.5 ºC/min) to an intermediate temperature (a constant-temperature plateau, Tplat), where the system was held for

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 49

two hours. The solution was observed to remain clear prior to reaching the plateau, but on the temperature plateau the reading from the FBRM and the IR changed due to nucleation. Supersaturation (Splat) at the plateau temperature induced primary nucleation, and the FBRM and IR were used to confirm that the crystal population was fully developed and that supersaturation had been depleted after the two-hour plateau period. Subsequently, another cooling stage was started and the solution was cooled to the final temperature of 0 ºC. Table 1: Experimental conditions for parameter estimation Initial concentration Name [g solute / kg solvent]

1st cooling rate [˚C/min]

Run 1 Run 2

370

0.5

Run 3

Tplat [˚C]

Splat after 1st cooling

40

1.4

44

1.3

46

1.2

2nd cooling [˚C/min]

0.25

When a crystallization run was completed, the crystals were recovered by filtration, washed with toluene, dried in an oven, and sieved to nine size fractions. The mass of crystals left on each sieve tray was weighed for estimating the volume density of the CSD. Observation under a microscope indicated the crystals were octahedral (kv = 1/3), which is typically associated with Form I of paracetamol 44. No significant agglomerates were found.

12


Page 13 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.5 Experimental data processing Three types of data could be used in the model-fitting process: the supersaturation profile (S vs. t), the CSDs estimated from the CLD histograms, and the CSD of the final product obtained by sieve analysis. The first two types of measurements were discussed in our earlier publications31, 32, and are summarized only briefly here. During batch cooling crystallization, the IR spectrum of the solution was monitored by ATR-FTIR. As the solution concentration and temperature changed, the heights and widths of peaks in the IR spectrum also changed. There are several chemometric methods available to quantify the solution concentration 45 from these data, including peak-height (or area) regression, principal-component regression, and partial least-squares regression. We found in our previous study 32 that correlating the ratio of two peak heights is simple and sufficiently accurate for our two-component system. As shown in Figure 3, two representative peaks were chosen for paracetamol (1667 cm-1) and ethanol (1048 cm-1), respectively. IR spectra were collected for various known solute concentrations and temperatures, and the concentrations were fit to a polynomial equation, in which paracetamol concentration cs is a function of both T and the peakheight ratio of paracetamol to ethanol. Solubility of paracetamol in ethanol is approximately by 32:

cs = 7.915×10−4T 3 − 6.439×10−1T 2 +1.765×102T −1617

(8)

where T is the solution temperature in K, and the unit of cs is in g-solute/kg-solvent.

13



1.2 Ethanol Paracetamol/ethanol

Ethanol 1 0.8

AU

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 49

0.6 0.4 Paracetamol

0.2 0 −0.2 800

1000

1200

1400

1600

1800

Wavenumber (cm−1)

Figure 3: Infrared spectrum of paracetamol dissolved in ethanol. Peaks at 1667 and 1048 cm-1 are used for estimation of the paracetamol concentration In our previous work 31, 32, we used CLD histograms generated by crystals of different sizes to develop a model that can translate CLD to CSD. The calibration procedure resulted in “fingerprint” CLD histograms by normalizing the CLD histogram with the number of crystals per unit volume of solution. The procedure was repeated for each of nine size fractions, as shown in Figure 4. The heights of the fingerprints suggest that the CLD of a crystal population is dominated by larger crystals, which leads to the ill-conditioning of the estimation problem given below. To deal with this challenge, we formulate a carefully tailored least-square problem. The fingerprints from each crystal size were used to construct the linear model:

b = Ax (9)

14


Page 15 of 49

where b ∈ R nb b is the CLD histogram measured by the FBRM with nb = 100, x ∈ R nx is the CSD histogram with the same bin partition as used in sieving, and A ∈ R nb ×nx is the matrix that maps x to b, consisting of nine fingerprint CLDs as its column vectors

Fingerprint [# of chords ⋅ mL/# of crystals]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


30 53−75 µm 75−106 µm 106−150 µm 150−212 µm 212−250 µm 250−300 µm 300−355 µm 355−425 µm 425−500 µm

25 20 15 10 5 0 200

400

600

800

Chord length [µm]

Figure 4: Fingerprints of paracetamol crystals of different sizes (nx = 9). To estimate x, the linear system is inverted using regularized regression. 2

min Ax - b 2 + λ f (y )

(10)

x

s.t.

xi ≥ 0, i =1,2,..., nx

yi = li3

(11)

xi , i = 1, 2..., nx Li − Li −1

(12)

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 49

2

 y − 4 yi+1 + 3yi  f (y) = ∑ i+2  , y0 = 0 (13) li+2 − li i =0   r −2

li = Li Li +1 , i = 1,2,..., nx (14) where Li and Li+1 are the boundaries of the ith size range, and li is the geometric mean of the boundary values in the ith size range. The first constraint, Equation (11), enforces non-negative values in the CSD histogram. Equation (12) defines the volume-weighted crystal size density. Equation (13) defines a penalty term that suppresses oscillations in the CSD estimate, which is the sum of squares of the first-order derivative by the forward difference approximation. The weight of the penalty term λ is chosen to be 1 × 10-10 [µm/mL] as used in our previous study 32.

2.6 Parameter estimation Kinetics parameters were estimated by minimizing the error between the experimental data and the model predictions from solutions of the PBE. Four indices are needed to define the estimation problem: (1) the type of measurements, using subscripts or superscripts “S” for supersaturation ratio and “sv” for sieving results; (2) the number of experiments (runs) Nr for each type of measurement; (3) the number of sample points Nd in each experiment; and (4) the number of measured variables Nm in each sample. The objective function includes errors from supersaturation ratio and sieving.

Φ(θ) = wSeS (θ) + wsvesv (θ) (15)

16


Page 17 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where N r ,s N d ,i

eS (θ) = ∑∑ ( Sîj − Sij ) 2 (16) i =1 j =1

Nr ,sv Nm ,sv

esv (θ) = ∑ ∑ (vîk − vik )2 (17) i =1 k =1

In the equations above, θ = {kb1,σ , kb2 ,α, β , kg , Ea , γ } is the parameter set, eS is the fitting error to supersaturation, esv is the fitting error to sieve analysis, wS , wsv are the weights for the error terms, Sîj is the supersaturation measured, vik is the volumetric density of the crystal population, and vîk is the volumetric density measured by sieving. Solution obtained from the CE/SE scheme is synchronized with experimental data by interpolation in the temporal domain. Since the CSD measurements from sieving only have nine bins, the simulated volumetric density of the CSD is averaged within the nine size ranges. The sum of the squared errors Φ(θ) is minimized by the fminsearchbnd function in MATLAB R2009a, which uses a derivative-free Nelder-Mead method 46. In this work, we use volumetric density to show CSDs instead of number density. These two representations have their own advantages and disadvantages. One of the advantages for volume density is that it is proportional to the crystal mass distribution, which was obtained directly from experimental sieve analysis. By employing volume density, we can avoid quantifying the number of crystals, which can be difficult for small crystals. 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 49

The confidence region of the parameters is calculated according to Rawlings et al. 47 and Bard 48. The confidence region around the estimated parameter θˆ is estimated by the following quadratic form. (θ − θ$ )T Vθ−1 (θ − θ$ ) ≤ χ N2 p ,α

(18)

in which Vθ is the convariance matrix of θ from different measured variables,

(Vθ )

−1

= (VθS ) + (Vθsv ) −1

−1

(19)

The sum of squared errors are assumed to follow a chi-square distribution with degree of freedom Np and α = 0.05 for 95% confidence, where Np is the number of parameters. For each type of measured variable, the covariance matrix Vθq is obtained from the model sensitivity Bj and the variance of the measurements Vq :

(V )

q −1

θ

= ∑ ( Bqj )T (V q ) −1 Bqj , q ∈ {S , sv}

(20)

j

The variance of measurement Vq is approximated by the fitting error:

Viiq =

1 N d ,q

N d ,q

∑e

2 i ,q

, q ∈ { S , sv} (21)

i =1

18


Page 19 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where Nd is the number of sample points considered in the experiments. The matrixBj represents the sensitivity of the jth sample to parameter set θ, and Bj ∈ RNm ×Npis obtained by a finite difference approximation.

B kj, q ≈

y kj, q (θ$ + hk e k ) − y kj, q (θ$ ) hk

, q ∈ { S , sv}

(22)

in which k = 1, 2, ..., Nm is the ith measured variable, j = 1, 2, ..., Nd is the jth sample points, and hk ek is the perturbation given to θˆ . In this study, Nm = 1 for supersaturation, Nm = 9 for sieving and in-situ CSD estimates, and hk ek is a 0.1% variation of each estimated parameter.

3.

Results

3.1 Two-stage cooling experiments and measurements The experimental conditions for Runs 1–3 are listed in Table 1, which are also visualized in Figure 5 (a). The only differences among the experiments are the plateau temperature Tplat, so that the importance of Splat can be captured from experiments. Other factors that may affect the process, such as cooling rates and total mass of solute, were not varied. Figure 5 shows the change of supersaturation during the course of each run and the final crystal size distribution obtained by sieving. The plateau temperature and corresponding supersaturation are given by the symbols Tplat and Splat. The in-situ supersaturation calculated from the ratio of measured concentration over the solubility is shown in Figure 5(b). Figure 5(c) shows

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 49

the sieve analysis of the recovered product crystals from each of the runs; it shows the volume density function v, which is related to the mass density function m by the expression m = ρkvv. The plateau temperatures were maintained for a period of two hours, during which time the crystals resulting from primary nucleation were observed and the supersaturation had been consumed (S = 1.0). As illustrated in Figure 5(a), a higher Splat was achieved by using a lower Tplat. The depletion of supersaturation also took a longer time for cases of lower Splat. When the cooling was resumed in the second cooling stage, supersaturation increased. Upon reaching 0 °C, the temperature was held constant and supersaturation gradually dropped back to 1.0.

20


Page 21 of 49

(a)

(b)

80

1.4

S

= 1.4

S

= 1.3

S

= 1.2

plat

Temperature [ °C]

60

plat plat

S

plat

= 1.4

Splat = 1.3

1.3

S

plat

40

= 1.2

S

1.2

20

1.1

0

1

−20 0

100

200 300 Time [min]

10

400

500

0

100

200 300 Time [min]

400

500

(c)

8

Volume density [ µm3/µm/mL]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


x 10

Splat = 1.4 8

S

plat

= 1.3

Splat = 1.2

6 4 2 0 0

100

200

300 Size [µm]

400

500

600

Figure 5: (a) Supersaturation in Runs 1–3; (b) volumetric density of final products in Runs 1-3 Crystals from each of the runs were recovered, processed, and sieved as described earlier. The data in Figure 5(b) show that the dominant size (peak) and spread of the distribution increased with a decrease in Splat. Since primary nucleation is usually more sensitive to supersaturation, we suspect that nucleation had a greater dominance over growth at a higher supersaturation. With more crystals produced at higher Splat, the second phase of cooling resulted in a fixed solute mass spread over a larger number of crystals and, thus, smaller product crystals. Also note the spread of the distributions increased as Splat decreased. This is thought to result from the increased

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 49

consumption of solute by crystal growth relative to nucleation at decreasing Splat, as well as by secondary nucleation. 3.2 Applicability of CSD Estimates from FBRM Figure 6 provides comparisons of crystal size distributions obtained from the final FBRM measurements and the subsequent sieve analyses for each run. The results show excellent agreement for Run 2, good agreement for Run 1, and relatively poor agreement for Run 3; we attribute the disparity in Run 3 to the crystal size distribution produced in that run, which produced significant numbers of crystals outside the range of sieve analysis (more than 18% of the final mass was greater than 500 µm), making the transformations from CLD to CSD less accurate. Given the mismatch, the in-situ CSD estimates were not used for kinetics parameters estimation but used for comparisons with the PBE model later.

22


Page 23 of 49

(a)

8

x 10

Volume density [ µm /µm/mL]

Sieving CSD estimates

3

8 6 4 2 0 0

100

200

300 400 Size [µm]

500

x 10

8 6 4 2 0 0

600


10

100

200

300 400 Size [µm]

500

600

(c)

8

12

(b)

8

12

10

3

Volume density [ µm /µm/mL]

12


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


x 10


10 8 6 4 2 0 0

100

200

300 400 Size [µm]

500

600

Figure 6 CSD of final product estimated by the CLDs obtained by FBRM, compared with sieving analysis: (a) Run 1; (b) Run 2; (c) Run 3 3.3 Parameter estimation The data set used in parameter estimation contains supersaturation and final volumetric density functions, from Runs 1 through 3. Given that each error source has a different order of magnitude, weighting coefficients were adjusted to balance the significance of each error term. They were chosen to be wS = 100 and wsv = 2 × 10-7. Table 2: Estimated kinetic parameters 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

kb1

σ

kb2

No./s/kg solvent

mJ/m2

No./s/kg solvent

0.192

4.25

1.00×105

α

2.08

β

0.713

Page 24 of 49

kg

Ea

γ

(m/s)(g/g)γ

kJ/mol

45.5

41.3

1.24

±0.5

±0.04

95% confidence interval obtained by 0.1% perturbation 2.28×10−4 –1.63×102

4.18–4.37

(0.82–1.23)×105 ±0.10 ±0.017 40.3-51.5

* In the parameter estimation, log kb1, σ3, log kb2, and ln kg were used in the minimization of the fitting error and the evaluation of the 95% confidence intervals. Since the logarithm and cube are nonlinear transformations, the estimated parameters may not center in their confidence intervals.

The parameters estimated from the experimental data are given in Table 2, which we can compare to previous work to confirm that our parameters are in the same order of magnitude as follows:

•

The estimated value of interfacial energy is 4.25 mJ/m2, while that predicted

by the method from Mersmann 49 is 7.60 mJ/m2.

•

The secondary nucleation parameters for the paracetamol-ethanol system

were estimated by Worlitschek and Mazzotti 15 with seeded experiments, but the equation was in a different form. Calculations with the parameters and mathematical expressions from the reference suggest that the average secondary nucleation rate in their seeded run was around 107/min/(kg solvent). Using our model, the secondary nucleation rate is predicted to be 105 –107 /min/(kg solvent) during the second phase of cooling.

•

Growth kinetics reported in other studies 15, 39 were estimated from seeded

paracetamol crystallizations. They are compared with the present work in Table 3, which 24


Page 25 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


shows the estimated activation energy to be very close to the literature values. The growth rate constant kg is also in the same order of magnitude. The exponents of crystal growth from the three studies are all between one and two, which is often considered as a reasonable range 10. Table 3: Comparison of growth parameters with other studies

kg [(m/s)(m3/kmol)γ

Ea [kJ/mol]

γ

Worlitschek and Mazzotti 15

21

41.6

1.9

Mitchell et al. 39

9.979

40.56

1.602

This study

5.95

41.3

1.24

Source

The 95% confidence intervals on the parameters also are shown in Table 2. The primary nucleation constant kb1 varies within six orders of magnitude at the 95% confidence level, which is much wider than the confidence interval for other parameters. This suggests that kb1 is insensitive to the measurements, since secondary nucleation is dominant over primary nucleation. It should be noted that these confidence intervals may contain error due to the finite difference approximation (Equation (22)), and may not represent the true sensitivities due to the quadratic approximation (Equation (18)); we use the confidence intervals to identify potentially insensitive parameters but not to evaluate sensitivities precisely.

25



Figure 7 and Figure 8 provide comparisons between the model and experimental data. The plots of supersaturation as a function of run time in Figure 7 show that the set of parameters provided a good fit to the three runs, including the onset of primary nucleation and the second phase of cooling. Plots of the final CSD in Figure 8 show close agreement with the sieve results. The model also captures the influence of Splat on the final CSD, which is that higher values of Splat lead to smaller mean sizes and narrower final CSDs.

(a)

(b)

1.5 Exp Fitted

1.4 1.3

1.3

1.2

S

S

1.2

1.1

1.1

1

1

0.9

0.9

0.8 0

Exp Fitted

1.4

100

200 300 Time [min]

400

0.8 0

500

100

200 300 Time [min]

400

500

(c) Exp Fitted

1.4 1.3 1.2

S

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 49

1.1 1 0.9 0.8 0

100

200 300 Time [min]

400

500

Figure 7: Supersaturation-time profile in experiments and simulations: (a) Run 1; (b) Run2; (c) Run 3

26


Page 27 of 49



x 10 Sieving Fitting

8 6 4 2 0 0

100

200

(b)

8

(a)

8

x 10

300 400 Size [ µm]

500

6 4 2 0 0

600

Sieving Fitting

8

100

200

300 400 Size [ µm]

500

600

(c)

8

x 10 Volume density [ µm3/µm/mL]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Sieving Fitting

8 6 4 2 0 0

100

200

300 400 Size [ µm]

500

600

Figure 8: Final CSD volume density in experiments and simulations: (a) Run 1; (b) Run 2; (c) Run 3 3.4 Verification of reproducibility A run identical to Run 2 was carried out to test the reproducibility of the experimental results. In Figure 9(a), supersaturation curves of these two identical runs are plotted. In these figures, desupersaturation (induction time) is observed at approximately 100 minutes, which is probably triggered by primary nucleation. A difference in the induction time of approximately 6 minutes is observed between the two identical experiments due to the stochastic nature of primary

27



nucleation. Nevertheless, as shown in Figure 9(b), the difference did not impact the CSD of the final crystal product. This reproducibility is the desired outcome of the temperature plateau strategy (Figure 2) . The rest of the crystallization course was almost the same, including the desupersaturation curve on the constant-temperature plateau and in the second phase of cooling.

1.2

1

0

100

200 300 Time [min]

400


Run 2 Repetition

1.4

(b)

8

x 10

(a)

S

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 49

6 4 2

0 0

500

Run 2 Repetition

8

100

200 300 Size [ µm]

400

500

Figure 9: Run 2 and its repetition: (a) S profile; (b) sieving results 3.5 Model validation with seeded crystallization A seeded crystallization run, which only involves secondary nucleation and growth, was performed to see if the model is able to predict the crystallization process. A mass of 0.03 g of crystals from sieve fraction 75 to 106 µm was added to an agitated supersaturated solution that had been generated by cooling a 370g/kg solvent to 48 °C (S=1.16). As can be seen in Figure 10(a), supersaturation decreased immediately once the seeds were added (t = 0). After seeding, the temperature was held at 48 °C until desupersaturation was no longer observed. A second cooling stage reduced the temperature by 0.40 °C/min until the temperature reached 0 °C (Figure 11(a)).

28


Page 29 of 49

The resulting supersaturation profile is shown in Figure 10(b), and the CSD of the final crystal product, in number density and volume density, is shown in Figure 10(c) and (d), respectively. (b)

Temperature [°C]

(a) 60

1.4

50

1.3

Experiment Prediction

40

1.2 S

30 20

1.1

10

1 0 50

3

150 200 Size [µm]

250

0.9 0

300

x 10

100

150 200 Time [min]

6 4 2 200 300 Size [µm]

300

(d)

8

100

250

1200

Prediction Experiment

10

0 0

50

(c)

8

12

100

Number density [#/µm/mL]

−10 0

Volume density [µm /µm/mL]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


400

800 600 400 200 0 50 100

500

Prediction Experiment

1000

200

300 Size [µm]

400

500

Figure 10 Seeded crystallization: (a) Temperature profile, (b) supersaturation profile; (c) final CSD in volume density; (d) final CSD in number density. Figure 10(b) shows discrepancies between the measured and predicted supersaturations. One reason could be that the initial CSD of seeds provided to the simulation was inaccurate; the actual seeds could contain fines adhering to their surfaces, known as initial breeding

50

. Furthermore, the

seeds, obtained from a separate crystallization run, may have acquired some surface features in seed 29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

preparation steps

51

Page 30 of 49

,such as washing, drying, and milling, which could affect the secondary

nucleation and growth rates. These effects are not included in the model and made the predicted value of S decrease more slowly on the temperature plateau and increase to higher values in the second cooling stage than the experimental observation. The initial breed and surface features, difficult to determine for the seeds, are the reason why seeded crystallization was excluded from parameter estimation. The CSD obtained from sieving and predicted by model are compared in Figure 10(b), which shows some overestimation of crystal growth. Again, this is consistent with the higher supersaturation in the simulation than the actual experiment. However, the shape of the simulated CSD based on number density is close to the sieving result, which can be seen in Figure 10(c).

4.

Discussion

4.1 Comparison of model prediction with in-situ CSD estimation The evolution of CSD estimates from FBRM for Runs 1 and 2 are shown in Figure 11(a) and Figure 12(a), respectively, which may be compared with the solution of the PBE in Figure 11(b) and Figure 12(b). The independent axes on these three-dimensional plots are crystal size and time, while the volume density is plotted as the dependent variable on the vertical axis. Notice that the time references t = 0 move to the beginning of the constant-temperature plateau, since no crystals were generated prior to this time. Both in-situ CSD estimates and model predictions show two major developments in the volumetric density functions. The first is a rapid change at approximately 50 min that corresponds 30


Page 31 of 49

to the onset of primary nucleation and subsequent crystal growth on the constant- temperature plateau. Then the CSD exhibits little change over the remainder of the plateau. The second major development starts with the second phase of cooling, which begins at approximately 200 min. It should be noted that the estimated CSDs were used only for comparison with the population balance model in a qualitative manner. We found that including the FBRM estimated CSDs in the parameter fitting does not strongly influence the parameter values or reduce the confidence regions from the ones in Table 1. This is because the accuracy of the CSD estimation by FBRM is not sufficient to improve the estimation; crystals smaller than 53 µm were ignored in the fingerprint approach (Figure 4), which led to mass balance errors.

8

x 10

8

x 10

8


6

3

Volume density [µm3/µm/mL]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4 2

0 600 400 200 0 Time [min]

0

8 6 4 2

0 600 400 200 0 Time [min]

100 200 300 400 500 Size [µm]

(a)

0

400 500 100 200 300 Size [µm]

(b)

Figure 11: Evolution of CSD volume density in Run 1: (a) estimates from CLD; (b) model predictions

31



8

x 10 8


x 10 8 6

8 6

3

Volume density [µm3/µm/mL]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 49

4 2

0 600 400 200 0 0 Time [min]

4 2

0 600 400 200 Time [min]

400 500 100 200 300 Size [µm]

(a)

0

0

500 100 200 300 400 Size [µm]

(b)

Figure 12: Evolution of CSD volume density in Run 2: (a) estimates from CLD; (b) model predictions 4.2 Analysis of nucleation Since the characteristics of products from batch crystallization are governed primarily by nucleation, determination of the mechanisms by which product crystals are generated is of great interest. Results from the present work provide an opportunity to examine whether primary or secondary nucleation was dominant in the unseeded batch runs of the study. Of course, primary nucleation is critical, as the process is unseeded; however, our concern is whether more crystals in the final population were produced by primary or secondary nucleation. This was determined by using the nucleation models given earlier (Equations 5 and 6) and the parameters in Table 2 to estimate numbers of nuclei generated by each mechanism as follows: 32


Page 33 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


t

N primary = ∫ B1 (t )dt

(23)

0

t

N secondary = ∫ B2 (t )dt

(24)

0

Values of parameters in Table 2 were used in preparing Figure 13, which compares the total number of nuclei (Nprimary + Nsecondary) and the number of nuclei generated by primary nucleation. The number of secondary nuclei is the difference between the two quantities. Both coordinates are in logarithmic scales to show the wide variations of the values, and the time periods corresponding to the temperature plateau are marked so that the first and second cooling stages can be seen. It should be noted that the rapid increase of primary nuclei prior to the plateau is caused by the generation of primary nucleation due to the supersaturation resulting from the first cooling stage. Little difference initially exists between the numbers of total nuclei and primary nuclei, but the gap between the two curves expands in the temperature plateau. In the constant-temperature plateau that follows, the number of primary nuclei increases only slowly. However, the number of total nuclei keeps increasing through secondary nucleation. The increase of total nuclei slows when supersaturation again starts to decrease.

33



log(N) [log(# of nuclei/kg solvent)]

10

5

0

T plateau

−5 Total nuclei Primary nuclei −10 32

64

128 Time [min]

256

512


(b)

(a)

10

5

0 T plateau


64

128 Time [min]

256

512

(c)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 49

10

5

0 T plateau


64

128 Time [min]

256

512

Figure 13: Total number of crystals calculated by the primary and secondary nucleation models: (a) Run 1 (Splat = 1.4); (b) Run 2 (Splat = 1.3); (c) Run 3 (Splat = 1.2). The time and nucleation rates are shown in logarithmic scales. Figure 13 shows that in all three runs far more crystals resulted from secondary nucleation than from primary nucleation. Because no seed crystals were present at the beginning of the runs, primary nucleation had to occur before secondary nucleation. This means that only a few crystals were generated by primary nucleation, which then initiated subsequent secondary nucleation.

34


Page 35 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Analyses of the experimental observations are consistent with the assessment in the preceding paragraph. For example, the data in Figure 9(a) show that at the beginning of the constant-temperature plateau (between 84 and 90 minutes), supersaturation remained constant. After a delay of approximately 6 minutes, the supersaturation ratio dropped suddenly to nearly 1.0, which is likely the result of rapid nucleation and subsequent growth. The nucleation model suggests this behavior is the result of the following two steps: (1) primary nucleation and subsequent crystal growth during the period between 40 and 90 minutes; and (2) when the crystal mass became sufficient, massive secondary nucleation was initiated between 90 and 110 minutes. In this self-catalyzed manner, the nucleation rate accelerates until a rapid desupersaturation occurred. The experimental observation described above cannot be explained by the primary nucleation model given in Equation (5), and the conventional assumption that primary nucleation dominates the formation of crystals in unseeded batch crystallization. The model for primary nucleation given in Equation (5) cannot describe a delay in the formation of crystals that is consistent with the experimental observations; i.e. the primary nucleation rate must be constant for a given value of S. Thus, the nucleation rate should increase immediately once the supersaturation S increases, and desupersaturation should also happen immediately when the supersaturation reaches Splat. However, this is inconsistent with our experimental observation in which supersaturation remained constant for approximately 6 minutes before rapid desupersaturation began.

35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 49

The experimental observations and numerical analyses in the present work are consistent with several reference showing that under certain conditions secondary nucleation can dominate primary nucleation under certain conditions. For example, Kadam et al., 8, 29 reported observing only one or very few large crystals in an unseeded stirred volume prior to the massive nucleation event. It is also consistent with work in which chiral symmetry in was broken in in unseeded crystallization of sodium chlorate 30. The simulation results in Figure 13 provide a mechanistic interpretation of the induction time for the unseeded experiments. The initial generation of a few primary nuclei results in few crystals of small size. However, these crystals can induce secondary nucleation once they reach sufficient mass, and the required period of growth is correlated with Splat, which makes the induction time shorter for the higher Splat run, as shown in Figure 5(a). The cooling strategy using a constant-temperature plateau to induce nucleation is a powerful internal seeding method for better reproducibility. Unlike our proposed method, internal seeding is usually achieved by continuous cooling, with linear or nonlinear cooling profiles, until nuclei are detected 52-55. The solution may be cooled to an undesired supersaturation level which results in too many fine crystals. On the other hand, if the proposed cooling strategy is implemented with a plateau temperature selected carefully, the initial nucleation can be manipulated more easily, and the consistency of final product can be greatly improved, as shown in Run 2 and its repetition in Figure 9. Moreover, with the kinetic parameters, model-based cooling strategy can be applied to

36


Page 37 of 49

optimize the temperature plateau strategy, which can be considered a special case of piecewise linear parameterization 43, 56, 57. 4.3 Mean size of final product from different Splat Figure 14 shows volume-weighted final mean sizes (fourth moment over third moment, µ4/µ3) determined from the model equations and from experimental sieve data. Even that there is approximately 30 µm difference between the two, both show that mean size decreases with increasing Splat. These results imply that the two-stage cooling strategy can be used to manipulate the mean size of the final product by adjusting Splat in an unseeded process.

Prediction Sieving

350 325 300

4

3

Volume−weighted mean size (µ /µ ), [µm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


275 250 1.1

1.2

1.3 S

1.4

1.5

plat

Figure 14: The effect of Tplat on the volume-weighted mean size

37



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

5.

Page 38 of 49

Conclusion

Experimental and numerical approaches are presented for estimating crystallization kinetics in unseeded processes, using in-situ and ex-situ process measurements and population balance modeling. Unique approaches have been demonstrated in our case study, including (i) a cooling strategy with temperature plateaus; (ii) simultaneous estimation of primary nucleation, secondary nucleation, and crystal growth in unseeded cooling processes. Crystal size distributions and measured supersaturations from unseeded two-stage cooling experiments were used to decouple primary nucleation, secondary nucleation, and growth. Analysis of product characteristics showed that this cooling profile has advantages over a linear cooling profile. These result from the crystallization process being divided into two stages: An initial stage on the temperature plateau involves all three crystallization mechanisms; both primary nucleation and crystal growth occur, but secondary nucleation also may participate. The second cooling stage, starting with the crystal population created on the temperature plateau, can be viewed as a seeded crystallization run in which the primary nucleation rate is insignificant. Acknowledgement The financial support from C. J. “Pete” Silas Chair, Georgia Research Alliance, and Consortium for Risk Evaluation with Stakeholder Participation (CRESP) is gratefully appreciated.

38


Page 39 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nomenclature Abbreviations ATR-FTIR

attenuated-total-reflectance-Fourier-transform infrared

CE/SE

conservation element/solution element

CLD

chord length distribution

CSD

crystal size distribution

CSH

crystal size histogram

CLH

chord length histogram

FBRM

focused beam reflectance measurement

PBE

population balance equation

Greek letters Ω

confidence region

Φ

fitting error

α

exponent of supersaturation in secondary nucleation

39



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

β

exponent of mass in secondary nucleation

γ

exponent in growth equation

λ

tuning parameter in CLD-CSD transformation

µi

ith moment of crystal population

ν

molecular volume

ρ

density of crystals

θ

parameter set

σ

crystal-solution interfacial tension

ω

weight of error term

English letters A

transformation matrix, x to b

B

nucleation rate

B1

primary nucleation rate

B2

secondary nucleation rate

Bj

sensitivity of simulation to parameters

40


Page 40 of 49

Page 41 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ea

activation energy of crystal growth

G

growth rate

L

crystal size

L1,0

0

number-averaged size of crystal

L4,3

volume-averaged size of crystal

N

number of experiments, sample points, measurements

N1

number of nuclei generated by primary nucleation

N2

number of nuclei generated by secondary nucleation R gas constant

S

supersaturation, c/cs

T

temperature

V

covariance matrix of measurements

Vθ

covariance matrix of parameters

b

CLH vector

c

concentration

∆c

supersaturation, c − cs 41



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

cs

solubility

ek

natural basis on kth dimension

e

fitting error

h

small perturbation for finite difference

ind

induction time

k

Boltzmann constant

kb1

coefficient of primary nucleation

kb2

coefficient of secondary nucleation

kg

coefficient of crystal growth

kv

shape factor

l

geometric mean of a size range

ms

mass of crystals

m

mass density of CSD

n

number density of CSD

rc

critical radius of nuclei

42


Page 42 of 49

Page 43 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


t

time

v

volumetric density function of CSD

x

CSH vector

x

size of crystals

∆x

discretization of size domain in CE/SE scheme

y

calculated volume density

Subscripts 0

initial condition

is

In-situ CSD estimates from FBRM

S

supersaturation

d

sample points

m

measured variables plat temperature plateau r experimental runs

sv

sieving

v

volume density

43



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 49

References (1)

Mullin, J. W. Crystallization. Butterworth-Heinemann: 2001.

(2) Chen, J.; Sarma, B.; Evans, J. M. B.; Myerson, A. S. Pharmaceutical crystallization. Cryst. Growth Des. 2011, 11, 887-895. (3) Lawton, S.; Steele, G.; Shering, P.; Zhao, L.; Laird, I.; Ni, X.-W. Continuous crystallization of pharmaceuticals using a continuous oscillatory baffled crystallizer. Org. Process Res. Dev. 2009, 13, 1357-1363. (4) Vekilov, P. G. Dense liquid precursor for the nucleation of ordered solid phases from solution. Cryst. Growth Des. 2004, 4, 671-685. (5) Nielsen, M. H.; Aloni, S.; De Yoreo, J. J. In situ TEM imaging of CaCO3 nucleation reveals coexistence of direct and indirect pathways. Science 2014, 345, (6201), 1158-1162. (6) Kulkarni, S. A.; Kadam, S. S.; Meekes, H.; Stankiewicz, A. I.; ter Horst, J. H. Crystal Nucleation Kinetics from Induction Times and Metastable Zone Widths. Cryst. Growth Des. 2013, 13, 2435-2440. (7) Jiang, S.; ter Horst, J. H. Crystal nucleation rates from probability distributions of induction times. Cryst. Growth Des. 2010, 11, (1), 256-261. (8) Kadam, S. S.; Kulkarni, S. A.; Coloma Ribera, R.; Stankiewicz, A. I.; ter Horst, J. H.; Kramer, H. J. M. A new view on the metastable zone width during cooling crystallization. Chem. Eng. Sci. 2012, 72, 10-19. (9) Knezic, D.; Zaccaro, J.; Myerson, A. S. Nucleation induction time in levitated droplets. J. Phys. Chem. B 2004, 108, (30), 10672-10677. (10)

Myerson, A. Handbook of Industrial Crystallization. Butterworth-Heinemann: 2002.

(11)

Paul, E. L.; Tung, H.-H.; Midler, M. Organic crystallization processes. Powder Technol. 2005, 150, 133-143.

(12) Aamir, E.; Nagy, Z. K.; Rielly, C. D.; Kleinert, T.; Judat, B. Combined Quadrature Method of Moments and Method of Characteristics Approach for Efficient Solution of Population Balance Models for Dynamic Modeling and Crystal Size Distribution Control of Crystallization Processes. Ind. Eng. Chem. Res. 2009, 48, 8575-8584. (13) Mitchell, N. A.; Frawley, P. J.; O'Ciardha, C. T. Nucleation kinetics of paracetamol-ethanol solutions from induction time experiments using Lasentec FBRM®. J. Cryst. Growth 2011, 321, 91-99. (14) O'Ciardha, C. T.; Hutton, K. W.; Mitchell, N. A.; Frawley, P. J. Simultaneous Parameter Estimation and Optimization of a Seeded Antisolvent Crystallization. Cryst. Growth Des. 2012, 12, 5247-5261. (15) Worlitschek, J.; Mazzotti, M. Model-based optimization of particle size distribution in batch-cooling crystallization of paracetamol. Cryst. Growth Des. 2004, 4, 891-903. (16) Togkalidou, T.; Tung, H.-H.; Sun, Y.; Andrews, A. T.; Braatz, R. D. Parameter Estimation and Optimization of a Loosely Bound Aggregating Pharmaceutical Crystallization Using in Situ Infrared and Laser Backscattering Measurements. Ind. Eng. Chem. Res. 2004, 43, (19), 6168-6181.

44


Page 45 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(17) Hu, Q.; Rohani, S.; Wang, D. X.; Jutan, A. Nonlinear kinetic parameter estimation for batch cooling seeded crystallization. AlChE J. 2004, 50, 1786-1794. (18) Gherras, N.; Fevotte, G. On the use of process analytical technologies and population balance equations for the estimation of crystallization kinetics. A case study. AlChE J. 2012, 58, (9), 2650-2664. (19) Larsen, P. A.; Rawlings, J. B.; Ferrier, N. J. An algorithm for analyzing noisy, in situ images of high-aspectratio crystals to monitor particle size distribution. Chem. Eng. Sci. 2006, 61, 5236-5248. (20) Kempkes, M.; Eggers, J.; Mazzotti, M. Measurement of particle size and shape by FBRM and in situ microscopy. Chem. Eng. Sci. 2008, 63, (19), 4656-4675. (21) Presles, B.; Févotte, G.; Pinoli, J.-C.; Debayle, J. Novel image analysis method for in situ monitoring the particle size distribution of batch crystallization processes. J. Electron. Imaging 2010, 19, 31207. (22) Wang, X. Z.; Roberts, K. J.; Ma, C. Crystal growth measurement using 2D and 3D imaging and the perspectives for shape control. Chem. Eng. Sci. 2008, 63, (5), 1173-1184. (23) Ochsenbein, D. R.; Schorsch, S.; Vetter, T.; Mazzotti, M.; Morari, M. Growth Rate Estimation of β l-Glutamic Acid from Online Measurements of Multidimensional Particle Size Distributions and Concentration. Ind. Eng. Chem. Res. 2013. (24) Nagy, Z. K.; Fujiwara, M.; Woo, X. Y.; Braatz, R. D. Determination of the kinetic parameters for the crystallization of paracetamol from water using metastable zone width experiments. Ind. Eng. Chem. Res. 2008, 47, 1245-1252. (25) Kail, N.; Briesen, H.; Marquardt, W. Analysis of FBRM measurements by means of a 3D optical model. Powder Technol. 2008, 185, (3), 211-222. (26) Worlitschek, J.; Hocker, T.; Mazzotti, M. Restoration of PSD from chord length distribution data using the method of projections onto convex sets. Part. Part. Syst. Charact. 2005, 22, 81-98. (27) Czapla, F.; Kail, N.; Oncul, A.; Lorenz, H.; Briesen, H.; Seidel-Morgenstern, A. Application of a recent FBRM-probe model to quantify preferential crystallization of DL-threonine. Chem. Eng. Res. Des. 2010, 88, 14941504. (28) Kee, N. C. S.; Arendt, P. D.; Goh, L. M.; Tan, R. B. H.; Braatz, R. D. Nucleation and growth kinetics estimation for L-phenylalanine hydrate and anhydrate crystallization. CrystEngComm 2011, 13, 1197-1209. (29) Kadam, S. S.; Kramer, H. J. M.; ter Horst, J. H. Combination of a Single Primary Nucleation Event and Secondary Nucleation in Crystallization Processes. Cryst. Growth Des. 2011, 11, 1271-1277. (30) Kondepudi, D. K.; Laudadio, J.; Asakura, K. Chiral symmetry breaking in stirred crystallization of 1, 1'binaphthyl melt. J. Am. Chem. Soc. 1999, 121, (7), 1448-1451. (31) Li, H.; Grover, M.; Kawajiri, Y.; Rousseau, R. W. Development of an empirical method relating crystal size distributions and FBRM measurements. Chem. Eng. Sci. 2013, 89, 142-151. (32) Li, H.; Kawajiri, Y.; Grover, M.; Rousseau, R. Application of an empirical FBRM model to estimate crystal size distributions in batch crystallization. Cryst. Growth Des. 2014, 14, 607-616.

45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 49

(33) Chang, S.-C. In Courant number insensitive CE/SE schemes, Proceedings of the 38th AIAA Joint Propulsion Conference, Indianapolis, Indiana, USA, 2002; AIAA-2002-3890. (34) Lim, Y. I.; Chang, S. C.; Jorgensen, S. B. A novel partial differential algebraic equation (PDAE) solver: iterative space-time conservation element/solution element (CE/SE) method. Comput. Chem. Eng. 2004, 28, 1309-1324. (35) Chang, S. C. The Method of Space-Time Conservation Element and Solution Element - A New Approach for Solving the Navier-Stokes and Euler Equations. J. Comput. Phys. 1995, 119, 295-324. (36) Galkin, O.; Chen, K.; Nagel, R. L.; Hirsch, R. E.; Vekilov, P. G. Liquid-liquid separation in solutions of normal and sickle cell hemoglobin. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, (13), 8479-8483. (37) Kashchiev, D.; Vekilov, P. G.; Kolomeisky, A. B. Kinetics of two-step nucleation of crystals. J. Chem. Phys. 2005, 122, (24), 244706. (38) Vekilov, P. G. Dense liquid precursor for the nucleation of ordered solid phases from solution. Crystal Growth & Design 2004, 4, (4), 671-685. (39) Mitchell, N. A.; O'Ciardha, C. T.; Frawley, P. J. Estimation of the growth kinetics for the cooling crystallisation of paracetamol and ethanol solutions. J. Cryst. Growth 2011, 328, 39-49. (40) Lim, Y. I.; Jorgensen, S. B. A fast and accurate numerical method for solving simulated moving bed (SMB) chromatographic separation problems. Chem. Eng. Sci. 2004, 59, (10), 1931-1947. (41) Motz, S.; Mitroviã, A.; Gilles, E. Comparison of numerical methods for the simulation of dispersed phase systems. Chem. Eng. Sci. 2002, 57, 4329-4344. (42) Qamar, S.; Elsner, M. P.; Angelov, I.; Warnecke, G.; Seidel-Morgenstern, A. A comparative study of high resolution schemes for solving population balances in crystallization. Comput. Chem. Eng. 2006, 30, (6-7), 1119-1131. (43) King, J. C.; Li, H.; Grover, M. A.; Kawajiri, Y.; Rousseau, R. W. In Optimization of two-stage cooling profile in unseeded batch crystallization, Proceedings of the 9th IFAC International Symposium on Advanced Control of Chemical Processes (ADCHEM), Whistler, British Columbia, Canada, 2015; 48 Elsevier: pp 297-302. (44) Beyer, T.; Day, G. M.; Price, S. L. The prediction, morphology, and mechanical properties of the polymorphs of paracetamol. J. Am. Chem. Soc. 2001, 123, 5086-5094. (45) Cornel, J.; Lindenberg, C.; Mazzotti, M. Quantitative application of in situ ATR-FTIR and raman spectroscopy in crystallization processes. Ind. Eng. Chem. Res. 2008, 47, (14), 4870-4882. (46)

Nelder, J. A.; Mead, R. A simplex method for function minimization. Comput. J. 1965, 7, (4), 308-313.

(47) Rawlings, J. B.; Miller, S. M.; Witkowski, W. R. Model Identification and Control of Solution Crystallization Processes - A Review. Ind. Eng. Chem. Res. 1993, 32, 1275-1296. (48)

Bard, Y. Nonlinear Parameter Estimation. Academic Press New York: 1974.

(49)

Mersmann, A. Calculation of interfacial-tensions. J. Cryst. Growth 1990, 102, 841-847.

46


Page 47 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(50) Girolami, M. W.; Rousseau, R. W. Initial breeding in seeded batch crystallizers. Ind. Eng. Chem. Process Des. Dev. 1986, 25, (1), 66-70. (51) Aamir, E.; Nagy, Z.; Rielly, C. Evaluation of the effect of seed preparation method on the product crystal size distribution for batch cooling crystallization processes. Cryst. Growth Des. 2010, 10, (11), 4728-4740. (52) Doki, N.; Seki, H.; Takano, K.; Asatani, H.; Yokota, M.; Kubota, N. Process control of seeded batch cooling crystallization of the metastable α-form glycine using an in-situ ATR-FTIR spectrometer and an in-situ FBRM particle counter. Cryst. Growth Des. 2004, 4, (5), 949-953. (53) Abu Bakar, M. R.; Nagy, Z. K.; Saleemi, A. N.; Rielly, C. D. The impact of direct nucleation control on crystal size distribution in pharmaceutical crystallization processes. Cryst. Growth Des. 2009, 9, (3), 1378-1384. (54) Saleemi, A.; Rielly, C.; Nagy, Z. K. Automated direct nucleation control for in situ dynamic fines removal in batch cooling crystallization. CrystEngComm 2012, 14, (6), 2196-2203. (55) Saleemi, A. N.; Rielly, C. D.; Nagy, Z. K. Comparative investigation of supersaturation and automated direct nucleation control of crystal size distributions using ATR-UV/vis spectroscopy and FBRM. Cryst. Growth Des. 2012, 12, (4), 1792-1807. (56) Shi, D.; El-Farra, N. H.; Li, M. H.; Mhaskar, P.; Christofides, P. D. Predictive control of particle size distribution in particulate processes. Chem. Eng. Sci. 2006, 61, 268-281. (57) Nagy, Z. K.; Chew, J. W.; Fujiwara, M.; Braatz, R. D. Comparative performance of concentration and temperature controlled batch crystallizations. J. Process Control 2008, 18, (3-4), 399-407.

47



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC graphic

48


Page 48 of 49

Page 49 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC graphic 44x23mm (600 x 600 DPI)


Modeling of Nucleation and Growth Kinetics for Unseeded Batch

Recommend Documents