Ind. Eng. Chem. Res. 2007, 46, 8009-8017
8009
Modeling of the Gas-Antisolvent (GAS) Process for Crystallization of Beclomethasone Dipropionate Using Carbon Dioxide Shawn Dodds,† Jeffery A. Wood,‡ and Paul A. Charpentier* Department of Chemical and Biochemical Engineering, UniVersity of Western Ontario, London, Ontario, Canada N6A 5B9
The purpose of this work was to develop a mathematical model to describe the crystal size distributions (CSDs) produced from the gas-antisolvent (GAS) technique on the crystallization of beclomethasone17,21-dipropionate (BDP), which is an anti-inflammatory corticosteroid commonly used to treat asthma. The solvent used was acetone, and the antisolvent was carbon dioxide (CO2). The GAS technique was chosen for its ability to produce micrometer-sized particles of uniform size. A better understanding of how the GAS process affects the CSDs of BDP is desirable to optimize the GAS experimental conditions for the production of inhalable powders for next-generation dry-powder portable inhalers (DPIs). To describe the pressurization during the GAS process, a mass balance and a phase equilibrium model were required. A predictive relative partial molar volume fraction (RPMVF) equilibrium model was used in the absence of existing phase data for the BDP-acetone-CO2 system. This model uses the binary two-phase solvent/CO2 phase equilibrium, and then relates it to the solid concentration. The model was tested successfully with the phenanthrenetoluene-CO2 model system, the naphthalene-toluene-CO2 model system, and the more complex cholesterolacetone-CO2 model system before being used to predict the BDP-acetone-CO2 system. A population balance was then used to model experimentally determined particle size distributions. Two models for secondary nucleation were used independently: (i) an empirical equation that is commonly used to model secondary nucleation, and (ii) a theory-based equation. The crystallization model was able to give good estimates of the cumulative volume (mass)-weighted size distributions metrics dp(10%), dp(50%), and dp(90%) of the experimental CSDs. 1. Introduction The production of particles with controlled size and variability has attracted significant interest in the scientific and industrial communities with applications for pharmaceuticals, food, nutraceuticals, chemical, paint/coating, and polymers.1-8 The important properties of these products are narrow particle size distribution and uniform morphology.9,10 The application of high-pressure gas and fluid techniques (such as near-critical and supercritical fluids) has attracted considerable interest as an emerging “green” technology for the formation of particles in these size ranges.11,12 Particle formation using near-critical and supercritical fluids can be performed according to several different techniques,7 including antisolvent techniques such as the gas-antisolvent (GAS) process. Antisolvent techniques exploit the low solubility of most compounds in the antisolvent (in particular, carbon dioxide (CO2)), which must be miscible with the organic solvent.13,14 In the GAS process, high-pressure CO2 is injected into the liquid-phase solution, which causes a sharp reduction of the solute solubility in the expanded liquid phase. As a result, precipitation of the dissolved compound occurs. The potential advantages of the GAS recrystallization process lies in the possibility of obtaining solvent-free, micrometer-sized, and submicrometer-sized particles with a narrow size distribution.15 By varying the process parameters such as pressure, temperature * To whom correspondence should be addressed. Tel.: (519) 6613466. Fax: (519) 661-3498. E-mail:
[email protected]. † Presently at Department of Chemical Engineering and Materials Science, University of Minnesota, 151 Amundson Hall, 421 Washington Avenue SE, Minneapolis, MN 55455. ‡ Presently at Department of Chemical Engineering, Queen’s University at Kingston, Kingston, Ontario, Canada K7L 3N6.
or antisolvent addition rate, the particle size, size distribution and morphology can be “tuned” to produce a product with desirable qualities. This makes the GAS technique attractive for the micronization of high-valued products, such as pharmaceuticals.6 Unfortunately, only a limited number of theoretical studies have been undertaken to develop the understanding, and to achieve a quantitative description of the GAS process. Most theoretical approaches to describe the GAS process rely on a coupling of the population balance of particles with the thermodynamic equilibrium. Muhrer et al.,16 to rationalize the precipitation kinetics in the GAS process, presented a model that accounts for solution thermodynamics and particle formation and growth. The model describes and explains the effect of the antisolvent addition rate on the average particle size, and on the particle size distribution. The developed model was constructed under the assumption of instantaneous phase equilibrium of the vapor and liquid phases on antisolvent addition (i.e., negligible mass-transfer resistance). Their results showed that the model correctly predicts the variation of particle size and particle size distribution with the main operating parameter, the antisolvent addition rate. Their findings also demonstrated the possibility of adjusting the antisolvent addition rate in accordance with the final product specifications. Elvassore et al.17 proposed a population balance model to study the kinetics involved in the GAS process. Their developed model accounts for particle nucleation, growth, aggregation, and settling with nucleation and growth being represented by the McCabe model. Bakhbakhi et al.18 used the approach of Muhrer et al. for phenanthrene, obtaining good quantitative agreement. The objective of this work is to provide a theoretical framework for our experimental results19 for the crystallization of beclomethasone-17,21-dipropionate (BDP), which is an antiinflammatory corticosteroid. BDP is an active pharmaceutical
10.1021/ie070758s CCC: $37.00 © 2007 American Chemical Society Published on Web 10/19/2007
8010
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
ingredient (API) of high value, because of its role in treating asthma.20 A mathematical model is developed to describe the elementary phenomena involved in the gas antisolvent crystallization process, i.e., the thermodynamics and the kinetics of crystal growth that governs the process. A model of the GAS crystallization process is desirable, as it will yield a better understanding of how the particles form, and so can also be used to optimize the API particle size distributions for inhalation therapy.
calculated for the binary solvent-antisolvent system. It is important to note that this method still requires an EOS to calculate the molar volumes for the solvent and CO2. The PMVF approach represents the contribution of the solvent to the total volume of the liquid phase and, thus, is related to the solvent power of that phase. Therefore, it has been proposed that the solubility of the solute is proportional to the relative partial molar volume fraction (RPMVF), which is the PMVF of the binary solution at a particular pressure/composition to the initial PMVF of the binary system.24
2. Theory The precipitation of solids in the GAS process is governed by two general phenomena: phase equilibrium and kinetics. A good understanding of the phase equilibrium underlying the GAS process can help one select the appropriate solvents, operating temperatures, pressures, etc. before performing any actual GAS experiments. Therefore, a thermodynamic model of the volume expansion process is first developed and compared to several model solute cases to determine applicability. In addition, the developed kinetic model requires information about the equilibrium concentration of solute in the liquid phase, as well as the volume of the liquid phase, as a function of antisolvent concentration. To understand the GAS crystallization procedure, and thus be better able to optimize process conditions to produce particles of a particular size and shape, it is important to develop a model for the crystallization kinetics. The ultimate goal of this model is to predict the final particle size distribution after GAS processing accurately, based on system-dependent physical constants and the initial experimental conditions. 2.1. Volume Expansion Model. In the GAS process, the supersaturation of the solute is created by dissolving a pressurized gas (i.e., CO2) as an antisolvent into the liquid-phase solvent. A volume expansion model derived from classical thermodynamic principles, such as a cubic equation of state (EOS) with mixing rules to describe fugacity and then phase equilibrium, is possible but is limited, because of the requirement of needing to know the critical constants of all three components (i.e., solute, solvent, and antisolvent). Although these constants are known for almost all solvents and antisolvents, they are rarely known for the solid pharmaceutical and were not available for BDP. In fact, quite often, the solid critical point has no physical meaning, because the compound will decompose before reaching the critical temperature. Although there are ways to evaluate/estimate physical constants, such as group contribution methods, these will not always provide physically realistic results for phase equilibrium. Precipitation in the GAS process occurs because the CO2 dissolves into the liquid phase, expanding this phase and, thus, causing a decrease of the solvents solvating power. Therefore, it seems intuitive that a description of the volume expansion process should be able to account for the drop in saturation of the solute in the liquid phase. Several attempts have been made to describe the relative volume expansion (in terms of the change in total volume of the liquid phase, divided by the initial volume of the liquid phase),21 the ratio of the current and initial molar volumes of the liquid phase,22 and the partial molar volume fraction (PMVF),23 the latter of which is described as
PMVF(T,P,Xs) )
X h sVL,s VL
(1)
where Xs is the mole fraction of solvent on a solute-free basis, VjL,s the partial molar volume of the solvent in the liquid phase, and VL the molar volume of the liquid phase, all of which are
xP )
PMVF(T,P,Xs) PMVF(T,P°,XSO)
x°P
(2)
where x°P is the initial mole fraction of solute dissolved in the liquid phase of the ternary system. This method does not incorporate the solid into the molar volume calculations; therefore, it does not require any of the physical constants for the solute that are needed for the expanded liquid-phase model. Also, because the solid is not considered strictly in the phase equilibrium, no special considerations must be taken with the mixing rules. This is especially advantageous in the case of BDP, for which there is no three-component phase equilibrium data with CO2-acetone mixtures, and which does not have welldefined critical properties (acentric factor, etc.). Implementation of the RPMVF model is a straightforward exercise in calculating the partial molar volumes, with the solute ignored in the equilibrium calculations. The fluid phase mixture was described using the van der Waals (vdW) mixing rules (quadratic for the a and b parameters), whereas the partial molar volume can be obtained from the Peng-Robinson EOS:
P)
a RT V - b V2 + 2bV - b2 Vji )
0)
( ) ∂(nV) ni
RT[(V - b) - (Vji - bhi)] (V - b)
2
(3) (4)
T,P,nj*ni
aji
+ (V + 2bV - b2) 2a[V(Vji + bhi) + b(Vji - bhi)]
-
2
(V2 + 2bV - b2)2
(5)
where
bhi ) aji )
( ) ∂(nb) ni
(6)
T,P,nj*ni
( )
2 1 ∂(n a) n ni
(7)
T,P,nj*ni
The RPMVF model was applied to phenanthrene and naphthalene in toluene, as shown in Figure 1. In both cases, the model was able to predict the mole fractions of all three components reasonably well. It is important to note that this method will predict negative solid mole fractions at high pressures, because the partial molar volume of the solvent becomes negative. Therefore, the “cutoff” for the model is the point at which the solute mole fraction becomes negative (see Figure 2). However, Figure 1 clearly shows that, up to this point, the model is able to handle the phase equilibrium reasonably well. The RPMVF model was then applied to the cholesterolacetone-CO2 system, because cholesterol is a large complex
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8011
Figure 1. Liquid-phase composition in three-component systems at 298 K, using the RPMVF model: (a) the phenanthrene-toluene-CO2 system and (b) the naphthalene-toluene-CO2 system. The points are experimental data (from ref 32), whereas the lines are from the model.
Figure 3. (a) Liquid-phase composition in the cholesterol-acetone-CO2 system at 308 K, using the RPMVF model. (b) Liquid cholesterol mole fraction at 308 K, using RPMVF model. The points are experimental data (from ref 33), whereas the lines are from the model.
Figure 4. Effect of CO2 addition rate on solute solubility at 298 K.
Figure 2. Acetone partial molar volume versus the CO2 mole fraction.
steroid molecule, similar in structure to the steroid BDP studied in this work, and there is existing data for this system. The results of the model for this system are shown in Figure 3. It is apparent that the RPMVF model was able to describe the equilibrium of all three components adequately. In particular, the solid predictions are very good, indicating that this method is able to work both with small and simple compounds as well as larger, more-complex compounds, and give physically meaningful behavior for equilibrium as the liquid phase expands. The RPMVF model is justified for application to the GAS process on the assumption that the volumetric expansion drives the increase in supersaturation. Therefore, this approach can be considered to be semiempirical, because it does not have a rigorous thermodynamic background but is based on the current understanding of the GAS process. However, there is no denying the usefulness of this approach or its apparent success in predicting three-phase equilibrium. Also, it is important to stress that this model does not consider the equilibrium of the solid phase and, therefore, is not a true equilibrium model for the
system, although it does consider the solvent-CO2 equilibrium. Also, it is interesting to note that eq 2 relates the solute mole fraction linearly to the solvent’s contribution to the molar volume. This implies that CO2’s role in the GAS process is not related directly to the antisolvent/solid interactions, but that it simply dilutes the liquid phase. After the phase equilibrium model was completed, it was possible to generate the profile of solute saturation versus time by coupling the equilibrium data to a material balance on the system.25 This allowed preliminary results on how the various physical process conditions would affect the particle size to be obtained. The effect of flow rate on the saturated solute concentration at 25 °C was examined, as shown in Figure 4. Note that the trend observed for a flow rate of 1 g/min is the same as that observed in the others, but that ∼5000 s were needed to approach zero, and so the graph was truncated. At a low flow rate, the decrease in mole fraction is quite gradual, thereby providing time for nuclei to grow. However, as the flow rate is increased, the mole fraction decreases much more quickly. In these cases, the supersaturation will reach higher levels, because the actual amount of solute dissolved cannot decrease fast enough. Therefore, under these conditions, one would expect much smaller particles, because they have a shorter time in which to grow, and this was the experimentally observed trend for BDP, as shown previously.26 2.2. Crystallization Model of the GAS Process. A population balance (i.e., a mass balance over an infinitely small subregion) is commonly used in crystallization kinetics.27 This type
8012
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
of model can be used as a predictive tool, as opposed to empirical models, which are often limited to a specific application. The general form of the one-dimensional population balance is as follows:
∂n ∂n n d(NLVL) +G + ) B(L) - D(L) ∂t ∂L NLVL dt
(8)
where n ) n(t,L) is the distribution function (number density) of the particles, L the size of the particles, G the growth rate of the particles (which is assumed to be independent of particle size), NL the liquid molar holdup, VL the molar volume of the liquid phase in the vessel, B(L) the crystal birth rate, and D(L) the crystal death rate. Note that the volume of the GAS crystallizer in this study was taken to be the liquid phase only (i.e., the vapor phase was ignored, because the liquid phase is where the precipitation and crystal formation occurs). This means that the volume of interest is not constant, despite the crystallizations being performed in a rigid vessel. The birth and death terms generally are used to account for any effect that creates unique particles or destroys existing particles, such as nucleation, aggregation, and breakage. However, if the nuclei are assumed to form with a zero radius, then nucleation can be incorporated into the boundary condition for eq 8.25 Agglomeration and breakage are somewhat more complicated, because they inherently involve both the birth and death of particles with a nonzero size (i.e., agglomeration involves the death of two small particles to create one larger particle). However, as described elsewhere,19 the CSDs measured by lowangle laser light diffraction were treated with ultrasonic dispersion to minimize aggregation (we defined this as weak interactions between particles), and the scanning electron microscopy (SEM) results showed that particle breakage did not seem to have a significant role. Therefore, we have ignored the birth and death terms, so the population balance takes on the following simplified form:
∂n n d(NLVL) ∂n +G + )0 ∂t ∂L NLVL dt
(9)
We are now able to start defining the different terms in the population balance. First, it has been shown that the growth rate can be considered to be size-independent, while maintaining good agreement with the experimental results for naphthalene.28 Based on this observation, we will use the same relatively simple expression for the crystal growth as that used in Bakhbakhi et al.18 and Muhrer et al.:25
G)
dL ) kg(S - 1)g dt
(for S > 1)
(10)
where kg and g both are system-dependent empirical constants, and S is the supersaturation. The volume derivative in eq 9 was calculated using the results from the volume expansion model previously discussed. There is no single strict definition for the supersaturation. Any form that is able to describe the amount of solute dissolved, relative to the amount of solute dissolved in a saturated solution under identical conditions, seems to suffice. As such, the supersaturation S is often defined as the difference between current and saturated concentration (S ) ∆C/Csat ) (C - Csat)/ Csat), as the ratio of fugacities between current and saturated conditions (S ) ˆf LP /fˆ L,sat P ), or simply as a ratio of mole fractions:27
S)
xP
(11)
xsat P
Equation 11 was used in this work primarily for the sake of simplicity of calculation, because the equilibrium solid mole fractions are calculated directly using the RPMVF model, as discussed previously. We now require one boundary and one initial condition to solve the population balance. First, we can assume that all of the solid in the system is initially dissolved, i.e., that
n(0,L) ) 0
(12)
For our boundary condition, we will make the assumption that the nuclei are formed with an initial diameter of zero, as stated previously. Therefore, the nucleation rate can be considered as the rate of change of the particle distribution, with respect to time, at a size of zero, i.e.,
B)
(dndt )
) Lf0
(dLdt dLdn)
(13)
Lf0
where B is the nucleation rate. However, the first derivative in the brackets is defined as the growth rate G (see eq 10), which is independent of size, so that the limit L f 0 does not require special consideration. The second derivative represents the particles that are formed at zero length and, therefore, is the distribution of the nuclei, n(t,0). When these observations are combined, we get our boundary condition:27
n(t,0) )
B G
(14)
where B is defined by the following set of equations:25
B) B′ ) 1.5DAB(cpNA)7/3
{
B′ + B′′ (if S > 1) (if S e 1) 0
() [
(15)
( )(
16π γ 3 Vp γ Vp exp kT NA 3 kT NA ln S
x
)] 2
(16) B′ and B′′ are the rates of primary and secondary nucleation, respectively. In this work, two different equations for secondary nucleation will be evaluated. The first (eq 17) has a theoretical basis but has not seen widespread use, whereas the second (eq 18) is entirely empirical but is commonly used to describe secondary nucleation:27
B′′ )
R′′aVDAB dM4
[ ( )]
2 -π γdM exp ln S kT
B′′ ) ksMTj(S - 1)i
2
(17) (18)
where DAB is the diffusion coefficient, cp is the particle concentration in the liquid phase, NA is Avogadro’s number, γ is the solid-liquid interfacial tension, k is the Boltzmann constant, T is the temperature, Vp is the molar volume of the solid particles, aV is a volumetric shape factor, dM is the molecular diameter, R′′ is an empirical constant (which determines the order of magnitude of secondary nucleation in eq 17; ks does the same in eq 18), MT is the suspension density, and i and j are empirical exponents. Several of these variables can be further defined as follows:25
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8013
DAB )
kT 2πηdM
(19)
be obtained by a simple mass balance:
xp )
where
dM )
() Vp NA
(20)
∫0∞ nLi dL
(21)
xp VL
(22)
cLp )
where η is the dynamic viscosity of the liquid phase, ka a surface shape factor (ka ) π for spherical molecules), mi the ith moment of the population density function n(t,L), and xp the mole fraction of solid dissolved in the liquid phase. The solid-liquid interfacial tension is an important parameter in crystallization, because the particle that is forming must overcome this energetic barrier. However, the experimental determination of the interfacial tension, while possible, presents a significant challenge, because the solid is soluble in the liquid, so traditional methods are often not applicable. Mersmann29 developed a theoretical approach to calculate the interfacial tension of binary solid-liquid systems based on a free energy balance over the interface, which leads to the equation
γ ) 0.414kT(cSp NA)2/3 ln
() cSp
cLp
(23)
where k is the Boltzmann constant, T is the temperature, NA is Avogadro’s number, and cSp and cLp are the (molar) solute concentration in the solid and liquid phases, respectively. It is important to note that this equation accounts for the change in interfacial tension as the system supersaturation increases, whereas the approach taken by Muhrer25 does not. Because eq 23 was developed purely from theory, it should be able to represent the interfacial tension in most systems with reasonable accuracy; it was tested in 58 different binary systems, showing good agreement in most cases.29 Therefore, this method was used to estimate the interfacial tension of BDP. Finally, to determine the supersaturation, the amount of solute currently in solution is required. To determine this, we must know how many moles of particles have precipitated in a given time. The volume of the precipitated particles can be determined using the third moment of the particle density distribution (n), along with a shape factor. However, the third moment describes the volume of particles per unit volume of the liquid phase. Therefore, this moment must be multiplied by the total liquid volume to obtain the total crystal volume:30
Vcrystal ) kVm3NLVL
(24)
where kV is a volumetric shape factor (kV ) π/6 for spherical crystals), and Vcrystal is the volume of the crystals. This can then be converted to the number of moles of particles by simply dividing by the molar volume of the solid, i.e.,
Np )
(26)
1/3
aV ) kam2 mi )
N°p - Np NL
kVm3NLVL Vp
(25)
Therefore, the mole fraction of solute remaining in solution can
As described below, the commercial software package FEMLAB was used to solve the population balance. However, to solve the population balance, the volume of the liquid phase and supersaturation must be known as a function of time. Therefore, a thermodynamic model of the equilibrium must also be performed. Although the crystallization model and the thermodynamic model are both inherently linked, it has been previously shown that the thermodynamic model with the material balance can be decoupled from the crystallization model without significant error.25 Therefore, the mass balance was solved independently of the population balance in this work. 3. Results and Discussion 3.1. Model Implementation. The population balance (eq 9) was implemented using the commercial software package FEMLAB, which uses the finite-element method to solve partial differential equations. The population balance was solved for particle sizes ranging from 1 nm to 1 mm. A dirichlet boundary condition was used, based on eq 14, at the 1 nm boundary, which implied that particles are formed with “zero” size, considering that the sizes of the actual particles formed were in the range of 1-100 µm. A cutoff point of 1 mm was chosen arbitrarily, because only a diameter that is significantly larger than the maximum observed particle size in the system is required. The phase equilibrium model/material balance was, as previously mentioned, solved independently from the crystallization, and the results were used as input variables for the crystallization model at a given time. The saturated mole fraction, the liquid volume in the precipitation vessel, the number of moles in the liquid phase, and the volume derivative in eq 9, were all input in this manner. The size of the mesh was chosen to balance computational time with numerical accuracy. This was done in two ways. First, particle size distributions normally scale logarithmically, with a large fraction occurring at smaller sizes (i.e., 1 µm) as opposed to bigger sizes (i.e., 1 mm). Therefore, the mesh spacing was scaled to be very fine at the left end (1 nm), and much larger at the right end (1 mm). This allowed good numerical accuracy to be achieved without an excessive number of elements. Second, the simulation was run using successively finer meshes to determine acceptable values for the mesh growth rate and minimum mesh spacing. The solution time was normally 5-20 min, depending on the mesh spacing used, and little change in the predicted particle size distribution was observed when the mesh spacing/growth rate were further decreased. Therefore, a minimum mesh spacing of 1 nm and a growth rate of 1.05 were used, giving 226 mesh elements over the region from 1 nm to 1 mm. 3.2. Model Results. In this work, the crystallization kinetics were described by three kinetic terms: primary nucleation, secondary nucleation, and size-independent growth. Primary nucleation is described by eq 16, which is a standard equation based on the free-energy change when a crystal forms out of solution and is fully described by the physical properties of the system being studied. The secondary nucleation was described by two equations separately (eqs 17 and 18). Equation 17 was developed based on the growth of a crystal on a smooth surface multiplied by the probability that this second crystal will
8014
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
Table 1. Summary of Experimental Conditions Studied in Crystallization Modeling
Table 3. Summary of Model Results Using Eq 17 dp(10%) (µm)
run
temperature (°C)
agitation (rpm)
flow rate (mL/min)
F1 F2 F3
25 25 25
1000 1000 1000
25 50 75
run
ln(R′′)
kg
g
AAD
F1 F2
-33.79 -33.70
2.645 × 10-6 2.574 × 10-6
2.184 2.170
0.476 0.500
R′′),31
dislodge (i.e., the parameter whereas eq 18 is empirical in nature, but is commonly used.27 Table 1 provides the experimental conditions used in the crystallization modeling of the GAS process. Equation 17 was used to model secondary nucleation, and the three constants (R′′, kg, g) were regressed from experimental data using Matlab; these are provided in Table 2. The average absolute deviation (AAD) was used instead of the AARD, because the model could not accurately predict the volume percentage when it approached zero, which inflated the AARD.
1 N
N
∑ | ycalc - yexp|
dp(90%) (µm)
F1 data model
18.97 26.47
47.28 46.63
86.64 79.64
% difference
39.54%
1.37%
8.09%
data model
17.22 22.18
36.51 39.72
68.14 68.18
% difference
28.76%
8.77%
0.05%
F2
Table 2. Regressed Parameters for Growth and Secondary Nucleation Using Eq 17
AAD )
dp(50%) (µm)
(27)
Therefore, the AAD was given so the error on small volume fractions would not skew the final result significantly. Also, note that the average was taken with respect to only the data points that had a moderate volume percentage (i.e., >0.01%). This method provided a more-accurate representation of the error. For simplicity, the effects of agitation rate and temperature were ignored in the crystallization model, and only the effect of flow rate was considered. As was shown elsewhere,19 at 25 °C, the dp(10%), dp(50%), and dp(90%) metrics were quite
Table 4. Regressed Parameters for Growth and Secondary Nucleation Using Eq 18 run F1 F2 F3
ln(ks) 21.32 21.99 24.13
i 3.552 3.370 3.713
j 1.157 1.201 1.199
kg 10-6
5.651 × 3.836 × 10-6 1.5734 × 10-5
g
AAD
1.910 1.855 1.594
0.224 0.402 0.374
similar to the point where no real difference could be inferred in most cases at 1000 RPM. The effect of temperature was also small when compared to that of flow rate, when considering the variation in experimental values, although there is possibly a more pronounced effect for dp(10%) at lower flow rates. The effect of flow rate was the most important for determining the size distribution of particles under the conditions studied in the previous experimental work, and, as such, only runs that varied the flow rate at constant temperature and an agitation rate of 1000 RPM were used for regressions. When eq 17 was used, the developed GAS model was able to simulate the experimental CSD curves very well, as shown in Figure 5. However, it overestimated the height of the primary mode by more than 2 vol % in each case for F1 and F2, which amounts to a relative error of ∼20%. Table 3 provides a summary of the model results for the dp(10%), dp(50%), and dp(90%) metrics of the experimental CSDs. These metrics, which represent the 10th [dp(10%)], 50th [dp(50%)], and 90th [dp(90%)] percentiles of the cumulative volume (mass)-weighted
Figure 5. Experimental results versus model results for the crystallization kinetics using eq 17 for flow rates of (a) 25 and (b) 50 mL/min. The points are experimental data (from ref 19), whereas the lines are from the model.
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8015
Figure 6. Comparison of experimental and predicted particle size distributions at 25 °C, 1000 RPM, and flow rates of (a) 25, (b) 50, and (c) 75 mL/min. The points are experimental data (from ref 19), whereas the lines are from the model.
size distributions, were chosen to represent the fine, central, and coarse portions, respectively, of each CSD. Increased nucleation would be expected to increase the proportion of smaller particles, therefore primarily affecting dp(10%), while increased agglomeration would be anticipated to affect the magnitude of dp(90%), which represents the coarse end of the distribution. The crystallization model was able to give excellent estimates of the larger particles, as shown by low-percentage differences between the data and model for the dp(50%) and dp(90%) metrics. However, it gave a reasonably high difference for dp(10%). Equation 18 was also used to model secondary nucleation, and the five constants (ks, i, j, kg, g) were regressed from experimental data using Matlab. The constants are given in Table 4, along with the associated AAD values. As shown in Figure 6, the developed GAS model was able to represent the experimental CSDs reasonably well. Table 5 provides a summary of the model results using eq 18 for the dp(10%), dp(50%), and dp(90%) metrics of the experimental CSDs. This crystallization model also gave excellent estimates of the larger particles, as shown by low-percent differences between the data and model for the dp(50%) and dp(90%) metrics. For runs
Table 5. Summary of Model Results Using Eq 18 dp(10%) (µm)
dp(50%) (µm)
dp(90%) (µm)
F1 data model
18.97 23.56
47.28 47.62
86.64 88.16
% difference
24.19%
0.72%
1.76%
data model
17.22 17.20
36.51 40.59
68.14 73.05
% difference
0.12%
11.17%
7.19%
data model
9.76 15.05
32.84 33.14
63.35 63.82
% difference
54.24%
0.93%
0.75%
F2
F3
F1 and F2, it also gave a better representation of the dp(10%) metric than the GAS model from eq 17, which is not surprising, considering the increase in regression parameters (five versus three). The effect of agitation, while important, was not considered at this time. The current model may be able to handle this effect implicitly through the fitting parameters, and the suitability of the model for this task will be evaluated in the future.
8016
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
While accounting for secondary nucleation with either model studied allowed for a reasonably accurate description of the upper two metrics of the CSD, it provided a poorer representation of dp(10%), which may be related to the sphericity of the particles. Laser diffraction calculates the spherical volume equivalent diameter and, thus, is able to construct a volumeweighted distribution. However, this can be problematic when the particles are not spherical, as was the case with the formed BDP particles,19,26 where particles occasionally grew from a central nucleus. 4. Conclusions To describe pressurization during the gas-antisolvent (GAS) process, a mass balance and a phase equilibrium model were required. A semiempirical phase equilibrium model (the relative partial molar volume fraction (RPMVF) model) was evaluated that uses the binary two-phase solvent/CO2 phase equilibrium and then relates it to the solid concentration in the three-phase system (BDP-acetone-CO2). The RPMVF model was tested successfully with the phenanthrene-toluene-CO2 model system, the naphthalene-toluene-CO2 model system, and the more-complex cholesterol-acetone-CO2 model systems before predictions were made for the system of interest in this work. A population balance was used to model the experimentally determined crystal size distributions (CSDs). Two models were used independently: an empirical equation, which is commonly used to model secondary nucleation, and a theory-based equation. Both models were determined to give reasonable estimates of the studied metrics of the crystal size distributions dp(10%), dp(50%), and dp(90%), although they predicted the higher end of the CSDs (i.e., dp(50%) and dp(90%)) better. Notation A ) dimensionless parameter in the Peng-Robinson equation of state (PR-EOS) a ) parameter in the PR-EOS (J m3/mol2) B ) dimensionless parameter in the PR-EOS b ) parameter in the PR-EOS (m3/mol) ˆf Ri ) fugacity of component i in phase R (Pa) Hfusion ) heat of fusion for the solute (J/g) kij, lij ) interaction parameters in the PR-EOS NA ) Avogadro’s number (mol-1) Ni ) molar holdup in phase i (mol) P ) pressure (bar) Q˙ a ) CO2 flow rate (mol/s) R ) ideal gas constant (m3 Pa/(mol K)) T ) temperature (K) t ) time (s) Vi ) molar volume of phase i (mol/m3) xi ) liquid-phase mole fraction of component i yi ) vapor-phase mole fraction of component i Z ) compressibility factor Greek Letters R ) parameter in the PR-EOS γ°3 ) activity coefficient for the solute under the initial conditions δ ) solubility parameter ((cal/cm3)1/2) φRi ) fugacity coefficient of component i in phase R ω ) acentric factor Superscripts and Subscripts ° ) initial condition
A ) antisolvent atm ) atmospheric pressure c ) critical property i ) dummy variable L ) liquid phase P ) solute S ) solvent V ) vapor phase R ) dummy variable Acknowledgment The authors acknowledge the financial support of the Natural Sciences and Engineering Council (NSERC) of Canada, and the Ontario Premiers Research Excellence Award (PREA). Literature Cited (1) Cansell, F.; Chevalier, B.; Demourgues, A.; Etourneau, J.; Even, C.; Garrabos, Y.; Pessey, V.; Petit, S.; Tressaud, A.; Weill, F. Supercritical Fluid Processing: a new route for materials synthesis. J. Mater. Chem. 1999, 9, 67-75. (2) Cooper, A. I. Recent Developments in Materials Synthesis and Processing Using Supercritical CO2. AdV. Mater. 2001, 13 (14), 11111114. (3) Subramaniam, B.; Rajewski, R. A.; Snavely, W. K. Pharmaceutical Processing with Supercritical Carbon Dioxide. J. Pharm. Sci. 1997, 86 (8), 885-890. (4) Woods, H. M.; Silva, M. C.; Nouvel, C.; Shakesheff, K. M.; Howdle, S. M. Materials processing in supercritical carbon dioxide: surfactants, polymers and biomaterials. J. Mater. Chem. 2004, 14, 1663-1678. (5) Ye, X.; Wai, C. M. Making Nanomaterials in Supercritical Fluids: A Review. J. Chem. Ed. 2003, 80 (2), 198-203. (6) Tan, H. S.; Borsadia, S. Particle Formation Using Supercritical Fluids: Pharmaceutical Applications. Expert Opin. Ther. Pat. 2001, 11 (5), 861-872. (7) Chow, A. H. L.; Tong, H. Y.; Chattopadhyay, P.; Shekunov, B. Y. Particle Engineering for Pulmonary Drug Delivery. Pharm. Res. 2007, 24 (3), 411-437. (8) DeSimone, J. M.; Tumas, W. Green Chemistry Using Liquid and Supercritical Carbon Dioxide; Oxford University Press: New York, 2003; p 259. (9) Yu, B.; Shekunov, P.; York, P. Crystallization processes in pharmaceutical technology and drug delivery design. J. Cryst. Growth 2000, 211, 122-136. (10) York, P. Strategies for particle design using supercritical fluid technologies. Pharm. Sci. Technol. Today 1999, 2 (11), 430-440. (11) Jung, J.; Perrut, M. Particle Design Using Supercritical Fluids: Literature and Patent survey. J. Supercrit. Fluids 2001, 20, 179-219. (12) Perrut, M. Supercritical Fluid Applications: Industrial Developments and Economic Issues. Ind. Eng. Chem. Res. 2000, 39, 4532-4535. (13) Gallagher, P. M. C., M. P.; Krukonis, V. J.; Klasutis, N. Gas antisolvent recrystallization: new process to recrystallize compounds insoluble in supercritical fluids. In Supercritical Fluid Science and Technology; Johnson, K. P., Penninger, J. M. L., Eds.; ACS Symposium Series 406; American Chemical Society: Washington, DC, 1989; pp 334354. (14) Gallagher, P. M.; Krukonis, V. J.; Botsaris, G. D. Gas anti-solvent recrystallization: Application to particle design. In Particle Design Via Crystallization; Ramanarayanan, R., Kern, W., Larson, M., Sidkar, S., Eds.; AIChE Symposium Series 284; American Institute of Chemical Engineers: New York, 1991; pp 96-103. (15) Muller, M.; Meier, U.; Kessler, A.; Mazzotti, M. Experimental study of the effect of process parameters in the recrystallization of an organic compound using compressed carbon dioxide as anti-solvent. Ind. Eng. Chem. Res. 2000, 39, 2260-2268. (16) Muhrer, G.; Lin, C.; Mazzotti, M. Modeling the gas antisolvent recrystallization process. Ind. Eng. Chem. Res. 2002, 41, 3566-3579. (17) Elvassore, N.; Parton, T.; Bertucco, A. Kinetics of Particle Formation in the Gas Antisolvent Precipitation Process. AIChE J. 2003, 49, 859-868. (18) 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 (19), 7345-7351.
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8017 (19) Dodds, S. Gas Antisolvent Micronization of Pharmaceutical Powders. MESc. Thesis, University of Western Ontario, London, Ontario, Canada, 2006. (20) Expert Panel Report: Guidelines for the Diagnosis and Management of Asthma, Update on Selected Topics 2002; National Asthma Education and Prevention Program, Bethesda, MD, 2003. (21) Kordikowski, A.; Schenk, A. P.; Van, Nielen, R. M.; Peters, C. Volume Expansions and Vapor-Liquid Equilibria of Binary Mixtures of a Variety of Polar Solvents and Certain Near-Critical Solvents. J. Supercrit. Fluids 1995, 8, 205-216. (22) de la Fuente, J.; Peters, C.; de Swaan Arons, J. Volume Expansion in Relation to the Gas-Antisolvent Process. J. Supercrit. Fluids 2000, 17, 13-23. (23) Mukhopadhyay, M. Partial Molar Volume Reduction of Solvent for Solute Crystallization Using Carbon Dioxide as Antisolvent. J. Supercrit. Fluids 2003, 25, 213-223. (24) Mukhopadhyay, M.; Dalvi, S. Partial Molar Volume Fraction of Solvent in Binary (CO2-Solvent) Solution for Solid Solubility Predictions. J. Supercrit. Fluids 2004, 29, 221-230. (25) Muhrer, G.; Lin, C.; Mazzotti, M. Modeling the Gas Antisolvent Recrystallization Process. Ind. Eng. Chem. Res. 2002, 41, 3566-3579. (26) Bakhbakhi, Y.; Charpentier, P. A.; Rohani, S. Experimental study of the GAS process for producing microparticles of beclomethasone17,21-dipropionate suitable for pulmonary delivery. Int. J. Pharm. 2006, 309 (1-2), 71.
(27) Randolph, A.; Larson, M. Theory of Particulate Processes, Second Edition; Academic Press: San Diego, CA, 1988. (28) Bakhbakhi, Y. Supercritical Crystallization of Pharmaceuticals: The Gas Antisolvent Process; Ph.D. Thesis, University of Western Ontario, London, Ontario, Canada, 2004. (29) Mersmann, A. Calculation of Interfacial Tensions. J. Cryst. Growth 1990, 102, 841-847. (30) Grosch, R.; Briesen, H. Getting Started with ParsiVal; RWTHAachen University: Aachen, Germany, 2004. (31) Worlitschek, J.; Mazzotti, M. Model-Based Optimization of Particle Size Distribution in Batch-Cooling Crystallization of Paracetamol. Cryst. Growth Des. 2004, 4, 891-903. (32) Dixon, D.; Johnston, K. Molecular Thermodynamics of Solubilities in Gas Antisolvent Crystallization. AIChE J. 1991, 37, 14411449. (33) Liu, Z.; Wang, J.; Song, L.; Yang, G.; Han, B. Study on the Phase Behavior of Cholesterol-Aceton-CO2 System and Recrystallization of Cholesterol by Antisolvent CO2. J. Supercrit. Fluids 2002, 24, 1-6.
ReceiVed for reView May 30, 2007 ReVised manuscript receiVed August 23, 2007 Accepted August 29, 2007 IE070758S
