A complete 3D simulation of a crystallization process induced by supercritical CO2 to predict particle size

Computers and Chemical Engineering 52 (2013) 1–9

Contents lists available at SciVerse ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

A complete 3D simulation of a crystallization process induced by supercritical CO2 to predict particle size Arnaud Erriguible a,∗ , Tarik Fadli a , Pascale Subra-Paternault a,b a b

Université de Bordeaux, Institut Polytechnique de Bordeaux (IPB), I2M UMR 5295, Site ENSCBP, 16 Avenue Pey-Berland, 33607 Pessac Cedex, France Université de Bordeaux, Institut Polytechnique de Bordeaux (IPB), CBMN UMR 5248, Allée Geoffroy Saint Hilaire, Bat B14bis, 33600 Pessac Cedex, France

a r t i c l e

i n f o

Article history: Received 19 July 2012 Received in revised form 11 December 2012 Accepted 11 December 2012 Available online 10 January 2013 Keywords: CFD Crystallization Supercritical fluids Large Eddy Simulation

a b s t r a c t Crystallization induced by compressed CO2 is a process that operates under several MPa of pressure. By rendering on-line measurements very difficult to perform, simulation appears as a suitable tool to better identify the important parameters of the process. A mathematical model is developed in the case of a spray-crystallization process in which a solution of minocycline-ethanol is injected into carbon dioxide as antisolvent. The model accounts for the main physical phenomena involved, i.e. hydrodynamics, mass transfer, phase equilibrium, crystallization kinetics. Simulations are performed in 3D with a special insight in turbulence modeling. Numerical results are compared with experimental data from literature. Although experimental and simulated PSD fit satisfactorily, results emphasize the major role of the crystal–fluid interfacial tension on the accuracy. Numerical investigations are further performed to highlight the effects of injection velocity and solution concentration on the spatial distribution of the important variables in crystallization. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction Green processes are nowadays more and more considered in industries particularly due to attitudes change and legislation that are increasingly demanding for environmental friendliness. The use of technology “supercritical fluids” in industry is directly linked to these concepts. Crystallization using supercritical fluids, and more specially CO2 was mostly developed as an alternative approach to conventional liquid crystallization, challenging for a better control of the particle size and a lower concentration of residual solvent in drugs. Beside the rapid expansion of a supercritical solution (RESS process) that is limited in terms of market prospectives, the use of CO2 as antisolvent offers a solution for recrystallizing a large spectra of molecules. In this technique, solid particles are produced from a mother liquor upon addition of CO2 . Despite the use of organic solvent, a major environmental gain for CO2 versus liquid antisolvent resides in the design of a one-stage process that minimizes energy and environmental impacts by saving the post-processing of produced powders. In the specific antisolvent variant investigated in this work, the so-called SAS process, the compound to recrystallize is initially dissolved in an organic solvent before being injected through a nozzle into a stream of compressed CO2 . Upon the mixing of the two fluids,

∗ Corresponding author. E-mail address: [email protected] (A. Erriguible). 0098-1354/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compchemeng.2012.12.002

the solubility of the solute decreases and the precipitation is initiated. The final particles characteristics like the mean size and the size distribution are influenced by both thermodynamic and kinetic factors, and when precipitation is induced by the addition of a third species like in this process, the time and space scales at which the mixing occurs become important for the control of the particle size. Literature provides numerous examples of recrystallization induced by the addition of CO2 to an organic liquor. Although the effect of parameters on particle characteristics is evidenced for each compound, contradictory trends are obtained (Fadli, Erriguible, Laugier, & Subra-Paternault, 2010) when the whole set of systems is considered, especially with pressure and temperature as operating variable. The challenge for optimizing a SAS recrystallization comes from the large number of operating parameters and their action on the several phenomena involved in the process, as mass transfer, phase equilibria and nucleation/growth of the species to be precipitated. Although an extensive investigation carried out for given species could provide experimental trends, it is highly difficult to rationalize the results and to identify the main mechanisms because of the interactions between parameters. A numerical approach hence becomes fundamental to clarify the contribution of the phenomena and to predict the effect of a specific parameter onto the process performance. A complete description of the Supercritical Antisolvent Process should be able to account for mass transfer related to the mixing of the solution and CO2 , for hydrodynamics related to the injection of the solution, for phase equilibria related to the evolution of the solute solubility with the

2

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9

composition of the forming CO2 –solvent mixture, and finally, for crystallization kinetics related to the nucleation of the solute and its growth as particles or crystals. Several complete simulation of the SAS process has been yet published. In their pioneer work, Martín and Cocero (2004) developed the simulation in a 2D axisymetrical configuration and used the method of moments to solve the population balance equation. The turbulence was taken into account in the hydrodynamic model via the k−ε model. Lately, Henczka, Bałdyga, and Shekunov (2005) proposed a similar approach to simulate the process in a 2D configuration, coupling the commercial CFD code FLUENT and the standard method of moments. The interaction between the mixing and precipitation was carried through the ␤-PDF approach coupled with a turbulent mixer model (Baldyga, 1989). The first 3D simulation was developed by Cardoso, Cabral, Palavra, and Geraldes (2008), using the FLUENT CFD code, who focused on the hydrodynamics of the reactor but did not describe the kinetics of the crystallization. More recently, Sierra-Pallares et al. (2012) carried out a complete numerical study by using a 2D simulation with a RANS turbulent model. A micromixing model was coupled with the hydrodynamics one in order to emphasize the role of the mixing scale on the particle size distribution. It is worthwhile to notice that all the previous studies have highlighted the importance of an accurate modeling of the micromixing in the reactor to correctly predict the sizes of the particles. Indeed, the specificity of this process is that it is intrinsically unsteady. Therefore, from a modeling point of view, an unsteady turbulent model and a 3D simulation appear to be a good combination to get an accurate description of the process. As a second approach, the 2D RANS model can be envisaged but high level micromixing models have to be added as proposed by SierraPallares et al. (2012) to describe the turbulent mixing phenomena. As described previously, the majority of the literature focused on this second approach whereas we favored the unsteady 3D-LES model to provide local and instantaneous description of the flow. Therefore, in this work, the effects of turbulence on the mixing are described by implementing a Large-Eddy Simulation model (LES) in a 3D simulation. This deterministic approach (Sagaut, 1998) relies on a scale separation of the resolved large turbulent scale and of the modeled small turbulent structures and allows for describing turbulent mixing with a better accuracy than the statistical RANS methods. DNS turbulence method would have certainly provided better results but the CPU time would have been prohibitive. The method proposed here was selected as a good compromise between precision and CPU cost. The 3D hydrodynamics model was further coupled to the crystallization kinetics in order to predict the particles size distribution. To author’s best knowledge, this work hence describes for the first time the complete simulation of the SAS process in 3D configuration with LES turbulence modeling. To model the crystallization kinetics, the population balance equation was solved by the standard method of moments assuming a monodisperse distribution. The simulation was validated by comparing the predicted sizes to experimental investigations. Although many examples of CO2 assisted recrystallization can be found in literature, the validation of the numerical model requires not only the particles size but the solubility data as well, a combination that is rarely found. Both solubility data and SAS experiments were reported by Cardoso et al. (2008) in the case of minocycline recrystallized from ethanol. Hence, this system and the related experimental conditions of the process were selected for the simulations. Once being validated, the numerical investigations were performed to highlight the effects of injection velocity and solution concentration on the spatial distribution of the important variables in crystallization, like the solute concentration, the supersaturation ratio and the number of particles.

2. Modeling 2.1. Thermophysical properties The simulation of the injection of a solution into a stream of CO2 will require the estimation of several properties of the newly formed mixtures. 2.1.1. Density The effect of pressure and composition on the fluid density is described through the Peng–Robinson equation of state with quadratic mixing rules, owing to its well-known ability to represent binary systems at high pressure and its ease of application. P=

RT

am =

 i

am v(v + bm ) + bm (v − bm )

−

v − bm

(1)

xi xj aij

j

aij = (aii ajj )0.5 (1 − kij ) bm =

 i

bij =

xi xj bij

(2)

j

(bii + bjj ) 2

(1 − lij )

with T the temperature, P, the pressure, v the molar volume, R the ideal gas constant, a and b, the attraction parameter and the Van Der Walls co-volume, respectively and kij and lij the binary interaction parameters. For the CO2 –ethanol system, kij = 0.066, lij = 0.005 (Wu, Ke, & Poliakoff, 2006). 2.1.2. Diffusion coefficient The diffusion coefficient of Ethanol, the solvent, in CO2 , the antisolvent, is calculated by the Hayduck–Minhas equation whereas the diffusion coefficient of the solute in the CO2 –EtOH mixture is estimated by the Wilke and Chang correlation (Fadli et al., 2010). 2.1.3. Viscosity The viscosity of the mixture is calculated by a simple mixing rule as proposed by Reid, Prausnitz, and Poling (1987): =



xi i

(3)

i

whilst the viscosiy of each fluid, i , is obtained by the Chung correlation (Reid et al., 1987). This ideal mixing rule for viscosity might appear as very simple, but rapid calculations showed that the deviation from more sophisticated correlations were of 10–30%, a range that is insufficient to modify significantly hydrodynamic results. 2.2. Hydrodynamics In the conditions investigated in this work, the pressure is above the critical pressure of the CO2 –EtOH mixture so the two fluids are fully miscible and merge as a monophasic flow. The flow is considered as incompressible since the variation of density is only due to the variation of the mixture composition. The equation of continuity and conservation of momentum is described by the Navier–Stokes model:

∇ · u¯ = 0





¯ ∂u ¯ ¯ · ∇u +u ∂t

(4)

 ¯ + ∇ t u] ¯ = −∇ p + g + ∇ · ( + t )[∇ u

(5)

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9

in which p is the pressure,  the density of the fluid, g the gravity vector,  the dynamic viscosity, t the turbulent viscosity, t the ¯ the filtered velocity vector. time and u The species continuity equations are expressed by taking into account the diffusion of species according to the Fick’s law: 

∂xa ¯ − (Da + Dt )∇ xa ) = 0 + ∇ · (xa u ∂t

(6)



∂xb ¯ − (Db + Dt )∇ xb ) = Sb + ∇ · (xb u ∂t

(7)

With xa and xb the mass fractions of solvent and solute, respectively, Da and Db the diffusion coefficients of solvent and solute in solvent/CO2 , respectively, and Sb a source term that represents the quantity of solute produced by nucleation and growth, as detailed hereafter. Dt is the turbulent diffusion coefficient calculated by the following relationship: Dt =

 Sct

(8)

with Sct the turbulent Schmidt number set to 0.9. The turbulence is modeled through a deterministic Large Eddy Simulation approach (Sagaut, 1998) that relies on a scale separation of the resolved large turbulent scale and of the modeled small turbulent structures. A structural mixed scale model is used for the turbulent viscosity: 1−˛ ¯ 1−˛ (2Sij Sij )˛/2  t = CS2˛ CTKE

 1 1−˛/2 2u u

2.3. Solute solubility The process under investigation involves the mixing of a liquid solvent, ethanol, with compressed CO2 . At the pressure of 13 MPa, ethanol and CO2 are fully miscible since the selected pressure is above the mixture critical pressure. Hence, there is no coexisting liquid and vapor phases, but a monophasic fluid whose composition evolves during the injection. The solubility of a solute in a mixture of CO2 and a solvent depends of temperature, pressure and of the fluid composition. In monophasic conditions and fixed T and P, the solubility only changes with the fluid composition. For the minocycline–ethanol–CO2 system considered here in conditions of 13 MPa and 313 K (Cardoso et al., 2008), the exponential decrease of the solubility with the decreasing ethanol content was satisfactorily correlated by the empirical fit for 0 < xa < 0.17: xsat = 10

−5 25,8xa

e

In addition, considering that the particles are spherical and have the same velocity than the flow, the population balance equations then write as:  

∂N ¯ = Bhom + ∇ · (N u) ∂t ∂m ¯ = + ∇ · (mu) ∂t

 6

(10)

2.4. Particle size distribution The general dynamic equation for nucleation and condensation is solved with the standard method of moments assuming that the particle size distribution follows a monodisperse distribution (Pratsinis, 1988). In this particular case, solving the original partial differential equation remains to solve only two ordinary differential equations while preserving a good solution. The number of particles and the volume of particles are deduced from only two moments of the distribution. Such method is well suited when coupling phenomena of nucleation-growth with a CFD code, because of its ease of application and its relative accuracy on the particle size prediction (Schwarzer, Schwertfirm, Manhart, Schmid, & Peukert, 2006).

(11) Bhom p xc3 +

 NGp dp2 2

 (12)

where N is the number of particles and m is the total mass of the particles. Bhom and G represent the nucleation and the growth rates, respectively. p , xc and dp are respectively the solid particle density, the critical nuclei size and the particle size. The quantity of solid solute created by nucleation and growth, that represents the source term Sb in Eq. (7), reads to: Sb = −

 6

Bhom p xc3 +

 NGp dp2 2



(13)

The supersaturation is the driving force of the precipitation process, since any value higher than unity creates potential conditions for nucleation to occur. This key parameter is defined by: S=

xb xsat

(14)

Assuming that the particles are formed by primary nucleation, the number of critical nuclei formed per unit time in a unit of volume can be calculated with the following rate expression, according to the classical nucleation theory (Mersmann, 2001):



(9)

The Smagorinski constant CS is chosen equal to 0.2, the Turbulent Kinetic Energy (TKE) constant CTKE is 0.06 and the turbulent Prandtl number is 0.85. In the mixed scale model, the ˛ parameter ¯ is the size of the LES filter, Sij is the deformation rate tensor is 0.5,  and u is the velocity fluctuation.

3

Bhom = 1.5Db (Csat SNa )

7/3

Vm exp kb T

 −



16 3 kb T

3 V 2  m ln2 (S) (15)

with Db the solute diffusion coefficient (see Section 2.1.2), Csat the solute concentration, Vm the solute molecular volume, kb the Boltzmann constant and NA the Avogadro number. The solid–fluid interfacial tension is specific to each fluid–solute system and although its range of order can be found for few compounds growing in liquid, there is no data reported for molecules in CO2 –solvent environment. Therefore, the value will be found by fitting numerical and experimental results of mean size, whilst the influence on this term will be pointed out by performing simulations with values in the range encountered in liquid crystallization. In presence of high supersaturation levels, the growth rate is defined by the classical relation: G = 2Sh

MDb Csat (S − 1) d dp

(16)

with M the molar weight of the solute and Sh the Sherwood number estimated by the Froessling equation: 1/2

Sh = 2 + 0.6Rep Sc 1/3

(17)

where Rep and Sc are the particle Reynolds number and the Schmidt number, respectively. The critical nuclei size xc and the particle size dp are deduced by the following relationship: xc =

4 Vm kb T ln S

(18)

6m Np

(19)



dp =

3

3. Results and discussion 3.1. Numerical procedure The set of partial differential equations (4)–(7), (11), and (12) is numerically solved to estimate the velocities, pressure,

4

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9

composition, number of particles and the total mass of the particles at each time step in the reactor. Other physical parameters are deduced from the solved variables at each time step. The numerical tool is the home-made CFD code “Thétis” developed at the I2M/TREFLE Department. The approximation of the model described in Section 2.1 is realized by the finite volumes method on a fixed staggered grid. The space derivatives for the inertia, conduction and viscous terms are discretized with a centered scheme (Patankar, 1990). The time discretization is implicit for all terms in the Navier–Stokes and energy equations, except for the inertial ¯ n+1 , ¯ n · ∇u term in the momentum equation that is linearized as u with n the time index corresponding to time n t. The coupling between pressure and velocity unknowns is satisfied through the following formulation according to the time-splitting algorithm of Goda (1979) used to perform incompressible simulations:





∇·

¯∗ −u ¯n u ¯∗ + u¯ n · ∇ u t

 ∇t 

¯ ∗ + ∇tu ¯ ∗] = −∇ pn + g + ∇ · ( + t )[∇ u



∇ p = ∇ · u¯ ∗

¯ n+1 = u ¯∗ − u



(20)

t ∇ p 

Fig. 1. Geometric configuration and boundary conditions.

pn+1 = pn + p

3.2. Validation and influence of the solid–fluid interfacial tension

For transport equations, a TVD (Total Variation Diminishing) scheme is used since it is particularly well suited for diffusion problems (Vincent & Caltagirone, 1999). The velocity–pressure coupling is solved by a solver BiCGSTAB II (Falgout & Yang, 2002), preconditioned with a Jacobi method for the prediction step and MCG approach [available in the bookstore HYPRE (Falgout & Yang, 2002)] for the correction step. The geometrical configuration and boundary conditions are shown in Fig. 1. The dimensions of the domain were fixed to 5 cm × 3 cm × 3 cm by consideration of the experimental merging zone between CO2 and the solution reported by Badens, Boutin, and Charbit (2005). The mesh is constituted by 10 580 000 nodes (800 × 115 × 115) and is refined near the injector (Fig. 2) to improve the description of the injection phenomenon. The injector diameter was of 125 ␮m. The refinement of mesh allows for describing finely the instantaneous instabilities due to the mixing and can be visualized through the velocity vector in Fig. 2. Because the high number of nodes required for describing with accuracy the micromixing, calculations are performed by MPI parallel programming on 512 processors.

The numerical model is validated by comparing the predicted sizes to an experimental result from literature, i.e. the recrystallization of minocycline in a CO2 /ethanol mixture (Cardoso et al., 2008). The system was selected owing to the existence of solubility data required to compute the nucleation and growth rates for every composition encountered when the solvent and CO2 are mixing. The experimental conditions were of 313 K for temperature, 13 MPa for pressure, 10 mg/ml for minocycline concentration in the initial solution and 3 ml/min for the solution flow rate. In these conditions, spherical particles of 203 nm as mean size were produced. As mentioned previously, the interfacial tension is the only unknown parameter in the equations. In liquid crystallizations performed at atmospheric pressure, typical values are ranging between 5 and 40 mN/m according to the correlation of Mersmann (1990). Simulations were carried out varying the interfacial tension within this domain, and the effect on the mean size of minocycline is shown in Fig. 3. The interfacial tension has a large influence on the particle size, since an increase of 10–20 mN/m provokes an increase of the size of more than one order of magnitude. In fact, a high

Fig. 2. Mesh refinement and velocity vector near the injector in the meridian plane of the reactor.

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9

5

Fig. 3. Evolution of particle size (nm) in function of the interfacial tension (mN/m).

interfacial tension decreases strongly the nucleation frequency so growth becomes the main contributor to the crystallization process with an increase of size as a consequence. This influence was also observed by Sierra-Pallares et al. (2012) for the ␤-carotene in a dichloromethane/CO2 mixture. Those works emphasize the major role played by the interfacial tension on the particle size, and therefore, simulations have first to be run on a fully described system to provide an estimate of this parameter before predicting the behavior for extra conditions. For the future, the estimation of the surface tension will become of major interest to consider the simulation as a predictive tool. As mentioned previously, the particle size of minocycline produced experimentally was of 203 nm. From Fig. 3, it can be seen that a value of 15 mN/m ensures a good fit with this experimental size, so that value that will be kept for the next calculations. In addition, the simulation allows for recovering the size distribution as well, that is given in Fig. 4. The main important variables for a crystallization process are shown in Fig. 5 as instantaneous spatial fields. The mixing behavior of the injected solution and the ambient CO2 can be viewed in Fig. 5a that displays the ethanol mass fraction. In those monophasic conditions, the solution and the CO2 merge without forming

Fig. 5. Instantaneous fields in the meridian plane of the reactor of several variables: (a) mass fraction of ethanol; (b) mass fraction of minocycline; (c) supersaturation ratio; (d) number of particles; and (e) mass of particles. T = 313 K, P = 13 MPa, vinj = 4 m/s.

Fig. 4. Distribution of the mean particle size in the transversal plane (x = 0.008 m).

droplets. Closed to the nozzle exit, the inlet core of the jet is mostly composed of ethanol, and as the merging between the solution and the ambient CO2 proceeds in time and in space through convection and species diffusion, the ethanol content decreases in favor of CO2 . The minocycline evolution (Fig. 5b) closely matches the ethanol field at least at proximity of the nozzle tip, but values are much lower since minocycline is in minority compared to the solvent. Moreover, since the solute is consumed by nucleation/growth,

6

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9

Fig. 6. Evolution of particle size in function of solute concentration.

the minocycline fraction that remains in solution decreases quickly from the nozzle exit so the jet of the solute appears to be much shorter than the ethanol one. Let us now analyze the process from a crystallization standpoint. The supersaturation ratio, that is the ratio between the solute concentration over its solubility in the local CO2 –ethanol composition, intervenes in both the nucleation and growth rates. The instantaneous supersaturation profile in the meridian plane of the domain is shown in Fig. 5c. In the conditions investigated here, the maximum attainable supersaturation ratio, calculated from the sprayed concentration of 10 mg/ml and of the solubility of minocycline in CO2 , is on the range of 17–18. In antisolvent process like the one described here, the mixing of the solution and of the antisolvent creates the conditions of supersaturation since the enrichment with the antisolvent provokes the decrease of the solute solubility. Therefore, one can expect areas of supersaturation to be localized at periphery of the jet where the solvent content is lower, which is indeed viewed in Fig. 5c. As indicated above, the solute is consumed by the precipitation events. The jet upstream thus appears to be depleted in high supersaturation ratio because no high concentration of minocycline is available as shown previously in Fig. 5b. The

Fig. 7. Influence of the minocycline concentration in the injected solution on the number of particles in the meridian plane of the reactor (a = 5 mg/ml; b = 10 mg/ml; c = 15 mg/ml) and on the volume distribution of the supersaturation levels in the whole domain (d). T = 313 K; P = 13 MPa; vinj = 4 m/s.

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9

7

consequence of supersaturation is the creation of particles. Fig. 5d and e shows the spatial distributions of the number of created particles and of the mass of solids they represent. Mass, number and diameter of particles are linked by Eq. (19). According to Fig. 5d, a high number of particles (0.5 × 1015 to 2 × 1015 ) are created at the periphery of the jet and within the 0.1–0.2 cm from the nozzle exit. They represent however a low mass (<0.7 × 10−18 , Fig. 5e) which indicates that particles are small. A high number of small particles is suggesting that this area is the center of nucleation events rather than growth ones. Above those 0.1 cm, the particles are less numerous (<0.46 × 1015 ) but their mass can reach 1.1 to 1.6 × 10−18 which indicates the presence of larger particles and localizes thus an area where growth is promoted over nucleation. 3.3. Influence of injection velocity and solute concentration

Fig. 8. Evolution of minocycline mean size as function of the injection velocity; C = 10 mg/ml, P = 13 MPa, T = 313 K.

One objective of the simulation is to assist the process optimization by estimating the evolution of particle sizes with operating parameters, and eventually, for conditions not investigated experimentally. However, this can be done only if the solubility data

Fig. 9. Mass fraction of ethanol in the meridian plane of the reactor (v = 4, 16, and 32 m/s) and volume distribution of supersaturation levels as function of the injection velocity. T = 313 K, P = 13 MPa, and C = 10 mg/ml.

8

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9

are available. In the selected case of minocycline, no data except at 313 K and 13 MPa were reported, so the effect of pressure or temperature could not be investigated. We therefore focused the investigations on two parameters susceptible to modify the spatial distribution of the solution and of the crystallization rate, namely the minocycline concentration and the injection velocity. 3.3.1. Influence of solute concentration The concentration of minocycline in the injected solution was varied from 5 to 15 mg/ml that according to Eq. (14) induces a rise in the overall attainable supersaturation ratio. As seen in Fig. 6, the particles become smaller as the concentration raises, a behavior that comes mostly from the stronger impact of S on the nucleation frequency than on the growth rate. Besides this overall behavior, it is interesting to see if the concentration impacts the spatial and instantaneous distribution of particles number in the mixing region (Fig. 7) and the supersaturation profiles as well. To quantify the supersaturation distribution, we calculated the volumes occupied by each supersaturation level, which sort out as the histogram displayed in Fig. 7d. According to the histogram, concentrations of 10 and 15 mg/ml are both characterized by the existence of supersaturation levels higher than 10. Their volume profiles exhibit a similar distribution although higher levels of saturation (above 20) are obtained for the 15 mg/ml data that can explain the overall smaller mean size observed in Fig. 6. At 5 mg/ml, the maximal supersaturation ratio is only of 10 and the levels occupy larger volumes than for previous cases. According to Eqs. (15) and (16), smaller supersaturations induce smaller nucleation frequencies, which in turn, induce the creation of less numerous particles as viewed in Fig. 7 (1011 versus 1014 and 1015 ). Less numerous particles in a larger volume seems then to favor growth other nucleation since particles of larger size were obtained in this condition (Fig. 6). 3.3.2. Influence of injection velocity The velocity was varied between 4 and 32 m/s. Since this parameter is susceptible to affect the mixing profile, we imaged the mass fraction of ethanol (Fig. 9a–c), the particle mean size (Fig. 8) and the supersaturation histogram (Fig. 9d). The solution concentration was kept at 10 mg/ml, so the overall attainable surpersaturation remains unchanged. Fig. 8 shows that the higher is the speed of injection, the smaller is the particles mean size, although the effect is much less significant than the effect of concentration since sizes evolved only from 205 to 140 ␮m. The histogram provided in Fig. 9d shows that the velocity does not affect strongly the volume distribution of the supersaturation levels nor the supersaturation values. Slight differences can be however seen at 4 m/s where the volume of fluid occupied by each level is smaller than at 16 or 32 m/s. This low impact can be understood by looking at the distribution of the solvent (Fig. 9a–c) that shows similar profiles as well. It has to be noticed that the increase of velocity tends to increase the length of the zone that remains rich in ethanol, so the area of precipitation is possibly displaced upstream but this does not impact notably the overall supersaturation behavior. 4. Conclusion The objective of the simulation developed in this work is to describe the physical phenomena involved in a spray crystallization induced by CO2 , i.e. hydrodynamics, mass transfer, phase equilibrium and crystallization kinetics, in order to rationalize the effect of a given parameter on the mean size of produced particles. A 3D model has been developed in order to get a close picture of the mixing between the injected solution and the supercritical antisolvent, which is the necessary but not sufficient condition for crystallization. To predict the particle size or at least their evolution

as function of the operating variables, the model requires thermodynamic data of the solute solubility in CO2 + solvent mixtures, that, for this work, were issued from literature. The simulation was validated by comparison with a literature result as well, and although the experimental and simulated size fit satisfactorily, the major role of the crystal–fluid interfacial tension on the accuracy was emphasized. The simulation also evidenced that the mixing area between the injected solution and the ambient CO2 is not homogeneous and distributes two different localizations of the precipitation events: nucleation dominated at the proximity of the injector exit, i.e. below 10 times the capillary diameter, whereas the growth process was significantly developed upstream this zone and within the full cone of the mixing area. The concentration of the solution was found to have a pronounced effect on the particle size, a behavior that fit the trends classically encountered in crystallization from solutions. By describing more finely the mixing area, the simulation showed that the diluted solution was characterized by lower supersaturation levels distributed over larger volumes of fluid than more concentrated solutions, a situation that, by producing larger particles, seemed to favor growth other nucleation. The injection velocity had the potentiality to modify the mixing between the solution and the ambient CO2 . However, within the investigated range of 4–32 m/s, only a decrease of size of 30% was obtained, due mostly to an only slight impact of velocity on the supersaturation volume distribution. In terms of process optimization, high velocity combined to high concentration would hence favor the production of smaller particle sizes. Future directions are to investigate the effect of pressure that is a specific parameter to a crystallization assisted by supercritical fluids (it requires the measure of solubility in CO2 + solvent) and the effect of the injector design, both in progress. Acknowledgements This work was granted an access to the HPC resources of CINES by GENCI (Grand Equipement National de Calcul Intensif) under reference number x2011026115. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compchemeng. 2012.12.002. References Badens, E., Boutin, O., & Charbit, G. (2005). Laminar jet dispersion and jet atomization in pressurized carbon dioxide. Journal of Supercritical Fluids, 36(1), 81–90. Baldyga, J. (1989). Turbulent mixer model with application to homogeneous, instantaneous chemical reactions. Chemical Engineering Science, 44(5), 1175–1182. Cardoso, M. A. T., Cabral, J. M. S., Palavra, A. M. F., & Geraldes, V. (2008). CFD analysis of supercritical antisolvent (SAS) micronization of minocycline hydrochloride. Journal of Supercritical Fluids, 47(2), 247–258. Fadli, T., Erriguible, A., Laugier, S., & Subra-Paternault, P. (2010). Simulation of heat and mass transfer of CO2 –solvent mixtures in miscible conditions: Isothermal and non-isothermal mixing. Journal of Supercritical Fluids, 52(2), 193–202. Falgout, R. D., & Yang, U. M. (2002). Hypre: A library of high performance preconditioners. Lecture Notes in Computer Science, 2331, 632–641; Sloot, P. M. A., Tan, C. J. K., Dongarra, J. J., & Hoekstra, A. G. (Eds.). (2002). Computational science – CARS Part III. Springer-Verlag. Goda, K. (1979). A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows. Journal of Computational Physics, 30(1), 76–95. Henczka, M., Bałdyga, J., & Shekunov, B. Y. (2005). Particle formation by turbulent mixing with supercritical antisolvent. Chemical Engineering Science, 60(8–9 special issues), 2193–2201. Martín, A., & Cocero, M. J. (2004). Numerical modeling of jet hydrodynamics, mass transfer, and crystallization kinetics in the supercritical antisolvent (SAS) process. Journal of Supercritical Fluids, 32(1–3), 203–219.

A. Erriguible et al. / Computers and Chemical Engineering 52 (2013) 1–9 Mersmann, A. (1990). Calculation of interfacial tensions. Journal of Crystal Growth, 102(4), 841–847. Mersmann, A. (2001). Crystallization technology handbook. New York: Marcel Dekker Inc. Patankar, S. V. (1990). Numerical heat transfer and fluid flow. New York: Hemisphere Publishing Corporation. Pratsinis, S. E. (1988). Simultaneous nucleation, condensation, and coagulation in aerosol reactors. Journal of Colloid and Interface Science, 124(2), 416–427. Reid, R., Prausnitz, J., & Poling, B. (1987). The properties of gases & liquids (5th ed.). McGraw-Hill. Sagaut, P. (1998). Large Eddy Simulation for incompressible flows – An introduction. Berlin: Springer-Verlag.

9

Schwarzer, H. C., Schwertfirm, F., Manhart, M., Schmid, H. J., & Peukert, W. (2006). Predictive simulation of nanoparticle precipitation based on the population balance equation. Chemical Engineering Science, 61(1), 167–181. Sierra-Pallares, J., Marchisio, D. L., Parra-Santos, M. T., García-Serna, J., Castro, F., & Cocero, M. J. (2012). A computational fluid dynamics study of supercritical antisolvent precipitation: Mixing effects on particle size. AIChE Journal, 58(2), 385–398. Vincent, S., & Caltagirone, J. (1999). Efficient solving method for unsteady incompressible interfacial flow problems. International Journal for Numerical Methods in Fluids, 30(6), 795–811. Wu, W., Ke, J., & Poliakoff, M. (2006). Phase boundaries of CO2 + toluene, CO2 + acetone, and CO2 + ethanol at high temperatures and high pressures. Journal of Chemical and Engineering Data, 51(4), 1398–1403.