Idealized Mathematical Model for Pharmaceutical Powder

Article pubs.acs.org/IECR

Idealized Mathematical Model for Pharmaceutical Powder “Micronization” Using Compressed Gas Antisolvent (Re‑)Precipitation (GASP): Predicted Performance for the Model Ternary System: Phenanthrene/Toluene/CO2 Daniel E. Rosner*,†,‡ and Manuel Arias-Zugasti†,¶ †

Department of Chemical & Environmental Engineering, Sol Reaction Engineering Group, Yale University, New Haven, Connecticut 06520-8286, United States ABSTRACT: We consider pharma-motivated processes in which compressed gas (here, CO2) is employed as an “anti-solvent” (AS), causing drug particle (re)precipitation from an injected presaturated organic solvent microdroplet assumed to be molecularly well-mixed at each instant. It is shown that conditions leading to appreciable homogeneous particle nucleation in this “expanded supersaturated liquid” environment can be sufficiently short-lived that they preclude appreciable particle growth because of encounters with solute molecules or other precipitate particles. This leads to relatively narrow predicted precipitated particle size distributions (PSDs), with a characteristic particle size determined by the typical critical nucleation size. Our mathematical techniques, which exploit the method of characteristics (MOC), make no presumption about PSD-shape and our numerical simulations employ realistic thermophysical properties for the previously studied model system: phenanthrene (surrogate “active pharmaceutical ingredient” [solute], toluene [solvent, S], and compressed CO2 [AS]). We explicitly consider initial solvent droplet diameters in the micrometer range and CO2 pressures of 52−64 bar and demonstrate the necessity of relaxing frequently made approximations (e.g., nucleation with constant surface energy, negligible Gibbs−Kelvin−Ostwald solubility corrections, solute diluteness, ...). Our new methods/parametrizations/performance results for this mathematical model should already help select optimal GASP/SASP operating conditions, and, with suitable extensions, ultimately lead to the development of nearly equally tractable yet sufficiently complete process models.

1. INTRODUCTION 1.1. Motivation. It is a pleasure and honor for us to contribute to this special issue of Industrial & Engineering Chemistry Research acknowledging the impact that Prof. Doraiswamy Ramkrishna (i.e., “Ramki”) continues to have on our profession in general and on our ability to mathematically model particle processes in particular. As will be illustrated in this paper, the generality of the “population balance” (PB) approach expounded in Ramki’s timely monograph1 enables its instructive extension/application to pharma processing situations such as the (re-)precipitation of drugs from compressed CO2-expanded organic solventsi.e., so-called GASP processes. We show that this class of applications is of special interest/importance, because of the many novel and rather fundamental PB aspects of the highly nonideal GASP application/environmentincluding the need for (i) systematic modifications of the solute molecular impingement frequency for both nucleation and growth kinetics; (ii) what we have termed “surface energy evolution” (SEE); (iii) highly variable mixture volume due to dilation; (iv) high particle mass loading; and (v) the need for a tractable process model capable of tracking/anticipating the resulting API powder size distribution functionsso important in determining powdered drug availability when administered to humans (orally, by injection, or inhalation). Our present (single, well-mixed solvent droplet) analysis (also see our earlier work2) will be seen to exploit the mathematical method of characteristics (MOC) to track the © 2015 American Chemical Society

evolving particle number density distribution function (NDDF) via particle growth along individual particle “trajectories” (sketched in Figure 1 and described more thoroughly in Section 4.3) using the particle volume vs time plane. While our present analysis is demonstrated to be quite instructive and selfconsistent for sufficiently small initial solvent droplet diameters

Figure 1. Qualitative behavior of the characteristics curves “emerging” from the descending v*(t) locus on the v, t-plane. (See Section 4.3 for relevant discussion.) Special Issue: Doraiswami Ramkrishna Festschrift Received: Revised: Accepted: Published: 10383

March 28, 2015 July 9, 2015 July 10, 2015 July 10, 2015 DOI: 10.1021/acs.iecr.5b01155 Ind. Eng. Chem. Res. 2015, 54, 10383−10400

Article

Industrial & Engineering Chemistry Research

here is to identify (Section 4) and address several of them, first under almost-isothermal conditions. 1.4. Reaction Engineering Approach to GASP Processes and Relevant Characteristic Time Ratios. Specifically, we adopt a “sol reaction engineering” perspective and focus attention here on tractable situations governed by the following three characteristic times: • tAS, which is associated with the rate of addition of antisolvent (CO2), and, hence the associated rates of dilation and supersaturation growth; • tN, which describes the homogeneous nucleation (N) (or “birth rate”) of new API particles; and • tG, which is defined by the growth rate of particles due to solute molecule attachment kinetics in the prevailing supersaturated “expanded” solvent mixture. Thus, for the present, we explicitly consider the instructive “well-mixed droplet” (WMD) limit (i.e., we deliberately postpone dealing with reality R3)focusing instead on the predicted effects of homogeneous nucleation (N) and (to a lesser extent) interface-curvature-dependent molecular growth (G) in the absence of appreciable particle−particle coagulation and/or deliberate8 or self-induced nonisothermal effects. This will be seen to enable formally exact numerical evolution solutions of the resulting particle population balance equation for N/G rate laws containing realistic features not previously taken into account in the field of GASP/SASP process modeling.7 For homogeneous nucleation kinetics, this includes the consequences of an environment-dependent (hence, variable) ef fective surface energy.2 For growth, we have introduced a growth species incorporation probability9 with a nonlinear thermodynamic activity-based driving force, along with self-consistent environment- and particle-size-dependent interface curvature (Gibbs−Kelvin) “corrections”. In all cases, we make no presumption about the shape of the resulting precipitate PSDsindeed, as will be seen (Section 5), our present methods often lead to unusual/complex PSD shapes. The present class of three-phase model systems is shown to be dictated by corresponding (Damköhler-type) ratios of the above-mentioned characteristic times (see, e.g., Chapter 7 in ref 10 ), as illustrated in Section 4 using typical property values for the specific but representative ternary model system: phenanthrene (surrogate “API” or solute), toluene (solvent, S), and CO2 (AS). By deliberately comparing output PSD characteristics with those predicted by neglecting surface energy evolution and Gibbs-Kelvin growth rate corrections, we highlight (in Section 4) the importance/consequences of systematically including these phenomena in modeling GASP processes. We are also now able to address (Section 5) the PSD consequences of solvent spray polydispersityperhaps for the first time. In Section 6, we summarize our principal conclusions to date and comment on the above-mentioned additional fundamental phenomena (reality R3) that now remain to be examined in an analogous fashion. Indeed, tractable yet hopefully instructive theoretical studies, initiated to examine the effects of CO2 diffusion limitations in “large droplets” (Section 6.5), will be reported in subsequent papers, based on this same ternary model system.

(Section 3.2), it will not be the “last word” in developing a much-needed GASP process model, because of the likely importance of intradroplet CO2 diffusion limitations for solvent droplets in the (supermicrometer) size range of maximum practical interest. Our approach to dealing with this remaining generalization, while retaining much of the simplicity of our present well-mixed droplet process model, is outlined in Section 6.5which completes this paper. 1.2. Relevant Earlier GASP Process Studies (e.g., Coagulation Effects). Before proceeding, it should be mentioned here that this unusual pharma processing environment has also motivated some of our earlier PB-oriented worki.e., on the dilation effect on coagulation rate constants3 and the role of local temperature gradients in the Brownian coagulation of particle populations.4−6 Here, however, we focus on the almost-isothermal CO2 uptake stage responsible for the birth of the solid phase from the supersaturated solutecontaining “expanded” solvent. Even growth due to monomer impingement is found to be of secondary importance during the impingement rate-controlled CO2 antisolvent uptake stage (Section 3.1). For this reason, we expect growth due to Brownian coagulation-coalescence to be important only in postCO2 uptake stages, which are not considered further here. In bench-scale experiments, near-critical (compressed) CO2 has been shown to be a convenient antisolvent to reprecipitate valuable organic compounds, including active pharmaceutical ingredients (“APIs”). Because this strategy can lead to relatively solvent-free submicrometer powders under otherwise mild processing conditions (thermally, mechanically), there has beenand remainsa considerable incentive to develop tractable gas antisolvent precipitation (“GASP”) process models (see ref 7 for a review) that could be used to economically select, optimize, control, and/or scale-up such batch, semibatch, or continuous processes. Our progress toward this goal is summarized here. 1.3. Summary of Rate Process in Spray-Based GASP. To fix ideas, if one considers, for example, injecting a spray of liquid solvent (S) containing a dissolved API into compressed antisolvent (AS) and one’s goal is to predict not only the precipitate yield but also the output solid particle size distribution (PSD), then it becomes clear that a comprehensive process model should at least deal with the following “realities”: R1: Nonideal, at least three-component, three-phase chemical thermodynamics and the following rate processes; R2: Molecular impingement-controlled AS uptake by the resulting droplets; R3: Intradroplet molecular mixing (dif f usion, ...) of antisolvent R4: Solid particle “birth” (nucleation) in the resulting gasexpanded, ultimately evaporating solvent droplets; R5: Individual particle growth, either by solute molecule attachment, and/or by the acquisition of previously nucleated particles; and R6: Biases associated with particle “harvesting” (collection). This formidable combination, which clearly transcends thermodynamics alone, explains the fact that, while considerable progress has been made in several of these areas (loc. cit.) and some of their interactions, much still remains to be accomplished, even after the passage of approximately two decades since the earliest such exploratory experiments. Indeed, rather fundamental issues evidently remain, and our objective

2. THERMOPHYSICAL MODEL 2.1. Present VLSE Model for This Three-Component, Three-Phase System. We have exploited the earlier studies of 10384

DOI: 10.1021/acs.iecr.5b01155 Ind. Eng. Chem. Res. 2015, 54, 10383−10400

Article

Industrial & Engineering Chemistry Research Dixon and Johnston,11 Kikic et al.,12 Werling and Debenedetti,13 Badilla et al.,14 and Muhrer et al..15 Evidently, reasonable vapor−liquid−solid equilibrium (VLSE) predictions for the nonideal three-phase, three-component S + API (surrogate) + AS system that is composed of toluene, phenanthrene, and carbon dioxide (T-Ph-CO2) can be made using a Peng−Robinson equation of state (PR-EOS)16 formulation for the fluid phases, and a triple-point based formulation for the fugacity of pure solid phenanthrene, leading to ⎛ p sat ⎞ VPh,s sat ln ϕPh = ln⎜⎜ϕPh,0 Ph ⎟⎟ + (p − pPh ) p ⎠ 9T ⎝

condition, which, in the present case, leads to a local vapor phase whose composition is rich in CO2, poor in toluene, and “very poor” in phenanthrene.18 In Section 3, we exploit the fact that this disparity simplifies the prediction of the initial rate of antisolvent influx, which is likely to be CO2 “impingement-ratecontrolled”. 2.2. Supersaturation Produced by Antisolvent Uptake. In Section 4 (below), we focus on the homogeneous nucleation rate of Ph from the liquid phase solution at an ef fective supersaturation (S), which can be expressed in terms of Ph mole fractions and corresponding Ph fugacity coefficients (ϕPh), i.e.,

(1)

⎛ x ⎞ ϕPh(x Ph ; p , T ) S = ⎜ Ph sat ⎟ sat ⎝ x Ph ⎠ ϕPh(x Ph ; p , T )

where sat pPh (T )

=

sat pPh (TPh,tr)

⎡ ΔH sub ⎛ TPh,tr ⎞⎤ Ph exp⎢ − − 1⎟⎥ ⎜ ⎢⎣ 9TPh,tr ⎝ T ⎠⎥⎦

(3)

Because of expanded liquid-phase solution nonideality, S is expected and found to be systematically different from the often-employed mole fraction ratio: xPh/xsat Ph ≡ Sapp. The maximum possible AS-induced supersaturations, written as Smax below, were calculated from the values of xsat Ph(1 bar, 298 sat sat (p, 298 K), and xCO (p, 298 K), assuming the K), xPh 2 (hypothetical) absence of precipitation despite saturation with the antisolvent CO2, and negligible solvent loss. In such cases, the ratio (xAPI/xS) would be the same in the AS-expanded case, as in the initial (p = 1 bar) mixture. This condition, combined with xAPI,max + xS + xAS = 1, fixes xAPI,max and, hence, Smax (see Figure 3). Our calculations of evolving liquid-phase supersaturations also influence the prediction of precipitated particle growth rates, as a result of Ph “monomer” acquisition. However, as discussed in Section 4, the appropriate thermodynamic driving force for growth must be reduced by a Gibbs−Kelvin (GK) correction factor that is explicitly dependent on the effective fluid/solid surface energy and particle diameter. We remark that most earlier modeling studies are comprehensive enough to include a particle population balance, formally applied “classical nucleation theory” in the (implicit) absence of such GK corrections to the particle subsequent growth ratea combination of assumptions considered here to be thermodynamically inconsistent. Whether this formerly implicit simplification is likely to be acceptable in GASP process models is among the issues explicitly addressed in Section 4 (below). Another fluid property that must be considered when particle populations are evolving in a compressible carrier fluid is the extent of specific (or molar) volume change (see, e.g., ref 15 and its “convective” consequences3). As implied above, in the present case, all ternary mixture molar volumes are predicted from the PR-EOS16 and aforementioned parameters2 are predicted via selection of the smallest real root of the dimensionless cubic equation governing the mixture compression factor Z ≡ pV /(9T ). 2.3. Precipitate Volume Fraction and GASP Yield. Of particular interest is the relative amount of API precipitate obtained in the spray-mode GASP process. This defines the GASP yield, ξ ≡ (NAPI,0 −NAPI)/NAPI,0 (where NAPI is the number of moles of API in the liquid phase and the subscript “0” refers to initial conditions). The relative amount of API precipitate can be computed in terms of the nucleation-ratecontrolled precipitate volume f raction, ϕp, by means of

(2)

In the present case (API surrogate = phenanthrene (C14H10)), we used the following values:14,17 ΔHsub Ph = 90.34 MJ/kmol, TPh,tr = 372.4 K, psat (T ) = 29.04 Pa, and VPh,s = 0.1549 m3/ Ph Ph,tr kmol. On the other hand, we have introduced ϕPh,0 = 0.1915 as a semiempirical correction factor11 to recover experimentally observed results in the low-pressure limit for the present threephase, three-component system. The equilibrium composition based on this PR-EOS (for the fluid phases) and eqs 1 and 2 (for the solid phase) VLSE formulation is in good agreement with the available experimental results11 in the extended pressure range between 1 bar and the critical pressure (ca. 64 bar) of this mixture (see Figure 2).

Figure 2. Predicted pressure dependence of liquid phase compositions (mole fractions, xi) in the T−Ph−CO2 system at T = 298.15 K in equilibrium with Ph(s) (1 < p < 64 bar). Solid curves represent the model, dots points represent experimental results.11

The component critical parameters used for our present calculations, as well as the six relevant dimensionless “interaction parameters”kij and Sij for predicting the PREOS parameters, aij and bij, respectively, for the “participating” ij (i ≠ j) binariescan be found in ref 2. (A summary of the parameters used in this work is given in Table 1.) At, for example, T = 298.15 K and pressures below about ca. 64 bar, a ternary liquid and dense vapor phase will coexist, with the pressure dependence of the saturation mole fraction of Ph in the liquid phase exhibited in Figure 2. Compositions across any VLSE tieline are related by the Gibbs equal species activity 10385


Article

Industrial & Engineering Chemistry Research Table 1. Numerical Values of Parameters Used in the Present Work Pure Component Properties (AS, S, API Surrogate) M [kg/kmol]

component CO2 (AS) toluene (S) phenanthrene (API)

Tc [K]

44.01 304.2 73.8 0.0940 92.14 591.8 41.1 0.3155 178.23 869.3 29.0 0.554 Binary Interaction Parameters and Mixing Rules for Peng−Robinson EOS

unlike pair (kii = Sii = 0) CO2−toluene CO2−phenanthrene toluene−phenanthrene binary coefficient

Zc 0.274 0.263 0.222

kij

Sij

0.090 0.078 −0.004

0.0 −0.030 0.059

aij = (1 − kij)(aiaj)1/2

(

)

bij = 1 − Sij

a= ∑i ∑j xi xj aij Phenanthrene Sublimation at Triple Point Data (TAPI,tp = 372.4 K)

mixing rules

ω 0.225 0.257 0.495

Vc [m3/kmol]

pc [bar]

bi + bj 2

b = ∑i ∑j xi xj bij

sublimation enthalpy ΔHsublim = 90.34 MJ/kmol saturation pressure ptp = 24.04 Pa solid condensate volume VPh,s = 0.1549 m3/kmol Surface Energy and Nucleation/Growth Rate Law-Related Parameters surface energy

nucleation/growth rate law

⎛ nS ⎞1/3 ⎟ Γ = 1.06 ln⎜ API L ⎠ ⎝ nAPI n=2

kG = 10−4

Γrepl = 106

N = NAPI + NS + NAS = (1 − ϕp)

= V

Ni = xiN

(6)

(where V is the liquid phase mixture molar volume and xi is the liquid phase mole fraction of component i), and assuming no appreciable loss of solvent or API during the isobaric, isothermal AS-uptake stage, we derive the following useful inter-relations: (1 − ϕp)

(8)

⎛ x ⎞ 1 − xAPI,0 − xAS,0 ξ = 1 − ⎜⎜ API ⎟⎟ ⎝ xAPI,0 ⎠ 1 − xAPI − xAS

NAvm =

(9)

where the molar volume of solid API (VAPI,S = NAvm) is assumed to be pressure-independent, and the “CO2-expanded” liquid molar volume (V) is computed according to the PREOS. These relations are exact and can be used to estimate bounds to the overall yield of the GASP process. For instance, reference “final” values of =/=0 , ϕp and ξ can be calculated assuming complete precipitation of API above xsat Ph after the ASuptake stage. Hence, assuming that, in the long time limit, xi approaches its VLSE saturation value, we find the following reference quantities (denoted by the subscript “∞”):

(4)

where the precipitate volume fraction ϕp is defined as ϕp ≡ (NAPI,0 − NAPI)

(7)

⎡ ⎤−1 ⎛ V ⎞ 1 − xAPI,0 − xAS,0 ⎢ ⎥ ⎟⎟ ϕp = 1 + ⎜⎜ ⎢⎣ ⎝ VAPI,s ⎠ xAPI,0(1 − xAS) − xAPI(1 − xAS,0) ⎥⎦

Figure 3. Maximal reference supersaturation reached as a function of pressure for an initially VLE CO2−T−Ph mixture at p = 1 bar and T = 298.15 K, with (initial) xPh/xsat Ph values between 0.1 and 1 at steps of 0.1. Solid curves represent the maximal reference supersaturation results using eq 3; dashed curves correspond to the apparent maximal max = supersaturation, given by xPh/xsat Ph. A thick horizontal gray line at S 1 has been added for reference.

ϕp /(NAvm) ⎛ = ⎞ ξ= ⎜ ⎟ xAPI,0/V0 ⎝ =0 ⎠

⎛ V ⎞ 1 − xAPI,0 − xAS,0 = =⎜ ⎟ =0 ⎝ V0 ⎠ 1 − xAPI − xAS

(5)

(Here, = is the total droplet volume, i.e., (1 − ϕp)= is the volume occupied by the liquid phase; note that this definition implies that ϕp(0) = 0.) Taking into account our definitions, together with the overall mass balance equations

(1 − ϕp, ∞) 10386

=∞ V sat ⎛ 1 − xAPI,0 − xAS,0 ⎞ = ⎜ sat sat ⎟ V0 ⎝ 1 − xAPI − xAS =0 ⎠

(10)


Article

Industrial & Engineering Chemistry Research ϕp, ∞

⎡ ⎤−1 ⎛ V sat ⎞ 1 − xAPI,0 − xAS,0 ⎢ ⎥ ⎜ ⎟ = 1+⎜ ⎟ sat sat ⎢⎣ ⎝ VAPI,s ⎠ xAPI,0(1 − xAS ) − xAPI(1 − xAS,0) ⎥⎦ (11)

⎛ x sat ⎞ 1 − xAPI,0 − xAS,0 ξ∞ = 1 − ⎜⎜ API ⎟⎟ sat sat ⎝ xAPI,0 ⎠ 1 − xAPI − xAS

(12)

which can be calculated in terms of thermodynamic (VLSE) quantities for any particular set of initial conditions. The behavior of =∞/=0 , ϕp, ∞ and ξ∞, as a function of pressure, is shown in Figures 4, 5, and 6, respectively, corresponding to an

Figure 6. Reference fraction of API precipitated (eq 12) reached vs p for an initial CO2−T−Ph VLE mixture at p = 1 bar, T = 298.15 K, with an (initial)xPh,0/xsat Ph value between 0.1 and 1 in steps of 0.1.

our nonequilibrium treatment is our principal current interest. This question is answered below (see Section 5).

3. ANTISOLVENT ADDITION RATE AND BALANCE EQUATIONS 3.1. Characteristic Times for Antisolvent Uptake. Of decisive importance in our present GASP process model is the rate of antisolvent uptake and its associated characteristic time (tAS). Because the liquid solvent containing API is often injected into an almost-pure AS vapor, often CO2, we consider that the initial rate of uptake per unit liquid/vapor interface is controlled by the molecular impingement flux Ż AS ″ calculated at the prevailing vapor state: yAS, p, T. The order of magnitude19 of Ż ″AS (kmol/(m2 s)) was previously expected to be set by the ideal gas kinetic theory result: 1 ″̇ ZAS,ref = nAS,w cAS,w ̅ (13) 4 where the antisolvent molecular number density can be written as nAS,w = yAS,wNA/VV,w, in terms of yAS,w (the local AS-mole fraction in the compressed gas) and the vapor mixture molar volume (VV,w, evaluated at the solvent droplet/vapor interface) (denoted by the subscript “w”). On this basis, we define the characteristic anti-solvent (AS) introduction time tAS as the time required to saturate a fixed inventory (e.g., solvent droplet of initial diameter d0) fluid of S + API when the rate of AS introduction is fixed by the initial AS impingement flux, e.g., Ż AS,w ″ in the vapor at the vapor/liquid interface. This leads us to the relation

Figure 4. Reference final droplet volume =∞ divided by its initial value =0 vs p for an initial CO2−T−Ph VLE mixture at p = 1 bar, T = 298.15 K, with (initial)xPh,0/xsat Ph = 0.1 (upper curve) and 1 (lower curve). (See eqs 10 and 11.) As can be seen, at high pressures, =∞ ≫ =0 .

Figure 5. Reference solid volume fraction (eq 11) reached vs p for an initial CO2−T−Ph VLE mixture at p = 1 bar, T = 298.15 K, with an (initial) xPh,0/xsat Ph value between 0.1 and 1 in steps of 0.1.

tAS ≡

sat ⎛ =∞ ⎞ d0xAS ⎜ ⎟ sat ″ ̇ 6V L ZAS,w ⎝ =0 ⎠

(14)

which is seen to scale linearly with the initial solvent droplet diameter (d0). In the former relation, =∞/=0 is the ratio between the reference value of the AS-saturated droplet volume divided by the initial droplet volume, which is given according to our previous estimations (eqs 10 and 11). In the range of pressures that leads to higher expected values of ξ (i.e., above p ≃ 50 bar), the characteristic AS uptake time tAS increases dramatically with p. This is a consequence of the high values reached by the volume ratio =∞/=0 in that pressure range, as xsat AS becomes closer to unity (see Figures 2 and 4). In that case, a different characteristic time (tAS,0),

initially 1 bar VLE mixture of toluene and CO2, with several values of initial API undersaturation (xPh,0/xsat Ph(1 bar, 298 K) between 0.1 and 1), suddenly compressed up to a final pressure level. As may be seen, the thermodynamically calculable estimations ϕp,∞ and ξ∞ drop for pressures below, ca., 54 bar and also for high pressures (close to the mixture critical pressure, ca. 64 bar), as xsat CO2 approaches 1, showing a maximum for pressures between the former limits. Whether or not these thermodynamically expected general trends are indeed found in 10387


Article


the characteristic time for AS (here, AS is CO2) diffusion across the liquid-like solvent + API droplet. The characteristic diffusion time, estimated as tmix = d02/(4DAS−mix), is quadratic in d0, whereas the AS-uptake time (eq 14, above) scales linearly with d0. This estimation (see Figures 7 and 8) shows that AS diffusion limitations would probably be non-negligible for any droplet size much above the nanometric level. However, at larger droplet sizes, internal circulation and/or turbulence could “come to the rescue”, broadening the applicability of WMDmodel results. In any case, we adopt this convenient well-mixed droplet model and hereafter consider the simplest case of a constant rate of AS introduction at constant pressure and temperature up to a reasonable multiple of the aforementioned time, tAS when the initial inventory of S + API becomes saturated with the antisolvent. For this purpose, we also exploit the useful earlier finding that, during the stage of rapid AS uptake, there is negligible evaporative loss of solvent S or solute API. These assumptions underlie the molar balance equations developed in ref 2. Because it is likely that diffusion limitations will lead to higher local supersaturations, we conjecture that the present results provide a lower bound to the amount of precipitated API that would be found if intradroplet diffusion were taken into account. Thus, in cases that lead to a significant API nucleation, the present results may provide an approximation to the corresponding volume-average results, based on finite rate intradroplet mixing. Analysis of this assumption will be the subject of a forthcoming publication. Adopting these spray-motivated simplifications, we might expect the resulting population of precipitated API particles to respond systematically to the following two dimensionless ratios of the three relevant characteristic times, i.e.,

defined also by the AS uptake rate, becomes relevant. Namely, we define tAS,0 as the time scale related to droplet composition evolution under initial conditions. In the present work, the AS uptake rate is assumed to be constant, relative to time: ̇ = πd 2ZAS ″̇ = πd0 2ZAS,ref ″̇ NAS

(15)

Hence, the time needed to achieve an appreciable variation in droplet composition is given as tAS,0 ≡

=0 , thus ̇ V0NAS

⎛ ⎞ tAS sat V0 =∞ = xAS ⎜ sat ⎟ tAS,0 ⎝ V =0 ⎠

(16)

xsat AS,

As a consequence of the behavior of as a function of pressure, these two characteristic times have the same order of magnitude at pressures in the range of p ≈ 40−50 bar. However, as higher values of pressure are considered, tAS becomes considerably higher than tAS,0. For instance, at p = 62 bar, tAS is ca. 2 orders of magnitude larger than tAS,0 (see Figures 7 and 8). In that case, a short initial transient

DamN ≡ Figure 7. Characteristic time scales related to particle nucleation (solid gray curve), particle growth (short-dashed curve), intradroplet mixing (long-dashed line), and AS uptake (solid curves), as a function of pressure for a 1 μm initial diameter solvent droplet.

tAS tN

and

DamG ≡

tAS tG

(17)

where tN and tG are, respectively, characteristic times associated with the rate processes of homogeneous nucleation (Section 4.1) and precipitate particle growth (Section 4.2). However, while particle nucleation and growth are, in principle, independent rate processes, they are both tied to the evolution of API supersaturation, with nucleation being far more sensitive. Assuming, for simplicity, a constant rate of antisolvent AS uptake (given by eq 15), and neglecting any S or API losses during the isobaric, isothermal, AS uptake stage, the corresponding carrier fluid balance equations (see, e.g., Chapter 2 in ref 10) for each one of the droplet components (API, S, and AS) lead to the evolution equations for the spatial average droplet composition and volumetreating the hypothetical expanding liquid “droplet” as well-mixed at each instant. The corresponding set of carrier fluid balance equations is given by 7−11 in ref 2. This set of ODEs, along with the PREOS for the liquid phase molar volume V(x;p,T), as a function of droplet composition (at the current pressure and temperature), and the population balance equation (PBE) governing the precipitated API particle size distribution n(v, t), to be considered in Section 4.3 below, comprise our simplified WMD-mathematical model for this nucleation rate-controlled multiphase system.

characterized by fast droplet composition variation is observed, with a duration much shorter than the time needed to achieve AS saturation. Hence, in the range of high pressures of interest, droplet composition evolution alone is enough to define two different time scales, making this a multiple scale problem. Of course, particle nucleation and growth define their own characteristic time scales, and the relation between these scales determines the evolution of the system. 3.2. Role of Antisolvent Intradroplet Diffusion Limitations? We note that if the characteristic AS uptake time, tAS, is much longer than the characteristic liquid phase dif f usion time (tmix), then we can presume the dilating liquid phase (in which precipitation commences above some threshold “supersaturation”) will be reasonably “well-mixed” at each instant, thereby simplifying the population balance considerably and leading to a situation having features in common with the semibatch (but variable pressure) model developed by the ETH group of ref 15 and also exploited by refs 20 and 21. Indeed, assuming a “worst-case” scenarioi.e., that intradroplet molecular diffusion is the only mechanism for mixing it is possible to roughly estimate what values of the initial solvent droplet diameter (d0) would be small enough at any particular pressure such that tAS (see above) will be larger than

4. MATHEMATICAL MODEL FOR THE NUCLEATION AND GROWTH RATE LAWS 4.1. Particle Nucleation Rate Law. Apart from our evaluation of the relevant (effective) surface energy,2 we invoke 10388


Article


Figure 8. Characteristic time scales related to particle nucleation (solid gray lines), particle growth (short-dashed lines), intradroplet mixing (longdashed lines), and AS uptake (solid lines, upper line = tAS, lower line = tAS,0), as a function of initial solvent droplet diameter (d0) and pressure level.

the formalism of classical nucleation theory (CNT), as modified to account for solution nonideality22−24 (and solute nondiluteness). All present nucleation rate calculations reported here (Section 5) are therefore based on the rate law: 3⎞ ⎛ ̇ = B0̇ ‴ exp⎜ −16π Γ ⎟ B‴ ⎝ 3 ln 2 S ⎠

where the liquid-phase molar volume is given by the PR-EOS at the (time-dependent) liquid composition (and constant p, T). An interesting corollary of eq 20 is that classical nucleation theory (CNT) will necessarily fail when the local supersaturation (S) exceeds the “limiting value” ⎡⎛ 32π ⎞1/3 ⎤ ⎟ Γ⎥ S lim CNT = exp⎢⎜ ⎢⎣⎝ 3 ⎠ ⎥⎦

(18)

where we write the pre-exponential factor in the factored form: ″̇ nAPI vm 2/3 Γ1/2 Γrepl B0̇ ‴ = 2ε ZAPI *

(because v*/vm, the number of API monomers in the critical cluster, must exceed unity). Finally, Γrepl is the so-called “replacement factor”i.e., a correction factor (probably, at most, 106, see ref 28) introduced to account for fluctuations of the center of mass of the liquid droplet representation of a molecular embryo. Values of the model parameters kG, n, and Γrepl considered in the present illustrative study are kG = 10−4 (at T = 298 K; Section 4.2), n = 2, and Γrepl = 106. While the naphthalene/CO2 crystal growth system has been studied experimentally34 more-specific growth kinetics data for the phenanthrene−CO2 system will be needed to refine these preliminary but plausible estimates (also see Section 4.2). 4.1.1. Characteristic Nucleation Time. It remains for us to establish an appropriate choice of characteristic time (tN) associated with the rate process of homogeneous nucleation, completing the definition of the dimensionless time ratio, DamN ≡ tAS/tN, used below. While there are clearly other possibilities, we can imagine that tN is of the form np/Ḃ ‴ ref, where np is a characteristic particle number density and Ḃ ‴ ref is a characteristic value of the particle birth rate per unit volume. This can be rewritten as

(19)

In the latter expression, ε* = ε ln S = kG(T, pCO2) lnn S is an effective API incorporation probability (see eq 27 in Section 4.2), when monomers impact on a particle of the singular size, d*, given by 4γvm 4Γvm1/3 d = = * k T ln S ln S B

(20)

Ż ″API is the impingement frequency of API monomers per unit area on a API particle of diameter d* (see eq 28 in Section 4.2), nAPI is the number density of API monomers, and Γ is the dimensionless surface energy, i.e., the ratio

Γ≡

γvm 2/3 kBT

(21)

where the evolving surface energy γ is conveniently estimated using the Nielsen−Sohnel−Mersmann (NSM) correlation of the following form:2,25,26 ⎛ ns ⎞1/3 ⎟ γ ∝ ln⎜ API L ⎝ nAPI ⎠

(23)

tN =

(22)

ϕp ,ref ̇‴ vp̅ ,ref Bref

(24)

where vp̅ is the number-mean precipitated particle volume. Thus, for ϕp,ref we choose the thermodynamically calculable terminal saturation value of the API particle volume fraction:

where the proportionality constant implicit in eq 22 has been determined based on previous estimates of γ for the binary Ph/ CO2 system2,27 (obtaining the value 1.06; see Table 1), and 10389


Article

Industrial & Engineering Chemistry Research ϕp,∞ (see eq 11 and Figure 5). For the reference particle size, we select a volume corresponding to the reference value of the CNT critical nucleation diameter: vp̅ ,ref ≡

π d ,ref 3 , with 6 *

4Γref 1/3 d ,ref ≡ vm * ln Sref

(25)

where Γref and Sref are reference values of dimensionless surface energy and supersaturation, respectively (see below). Finally, the reference value of particle birth rate (Ḃ ref ‴ ) is defined as its corresponding value at Γref and Sref. As can be expected that the value of Ḃ ‴ ref is strongly dependent on the values chosen for Γref and Sref, especially on the former, which is cubed inside the exponential, while S enters through its natural logarithm. On the other hand, it is not entirely trivial to select appropriate reference values for Γ and S. This is because each one of these variables undergoes a substantial variation with time from its starting value. While Γ increases by a factor of ∼4 (at high pressures), S first increases as the AS enters the droplet, and finally decreased as the API precipitates. Since the dependence of Ḃ ‴ ref on Sref is less severe than that on Γref, the reference value of supersaturation has been defined directly as its corresponding hypothetical value after CO2 addition (up to the saturation level) neglecting phenanthrene condensation (see Figure 3), which is a reasonable estimate of S, at least at values of time corresponding to maximal API precipitation. However, it has been observed that if we estimate Γref accordingly (as its corresponding hypothetical value after CO2 addition, neglecting phenanthrene condensation, which is defined as Γ∞), then the reference value obtained for Ḃ ref ‴ at high pressures is unrealistically low, by many orders of magnitude, which is a direct consequence of the strong dependence of Ḃ ‴ on surface energy.2 The numerical results show that, in the high-pressure limit, API condensation is limited to a short transient, after which condensation cannot take place, because of the increased value of Γ. Thus, a meaningful estimation of Γref should take into account not only its value in the long time limit (which reference value is Γ∞), but also its initial value (Γ0). Hence, in the present work, the reference value Γref has been defined as the following linear combination of initial and (estimated) final values of Γ: Γref ≡ (1 − α)Γ0 + αΓ∞, where the value of α has been chosen as 1/3 or 0, depending on whether SEE has been taken into account or ignored (respectively). The reference value adopted here for supersaturation is given by its theoretical maximal value, and, accordingly, the reference value adopted for particle size (vp̅ ,ref) is a lower bound, rather than a typical value. However, as will be shown below, the volume of most nucleated particles lies in a range between v / vp,ref ̅ = 6(10), which makes the former ̅ = 6(1) and v / vp,ref estimation a convenient choice. On the other hand, we note that this estimation of vp̅ ,ref is entirely based on particle nucleation, and neglects particle growth. The reason for this is because, under the conditions considered here, typical particle growth times are orders of magnitude longer than the typical nucleation times (see Figures 7 and 8), making particle growth almost negligible on the present (tAS) time scale. In Figure 9, we show the results for the reference particle nucleation volume (eq 25) divided by the corresponding API molecular volume as a function of pressure for several values of initial API under-saturation level (at p = 1 bar). We recall that this estimation of v* is a lower bound, since it is based on Smax.

Figure 9. Reference particle volume v ̅ p,ref divided by API molecular −1/2 (upper line) and volume vm versus pressure level for xPh/xsat Ph = 10 = 1 (lower line). xPh/xsat Ph

As a consequence, this estimation shows that the CNT continuum approximation used in the present work is justified. 4.2. Rate of Particle Growth from Gas-Expanded Solvents in the Presence of Interface Curvature. Previous studies of organic particle crystal growth have employed particle growth laws, which are not well-suited to near-critical ASP applications of concern here, in part because they were developed for either ideal or liquid-like crystal growth environments,29 complicating the intercomparison of such data, and extension to describe arbitrary density environments. This worthwhile goal is simplified by explicitly factoring out the API molecular impingement frequency (itself dependent on both temperature and molecular volume fraction; see ref 9 and below), thereby isolating the kinetically significant dimensionless effective “incorporation rate constant”, written as ε below. For these reasons, we write the linear growth rate in the factored form:9 ⎛ F(GK) ⎞ ⎛ S ⎞ ″̇ k G(T , pCO , ...) (ln n − 1 S)⎜ ⎟ ln⎜ Ġ = vm ,API ZAPI ⎟ 2 ⎝ S ⎠ ⎝ F(GK) ⎠ (26)

where the last factor (also dimensionless) is proportional to the chemical potential difference “driving” the condensation process, vm, API is again the volume/molecule of the API in the solid state and Ż API ″ is the API monomer impingement frequency on a unit area of the growing solid surface under curved interface saturation conditionsi.e., accounting for the Gibbs−Kelvin(−Ostwald) effect of small particle interface curvature using, for example, F(GK) ≃ Sd*/dp. Much of the experimental growth rate data in the literature appears to be based on observations on particles large enough so that interface curvature does not affect their solubility, yet not so large that external molecular diffusion limitations become appreciable (thereby avoiding the evaluation of the effective supersaturation, S, existing at the fluid/solid interface) (see, e.g., ref 29). For example, a commonly employed phenomenological separable form for growth in ordinary liquids is Ġ ∝ K(T) lnn S, where n is the so-called growth “order”. Power laws of this type even follow from simple terrace−ledge−kink (TLK) crystal growth rate models.30 For these reasons, we write the individual particle growth rate law for S > 1 as above, but with the explicit choice:9 ε(T , pCO , S) = k G(T , pCO ) ln n − 1 S 2

10390

2

(27)


Article

Industrial & Engineering Chemistry Research where ε and kG will often be of the same order of magnitude (especially if n < 2). In the absence of appreciable fluid phase diffusional limitations (as will be assumed in what follows), S is approximately the value of the effective API supersaturation in the well-mixed “bulk” (the prevailing AS-expanded solvent + API mixture). It remains to estimate the API impingement frequency per unit area, Ż ″API above, in an environment of intermediate molecular volume fraction (between zero [the ideal gas kinetic theory limit] and liquid-like ca. ϕL = 0.64). For our present purposes, the following simple interpolation method is expected to be adequate, irrespective of fluid phase (liquid or vapor):31 Zi̇ ″ =

(L) ni ci̅ ⎛⎜ 1.6Di ,mix ⎞⎟ 4 ⎜⎝ (1/4) ci̅σmix ⎟⎠

environments (T, S, p, ... combinations) such that kG ln S exceeds unity. 4.2.1. Characteristic Growth Time. We define a characteristic growth rate Ġ ref as the value of Ġ for a hypothetical planar surface after CO2 addition but neglecting phenanthrene nucleation, i.e., ⎛ ln n Sref ⎞ ̇ ≡ vm ,API ZAPI,ref ″̇ Gref k G(T , pCO )⎜ ⎟ 2 ⎝ Sref ⎠

(29)

The characteristic evolution time related to particle growth, tG, is then defined as the time needed to achieve a particle radial growth of d*,ref/2 under a growth rate Ġ ref: d ,ref tG ≡ * ̇ 2Gref

ϕ / ϕL

(30)

where the reference particle diameter is defined in eq 25. If kG is fixed at 10−4, for example, the dependence of tG on pressure is shown in Figures 7 and 8, which have been computed for an initially saturated API + S + AS mixture at 1 bar suddenly but isothermally compressed to the final pressure level p. These figures clearly show that, under the conditions considered here, particle growth is negligible in the time scale dictated by the AS uptake stage. Tables 2 and 3 (presented later in this work) summarize the general trends that can be derived from our mathematical model, for the case of a solvent droplet with an initial diameter of 1 μm, as a function of pressure and initial phenanthrene undersaturation. Figures 10 and 11, on the other hand, show the

(28)

where i = API, c ̅ i is the mean molecular thermal speed (c ̅ i ≡ [8kBT/(πmi)]1/2), D(L) i,mix is the Fick diffusion coefficient when ϕ = ϕL = 0.64 (random close-packed hard spheres), and σmix is the effective molecular hard sphere diameter of the prevailing carrier fluid at temperature T. As mentioned earlier, this method is also invoked to calculate the AS impingement frequency on the initial solvent droplet2 and is also used for the API monomer impingement frequency in implementing “classical nucleation theory” (Section 4.1). For the numerical simulations to be discussed in Section 5, we seek to isolate/display the effects of various growth rate-law simplifications, as well as quantify the sensitivity of calculated API particle population characteristics to growth-law parameters (e.g., kG and n). For this reason, we deliberately selected an efficient computational method (Section5), which imposes no restriction on the complexity of the single-particle growth law. This information complements earlier studies that either claimed an insensitivity to the functional form of growth law (e.g., ref 15, which simply assumed Ġ ∼ (S − 1)n with n = 2 or 3/2), or ref 32, which (in a dilute solute RESS situation) actually presumed an API diffusion-controlled rate law irrespective of particle size, in which case Ġ ∼ (S − 1)/dp, enabling small particles to grow faster than large ones. The assignment of realistic growth rate parameters is complicated by the sparsity of phenanthrene crystal growth rate data in the physical chemistry literature and the fact that the parametrization of earlier GASP workers has not taken advantage of a “collision theory” perspective. Based on PSD measurements carried out in a variable-pressure batch device over probably a two-decade range of Ph impingement rates, in refs 15 and 20, we find best-fit values of the ratio: Ġ /(S − 1)n, with n near (or equal to) 2. From our perspective, this particular ratio should not be expected to be a constant. If we, instead, make use of eqs 26−28 above, we can only infer from earlier work that our kG (at T = 298 K) is probably of the order of ca. 100 ppm, used as our “baseline” value (Section 5), along with n = 2. Other data in the literature of polycyclic hydrocarbon crystal growth rate kinetics33 also appear to point to kG(298 K) values in the range of 1 to 100 ppm, with some papers reporting n = 1 (as in the case of naphthalene single crystal growth from supercritical CO2 solutions near 308−318 K34). We remark that these kG assignments correspond to modest effective incorporation probabilities.9 Clearly no growth law should be formally applied in

Figure 10. Contours of constant log DamN (eq 17) vs pressure and initial API mole fraction.

dependence on pressure and initial API undersaturation of the nucleation and growth Damköhler numbers (eq 17). As can be seen, the nucleation Damköhler number is very strongly dependent on the initial API saturation level. This is due in part to the dependence of supersaturation on xAPI, 0, as could be expected, but also the strong dependence of the nucleation rate on the evolving surface energy, which increases as nAPI decreases. Because of SEE (Section 4.1 and ref 2), nucleation 10391


Article


The second term is seen to play the role of a homogeneous sink for particles of any size, with a first-order rate constant given by (d/dt ) ln[(1 − ϕp)=], where = is the instantaneous total droplet volume and ϕp is the present particle volume fractioni.e., the integral of n(v,t)v dv over all particle volumes. This second term is quite important, because of the well-known volumetric expansion of the initially liquid solvent upon CO2 uptake. In thinking about qualitative behavior, as well as developing a practical numerical scheme to solve our PBE, it is useful to step back and briefly consider the more general case in which monomer growth is not negligiblei.e., which would correspond to the addition of a term (∂/∂ v)(v̇n) to the lefthand side (LHS) of eq 31. Then, using the full population balance PDE, on the v vs t plane (i.e., the plane of independent variables), we would be able to construct the “trajectories” of particles born along the v*(t) nucleation locus (see Figure 1). Each of these trajectories would have a slope dv/dt, dictated by the prevailing monomer growth rate law (eq 26), with dv/dt = Ġ a(v) (see Figure 1), and along each such trajectory, n would change in accordance with the second term of eq 31, except when the trajectory crosses the v*(t) locus, when n would “jump” by the amount Ḃ ‴(t*)/|v̇ − v̇*|t*, as shown in ref 2. In our present situation, these trajectories (which are called “characteristic curves” of the first-order PDE) are virtually horizontal lines and most particles will be born in the immediate vicinity of the minimum value of v*. Early studies attempting to predict the evolution of the PSD in transcritical ASP processes considered the relevant surface energy appearing in the nucleation rate expression to be a property of the embryonic solid phase alone. However, because of nondiluteness and highly time-dependent fluid phase solute concentration, this is far from the case in ASP-based micronization processes.2 Employing an independently tested, yet tractable quasi-steady estimation method due to Nielsen, Sohnel, and Mersmann25,26 (NSM, see eq 22), we showed that, in this nucleation-dominated situation, the surface energy evolution (SEE) dramatically reduced the precipitate particle yield and altered the predicted shape of the resulting API powder size distribution (for the model ternary system: phenanthrene + toluene + CO2).

Figure 11. Contours of constant log DamG for kG = 10−4 (eq 17) versus pressure and initial API mole fraction.

can be too slow to be observable, even when the supersaturation appears to be more than adequate. 4.3. Notable Explicit and Implicit Features of the Governing Population Balance Equation. Because the purpose of the GASP process is the gentle “micronization” of a previously coarse API powder, the sine qua non of a sufficiently complete and successful GASP process model is the ability to predict the output particle size distribution function information that transcends the total precipitate yield. This requirement brings us face to face with the relevant population balance equationwhich, as shown below, has several unusual features that become particularly prominent in this class of pharma-motivated “transcritical”, nondilute environments. Taking advantage of our expectation that particle nucleation rates would often dominate monomer-limited particle growth rates (during the rapid CO2 uptake period) and the absence of a spatial dependence in the assumed well-mixed droplet (WMD) limit, the evolution equation for the particle size distribution function simplifies considerably. Using particle volume v to quantify particle size and defining n(v,t) such that its integral over all v would be the instantaneous total particle number density, then, in our previous notation, we would expect n(v,t) to satisfy the differential equation ∂n ̇ δ(v − v ) − n d ln[(1 − ϕ )=] = B‴ p * ∂t dt

5. RESULTS AND DISCUSSION The results shown in the present work have been calculated assuming an initially saturated API + S + AS mixture at 1 bar, instantaneously compressed to the final (constant) pressure level, p, of 48−64 bar. This information, together with the PREOS, fixes the initial value of VL in this isobaric/isothermal system. On the other hand, the initial condition considered here for the NDDF is n(v,0) = 0, which implies ϕp(0) = 0. The values considered for the initial droplet diameter d0 are given by d0 (μm) = m × 10n, where m takes the values 1, 2, and 5 and n goes from −3 to +1 in steps of 1. In all cases, the results have been computed with and without considering the expected effects of surface energy evolution.2 The MOC-based numerical method used to integrate the present simplified WMD model is described in ref.2 For convenience, all of our present parameter assignments (thermodynamic + kinetic) are collected in Table 1 (presented earlier in this work). In Table 2, we give the PR-EOS-predicted numerical values sat sat max of xsat , =ref /=0, ϕp, ∞ and ξref for each choice of Ph, xCO2, VL , S pressure level p (bar) and initial API undersaturation (xPh,0/

(31)

where the first term on the right-hand side (RHS) accounts for the birth rate by homogeneous nucleation of the API solute, and the second term is associated with the rate of increase of the volume accessible to new particles. In the first term on the RHS, Ḃ ‴ is the previously discussed CNT-based rate of homogeneous nucleation (eqs 18 and 19), which is sensitive to the instantaneous API solute concentration primarily via the associated supersaturation S, the corresponding effective surface energy of the crystal/fluid interface (see below), and, to a lesser extent, the impingement frequency appearing in the pre-exponential term Ḃ ‴ 0 . δ(v −v*) is the Dirac function centered about the CNT critical volume v* (vanishing elsewhere), and v* is assumed here to be given by eq 20. 10392


Article


Table 2. Saturation Ph and CO2 Mole Fraction (xisat) and Molar Volume (VLsat), as a Function of Pressure (p), Together with Corresponding Reference Values of Supersaturation (Smax), Droplet Volume Ratio (=/=0 ), Volume Fraction of Solid Deposit (ϕp,∞) and Fraction of API Precipitated (ξref), as a Function of Initial Droplet Undersaturation (xPh,0/xsat Ph) after CO2 Uptake up a to the Saturation Level, with xsat (1 bar, 298 K) = 0.2151 Ph p [bar]

a

−3 xsat Ph [× 10 ]

46 46 46

71.31

48 48 48

61.44

50 50 50

50.24

52 52 52

36.61

54 54 54

17.30

xCO2sat 0.4995

0.5392

0.5862

0.6479

0.7531

VLsat [m3/kmol]

xPh,0/xsat Ph

Smax

ϕp,∞ [%]

ξref [%]

0.08197

0.1 0.3162 1

0.1499 0.4741 1.566

1.715 1.581 1.307

8.222

39.93

0.1 0.3162 1

0.1642 0.5145 1.653

1.773 1.634 1.359

8.791

44.39

0.1 0.3162 1

0.1849 0.5728 1.778

1.858 1.712 1.434

9.394

50.04

0.1 0.3162 1

0.2212 0.6747 1.998

2.012 1.854 1.563

9.999

58.05

0.1 0.3162 1

0.3312 0.9823 2.671

2.500 2.305 1.945

0.07881

0.07517

0.07064

0.06388

=ref /=0

10.08

72.76

56 56 56

3.545

0.8830

0.05890

0.1 0.3162 1

0.7303 2.111 5.258

4.666 4.355 3.502

1.254 6.824

57.67 88.71

58 58 58

1.224

0.9365

0.05933

0.1 0.3162 1

1.123 3.262 8.153

8.568 7.963 6.239

4.193 × 10−2 8.724 × 10−1 4.011

11.60 73.35 92.89

60 60 60

0.5434

0.9662

0.06136

0.1 0.3162 1

1.341 3.950 10.23

16.61 15.38 11.86

4.939 × 10−2 4.795 × 10−1 2.137

26.48 77.84 94.09

62 62 62

0.2653

0.9857

0.06456

0.1 0.3162 1

1.169 3.521 9.684

41.38 38.20 29.18

1.124 × 10−2 1.845 × 10−1 8.603 × 10−1

15.01 74.38 93.17

64 64 64

0.1318

0.9990

0.06922

0.1 0.3162 1

0.1630 0.5110 1.581

2.250 × 10−2

43.28

741.2 683.2 518.4

Cases corresponding to maximal (reference) supersaturations below unity have been discarded.

xsat Ph) considered here. Corresponding values of DamG and DamN follow immediately from these thermodynamically dictated quantities once we specify the particle growth rate parameters kG and n and the initial solvent droplet diameter (d0). The resulting values of DamG and DamN for the case of a d0 = 10−6 m droplet, are summarized in Table 3 and Figures 10 and 11 in the parameter “space” explicitly considered here. The estimations shown in Table 3 clearly show the strong dependency of DamG and (more remarkably) DamN on pressure and initial API saturation level. As a consequence, we conclude that the possibility of API nucleation is nonnegligible only at high initial API saturation. 5.1. Surface Energy Evolution Effects on the Nucleation Damköhler Number and Precipitate Yield. Results for the nucleation Damköhler number (DamN) for all initial droplet sizes considered are shown in Figure 12, as a

function of pressure. We see that DamN increases with d0, because of the linear dependence of tAS with d0. On the other hand, Figure 12 also shows that, if surface energy evolution is neglected, DamN continues to increase monotonically with p, whereas DamN decreases at pressures near 64 bar, after reaching a maximum value. This decrease is produced by a sharp increase in the characteristic nucleation time, as a consequence of the higher values of surface energy reached when the saturation level of CO2 approaches unity, similar to that which occurs at high pressures (see Figure 2). However, a corresponding decrease in the GASP process yield (ξ) is not observed at high pressures (see Figure 13). This is because, at high CO2 pressures, nucleation takes place as a short and intense “burst”, which takes place during a short time interval close to the initial time, before the increase in effective surface energy γ produced by the high level of CO2 in the droplet 10393


Article


Table 3. Summary of Parametric Cases (Combinations) and Their Associated Damköhler Numbers (DamG and DamN), Corresponding to an Initial Droplet Diameter d0 = 10−6 ma p [bar]

xPh,0/xsat Ph

tAS [ms]

v*,ref [nm3]

v*,ref/vm

DamG

DamN

11.91

46.31

2.220 × 10−5

3.767 × 10−1

8.743

33.99

3.282 × 10−5

2.210

46 46 46

0.1 0.3162 1

1.483 × 10 1.367 × 10 1.130 × 10

−5

48 48 48

0.1 0.3162 1

1.641 × 10 1.513 × 10 1.258 × 10

−5

50 50 50

0.1 0.3162 1

1.872 × 10−5 1.725 × 10−5 1.445 × 10−5

6.074

23.62

5.217 × 10−5

1.176 × 101

52 52 52

0.1 0.3162 1

2.278 × 10−5 2.100 × 10−5 1.770 × 10−5

3.738

14.53

9.672 × 10−5

6.440 × 101

54 54 54

0.1 0.3162 1

3.480 × 10−5 3.208 × 10−5 2.707 × 10−5

1.539

5.982

2.890 × 10−4

5.066 × 102

56 56 56

0.1 0.3162 1

7.902 × 10−5 7.376 × 10−5 5.931 × 10−5

21.11 4.659 × 10−1

82.06 1.811

6.567 × 10−5 1.088 × 10−3

4.996 × 10−9 2.675 × 103

58 58 58

0.1 0.3162 1

1.463 × 10−4 1.359 × 10−4 1.065 × 10−4

1.774 × 104 6.559 3.169 × 10−1

6.897 × 104 25.50 1.232

1.429 × 10−7 2.043 × 10−4 1.650 × 10−3

0.0 2.401 × 10−2 4.586 × 103

60 60 60

0.1 0.3162 1

2.707 × 10−4 2.506 × 10−4 1.934 × 10−4

1.288 × 103 5.154 3.188 × 10−1

5.007 × 103 20.04 1.240

2.123 × 10−6 2.875 × 10−4 1.883 × 10−3

0.0 1.070 × 10−1 5.857 × 103

62 62 62

0.1 0.3162 1

6.255 × 10−4 5.775 × 10−4 4.411 × 10−4

1.044 × 104 8.734 5.050 × 10−1

4.057 × 104 33.96 1.963

3.801 × 10−7 2.524 × 10−4 1.821 × 10−3

0.0 9.166 × 10−5 4.798 × 103

64 64 64

0.1 0.3162 1

10.12 × 10−3 9.323 × 10−3 7.074 × 10−3

1.581 × 102

6.145 × 102

8.254 × 10−5

0.0

−5 −5

−5 −5

a

The nucleation Damköhler number is strongly dependent on the initial undersaturation and pressure. However, there is a window of pressures that lead to nucleation times much shorter than the corresponding AS uptake times for initial apparent saturation close to unity (highlighted DamN values). Cases corresponding to maximal (reference) supersaturations below unity have been discarded.

consequence of the higher values of supersaturation reached in that case in the small time limit, which produces a very intense nucleation burst during a short transient, before the saturation level of CO2 is approached. The numerical results for ξ obtained when surface energy evolution is taken into account, on the other hand, are considerably lower than the corresponding results obtained for constant γ. The SEEinduced reduction in ξ is at least a factor of ca. 2 for d0-values in the micrometer range, and reductions by even higher factors are observed for smaller values of d0. Clearly, the often overlooked phenomenon of surface energy evolution introduces a serious limitation on the ef f iciency of the GASP process. Figure 13 also suggests that one way to minimize this reduction in GASP yield is to employ larger droplets (above the micrometer range) and high pressures.

totally inhibits particle nucleation. As already mentioned, of particular interest are our predictions for the overall GASP process yield (ξ), shown as a function of pressure and initial droplet diameter in Figure 13. In this figure, we see that ξ increases monotonically with pressure and initial droplet diameter in all of the cases. The increase with p takes place because higher values of supersaturation are reached in that case, and the increase with d0 occurs because there is more time available for nucleation. In Figure 13, we see that our predictions based on a constant value of γ have a tendency to approach the reference value of ξ(ξ∞) as higher values of d0 are considered (see gray curve in Figure 13; also see Figure 6), although the dependence of ξ with d0 becomes weak when d0 is >1 μm. On the other hand, the decrease in the reference quantity ξ∞ predicted for pressures close to 64 bar is not observed in our numerical results. This is interpreted as being a 10394


Article


Figure 12. Characteristic nucleation Damköhler number (DamN) values vs p for an initially (1 bar) saturated API + S + AS droplet, suddenly exposed to a GASP process pressure level p. Solid curves have been computed taking into account surface energy evolution (SEE), and broken curves have been computed neglecting SEE. Results shown correspond to different values of the initial droplet diameter d0, given by d0 = m × 10n μm (where m = 1, 2, 5 and n = −3, −2, −1, 0, 1).

Figure 14. Sauter mean diameter [nm] vs p for an initially (1 bar) saturated API + S + AS droplet, relative to the initial droplet diameter d0, as described in Figure 12. Solid curves have been computed by considering nucleation control in the presence of SEE; dashed curves have been computed by neglecting the inevitable SEE.

supersaturations. However, this effect becomes remarkably weak as soon as the pressure level that leads to high nucleation rates is reached (at pressures p higher than ca. 55 bar), because, in that case, the high nucleation rates inhibit further increases of supersaturation. On the other hand, the higher values of supersaturation reached when larger droplets are considered (larger tAS) make SMD a monotonically increasing function of initial droplet size, regardless of pressure level. Our results, which have been obtained with a more plausible environmentdependent surface energy, show the same trends. The main difference is that the conditions leading to higher values of supersaturation also lead to higher values of surface energy, producing a partial compensation of its effects on precipitate mean particle size. As a consequence, the sensitivity of the SMD results shown in Figure 14 with SEE (solid curves) is smaller than the corresponding results neglecting SEE (broken curves). This result can also be seen in Figure 15, where the dimensionless precipitate population spread (diameter-based standard deviation σ divided by the SMD value) is shown as a function of pressure and d0. In this case, we find that, when SEE is neglected, the marked increase of SMD with d0 compensates

Figure 13. Process “yield” (ratio of precipitated API/initially available API) vs p for an initially (1 bar) saturated API + S + AS droplet, with an initial droplet diameter d0 as in Figure 12. Solid curves have been computed by considering surface energy evolution (SEE); dashed lines have been computed by neglecting SEE. The thick gray curve corresponds to the reference ξ∞, given by eq 12. [Note that the resulting nucleation rate in each swelling droplet limits the amount of API precipitate that can be obtained using the spray-mode GASP process, in part because SEE decreases the nucleation rate (via increased surface energy) as the API mole fraction decreases.]

5.2. Factors Governing the Mean Precipitate Particle Diameter and PSD-Spread; Trends with Operating Pressure and Solvent Spray Droplet Diameter. Results for the typical precipitate particle size (Sauter mean diameter (SMD [nm])) are shown in Figure 14. Since particle aggregation is neglected (during the CO2 uptake stage) and particle growth remains a small correction during the time interval considered here, the sizes of the precipitated particles remain very close to the critical nucleation size at the corresponding nucleation time. Thus, time-dependent supersaturation is the principal variable that affects the critical nucleation size in the results for SMD based on neglecting SEE. In this case, Figure 14 shows that the SMD decreases as higher pressures are considered, because those cases lead to higher

Figure 15. Dimensionless diameter-based population standard deviation σ/SMD for an initially (1 bar) saturated API + S + AS droplet, relative to the initial droplet diameter d0, as described in Figure 12. Solid curves have been computed by considering nucleation control in the presence of SEE; dashed curves have been computed by neglecting the inevitable SEE. 10395


Article


Figure 16. Dimensionless number density distribution functions predicted by the present WMD-GASP process model at 10 equispaced values of time between t = 0 and AS saturation time, for several values of pressure and d0 = 1 μm. Left column has been computed assuming nucleation control but neglecting SEE; right column has been computed assuming nucleation control but including SEE (log−log plots).

that of σ, making the dimensionless spread a monotonically decreasing function of d0 for all pressures. The same result is obtained when SEE is taken into account but only for low pressures, and a reversed trend is actually observed in the highpressure limit. In any case, we also find that the precipitate population is quite narrow, with a dimensionless spread close to ca. 10% in all the cases.

5.3. Existence of Optimal Solvent Droplet Spray Diameter? An interesting conclusion derived from our present theoretical model is that, in order to obtain a high-yield ultrafine (almost-)monodisperse precipitate via GASP, one must operate under high CO2 pressures with solvent droplets large enough to allow for a high yield, but still small enough to avoid the increased values of SMD and dimensionless population spread found when higher values of d0 are 10396


Article

Industrial & Engineering Chemistry Research considered. Hence, we anticipate the existence of an optimal d0, the value of which (in the ca. micrometer size range) is dependent on the relative importance of the conflicting goals of high yield, on one hand, and small SMD and dimensionless spread, on the other hand. 5.4. Predicted Precipitate PSD Shapes via Present MOC. Figure 16 shows the corresponding number density distribution f unctions obtained for d0 = 1 μm and several values of pressure, artificially neglecting SEE (left column) and realistically including SEE (right column). It is immediately clear that in no case is the NDDF a simple log-normal (or “Gaussian”) shape, suggesting that approximate moment-based PBE methods that presume a constant PDF shape will not be useful to model GASP processes. In all cases considered here, the NDDF essentially tracks the time history of the critical nucleation volume (v*). The critical nucleation size starts with a large value at the initial time, when supersaturations (and nucleation rates) are very low, and decreases with time as supersaturation increases. The sharp peak observed in w·f for small values of w and final values of time in the high pressure cases shown (p = 58 and 62 bar) corresponds to the minimum value reached by v*, as a function of time under those conditions (as sketched in Figure 1). On the other hand, Figure 16 also shows the effects on the NDDF of the available droplet volume ((1 − ϕp)= ), which acts as a homogeneous sink of particles as it increases with time, producing the corresponding decrease in the dimensionless NDDF ( f). This effect is especially remarkable in the cases corresponding to high CO2 pressures (p = 58 and 62 bar), because, in those cases, nucleation takes place mainly during a short time interval close to the initial time. On the other hand, the (mild) nucleation observed in the “low” pressure cases considered (see results for p = 50 and 54 bar) continues to increase during the entire time interval, as the supersaturation rises, thus compensating for the sink effect produced by the increase in available droplet volume.

based on a saturated (1 bar, 298.15 K) ternary mixture formed by phenanthrene (API surrogate), toluene (representative organic solvent S), and CO2 (antisolvent AS), suddenly sprayed into a chamber at the final (constant) pressure level, in the range of 45−63 bar. Thermodynamic and VLSE relations have been described using a Peng−Robinson EOS formulation for the fluid phases, together with a triple-point based formulation for the fugacity of pure solid phenanthrene. The values of the initial droplet diameter formally considered cover the range of 10−3−10+2 μm (see Section 6.4). We have deliberately focused on the evaluation of the characteristic times expected to govern the dynamics of this system, namely, homogeneous nucleation time, growth time, AS uptake times (including the initial AS uptake time and AS uptake time needed to reach AS saturation), coagulation time, and intradroplet diffusion time. In this regard, it is shown that, for droplets injected into pure CO2, during the impingementcontrolled CO2 AS uptake time particle nucleation is the dominant process, with negligible particle growth or aggregation. While the simplified mathematical model adopted in our present paper is based on the comparatively tractable well-mixed droplet limit (which neglects intradroplet CO2 diffusion-rate limitations), our estimation of characteristic times shows that diffusion limitations will be important in the cases previously tested in the laboratory (i.e., solvent droplet size well above the micrometer range). Nevertheless, pending our current extensions to anticipate the principal PSD consequences of intradroplet AS-diffusion limitations (see refs 9 and 35, and Section 6.5), because this phenomenon is expected to produce even higher local supersaturations, our present results would be expected to produce a lower bound to the (volume-averaged) amount of precipitated API. Hence, besides revealing potentially valuable overall trends and qualitative behavior, our present mathematical model may provide a useful approximation to droplet volume-averaged results, at least in those cases leading to intense API nucleation. Since, for the conditions and property values (Table 1) considered here, the rate processes of particle growth or Brownian aggregation have been shown to be negligible, compared to particle nucleation, the precipitated API population essentially follows the time history of the critical nucleation size, showing a marked peak whenever the critical nucleation size reaches a local extremum. In this regard, very narrow precipitate size distributions are expected, with a diameter-based standard deviation close to 10% of the Sauter mean diameter (SMD) of the population. 6.3. Features of the Present PBE Strategy, as Applied to the GASP Environment. As outlined in Section 4.3, the requirement of a faithful representation/solution of the governing population balance equation for this highly coupled class of pharma-motivated applications led us to consider, at least as a first approximation for sufficiently small solvent spray droplets, the CO2-antisolvent uptake stage in the molecularly well-mixed droplet (WMD) limit. This proves to be an instructive exercise, enabling exploitation of the power of the method of characteristics (MOC) in a setting requiring basic changes in the way the particle “classical” homogeneous nucleation rate and growth rate laws are formulated and treated. Thus, our present calculations incorporate what is presently known about (i) the nonideality of these nondilute ternary solutions, (ii) solute impingement rates from fluid mixtures that are neither gaslike nor liquidlike, and (iii) the dependence of the solid/fluid interfacial energy on the evolving

6. CONCLUSIONS AND IMPLICATIONS 6.1. Goals and Scope of the Present Paper. A solvent spray-based variant of the Gas Anti-Solvent Precipitation (GASP) process being explored by the pharmaceutical industry to grow narrow distributions of ultrafine particles has been formulated and analyzed here, with the ultimate goal of providing a tractable process model that can be used to optimize/control such “micronization” processes in the future. As noted in the Introduction, this ambitious goal, which requires explicitly considering the population balance equation and many interacting “physical” rate processes (along with three-component, three-phase nonideal fluid thermodynamics) has understandably remained elusive for more than two decades. However, we have shown that the relatively tractable asymptotic limit explicitly considered here and in ref 2 (i.e., molecularly well-mixed droplet) enables access to product PSD information as well as the nonequilibrium process yield, without making previously used but overly restrictive physical or numerical assumptions. Moreover, we now foreseeand have, in fact, already initiatedtractable extensions (Section 6.5) which are expected to embrace the principal effects of systematic departures from the above-mentioned molecularly well-mixed droplet limit. 6.2. Model Ternary System and Assessment of Relevant Characteristic Times for Participating Rate Processes. Our present numerical results (Figures 12−16) are 10397


Article


predicting both rate-based precipitate product yield and resulting powder PSD yet exists. Our present tractable “wellmixed droplet” (WMD) model nearly attains this ambitious goal, but as seen in Section 4.1.1, only for injected solvent droplets much smaller than those actually used in previous demonstration experiments. Indeed, for solvent droplet sizes in the usual (supermicrometer-)size range associated with highthroughput requirements, the above-mentioned WMD idealization evidently fails, compromising its predictive value in this important domain. However, considering the computational demands of a practical process model, it may be possible to retain the most attractive features of the present WM limit for the all-important precipitate population dynamics, and yet anticipate the systematic effects of CO2 diffusion limitations for supermicrometer-diameter solvent droplets. On this basis, we have initiated work on a plausible yet tractable extension by simply imagining/treating a WM nucleation/growth zone (annulus) that propagates into the “fresh” liquid solvent core at a rate that is controlled by CO2 molecular diffusion. In effect, in this next-generation model, we view the moving interface between the well-mixed outer annulus and the inner solventrich core as a semipermeable expandable “membrane”, impervious to recently born API nanoparticles but permeable to the molecular antisolvent CO2. After accounting for the polydispersed nature of such supermicrometer sprays, comparison of such results with the observed performance of prototypical systems may then be sufficient to validate this tractable “hybrid” approach, opening the way for future pharma process design/control applications.

solute concentration. A further and somewhat surprising simplification in this mathematical model has been the smallness of the effects of monomer-based precipitate growth during this antisolvent uptake stage. While this may not be general for systems with much larger molecular incorporation probabilities (see, e.g., the 2014 work of Rosner, as outlined in greater detail in ref 9, we have exploited the presently expected disparity in characteristic times (Section4) via a perturbation schemein the first (lowest-order) approximation for which monomer-based growth is neglected completely (i.e., the characteristic shown in Figure 1 become horizontal lines). Our present PSD expectations (see Figures 13−16, as well as Section 6.4 below) are “simply” the results of forward integrating the coupled system of transient WMD-species balances and the PBE (Section 4.3) via a sufficient number of selected characteristic “trajectories”, which are initiated along the critical volume v* (vs time) locus. Our provisional suggestions for addressing the effects of systematic departures from the above-mentioned WMD limit, which are expected to be important to address with the common use of supermicrometer-sized solvent spray droplets, are briefly discussed in Section 6.5, which concludes the present paper. 6.4. Results of Well-Mixed Droplet Mathematical Model. It is found that smallest nominal precipitate size (e.g., SMD) and population spread (diameter-based standard deviation) are obtained at high pressures, because these are the conditions in which homogeneous nucleation takes place during a short transient, as a nucleation “burst”. Regarding the effect of initial solvent droplet size, it has been observed that larger droplets lead to higher values of the AS-uptake time and, hence, higher relative amounts of precipitate during AS uptake. On the other hand, initial solvent droplet size has only a modest impact on the precipitate SMD, with larger particle diameters corresponding to larger droplets. Initial solvent droplet size also has a relatively weak effect on the precipitate population spread. In this regard, we find that, in the range of pressures leading to intense API nucleation (p > 55 bar), when surface energy evolution (SEE) is taken into account, larger droplets produce API populations with larger dimensionless standard deviations (diameter-based standard deviation divided by the SMD). Ironically, the opposite result is obtained for lower pressures, and over the entire pressure range when SEE is (incorrectly) neglected. As an overall conclusion from the present tractable asymptotic spray mode-GASP-process model, it appears that high yield ultrafine/(almost-)monodisperse precipitates will require high CO2 pressures and injected solvent droplets with greater than micrometer-sized initial diameters. In principle, larger droplets will produce higher API precipitate yields, but at the risk of also producing a precipitate with larger SMD values and size dispersion. On this basis, an intermediate optimal initial solvent droplet diameter might be expected (perhaps in the ca. micrometer size range), depending on the relative importance of the conflicting goals of high yield on one hand and small precipitate population SMD and dimensionless spread on the other hand. We conjecture that these trends will remain valid in the presence of intradroplet antisolvent diffusion limitations, which is currently under investigation. 6.5. Instructive Next Step: A Molecularly Well-Mixed Annulus (WMA) Model? For process design/control situations in which the API-containing liquid solvent is sprayed into compressed CO2 antisolvent, it is still fair to say that no self-consistent theoretical/mathematical model capable of

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

́ Dept. Fisica Matemática y de Fluidos, UNED, Apdo 60141, 28080 Madrid, Spain. ¶

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Supported by Industrial Affiliates of the Sol Reaction Engineering (SRE) Research Group at Yale, Yale ChE graduate alumni, and the Yale School of Engineering and Applied Science. M.A.Z. also wishes to acknowledge the support of Ministerio de Ciencia e Innovación (Nos. CSD2010-00011 and ENE2011-26868) and Comunidad de Madrid (No. S2009/ ENE-1597) at UNED (Madrid). Short presentations describing our approach and earlier progress were made at 2014 International Aerosol Conference IAC 2014 (Korea, Sept 2014), AIChE 2013 Meeting (San Francisco, CA; Paper Nos. 300f, 394c, Sept. 5, 2013) and at two earlier aerosol conferences: 12th European Aerosol Conference (Sept. 6, 2012, Granada) and 13th European Aerosol Conference (Prague, Sept. 1, 2013). Professor Ramkrishna’s successful application of a population-balance formulation to mathematically model challenging chemical reaction engineering environments, (including even bubbling fluidized beds) has not only influenced the transport-oriented teaching of one of us (D.E.R.), but also encouraged our research group’s interest in applying/extending PB methods to spatially nonuniform flame synthesis (e.g., AIChE J. 2002, 48 (3), 476−491) and the 10398


Article

Industrial & Engineering Chemistry Research AS = antisolvent (here, CO2) API = active pharmaceutical ingredient crit = critical value diff = diffusion-controlled eff = effective value f = fluid phase G = growth i, j = component i, j L = liquid phase N = nucleation Ph = phenanthrene (API surrogate) C14H10 ref = reference value s = solid phase S = solvent (here, toluene) sat = saturation value (i.e., at VLE or SLE) T = toluene (C7H8) tr = triple point (VLSE condition) v = vapor phase w = at the fluid/solid interface

present class of three-phase transcritical nonequilibrium environments of current pharma interest.

■

NOMENCLATURE a(v) = surface area of spherical particle with volume v; a(v) = (36π)1/3v2/3 Ḃ ‴ = birth rate of particles per unit volume and time c ̅ i = mean thermal speed of molecules of component i (IGKT) d = diameter (of particle or droplet) D = Fick diffusion coefficient Dam = Damköhler number (ratio of characteristic times) F(GK) = correction factor due to interface curvature Ġ = linear growth rate of particle or crystal kB = Boltzmann constant K = temperature-dependent part of phenomenological particle growth law kG = dimensionless growth rate constant (see eq 27) mi = molecular mass of component i Mi = molar mass of component i (e.g., kg/kmol) n = kinetic “order” parameter in Ġ (S) expression n = total number (or molar) density of mixture n(v,t) = particle number density distribution function NA = Avogadro’s number; NA = 0.6023 × 1027 molecules/kg mol) ni = number density of component i (molar or molecular) Ni = number of moles of component i Ṅ i = molar flow rate of component i np = particle number density p = pressure (bar) 9 = universal gas constant S = supersaturation (see eq 3) t = time (reckoned from the onset of AS addition) tAS = characteristic time related to antisolvent (CO2) addition rate tG = characteristic time related to (solid) API growth rate tN = characteristic time related to API nucleation rate T = absolute temperature (K) v = volume of an individual precipitate particle V = molar volume of mixture vm = molecular volume (for API in solid state) = = total volume of droplet (containing S + API + AS) x = composition vector {x1, x2, x3}T xi = mole fraction of component i (liquid phase) yi = mole fraction of component i (vapor phase) Z = compression factor (Z ≡ pV /(9T )) Zc, i = critical value of compression factor for component i Ż ′i ′ = molar (or molecular) impingement flux of component i

Acronyms, Abbreviations

■

Greek Letters

γ = surface free energy Γ = dimensionless surface energy (γ v2/3 m /(kBT)) ε = incorporation probability; eq 27, see ref 9. μf = shear viscosity of prevailing fluid mixture ξ = fraction of API precipitated; eq 4 σ = molecular diameter τ = dimensionless time (τ ≡ t/tAS) ϕi = fugacity coefficient of component i in fluid mixture ϕp = API precipitate volume fraction (first volume moment of n(v, t))

AS = antisolvent (here, CO2) API = active pharmaceutical ingredient ASP = antisolvent precipitation CNT = classical nucleation theory (see, e.g., ref 24) EOS = equation of state G = growth GASP = gas(-induced) antisolvent precipitation GK = Gibbs-Kelvin IGKT = ideal gas kinetic theory MOC = method of characteristics N = nucleation NDDF = number density distribution function ODE = ordinary differential equation PB = population balance PBE = population balance equation PDE = partial differential equation PDF = probability density function Ph = phenanthrene (C14H10) PR = Peng−Robinson PSD = particle size distribution RESS = rapid expansion of a supercritical solution S = solvent SASP = supercritical (vapor) antisolvent precipitation SEE = surface energy evolution SMD = Sauter-mean diameter T = toluene (C7H8) VLSE = vapor−liquid−solid equilibrium

REFERENCES

(1) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: New York, 2000. (2) Rosner, D.; Arias-Zugasti, M. Surface Energy ‘Evolution’ (SEE) in Pharmaceutical Powder ‘Micronization’ Using Compressed Gas AntiSolvent (Re-) Precipitation (GASP). Ind. Eng. Chem. Res. 2014, 53, 4489−4498. (3) Rosner, D.; Arias-Zugasti, M. Coupling between homogeneous rate processes and fluid deformation rate: Brownian particle coagulation in a rapidly dilating solvent. AIChE J. 2011, 57, 307−318. (4) Rosner, D.; Arias-Zugasti, M. Thermophoretically Dominated Aerosol Coagulation. Phys. Rev. Lett. 2011, 106, 015502. (5) Arias-Zugasti, M.; Rosner, D. Thermophoretically-Modified Aerosol Brownian Coagulation. Phys. Rev. E 2011, 84, 021401.

Subscripts and Superscripts

0 = initial value (subscript) ∞ = final (reference) value 10399


Article


expansion of a supercritical solution under subsonic conditions. J. Supercrit. Fluids 2002, 23, 65−80. (28) McGraw, R. A. Corresponding States Correlation of the Homogeneous Nucleation Thresholds of Supercooled Vapors. J. Chem. Phys. 1981, 75, 5514−5521. (29) Mullin, J. W. Crystal Growth. In Crystallization, 4th Edition; Elsevier/Butterworth−Heinemann: Boston, MA, 2001; Chapter 6, pp 224−296. (30) Burton, W. K.; Cabrera, N.; Frank, F. C. The Growth of Crystals and the Equilibrium Structure of their Surfaces. Philos. Trans. R. Soc., A 1951, 243, 299−358. (31) Indeed, the factor correcting the well-known IGKT, ϕ ≪ 1 result: (1/4)nic ̅ i is sometimes not very far from unity, although the corrections are not negligible, especially in the liquid phase. For example, at 56 bar, near t = tAS we estimate ϕ = 0.49 and an impingement flux correction factor of 0.65 for CO2 uptake and a corresponding correction factor of 0.1 for Ph impingement rate in the liquid phase. As a consequence, the characteristic AS uptake times are corrected by a factor of ∼2, and the nucleation times by a factor of ∼10, producing an overall reduction in the nucleation Damköhler number (DamN). (32) Shaub, G. R.; Brennecke, J. F.; McCready, M. J. Radial model for particle formation from the rapid expansion of supercritical solutions. J. Supercrit. Fluids 1995, 8, 318−328. (33) Sloan, G. J. Kinetics of Growth of Anthracene from Vapor. Mol. Cryst. 1967, 2, 323−333. (34) Uchida, H.; Manaka, A.; Matsuoka, M.; Takiyama, H. Growth Phenomena of Single Crystals of Naphthalene in Supercritical Carbon Dioxide. Cryst. Growth Des. 2004, 4, 937−942. (35) Rosner, D.; Arias-Zugasti, M. Effects of Intra-droplet CO2 Diffusion Limitations on Spray-Based Anti-solvent Precipitation (ASP-) Process Performance. J. Supercrit. Fluids, Manuscript in preparation (Fall 2015).

(6) Rosner, D.; Arias-Zugasti, M. Novel Features of Aerosol Coagulation in Non-isothermal Environments. Ind. Eng. Chem. Res. 2011, 50, 8932−8940. (7) Martin, A.; Cocero, M. Micronization Processes with Supercritical Fluids: Fundamentals and Mechanisms. Adv. Drug Delivery Rev. 2008, 60, 339−350. (8) Exploitation of initial temperature contrast (between solvent droplets and the CO2 environment) to control/modify GASP-process performance will not be considered further here. Indeed, this attractive option motivated our recent study of Brownian particle coagulation in a local temperature gradient.4−6 (9) Rosner, D. Reconciliation of Kinetic Data for Organic Crystal Growth in Supercritical-CO2: Exploiting a Molecular Collision Theory Viewpoint, with Implications for Precipitator Process Design. Cryst. Growth Des. 2014, 14, 3783−3790. (10) Rosner, D. Transport Processes in Chemically Reacting Flow Systems; Dover Publications: Mineola, NY, 2000. (11) Dixon, D.; Johnston, K. Molecular Thermodynamics of Solubilities in Gas Antisolvent Crystallization. AIChE J. 1991, 37, 1441−1449. (12) Kikic, I.; Lora, M.; Bertucco, A. A Thermodynamic Analysis of Three-Phase Equilibria in Binary and Ternary Systems for Applications in Rapid Expansion of a Supercritical Solution (RESS), Particles from Gas-Saturated Solutions (PGSS), and Supercritical Antisolvent (SAS). Ind. Eng. Chem. Res. 1997, 36, 5507−5515. (13) Werling, J.; Debenedetti, P. Numerical modeling of mass transfer in the supercritical anti-solvent process. J. Supercrit. Fluids 1999, 16, 167−181. (14) de la Fuente Badilla, J.; Peters, C.; de Swaan Arons, J. Volume Expansion in Relation to the Gas-antisolvent Process. J. Supercrit. Fluids 2000, 17, 13−23. (15) Muhrer, G.; Lin, C.; Mazzotti, M. Modeling the Gas Antisolvent Recrystallization Process. Ind. Eng. Chem. Res. 2002, 41, 3566−3579. (16) Smith, J. M.; Van Ness, H.; Abbott, M. M. Introduction to Chemical Engineering Thermodynamics, Seventh Edition; McGraw− Hill: New York, 2000. (17) Roux, M. V.; Temprado, M.; Chickos, J. S.; Nagano, Y. Critically Evaluated Thermochemical Properties of Polycyclic Aromatic Hydrocarbons. J. Phys. Chem. Ref. Data 2008, 37, 1855−1996. (18) For example, at p = 45 bar and T = 298 K (not far from the condition of maximum dissolved Ph fugacity), when the liquid phase composition is approximately xCO2 = 0.48, xT = 0.44, xPh = 0.08, the corresponding vapor phase composition is estimated to be: yCO2 = 0.998,yT = 0.002 and yPh = 3.5 × 10−8. (19) But, more generally, we expect a correction to the RHS of eq 13, which will be dependent on the vapor phase molecular volume fraction, ϕ. For a rational estimate of this correction, see Section 4. (20) Bakhbakhi, Y.; Rohani, S.; Charpentier, P. A. Micronization of Phenanthrene Using the Gas Antisolvent Process. 1. Experimental Study and Use of FTIR. Ind. Eng. Chem. Res. 2005, 44, 7337−7344. (21) Bakhbakhi, Y.; Rohani, S.; Charpentier, P. A. Micronization of Phenanthrene Using the Gas Antisolvent Process: Part 2. Theoretical Study. Ind. Eng. Chem. Res. 2005, 44, 7345−7351. (22) Debenedetti, P. Homogeneous nucleation in supercritical fluids. AIChE J. 1990, 36, 1289−1298. (23) Kwauk, X.; Debenedetti, P. Mathematical modeling of aerosol formation by rapid expansion of supercritical solutions in a converging nozzle. J. Aerosol Sci. 1993, 24, 445−469. (24) Kashchiev, D. Nucleation: Basic Theory with Applications; Butterworth−Heinemann: Boston, MA, 2000. (25) Nielsen, A.; Sohnel, O. Interfacial Tensions Electrolyte CrystalAqueous Solution, from Nucleation Data. J. Cryst. Growth 1971, 11, 233−242. (26) Mersmann, A. Calculation of Interfacial-Tensions. J. Cryst. Growth 1990, 102, 841−847. (27) Weber, M.; Russell, L. M.; Debenedetti, P. Mathematical modeling of nucleation and growth of particles formed by the rapid 10400


Idealized Mathematical Model for Pharmaceutical Powder

Recommend Documents