Secondary Nucleation by Interparticle Energies. II. Kinetics.

This work presents a mathematical model that describes growth, homogeneous nucleation, and secondary nucleation that is caused by interparticle interactions between seed crystals and molecular clusters in suspension. The model is developed by incorporating the role of interparticle energies into a kinetic rate equation model, which yields the time evolution of nucleus and seed crystal populations, as in a population balance equation model, and additionally that of subcritical molecular clusters, thus revealing an important role of each population in crystallization. Seeded batch crystallization at a constant temperature has been simulated to demonstrate that the interparticle interactions increase the concentration of the critical clusters by several orders of magnitude, thus causing secondary nucleation. This explains how secondary nucleation can occur at a low supersaturation that is insufficient to trigger primary nucleation. Moreover, a sensitivity analysis has shown that the intensity of the interparticle energies has a major effect on secondary nucleation, while its effective distance has a minor effect. Finally, the simulation results are qualitatively compared with experimental observations in the literature, thus showing that the model can identify operating conditions at which primary or secondary nucleation is more prone to occur, which can be used as a useful tool for process design.

Introduction

Crystallization is an essential unit operation used as a purification technique in various industries for the production of chemicals, agrochemicals, nutritionals, and pharmaceuticals. Crystallization process is often initiated by the addition of seed crystals (seeds) of the target product in a supersaturated solution. While these seeds grow in the solution, they can generate new crystals (i.e., nuclei), through the phenomenon called secondary nucleation, which is widespread in industrial crystallization,[1] in particular, in the continuous crystallization of active pharmaceutical ingredients.[2] Hence, to operate and design continuous crystallization efficiently, we need a sound understanding of secondary nucleation mechanisms.

Traditionally, secondary nucleation has been largely attributed to microbreakages of seed crystals caused by crystal-impeller collisions,[3] that is, attrition.[1,4,5] However, recent research has noted that this mechanism cannot describe secondary nucleation completely,[3,6−13] based on two main reasons. First, secondary nucleation can occur even when the crystal microbreakages are prevented by precluding any crystal-impeller collision.[11−13] Second, the attrition mechanism cannot explain some crucial experimental observations according to which secondary nucleation can produce nuclei having a polymorphic[14] or chiral[15−17] form different from that of the seeds. In contrast with these observations, secondary nuclei generated by the attrition mechanism must exhibit the same polymorphic form and the same handedness in the case of chiral molecules of the seeds because these nuclei are the broken pieces of the seeds. These reasons clearly highlight that when describing secondary nucleation, we need to consider not only the attrition mechanism but also other secondary nucleation mechanisms.

Another perspective for rationalizing secondary nucleation suggests that secondary nucleation is induced by interparticle interactions,[18] in particular, the interaction energies between seed crystals and nearby molecular clusters. From this perspective, the interactions facilitate nucleation by reducing the energy barrier for nucleation in the vicinity of the seeds. This type of secondary nucleation can be called secondary nucleation caused by interparticle energies (SNIPE). As opposed to secondary nucleation by attrition, the SNIPE mechanism can explain why some secondary nuclei exhibit a polymorphic/chiral crystal structure different from that of the seeds; this is because the interparticle energies decrease the energy barrier for nucleation independent of the crystal structure of the relevant cluster,[15,18] thus causing secondary nucleation of any crystal structures.

A simple version of the SNIPE mechanism is the embryo-coagulation secondary nucleation (ECSN),[18] a key mechanism according to recent reviews on secondary nucleation.[1,4,5] However, despite its important contribution to our understanding of the SNIPE mechanism, the ECSN model has three critical limitations. First, the ECSN model assumes that subcritical molecular clusters are monodispersed. This assumption contradicts the classical nucleation theory[19] (CNT) whereby a cluster population develops in solution following a Boltzmann distribution. Second, the ECSN model further assumes that the subcritical molecular clusters grow by aggregation only, thus being inconsistent with the CNT where the clusters grow via a sequence of molecule attachments. Third, the ECSN model cannot describe the temporal evolution of any of the three key populations in crystallization—the (subcritical) molecular clusters, nuclei, and seed crystals—thus limiting a deeper understanding of how the SNIPE mechanism works.

The goal of this work is to develop a mathematical model that simultaneously describes homogeneous nucleation, growth, and the SNIPE phenomenon while overcoming the three limitations of the ECSN model. To this aim, the present work extends the conventional kinetic rate equation (KRE) model of nucleation by including a description of the SNIPE mechanism. The structure of this article is as follows. In Section , we provide an overview of the conventional KRE model that describes homogeneous nucleation and growth while being consistent with the CNT. In Section , we incorporate the description of the SNIPE mechanism into the conventional KRE model. In Section , using the developed model, we demonstrate how the SNIPE mechanism works, show the results of a sensitivity analysis, and compare the simulation results qualitatively with available experimental observations.

Conventional KRE Model of Nucleation

To describe the SNIPE mechanism, we extend the conventional KRE model of nucleation (i.e., the Szilard model[19,20]) by integrating the effect of the interparticle energies. Hence, knowing the working principle of the Szilard model is essential for understanding how our SNIPE model works. To give a brief overview, the basic principle of the model and its governing equations are presented in this section, followed by equations for the attachment and detachment rate constants that are necessary for solving the governing equations.

Principle

The Szilard model is a well-established theoretical framework that describes homogeneous nucleation and growth via a series of molecule attachments to and detachments from the molecular clusters,[19,21] that is, a description consistent with the CNT.[19] In this approach, a molecular cluster consisting of n solute molecules (i.e., a n-sized cluster) represents a solid particle (i.e., a crystal), when its size n is larger than the critical nucleus size n* (i.e., n > n*) defined in the CNT.[19] Otherwise (i.e., when n ≤ n*), the n-sized cluster is considered as part of the liquid phase; hence, its volume is excluded from the calculation of the solid-phase volume as in the literature.[19] The n-sized clusters can grow to (n+1)-sized clusters via molecule attachments, while (n+1)-sized clusters decay to the n-sized clusters through molecule detachments. This mechanism can be viewed as the following pseudo reaction[19,21]where W represents a n-sized cluster, ka is the rate constant of monomer attachment to the n-sized clusters, and kd is the rate constant of monomer detachment from the (n+1)-sized clusters. This formula suggests that the temporal evolution of the cluster concentrations can be described using rate equations once the rate constants ka and kd are known.

Governing Equations

On the basis of eq , we can formulate a system of KREs for a closed system, each characterizing the rate at which the concentration of the n-sized clusters evolves upon molecule attachments and detachments:[19]where Z(t) is the concentration of the n-sized clusters (i.e., their number per unit suspension volume) and Z(t) ≡ 0, with nmax being the upper bound of the cluster size, which must be sufficiently large to avoid any effect on the numerical solution.[22]Equation is the mass balance for the solute, ensuring that the total number of solute molecules is conserved. The initial condition is as followswhere Z is the initial cluster size distribution. The system of eqs and 2b can be solved with the initial condition (eq ) to obtain the temporal evolution of any initial cluster size distribution for a given set of rate constants (possibly evolving in time, for instance, due to temperature changes).

Attachment and Detachment Rate Constants

The attachment rate constant ka is determined by the rate-limiting mechanism of the monomer transport to the clusters. It is common that the monomer transfer is limited by the monomer integration step on a cluster surface, for example, in the crystallization of l-glutamic acid in an aqueous solution,[23] lactose in an aqueous solution,[24] and paracetamol in an ethanol solution.[25] Under this surface-integration limited mechanism, ka is given by[19,26]where α is the sticking coefficient, ks is the surface area shape factor (e.g., ks = π for a sphere), kv is the volume shape factor (e.g., kv = π/6 for a sphere), V1 is the volume of a solute molecule, and D is the diffusion coefficient of the solute molecule in the solution. Equation shows that ka is a purely kinetic quantity,[19] as it does not contain any thermodynamic parameter. The detachment rate constant kd can be derived based on thermodynamic considerations about the concentration of clusters.[19] For instance, at equilibrium, the monomer attachments to the n-sized clusters occur with the same frequency as the monomer detachments from (n+1)-sized clusters, thus kinetically ensuring the microscopic balance of the pseudo reaction (eq ); this implieswhere C is the equilibrium concentration of n-sized clusters at the actual level of Δμ, that is, the thermodynamic driving force for crystallization;[19] therefore, we write C(Δμ). The equilibrium cluster concentration C can be expressed as[19]where C0 is the concentration of nucleation sites in the system, ΔG(Δμ) is the Gibbs free energy for the formation of a n-sized cluster from n monomers at a given driving force Δμ, kB is Boltzmann constant, and T is the absolute temperature of the system. The nucleation site concentration C0 is typically defined as V1–1 in the literature,[19,26] but a different value of C0 calculated from a mass balance was used in this work to consider the bulk solubility ce of a system, as explained in Appendix A. It is worth highlighting that C0 is cancelled out in the nondimensional Szilard model (see eqs –A.8 with Est = 1), so the general behavior of the model is not affected by the choice of a specific value of C0. We further note that the probability C/C0 of finding a n-sized cluster at a given nucleation site is exponentially proportional to −ΔG/kBT,[19] thus indicating that the equilibrium cluster size distribution C is a Boltzmann-type distribution. The driving force Δμ is defined as[23]where C1,e≡ C1(0) is the monomer concentration at phase equilibrium (Δμ = 0) and s = Z1/C1,e is the supersaturation (ratio) based on the actual monomer concentrations. According to the CNT, ΔG is given by[19]where γ is the specific surface energy of a cluster and b = ks(V1/kv)2/3 is the surface area of a solute molecule. Combining eqs and 6, with a consideration that C1 is very close or equal to Z1,[19] leads toor with the help of eqs and 8where k0a = 0 and k1d = 0 by definition.[19] Note that eq shows that kd is determined by both kinetics (i.e., ka) and thermodynamics (i.e., ΔG).[19]

Modeling SNIPE

In this section, we extend the Szilard model to describe the SNIPE mechanism by incorporating the assumption that interparticle interactions stabilize cluster formation. This stabilization effect alters the thermodynamics of nucleation (i.e., ΔG) near seed crystals such that the effective detachment rate constants kd,eff(t) become smaller than those without stabilization, thus accelerating nucleation. For implementation, the identical governing equations of the Szilard model (eqs and 2b) can be reused after replacing the original detachment rate constants kd with kd,eff whose formulation is described in the following by considering that only part of a cluster population close to seed crystals gets stabilized by the interparticle energies.

Interparticle Energies

In part I of this series,[27] the thermodynamic aspects of the SNIPE mechanism are examined in detail by considering two classical interparticle forces, that is, van der Waals forces and Born repulsion forces. This contribution, as part II of this series, examines the kinetic aspect of the SNIPE mechanism by including the role of these interparticle interactions in the Szilard model under the following simplifying assumptions:

As illustrated in Scheme a, a seed crystal is surrounded by a stabilization volume Vst(t) consisting of a thin shell with thickness d1lst, where d1 is the characteristic length of a solute molecule and lst(≥0) is an adjustable nondimensional parameter.

Scheme 1

Schematic representation of a seed crystal in a solution when a crystal is spherically shaped; not drawn to scale

(a) A seed crystal consisting of n solute molecules is surrounded by the stabilization volume, consisting of a thin shell with thickness d1lst, where d1 is the characteristic length of a solute molecule and lst is an adjustable parameter. (b) A suspension volume consists of a solid phase volume Vs, a stabilization volume Vst, and a bulk solution volume Vb. Cluster formation occurring in Vb is governed by the Gibbs free energy ΔGb, whereas the cluster formation in Vst is controlled by the Gibbs free energy ΔGst in that volume. The concentrations Zb and Zst represent the cluster concentrations in Vb and Vst, respectively. A flux of molecular clusters fst→b from Vst to Vb through the interfacial area A is represented by an arrow.

A uniform field of interparticle energies is applied within Vst, which vanishes outside Vst; this external field stabilizes the cluster formation by decreasing ΔG.

Under these assumptions, we can define the Gibbs free energy ΔGst(Δμ) for the formation of a n-sized cluster in Vst at a given Δμ as follows (see also eq in part I of this series)where Est(≥1) is a nondimensional parameter that accounts for the intensity of the interparticle energies. By substituting ΔG in eq with ΔGst, the detachment rate constant kd,st in Vst can be written asHere, kd,b(=kd) is the detachment rate constant in the bulk solution volume Vb(t), with the corresponding Gibbs free energy ΔGb(=ΔG). Equation shows that when the interparticle energies are present (i.e., Est > 1), kd,st is smaller than kd,b, thus implying that clusters grow faster in Vst than they do in Vb. Obviously, when the interparticle energies are absent (i.e., Est = 1), kd,st is equivalent to kd, thus yielding no stabilization effect. To consider different kinetics in Vb and Vst and to describe the evolution of an entire cluster population, the system can be divided into multiple compartments, each described by a mass balance with a corresponding kinetic description. In this work, the system/suspension volume Vsusp is partitioned into three compartments, as illustrated in Scheme b: the bulk solution volume Vb, the stabilization volume Vst, and the volume of the solid phase Vs, with corresponding time-dependent volume fractions defined aswith ϕb(t) + ϕst(t) + ϕs(t) = 1. The solid volume fraction ϕs can be calculated by considering that the solid phase consists of supercritical clusters[19]where V1n approximates the volume of a n-sized cluster and n* = (2γb/3Δμ)3 is the critical size.[19] Likewise, ϕst can be calculated by considering all stabilization (shell) volumes in the suspensionwhere nst is the minimum cluster size that enables the existence of the stabilization volume and d is the characteristic length of a n-sized cluster.

For simplicity, nst is set to be the minimum size of the initial seed crystals, being on the order of 1013 in the present work (corresponding to 25–50 μm for small organic molecules). The approximation in eq is valid, provided that the size of seed crystals (i.e., d = d1n1/3) is much larger than the thickness of the stabilization volume d1lst, that is, n1/3 ≫ lst, where lst is varied from 1 to 100 in this study. Note that eq shows that ϕst is linearly proportional to the parameter lst and the total surface area of seed crystals. The concentrations of n-sized clusters in Vb and Vst, denoted by Zb and Zst, respectively, can be determined by considering that the total number of n-sized clusters is conserved, that is,and by assuming that their ratio Zb/Zst is well approximated by that at equilibrium, Cb/Cst, given bywhere Cst(Δμ) and Cb(Δμ), likewise in eq , are the equilibrium concentrations of the n-sized clusters in Vst and Vb, respectively. Combining eqs , 11, 13, 16, and 17 yields

The mass balances for Zb and Zst can be written aswhere fst→b is the flux of the n-sized clusters from Vst to Vb through the interfacial area A, illustrated by an arrow in Scheme b, and the monomer concentration in each volume is approximated by Z1, as in part I of this series.[27] Summing up these two equations and using eqs , 13, 16, 18, and 19 lead to a system of governing equations for Z (n = 2, 3, ..., nmax), analogous to eq 2a, that can be cast as followswhere the effective detachment rate constant kd,eff isEquation implicitly assumes that all clusters including the supercritical (n > n*) clusters experience the stabilization effect; it has been used in all simulations reported in this paper. This is because the influence of kd,eff on the simulation results Z(t) is negligible compared to that of kd,eff, and setting kd,eff ≈ kd increases the computational time by a factor of 20. Summarizing, a system of the governing equations (eqs and 22) can be solved with the initial condition (eq ) for simulating a crystallization process in the presence of the stabilization effect that is described by eqs , 10, 13–15, and 23. The developed model contains 11 parameters: αs, ks, kv, V1, D, T, γ, C1,e(or ce), Z, Est, and lst. As elaborated in Appendix B, the model can be made dimensionless by introducing two dimensionless variableswhere Y and τ are the nondimensional concentration and time, respectively. Consequently, the nondimensional model contains five dimensionless parameters: β, Ω, Y = Z/C1,e, Est, and lst. Here, β represents the volume fraction of monomers at phase equilibrium, defined asand Ω is the dimensionless surface energy, defined as

The nondimensional model has been solved numerically as explained in Appendix C.

Results

To demonstrate the SNIPE mechanism and the application of the developed model to experimental observations in the literature, seeded batch crystallization processes have been simulated at a constant temperature with varying initial supersaturation, initial seed size distribution, strength, and effective range of the stabilization effect. For simulating such conditions, two parameters (Ω and β) have been kept constant, while the other three [Est, lst, and Y0(n)(≡Y)] have been varied. The values of these five parameters have been specified by adopting the physicochemical properties of a reference case, that is, paracetamol crystallization from an ethanol solution at 20 °C; the physicochemical properties and the corresponding values of the model parameters are reported in Tables and 2, respectively.

Table 1

Physicochemical Properties of the Paracetamol–Ethanol System at 20 °C

properties	value	description
molecular weight of the solute, Mw	0.151 kg mol–1
solid density, ρs	1263 kg m–3
solvent density, ρsolv	789 kg m–3
mass-based bulk solubility, c̅e	0.180 kg kg–1	from ref (28).a
bulk solubility, ce	5.65 × 1026 m–3	≡c̅eNAρsolv/Mwb
specific surface energy, γ	0.0111 J m–2	assumed value.c
molecular volume of the solute, V1	1.99 × 10–28 m3	≡Mw/ρsNAb
molecular surface area of the solute, b	1.65 × 10–18 m2	≡ks(V1/kv)2/3d
monomer solubility, C1,e	4.82 × 1026 m–3	≡C1(0)

The solubility is based on solute mass per solvent mass.

NA is the Avogadro number.

The specific surface energy is within the range from an experimentally fitted value (0.0043 J m–2 at 40–46 °C by ref (29)) to a theoretically predicted value (0.0134 J m–2 at 20 °C by eq 12 in ref (26)).

Spherical shape is assumed: ks = π and kv = π/6.

Table 2

Model Parameters Used in the Simulations

symbol	value			description
β	0.0958			monomer volume fraction at phase equilibrium
Ω	4.5			nondimensional surface energy
Est	1, 1.1, 1.2			strength of the stabilization effect
lst	0, 1, 10, 100			effective range of the stabilization effect
Y0(n)	n = 1	1.52–2.34		initial supersaturation s0
	2 ≤ n ≤ n*	0		initial subcritical and critical clusters
	n* < n	ϕs,0	0.0016	initial solid volume fraction
		q3(L)a	normal(50, 0.5) or log–normal(5.14, 0.51)	normal dist. in Sections 4.1 and 4.2; log–normal dist. in Section 4.3

Volume-weighted probability density function of the crystal size in micrometer .

The parameters β and Ω have been calculated using eqs and 26, respectively, with the physicochemical properties in Table : β = 0.0958 and Ω = 4.5. The parameter Est = {1, 1.1, 1.2} reflects the strength of interparticle interactions, and its values are conservative estimates of Est, ranging between 1.04 and 1.26 (see Section in part I of this series[27]). The parameter lst represents the range of the stabilization effect. A wide range of {0, 1, 10, 100} has been investigated; a realistic value has been estimated to be around 4 in part I.[27] By setting either Est = 1 (zero strength) or lst = 0 (zero thickness) or both, the stabilization effect can be eliminated, thus reducing the new SNIPE model to the Szilard model.

The last parameter Y0(n) determines the initial condition: the initial supersaturation s0 via Y0(n = 1), the initial subcritical and critical cluster concentrations via Y0(2 ≤ n ≤ n*), and the initial seed size distribution via Y0(n > n*). The initial supersaturation s0 is equivalent to Y0(n = 1) (see eqs and 24b), and the studied range is s0 = 1.52–2.34. For the subcritical and critical clusters, it is typical to assume that they are initially absent,[21,22] thus setting Y0(2 ≤ n ≤ n*) = 0. To define the initial seed crystal distribution Y0(n > n*), it is required to set the initial amount of seeds and the initial probability density function (PDF) of the crystal size. The initial amount of seeds has been kept constant by fixing the initial solid volume fraction ϕs,0 = 0.0016. This condition corresponds to 1 g of seeds in a 500 mL solution, a reported experimental condition that enables secondary nucleation in the paracetamol–ethanol system at 20 °C.[30] For the initial PDF, two cases are considered. In Sections and 4.2, we use a narrow Gaussian distribution (characterized by the mean size 50 μm and the standard deviation 0.5 μm) to discuss the general behavior of the developed model. In Section , to compare our simulation results with experimental observations qualitatively, we use a log–normal distribution that adequately characterizes an experimental seed size distribution (explained in detail in Section ).

The simulation results are presented in three parts. In Section , we demonstrate the key features of the model in describing the time evolution of a cluster size distribution, including the formation of primary and secondary nuclei, thus leading to a bimodal crystal size distribution. Moreover, we verify that the stabilization effect can indeed cause secondary nucleation. In Section , we present the results of a sensitivity analysis that elucidates the effect of the simulation parameters on the final crystal size distribution. Finally, in Section , we make a qualitative comparison between simulation results obtained using the SNIPE model and experimental observations reported in the literature.

Demonstration of the Stabilization Effect

The developed model can provide a comprehensive description of batch crystallization processes in which the growth of seeds can occur together with secondary nucleation induced by the interparticle energies. These processes are simulated in the model such that the initial cluster size distribution evolves over time with or without the stabilization effect while experiencing a series of molecule attachments (i.e., growth) and detachments (i.e., dissolution). The model yields the time evolution of a cluster size distribution Y(τ) that reveals insightful dynamics of the three key populations in the crystallization processes, namely, the populations of subcritical clusters, nuclei, and seed crystals. We emphasize that this type of insightful information cannot be attained using a (classical) population balance model.[21]

In this section, we demonstrate how the stabilization effect can cause secondary nucleation by examining three simulation results. Two results are obtained from the simulations with no stabilization effect (i.e., Est = 1 and lst = 0) at two levels of the initial supersaturation s0: high level s0 = 2.22 (Figure a) and low level s0 = 1.93 (Figure b). The other simulation result is obtained when including the stabilization effect (corresponding to Est = 1.2 and lst = 10) and at the low supersaturation s0 = 1.93 (Figure c). All three simulations are illustrated in Figure in terms of the volume-weighted size distributions nY (solid lines) and the critical sizes n* (dashed lines) at three different times. Note that nY is called volume-weighted distribution in the literature,[21] though the distribution is weighed by the number of molecules because it is assumed quite rightly that the volume of a n-sized cluster is proportional to the number of molecules that comprise the cluster.[19]

Figure 1

Comparison of three simulations conducted with two levels of the initial supersaturation s0 and in the absence (a,b) or the presence (c) of the stabilization effect. For each simulation, the volume-weighted size distributions (nY(τ) vs n) at three times τ are depicted by the solid lines, and the corresponding critical sizes (n*(τ)) are depicted by the dashed lines. In all relevant subplots, the point where the solid line intersects the dashed vertical line corresponds to the critical clusters, and its vertical coordinate is an indication of their concentration. A set of parameters for each simulation is as follows: (a) Est = 1; lst = 0; s0 = 2.22. (b) Est = 1; lst = 0; s0 = 1.93. (c) Est = 1.2; lst = 10; s0 = 1.93.

The simulation without the stabilization effect at high supersaturation is analyzed first. The earliest snapshot [Figure (a1)] reveals the fast dynamics of a subcritical cluster population; a large number of subcritical clusters form rapidly. A small part of this population grows larger than the critical size, thus forming a population of primary nuclei through homogeneous nucleation. On the right-hand side of the figure, one can observe a seed crystal population in the form of a narrow Gaussian distribution. In the next snapshot [Figure (a2)], a growing nucleus population dramatically changes the cluster size distribution by increasing both the number and size of the nuclei. The third snapshot [Figure (a3)] is mainly characterized by a bimodal crystal size distribution consisting of both the seeds and the population of primary nuclei. All snapshots at the bottom row of Figure represent the results at the final process times τf at which the supersaturation is close to unity and the number of crystals is at its maximum (see Section S2 in the Supporting Information for more details).

The next simulation, illustrated in Figure b, is performed with no stabilization effect and a low level of the initial supersaturation s0 = 1.93. In this simulation, a nucleus population does not form (see Figure (b1,b2)); consequently, the seed crystals completely dominate the crystallization process through their growth, thus resulting in a unimodal crystal size distribution, as presented in Figure (b3). A comparison of the two simulations discussed so far (Figure a,b) indicates that the low initial supersaturation s0 = 1.93 is insufficient to trigger primary nucleation in the bulk solution for the system under consideration here.The third simulation is illustrated in Figure c. Although this simulation begins with the low initial supersaturation, a substantial amount of nuclei appear [see Figure (c1,c2)] because of the stabilization effect. This stabilization effect enhances the overall driving force for crystallization and thus increases the concentration of the critical clusters by several orders of magnitude, as shown in Figure (c1,c2). Their subsequent growth and that of the seeds yield a bimodal crystal size distribution [see Figure (c3)]. The observed nucleation can be regarded as secondary nucleation caused by the stabilization effect, in other words, SNIPE. This result, obtained using the SNIPE model developed here, can explain why secondary nucleation can occur even when the initial supersaturation is at such a low level that primary nucleation does not occur.

Sensitivity Analysis

A sensitivity analysis has been conducted in two steps. First, various simulations have been carried out by altering three model parameters: the intensity of the stabilization effect Est = {1, 1.1, 1.2}, the effective range of the stabilization effect lst = {0, 1, 10, 100}, and the initial supersaturation s0 = {1.52, 1.64, 1.76, 1.87, 1.93, 1.99, 2.11, 2.17, 2.22, 2.34}. Second, simulation results have been compared with respect to two characteristic quantities: the degree of nucleation, ψN, which indicates how many new crystals have formed with respect to their initial count, and the degree of growth, ψG, which specifies how much seed crystals have grown with respect to their initial size.

Degree of Nucleation and Growth

These two characteristic quantities can be readily calculated from a cluster size distribution, Y(τ), by using its moments. Hence, we first express the j-th moment m(p,τ) of the distribution Y(τ) in which the considered cluster size n is larger than a given size p

Based on this definition, the degree of nucleation ψN can be defined by dividing the zeroth moment of the final crystal size distribution, m0,f*, by the initial one, m0,i*where n*(τ) is the critical size at time τ. This equation gives ψN ≈ 1 when nucleation is negligible, simply because the final number of crystals is close to the initial one (i.e., m0,f* ≈ m0,i*). On the contrary, this definition yields ψN ≫ 1 when nucleation occurs because the final number of crystals is much larger than the initial one (i.e., m0,f* ≫ m0,i*). Similarly, the degree of growth ψG can be defined by dividing the final volume-weighted average size of the seed crystals, n̅seed,f, by the initial one, n̅seed,iwhere nseed,min(τ) is the smallest size of a seed crystal size distribution at time τ. This equation gives ψG ≈ 1 when the growth of seed crystals is negligible (n̅seed,f ≈ n̅seed,i), whereas it gives ψG > 1 when the seeds grow noticeably (n̅seed, f > n̅seed,i).

Results

The effect of the initial supersaturation s0 on ψN and ψG in the absence of the stabilization effect (Est = 1 and lst = 10) is illustrated in Figure by the circle markers. When the initial supersaturation s0 is below the threshold value sth(Est = 1) = 2.17 (the left side of the vertical dotted line in Figure ), ψN remains around 1 (see Figure a), while ψG increases linearly with s0 (see Figure b). By increasing s0 beyond the threshold value sth(Est = 1), ψN increases by orders of magnitude, whereas ψG drops to 1. This effect of s0 on ψN and ψG is also present in the simulations with the stabilization effect (Est = {1.1, 1.2}) but with the lower threshold values, sth(Est = 1.1) = 2.10 and sth(Est = 1.2) = 1.87, as illustrated in Figure by the square markers for Est = 1.1 and the triangle markers for Est = 1.2. In all three cases, if s0 < sth, the growth of the seeds is the only dominant crystallization mechanism, indicated by ψN ≈ 1 and ψG > 1. By increasing s0 beyond sth, the dominant crystallization mechanism continuously shifts from growth to primary nucleation (in the case of Est = 1) or secondary nucleation (in the case of Est > 1), indicated by ψN ≫ 1 and ψG → 1. Furthermore, the threshold initial supersaturation sth decreases with increasing strength of the stabilization effect Est. This suggests that secondary nucleation by interparticle energies can occur with a smaller value of s0 if a system exhibits strong interparticle interactions between seeds and molecular clusters.

Figure 2

Degree of nucleation ψN is shown in plot (a) and degree of growth ψG is presented in plot (b) corresponding to the simulations conducted with varying values of the initial supersaturation s0 and the stabilization effect strength Est and a fixed value of the stabilization effect range lst = 10. A vertical dotted line indicates the threshold initial supersaturation for Est = 1: sth(Est = 1). Other lines are just a guide for the eye.

Next, we shift our attention to another group of simulations for studying the influence of s0 and lst = {1, 10, 100} at a fixed value of Est = 1.2. In the following, the degree of nucleation ψN alone is sufficient to demonstrate the impact of the model parameters. A higher value of lst means a longer effective range of the interparticle energies. As illustrated in Figure , although we increase lst by 1 or 2 orders of magnitude, the threshold supersaturation barely decreases. This minor effect of lst on the degree of nucleation ψst sharply contrasts with the major effect of Est.

Figure 3

Degree of nucleation ψN of the simulations conducted with varying values of the initial supersaturation s0 and the stabilization effect range lst and a fixed value of the stabilization effect strength Est = 1.2. Lines are just a guide for the eye.

To further understand the relative importance of these two parameters, a set of simulations have been executed with varying values of Est and lst at a constant value of s0 = 2.22, as illustrated in Figure . This figure clearly shows that the degree of nucleation ψN changes more dramatically when we change Est than when lst is varied. This result can be explained by considering the functional form of eqs and 23, where Est determines the effective detachment rate constant kd,eff more dominantly than lst does.

Figure 4

Degree of nucleation ψN of the simulations conducted with varying values of the stabilization effect strength Est and the stabilization effect range lst and a fixed value of the initial supersaturation s0 = 2.22.

Comparison with Experimental Observations: Paracetamol Crystallization from an Ethanol Solution

In this section, we show the relevance of SNIPE-based simulations to experimental crystallization processes by comparing our simulation results with available experimental observations qualitatively. For the comparison, simulation results [i.e., Y(τ)] were converted into experimentally relevant quantities, namely, the bulk supersaturation S and the volume-weighted PDF of the crystal size q3(L), as explained in Appendix A and Appendix D, respectively. The conditions of the reference experiments were set up in the model as detailed in the following.

Setting up Experimental Conditions in the Model

References (30) and (31) are selected as benchmark experimental studies for two main reasons. One reason is that these experiments were conducted with the paracetamol–ethanol system, a widely studied system from various perspectives.[25,28,30−34] Another reason is that the chosen works were carried out systematically as a series; prior to their investigation of secondary nucleation,[30] the authors extensively studied primary nucleation[33,34] and growth[31] for the same system, thus offering various types of experimental data sets obtained by the same researchers in a consistent manner.

In the following, we summarize the experimental conditions in refs (30) and (31) and explain how these conditions are reflected in the model parameters summarized in Table . The reference experiments consist of seeded batch crystallization of paracetamol from an ethanol solution at 20 °C; these conditions were implemented in the model by setting β = 0.0958 and Ω = 4.5, calculated using eqs and 26, respectively, with the properties reported in Table . The experiments are started by seeding 1 g of crystals in a 500 mL supersaturated solution; the corresponding seed size distribution is shown in Figure in terms of the volume-weighted PDF q3(L) of the crystal size L; these conditions were reproduced in the model by setting the initial solid volume fraction ϕs,0 = 0.0016 and by approximating the experimental seed size distribution with a log–normal distribution (see the dash-dotted line in Figure ), characterized by the mean size 5.14 and standard deviation 0.51 of ln L.

Figure 5

Experimental data of the volume-weighted PDF of the initial crystal size (q3(L) vs L) shown by a bar graph and a fitted log–normal distribution (characterized by the mean 5.14 and standard deviation 0.51 of ln L) by a dash-dotted line. The displayed experimental data are adapted with permission from Figure 5c in ref (31). Copyright Elsevier (2011).

Comparison

Two sets of simulations were performed, namely, either in the absence or presence of the stabilization effect (Est = 1 and lst = 0 and Est = 1.2 and lst = 10, respectively), each for estimating a threshold initial bulk supersaturation for primary nucleation and secondary nucleation, respectively. These estimates were compared with the values inferred from the experiments.

The results of the first set are presented in Figure a in terms of q3(L). When S0 is larger than around 1.8, the simulation yields a bimodal size distribution, thus indicating the occurrence of primary nucleation. This suggests that the threshold initial bulk supersaturation is around Sth(Est = 1) ≈ 1.8. It is important to note that Sth(Est = 1) varied around ±5% when the specific surface energy γ in Table was varied around ±10%. A corresponding experimental value for the threshold supersaturation can be inferred from the experimental data in ref (31) that characterized the metastable zone width, reproduced in Figure c. This figure shows that, at 20 °C, primary nucleation was detected when the initial bulk supersaturation was larger than 1.84 (corresponding to the point B in the figure), which is close to the value estimated by the model. This comparison demonstrates that the model can estimate a threshold supersaturation for primary nucleation adequately.

Figure 6

(a) Simulations conducted with no stabilization effect (Est = 1 and lst = 0) at two values of the initial bulk supersaturation S0 = {1.78, 1.83}. (b) Simulations performed with the stabilization effect (given by Est = 1.2 and lst = 10) at two values of the initial bulk supersaturation S0 = {1.5, 1.6}. (c) Solubility and metastable zone width data. [Reproduced with permission from ref (31). Copyright Elsevier (2011)]. (d) SEM image of the final crystals showing agglomerates. [Reprinted with permission from ref (30). Copyright Elsevier (2012.)].

The outcome of the second simulation set is illustrated in Figure b, indicating that the threshold supersaturation for secondary nucleation falls in the range between 1.5 and 1.6. This range agrees well with the experimental observations in ref (30) where a value of the initial supersaturation in the range S0 = 1.42–1.71 indeed caused secondary nucleation. This again shows that the developed model could be used to support future experiments, for instance, when selecting the initial supersaturation for a crystallization process.

We do acknowledge that the final PDFs in our simulations may not quantitatively match the experimental ones because of two main reasons. One reason is that, in experiments, secondary nucleation is caused by a combination of multiple mechanisms (e.g., attrition and SNIPE), but the developed model considers only the SNIPE mechanism. The second reason is that the final PDF can be critically modified by other crystallization phenomena (e.g., crystal breakage[35−37] and agglomeration[38,39]), and these phenomena are not always predictable. In the reference experiment (ref (30)), agglomeration clearly occurred, as evident from the reported SEM image of the final crystal population (see Figure d). Since the developed model describes neither agglomeration nor breakage, the simulation outcome cannot quantitatively agree with the experiments where any of these phenomena may occur. Though modeling crystal agglomeration and/or crystal breakage in a KRE model is possible, this is nevertheless beyond the scope of this work.

Conclusions

We have developed a mathematical model that simultaneously describes growth, homogeneous nucleation, and secondary nucleation caused by the interparticle energies, that is, due to the stabilization effect caused by the interactions between seeds and clusters in suspension. The key feature of the model is the incorporation of the stabilization effect into the conventional KRE model of nucleation. By simulating seeded batch crystallization of paracetamol from an ethanol solution at 20 °C, this work demonstrates how the stabilization effect causes secondary nucleation at a low initial supersaturation by increasing the concentration of the critical clusters by orders of magnitude. Moreover, a sensitivity analysis has been carried out to understand the relative importance of three key simulation parameters: the initial supersaturation, s0, the strength of the stabilization effect, Est, and the effective range thereof, lst. The analysis shows that primary and secondary nucleation are enhanced to varying degrees on increasing any of the three parameters, and the parameter Est has a major effect on the stabilization effect, while the parameter lst has a minor effect. The simulation results are compared with the experimental observations in the literature qualitatively. The comparison verified a capability of the model to estimate a threshold initial bulk supersaturation at which primary or secondary nucleation occurs as the estimated values adequately agree with the values inferred from the reference experiments.

The model can be employed to obtain qualitative information on threshold supersaturation for nucleation, which is a useful tool for designing crystallization processes. Moreover, the developed model presents a new perspective on secondary nucleation mechanisms through which one can rationalize some experimental observations that cannot be explained by the attrition mechanism, for instance, the formation of secondary nuclei whose polymorphic form or chiral structure is different from that of a seed,[14−17] the occurrence of secondary nucleation caused by a seed crystal that was held stationary,[11−13] and the appearance of secondary nuclei in a no-impeller crystallizer (e.g., airlift crystallizer[40,41]).

Finally, demonstrating that the theory developed in part I and II of this series has potential for application and experimental validation remains a challenge. To address this, the complex KRE model used here needs to be transformed into a nucleation rate model, which can be used in a standard population balance modeling framework. This would prove that the SNIPE theory can describe experiments and would allow for a comparison of the SNIPE rate model with other well-recognized secondary nucleation models. These aspects will be addressed in part III of this series,[42] where the practical relevance of secondary nucleation rate by interparticle energies will be proved.