Modeling of particle formation processes using gas saturated solution atomization

J. of Supercritical Fluids 41 (2007) 115–125

Modeling of particle formation processes using gas saturated solution atomization Matteo Strumendo ∗ , Alberto Bertucco, Nicola Elvassore DIPIC “I. Sorgato”, University of Padova, Padova, Italy Received 6 December 2005; received in revised form 30 August 2006; accepted 8 September 2006

Abstract In the PGSS process, a gas saturated solution is atomized through a nozzle and a spray is formed; afterwards, the gas (carbon dioxide) evaporates and droplets solidify. In this work, the behaviour of a carbon dioxide super-saturated solution drop in a low pressure environment is investigated. A mathematical model, based on the general multi-component equations of change and on the Stefan conditions to account for the moving boundaries, is solved numerically providing the spatial and temporal profiles of temperature, composition, mass flow inside and outside the droplet. The time required to attain the solid–liquid equilibrium condition at the droplet boundary decreases monotonically when the initial carbon dioxide content is increased, and increases together with the initial solution temperature. These trends are consistent with experimental data, showing that the average final particle size decreases if the mixing vessel pressure is raised and increases together with the spray temperature. © 2006 Elsevier B.V. All rights reserved. Keywords: Gas saturated solutions; Moving boundary problem; Particle formation; Particle size; Solid–liquid equilibrium

1. Introduction In recent years, a growing interest has been paid to supercritical fluids and their industrial applications. Taking advantage of supercritical fluids peculiar properties [1], new techniques have been developed and existing ones have been improved in many fields, i.e. separation processes, polymerization/depolymerization processes and particle production. Particularly, the use of supercritical fluids has been attractive for particle formation [2–4], above all for applications in the pharmaceutical, cosmetic, food, particle coating and fine chemistry industries. Some advantages of supercritical fluids-based processes are: the possibility of avoiding (or at least minimizing) the use of organic solvents; the flexibility to treat a wide range of materials, including thermolabile ones; the ability to produce particles with a wide range of sizes, shapes and morphologies; the possibility of loading the particles with an active ingredient [2], which can be dispersed in a matrix (microspheres) or surrounded by a shell (microcapsules). ∗ Corresponding author. Present address: Department of Chemical and Environmental Engineering, Illinois Institute of Technology, Chicago, IL, USA. Tel.: +1 312 567 3066; fax: +1 312 567 8874. E-mail address: [email protected] (M. Strumendo).

0896-8446/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.supflu.2006.09.003

The more important processes for particle production using supercritical fluids are RESS, SAS (and its implementations) and PGSS (Particles from Gas Saturated Solutions); comparisons among the three processes are presented in the review papers [2–4]. In PGSS [5], the material to be processed is melted and mixed at high pressure with the supercritical fluid (typically carbon dioxide) in a mixing unit. In this way, a gas saturated solution is obtained, which is then flowed through a nozzle to an expansion unit (at low pressure and temperature) so that a spray is formed. Carbon dioxide evaporates very quickly from the droplets, which eventually solidify due to the temperature decrease and to the carbon dioxide loss. Finally, solid particles are separated from the carbon dioxide by means of filters, cyclones and electrostatic separators. One of the key features of PGSS is the strong reduction of the melted substance viscosity once it is mixed with the supercritical fluid. A second effect of mixing the melted substance with the supercritical fluid is the melting point temperature depression; this phenomenon prevents the solidification inside the nozzle as well. Another interesting feature concerns the decrease of the temperature and the subsequent droplets solidification. This is usually ascribed [5] to a combination of the Joule–Thompson

116

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

Nomenclature cP f H j K 2 PM q r R Rgas RF t T u v∞ Y

specific heat at constant pressure fugacity Henry constant and enthalpy mass flux (radial direction) thermal conductivity upper limit of the gas phase domain molecular weight heat flux (radial direction) radial coordinate droplet radius gas constant reduction factor in the grid definition time temperature gas velocity (radial direction) droplet/gas approach velocity mass fraction

Greek letters ℘ mass diffusivity Π pressure λ enthalpy variation between the CO2 enthalpy in the gas phase and the CO2 partial enthalpy in the solution μ fluid viscosity ρ density τ stress tensor ξ dimensionless radial coordinate Subscripts/superscripts C carbon dioxide eq equilibrium value i cell index in initial value int intercept of a linear relationship mol molar N total number of cells sl slope of a linear relationship trist tristearin 1 liquid phase (solution) 2 gas phase (environment) ˆ value in the mixture – partial molar

effect and of the cooling due to the supercritical fluid evaporation. Further advantages of the PGSS process [5] are: solvent-free products; low consumption of carbon dioxide; the possibility of treating thermo-labile materials; a wide range of potential applications, since the solubility of compressed gases in liquids and solids is often high. So far, PGSS has been used to produce glycerides [5], pharmaceutical compounds like nifedipine, felodipine, dihydropyri-

dine, fenofibrate [5–7] and polymers [8]. Recently, patents have been proposed (a review is presented in [4]) to obtain microcapsules with polymeric shells by PGSS and/or its variations. A fairly general experimental knowledge of the effects of the operative parameters on particle properties (size, shape, porosity, morphology) is not available yet, but some correlations have been proposed [5–10]. It seems that an increase of the mixing unit pressure results in smaller particle size, but sometimes [8] particle size showed no dependency on this parameter. Particle size was found to increase monotonically with respect to the spray temperature [5,9], and also with respect to the nozzle diameter [8]. Different kinds of particle morphologies have been recognized, such as porous structures, spheres, sponges, fibres [6,8–10] but explanations on which formation mechanism should be related to each specific structure are still hypothetical. Empirical correlations to determine the particle morphology and bulk density as a function of the process conditions have been obtained for polyethyleneglycol-carbon dioxide systems [9,10]. Theoretical investigation of PGSS has just begun and few papers have been published to date in this area. Nevertheless, it is clear that a rational understanding can be useful in improving the PGSS process. Fundamental modeling of PGSS physics is particularly essential to the design of particles with specific characteristics. However, an exhaustive description of the PGSS process demands the investigation of many different fields and phenomena, including supercritical fluid thermodynamics, jet/spray hydrodynamics, droplet fluid dynamics, crystallization kinetics, bubble formation and droplet coalescence. In a schematic way, the overall investigation can be split into the following parts/steps: 1. thermodynamics; 2. nozzle hydrodynamics and atomization; 3. isolated droplet fluid dynamics, evaporation and solidification; 4. behaviour of the droplets ensemble in the spray, accounting for the droplet interactions and coalescence phenomena. Thermodynamics modeling is required to calculate the “initial condition”, i.e. the mixture composition for a given temperature and pressure in the mixing vessel; it is necessary also in the next steps, for example, to estimate the driving force for the energy and mass transfer. Such thermodynamics information is often provided in terms of melting point temperature curve as a function of the carbon dioxide pressure or in terms of carbon dioxide solubility isotherms. Investigations on the PGSS thermodynamics have been performed by Knez [8] and by Elvassore et al. [11]. Both these authors measured solubility curves and melting point temperature curves, obtaining similar results. Particularly, the melting point temperature was found to decrease when increasing the carbon dioxide pressure, until a minimum in the melting point temperature was reached; afterwards, the melting point temperature was increasing together with the carbon dioxide pressure. Besides, Elvassore et al. [11] reproduced experimental data imposing thermodynamical equilibrium conditions and using the Perturbed Hard-Sphere Chain

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

(PHSC) theory. They were able to develop operating charts to define the initial conditions (mixing vessel pressure and temperature) which would lead to solid or solid/liquid or liquid final products. The only work on the nozzle hydrodynamics and the PGSS process is contained in the papers of Li et al. [12,13]. One of the main goals of these investigations was to compute the droplet diameter produced by atomization, or the particle size formed inside the nozzle. In the first paper [12], Li et al. assumed that two phases are flowing through the nozzle, one being the carbon dioxide rich phase, the other being the gas saturated solution. They solved continuity, momentum and energy balances of the two phases and obtained the pressure, composition, temperature, density and velocity profiles in the nozzle. Atomization was supposed to be the mechanism for particle formation; the Jasuja relation was solved [14,15], coupled with the solution of the mentioned balance equations, in order to compute the average droplet diameter produced at the nozzle exit. In a second paper [13], Li et al. considered crystallization inside the nozzle as the mechanism for particle formation; assuming one pseudo-phase flow, they solved conservation equations for mass, momentum, energy and moments of the particle size distribution. In both the papers recalled above [12,13], authors claim that their models are able to capture qualitatively some experimental trends, however the predicted particle size are much smaller than the experimental values. It is clear though that further investigation on this field will be extremely valuable in order to clarify the fluid dynamics inside the nozzle and to predict correctly the droplet or particle size at the nozzle exit (a brief discussion on this subject is reported in Section 4 below). To the best of our knowledge, no papers have been published so far concerning the single droplet fluid dynamics in the framework of the PGSS process nor about the behaviour of the population of droplets and particles in the spray, accounting for their interactions. The present paper is an attempt to investigate the fluid dynamics of an isolated droplet, considering mass and heat transfer phenomena and the carbon dioxide evaporation until the droplet solidification process begins. In a first approximation, the behaviour of the gas saturated solution droplet will be modelled from the time at which it is created in the spray (at the exit of the nozzle) until the solid–liquid equilibrium condition (SLEC) is achieved; interactions of any kind with other droplets will not be addressed. As an example, tristearin and carbon dioxide will be considered as the components of the gas saturated solution; the choice of the tristearin is determined by its relevant applications in the pharmaceutical industry [11]. The main physical phenomena can be summarized as follows. At the initial time, the gas saturated solution droplet is placed in a stagnant carbon dioxide atmosphere at low pressure (corresponding to the expansion unit conditions of the PGSS process). Temperature and composition of the droplet are actually those of a super-saturated solution since they are close to the conditions in the high pressure mixing unit, whereas the stagnant carbon dioxide environment is at low pressure. Therefore, carbon dioxide tends to evaporate from the saturated solution. The carbon dioxide change of phase can happen by surface evaporation or by bubble formation inside

117

the droplet itself. If the first mechanism is prevailing, the droplet will shrink because of the carbon dioxide mass loss, although the droplet radius reduction will be generally small. Carbon dioxide evaporation remarkably cools down the droplet surface. This, together with the droplet temperature decrease due to the heat conduction between the droplet and the environment, causes the droplet to reach the SLEC, which is given in terms of melting point temperature curve [11]. The description becomes more complex if also bubble formation inside the droplet is considered, a mechanism which is not taken into account in this paper. The modeling approach in the present paper follows the concepts of the Stefan problem formulation for the solution of moving boundary problems [16–18]. Typically, different phases coexist and are separated by a moving boundary; governing equations are solved in each phase while the behaviour at the moving boundary is described by the Stefan conditions, which represent surface balances at the moving boundary. The solution of the problem is given in terms of spatial-temporal profiles of density, velocity, temperature, compositions and in terms of the temporal evolution of the boundary displacement. Applications of these techniques in the area of supercritical fluids have been already presented by Werling and Debenedetti [19,20], as far as the SAS process and by Henczka et al. [21], as far as the spray-freeze drying process. In order to carry out a first, basic modeling of the droplet behaviour in the PGSS framework, some assumptions have been introduced. On the one hand, it is essential to take into account both the mass and heat transfer, because the SLEC is attained when the melting point temperature is reached and the melting point temperature is a function of the carbon dioxide content (see [8,11] and the discussion in Section 2). On the other hand, the solution density dependence with respect to the carbon dioxide content is neglected; simulations were stopped when the SLEC was reached on the droplet boundary (as it will be discussed in Section 2); both the evolution of the solid front inwards in the droplet and the crystallization kinetics are not simulated. Additionally, only the conductive contribution has been used in the energy balance within the droplet, neglecting the interdiffusion contribution. Besides, no relative motion is assumed between the droplet and the environment (the droplet is still in an initially stagnant atmosphere), so convective effects related to that have been neglected; this assumption will be discussed with more details in Section 4. Finally, gravity has not been considered. 2. Governing equations A spherical droplet of a solution of tristearin and carbon dioxide is placed in a stagnant carbon dioxide atmosphere; the initial droplet radius is R(0). Within the droplet the initial density, temperature and composition are uniform and the velocity is zero. Outside of the droplet, the temperature and pressure of the stagnant medium are also uniform at time zero; particularly the pressure is equal to 1 atmosphere. The initial gas velocity is also zero. Governing equations are derived from the general nonisothermal multicomponent form of the equations of change

118

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

[22]. Spherical symmetry is assumed. The pressure field around and within the droplet can be considered constant. Subscripts 1 and 2 denote liquid and gas phases, respectively. Since the droplet radius is a function of time R(t), the governing equations for the phase 1 will be solved in the domain 0 ≺ r ≺ R(t), while the governing equations for the phase 2 will be solved in the domain R(t) ≺ r ≺ 2 . Here 2 denotes the point in which boundary conditions for the phase 2 are set; this quantity would be equal to infinity if an analytical solution were sought for, however, it will be set (Section 3) equal to a finite value in order to solve the problem numerically. As discussed in Section 1, in phase 1 (solution) the density has assumed to be constant and the solution hydrodynamic velocity is zero (no convective motion inside the droplet). Accordingly, the only equations to be considered are the transport equations for the energy and for one of the species (carbon dioxide). Neglecting interdiffusion heat fluxes and gravity effects, the following balance equations result, which are valid inside the droplet (0 ≺ r ≺ R(t)): ρ1 cP,1 ρ1

∂T1 1 ∂r 2 q1 =− 2 ∂t r ∂r

(1)

∂Y1,C 1 ∂r 2 j1 =− 2 . ∂t r ∂r

(2)

In Eqs. (1) and (2), ρ1 is the solution density, cP,1 the solution specific heat (constant pressure), T1 the temperature, q1 the heat flux (conductive contribution) in the radial direction, Y1,C the carbon dioxide mass fraction and j1 is the mass flux (diffusive contribution) in the radial direction. In phase 2 (carbon dioxide gas), variations in the pressure field can be neglected and the carbon dioxide can be assumed incompressible [23]. Therefore, we will consider momentum and energy balance equations, which have to be solved in the domain R(t) ≺ r ≺ 2 : ρ2

1 ∂r 2 τrr,2 τθθ,2 + τφφ,2 ∂u2 ∂u2 + ρ 2 u2 =− 2 + ∂t ∂r r ∂r r

ρ2 cP,2

∂T2 ∂T2 1 ∂r 2 q2 + ρ2 cP,2 u2 =− 2 . ∂t ∂r r ∂r

(3) (4)

In Eqs. (3) and (4), ρ2 is the carbon dioxide density, u2 the carbon dioxide radial velocity, τ rr,2 , τ φφ,2 , τ θθ,2 the three normal components of the momentum flux viscous part, cP,2 the carbon dioxide specific heat (constant pressure), T2 the temperature and q2 is the heat flux (conductive contribution) in the radial direction. In Eq. (4), viscous dissipation was neglected. The equations of change (1)–(4) were solved with respect to the variables T1 , Y1,C , u2 , T2 . Appropriate boundary conditions are required for these variables. Symmetry boundary conditions at the droplet centre are set for T1 and Y1,C . At 2 , the following boundary conditions are prescribed: u2 (2 , t) = u2 (r, 0) = 0 T2 (2 , t) = T2 (r, 0) = T2,in

.

(5)

In Eq. (5), T2,in denotes the initial, uniform value of the environment temperature.

The remaining boundary conditions for T1 , Y1,C , u2 , T2 are prescribed at the droplet radius. These, together with the rate of displacement of the droplet radius dR(t)/dt, are obtained by the solution of the following equations, valid at the droplet radius r = R(t): q2 (R) − q1 (R) = λρ1 −ρ2 u2 = (ρ1 − ρ2 )

dR dt

(6)

dR dt

j1,C (R) = −ρ1 (1 − Y1,C )

(7) dR dt

T1 (R) = T2 (R) = Teq

(8) (9)

 f 1,C = f2,C .

(10)

The first three Eqs. (6)–(8) are the generalized Stefan conditions, representing surface balances of energy, mass and carbon dioxide on the droplet moving boundary. Eq. (6) essentially states that the heat flux coming from the interior of the droplet partly is used to vaporize the carbon dioxide and partly heats the surrounding phase 2 (the CO2 atmosphere). The coefficient λ is the enthalpy variation between the CO2 enthalpy in the gas phase and the CO2 partial enthalpy in the solution at the boundary conditions. Eq. (6) can be obtained from the general form of the energy balance at the surface by assuming that: mechanical energy contributions are negligible; the CO2 mass fraction is negligible at the boundary, because thermodynamical equilibrium is assumed at the boundary (according to Eqs. (9) and (10)) and the pressure is low (1 atmosphere); the hydrodynamic velocity in the phase 1 is zero. Eq. (7) shows that a convective flux (vaporization) arises in the gas phase because of the density difference between the two phases. Eq. (8) explains that a carbon dioxide diffusive flux in the liquid phase is necessary to sustain the carbon dioxide vaporization. Eqs. (9) and (10) represent the thermodynamical equilibrium conditions at the droplet surface. Since the pressure effects are not relevant and since there is no tristearin in the gas phase, only the thermal equilibrium condition (Eq. (9)) and the iso-fugacity equilibrium condition for the carbon dioxide (Eq. (10)) are considered. The previous Eqs. (1)–(4), (6) and (8) must be coupled with the constitutive equations for the heat, carbon dioxide and momentum fluxes. These are given by [22]: q1 = −K1

∂T1 ∂r

(11)

q2 = −K2

∂T2 ∂r

(12)

j1 = −ρ1 ℘1

∂Y1,C ∂r

  ∂u2 τrr,2 = −μ2 2 , ∂r  u  2 τϕϕ,2 = −μ2 2 . r

(13)  u  2 τθθ,2 = −μ2 2 , r (14)

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

In the above equations, K1 and K2 are the thermal conductivities of the phases 1 and 2, respectively, ℘1 is the binary (tristearin-carbon dioxide) mass diffusivity in the phase 1 and μ2 is the carbon dioxide viscosity in the gas phase. The heat flux is given by the conductive contribution, the carbon dioxide mass diffusional flux by the ordinary diffusion term and the momentum flux by the Newton contribution. The iso-fugacity relation (10) can be rewritten in terms of solubility of carbon dioxide in the tristearin: mol H(Teq )Y1,C = Π.

(15)

In the previous equation, Π is the total pressure, assumed mol is the carbon equal to the initial pressure of 1.01325 bar, Y1,C dioxide molar fraction in the solution and H(Teq ) is the Henry constant. The dependence of the Henry constant with respect to the temperature for this system has been evaluated from the data of Elvassore et al. [11]. A linear relationship holds between H and Teq : H(Teq ) = Hsl Teq + Hint . Hsl = 0.2723 bar/K Hint = −63.5797 bar

(16)

The carbon dioxide density in the phase 2 has been computed by the ideal gas law: ρ2 =

ΠPMC Rgas T2,in

(17)

where the initial uniform value of the temperature has been used; PMC is the carbon dioxide molecular weight and Rgas is the gas constant. The density of the saturated solution in the phase 1 has been estimated as 898.2 kg/m3 (value at 75 ◦ C from [24] for the tristearin; 75 ◦ C is a typical value used in simulations for phase 1). Constant values are set for the coefficients of mass diffusivity, for the thermal conductivities in both the phases, for the gas viscosity, for the specific heats (at constant pressure) in both the phases and for the enthalpy variation λ; in fact, their variations are not expected to affect the main trends of the process. The carbon dioxide specific heat (phase 2) has been estimated as 870 J/(kg K) (value at 40 ◦ C, typical value of the phase 2 temperature); the carbon dioxide viscosity is equal to 1.5526 × 10−5 kg/(m s) (at 40 ◦ C); the carbon dioxide thermal conductivity is equal to 1.71 × 10−2 J/(m s K) (at 40 ◦ C). The solution specific heat has been taken equal to the tristearin specific heat, which can be evaluated as 2208 J/(kg K) (at 75 ◦ C from [24]). The value used for the solution thermal conductivity is the stearic acid thermal conductivity (supposedly close to the tristearin one) and is equal to 0.1584 J/(m s K) (at 75 ◦ C from [24]). In the model presented, mass diffusivity values are critical above all at the boundary, because that affects the boundary conditions (Eqs. (6)–(10)); therefore, the infinite dilution value is used, equal to 0.92 × 10−9 m2 /s. According to the suggestion of Lockermann and Schlunder [25], mass diffusivity was computed using the formula proposed by McManamey and Woollen ([26]; tristearin viscosity at 75 ◦ C is evaluated from [24]). The

119

enthalpy variation λ has been estimated by: mol ∂H1 ¯ C = H1 + Y1,trist . H mol ∂Y1,C

(18)

¯ C is the partial molar enthalpy of the carbon dioxide In Eq. (18), H mol and Y mol in the solution, H1 is the solution molar enthalpy, Y1,trist 1,C are the liquid molar fractions of the tristearin and the carbon dioxide, respectively. Using the liquid enthalpy data (at 75 ◦ C and at infinite dilution for the carbon dioxide) by Elvassore et al. [11], the enthalpy variation λ results equal to 6.249 J/mol (the influence of pressure on the enthalpy data used is negligible). One of the goals of the present work is to compute the time required for the SLEC to occur in the saturated solution droplet. This can be considered as the beginning of the solidification process. The SLEC is given by the melting point temperature curve (Fig. 1), which has been obtained from the data of Elvassore et al. [11]. These data have been replotted in terms of melting point temperature versus carbon dioxide content and the results of the PHSC correlation, developed by Elvassore et al. [11], have been fitted using a second order polynomial. Under the hypotheses stated in Section 1, it is likely that the SLEC is reached first at the droplet boundary, because that is the point inside the droplet where both the temperature is lower and the carbon dioxide content is lower (i.e. the melting temperature is higher). Anyhow, assuming that the SLEC happens first at the droplet boundary is not really restrictive, because in the numerical solution (Section 3) it is possible to check at any point within the droplet whether the SLEC occurs or not. Finally, a brief discussion is necessary regarding the data used to draw the melting point temperature curve (Fig. 1). In fact, these data have been measured by Elvassore et al. [11] in solid–liquid–gas equilibrium conditions and, because of that, pressure variations and carbon dioxide content variations are linked (one degree of freedom). Instead, an ideal procedure

Fig. 1. Melting point temperature profile: carbon dioxide mass fraction vs. melting temperature. Results of the PHSC correlation (circles) developed by Elvassore et al. [11] are fitted by a second order polynomial (solid line). The PHSC model represents fairly well the experimental data (squares) of Elvassore et al. [11].

120

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

would be to consider a melting point temperature as function of the pressure and of the carbon dioxide content, where the pressure and the carbon dioxide content are independent variables; then, to set the pressure equal to the operating conditions (1 atmosphere in the simulations performed in this work), in order to obtain the dependence of the melting temperature with respect to the composition only. This is not possible with the available data because the pressure and composition variations are linked. However, data of Elvassore et al. [11] can be used because of the following reasons: (1) pressure has a weak effect on the melting point curve of the pure tristearin [27] and, anyway, the melting temperature increases together with the pressure; (2) therefore, the melting point temperature variation with respect to the joint variation of pressure and carbon dioxide can be essentially ascribed to the carbon dioxide influence; (3) this is especially true at low-medium carbon dioxide pressures, where the melting point temperature decreases (Fig. 1) while the pure tristearin melting temperature increases together with the pressure, and where the melting temperature decrease (Fig. 1) is much bigger than the variations induced by the pressure in the pure tristearin [27]; (4) carbon dioxide mass fraction values at the droplet boundary are always low (this will be clear in the Section 4). Therefore, a melting point temperature curve derived from the three-phase equilibrium data of Elvassore et al. [11] can be legitimately used. 3. Numerical technique and validation A numerical solution of the mathematical model outlined in the previous Section 2 was developed. Particularly, Eqs. (1)–(4), (6)–(9) and (15) were solved numerically. Moving boundary conditions (Eqs. (6)–(9) and (15)) have been dealt by a frontfixing method [16]. Specifically, the following coordinate transformations have been applied in the phase 1 domain: ξ1 =

r R(t)

(19)

and in the phase 2 domain: ξ2 =

r − 2 . R(t) − 2

(20)

In this way, both the coordinates ξ 1 and ξ 2 range from 0 to 1. Model equations were rewritten in dimensionless form using the following scale factors: 2 for r, R(t) and 2 itself; 22 (K1 /ρ1 cP,1 )−1 for the time; the initial, uniform value of  the environment temperature T2,in for the temperatures; Rgas T2,in /PMC for the gas velocity. Eqs. (1)–(4), after transformation in the new reference frame and in dimensionless form, were discretized using a finite volumes technique [28,29]. Note that a finite volumes discretization cannot be fully applied for the following reasons: 1. Eqs. (3) and (4) are formally expressed not in the conservation (or divergence) form, as far as concerns the convective terms. The reason is that the incompressible flow form is used for

Eqs. (3) and (4); eventually, a full conservative form for these terms could be implemented. 2. In Eq. (3) the terms of the normal not radial components of the stress tensor are not formally expressed in the conservation (or divergence) form. This is because a spherical reference frame is used. 3. The coordinate transformation (19) and (20) introduces additional terms, not expressed in the conservative (or divergence) form. Therefore, these terms, not expressed in conservative form, were discretized as follows: 1. In the cell integration of the convective terms in Eqs. (3) and (4), integrands were assumed constant, using the central value for the gas velocity and central differences for the derivative terms. 2. Central values were used in the cell integration of the normal not radial components of the stress tensor terms. 3. In the cell integration of the additional terms introduced by the coordinate transformation (19) and (20), integrands were assumed constant, using the central value for the ξ 1 or ξ 2 coordinates and central differences for the derivative terms. Derivatives in the expression of the fluxes in Eqs. (11)–(14) were evaluated at the cells boundaries; central differences with respect to the cell boundary point were used. The discretization of Eqs. (1)–(4) and (11)–(14) leads to an ordinary differential equations (ODE) system, which was solved numerically using a backward differentiation formulas (BDF) technique [30]. At the moving boundary (droplet radius), Eqs. (6)–(9) and (15) must be solved at each time step and the first derivatives were discretized by backward differences. After discretization, one of these equations (say Eq. (6)) provides the temporal evolution of the moving boundary position R(t), the remaining equations form an algebraic equation system whose solution provides the boundary values of the temperatures, composition and gas velocity. Centers of the cells were chosen according to the criterion of having more points close to the droplet radius where gradients are steeper (cell faces were located midway between the centers of the cells). Center points ξ i in each cell i were placed as follows: ξi =

1 + RF + · · · + RFi−2 1 + RF + · · · + RFN−2

i = 2, · · ·, N.

(21)

In Eq. (21), N is the total number of cells and RF is a reduction factor, equal to the constant ratio: ξi+1 − ξi . ξi − ξi−1 The center of the first cell is at ξ i = 0 and the center of the last one is at ξ i = 1. When RF is equal to 1, the grid is uniform. Lower values of RF can be chosen when gradients are steeper, particularly when the initial carbon dioxide mass fraction is higher. In the simulations performed, RF was equal to 0.97 (in both phases). Tests were performed to check the grid independence

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

Fig. 2. Comparison between the analytical solution [23] and the numerical results: moving boundary displacement vs. time profile. Both the quantities are dimensionless.

of the results. The value of 2 was set so that the ratio R(0)/2 was equal to 1/25. The numerical technique was validated by comparison with a number of analytical solutions [17,31,23], obtaining always results which follow closely the exact solution. As an example, the numerical solutions computed in the simulations of the Stefan problem with convection in spherical coordinates [23] are plotted in Figs. 2–4; the comparison with the respective analytical solutions is also shown. The Stefan problem with convection describes the behaviour of a solid phase surrounded by a fluid phase; the solid phase initial size is zero and its temperature is constant and equal to the solidification temperature. For this sample case the solution of the governing equations [23] describes the growth of the solid phase in terms of solidification front displacement (Fig. 2), temperature radial distribution in the fluid phase (Fig. 3) and fluid velocity radial distribution (Fig. 4).

Fig. 3. Comparison between the analytical solution [23] and the numerical results: temperature spatial distribution in the fluid phase at the initial and final times; the initial profile is the flat one. Temperature is expressed in dimensionless form, while the radial coordinate is in meters.

121

Fig. 4. Comparison between the analytical solution [23] and the numerical results: radial velocity spatial distribution in the fluid phase at the initial and final times; the initial profile is the flat one. Velocity is expressed in dimensionless form, while the radial coordinate is in meters.

4. Results and discussion The equations of change (1)–(4) were solved following the numerical technique in Section 3 obtaining, as a result, the spatial (radial) and temporal profiles of the variables T1 (solution temperature), Y1,C (CO2 mass fraction in the solution), u2 (gas radial velocity), T2 (gas temperature). The equations of change (1)–(4) were solved subject to the initial conditions (uniform values for all the variables, zero for u2 ) and to the boundary conditions. These are given by the symmetry conditions at the droplet center, by Eq. (5) at 2 and by Eqs. (6)–(9) and (15) at the droplet interface. One of the variables in Eqs. (6)–(9) and (15) is the interface position R(t), which is a function of time; therefore, the interface displacement versus time profile was also determined, given the initial value R(0). Particularly, it is interesting to study how the profiles of T1 (r,t), Y1,C (r,t), u2 (r,t), T2 (r,t), R(t) vary with respect to the initial values of the temperature and CO2 mass fraction within the droplet T1,in , Y1,C,in and with respect to the initial interface position R(0). Simulations were stopped when the SLEC was attained. As discussed in Section 2, the SLEC is assumed to be attained always at the droplet interface when the melting point temperature is reached (Fig. 1). The dependence of the SLEC time (the time required to reach the SLEC at the droplet interface) versus T1,in , Y1,C,in and R(0) was investigated. A typical example of the results obtained is shown in Figs. 5–7, where the radial profiles of the CO2 mass fraction (inside the droplet), of the temperature and of the convective radial gas velocity (outside the droplet) are shown, respectively. These results refer to Y1,C,in equal to 0.1, T1,in equal to 348.15 K, T2,in (initial environment temperature) equal to 313.15 K. The corresponding trend of the interface position versus time is given in Fig. 8. The droplet radius is 20 ␮m.

122

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

Fig. 5. Radial carbon dioxide mass fraction profile: dotted line represents the initial condition, solid line represents the final solution (at the SLEC time). Initial condition uniform values: the CO2 mass fraction is 0.1, the droplet temperature is 348.15 K, the environment temperature is 313.15 K.

The values of Y1,C,in and T1,in , in the PGSS process, can be somehow related to the values of the pressure in the mixing vessel before the atomization. This demands the solution of the fluid dynamics from the mixing vessel to the nozzle. If this, in a first approximation, is neglected, the mixing vessel pressure value corresponding to Y1,C,in = 0.1 and T1,in = 348.15 K can be obtained from the data of Elvassore et al. [11] and it is approximately equal to 6 MPa. As expected [18,32]), mass diffusion in the gas saturated solution is much slower than heat conduction. Comparing Fig. 5 with Fig. 6, at the final time about 1/3–1/2 of the droplet (the more internal part) has still the initial value of the CO2 , while in the remaining part of the droplet the composition has changed and

Fig. 6. Radial temperature profile: dotted line represents the initial condition, solid line represents the final solution (at the SLEC time). In the left side, there is the droplet, in the right side there is the environment. Initial condition uniform values: the CO2 mass fraction is 0.1, the droplet temperature is 348.15 K, the environment temperature is 313.15 K. Scale factor for the temperature is the environment temperature.

Fig. 7. Radial carbon dioxide radial velocity profile in the gas phase (environment): dotted line represents the initial condition, solid line represents the final solution (at the SLEC time). Initial condition uniform values: the CO2 mass fraction is 0.1, the droplet temperature is 348.15  K, the environment temperature is 313.15 K. Scale factor for the velocity is Rgas T2,in /PMC .

diffusion has occurred. Instead, temperature at the final time has changed in the entire droplet and the temperature profile within the droplet is relatively constant (compared with the CO2 profile). This is a general remark, that was observed in all the simulation results, at different values of T1,in , Y1,C,in . This behaviour is enhanced when the SLEC times are shorter (as it will be discussed in the following, SLEC times are shorter when Y1,C,in is higher and T1,in is lower). A simplified picture of the process is as follows. The low pressure (typically 1 atmosphere) of the gas phase surrounding the droplet determines, through the chemical equilibrium at the interface (Eqs. (10) or (15)), a low value of the CO2 content at the boundary. This sets the melting point temperature (Fig. 1), which

Fig. 8. Interface position vs. time profile: the numerical solution is the solid line, while the circles represent a sqrt (time)-law. The scale factor for the displacement is 2 , the scale factor for the time is 22 (K1 /ρ1 cP,1 )−1 .

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

is usually close to the pure tristearin melting point temperature at atmospheric pressure. Therefore, since the SLEC time can be seen as the time required for the droplet temperature (particularly at the boundary) to reach the melting temperature, the process basically depends on the cooling rate within the droplet and/or on the initial value of the droplet temperature (the initial distance of the temperature with respect to the melting temperature). For a given initial droplet temperature, the cooling rate is essentially given by the CO2 content within the droplet, because an increase of the CO2 content determines an increase of the evaporation rate (Eq. (8)), and this enhances the droplet cooling rate. The radial convective velocity of the evaporating gas (Fig. 7) is determined at the droplet boundary (through Eq. (7)) by the interface rate of displacement (Fig. 8). The latter has its maximum value at the beginning of the evaporation process, thereafter it decreases monotonically; the former follows the same trend. In all the simulations, the droplet shrinks during the evaporation, but the SLEC at the boundary is reached so quickly that the degree of shrinking is very small; besides, the trend of the interface position versus time seems to approximate closely a relationship of proportionality between the interface position (with respect to the initial interface position) and the square root of the time. This observation might be eventually useful in a more complex model, accounting also for the droplet interactions. In fact, in that framework, it is essential to reduce the number of parameters which characterize the specific state of each droplet and to have approximate but simple expressions describing the droplet state. In this respect, also the findings about the different relaxation times between the mass and heat diffusion within the liquid phase or between the diffusional processes in the liquid phase compared to the gas phase ones might be valuable pieces of information. Simulations were performed for different values of T1,in , Y1,C,in , particularly for T1,in ranging from 338.15 to 358.15 K and for Y1,C,in ranging from 0.02 to 0.3. As already mentioned for the case of Y1,C,in = 0.1 and T1,in = 348.15 K, the corresponding range of values in terms of mixing vessel pressure is between 0.1 and 20 MPa (data of Elvassore et al. [11]). In Figs. 9 and 10, simulation results representing the SLEC time as a function of the initial CO2 content (i.e. pressure) and of the droplet temperature are plotted, respectively. The SLEC time decreases monotonically when the initial CO2 content (mixing vessel pressure) is increased (Fig. 9). For example, considering the curve at 348.15 K, the SLEC time decreases of about one half, from 0.011 s (value at the lower pressure and CO2 content) to about 0.006 s (value at the higher pressure and CO2 content). This trend would be even more remarkable if the mass diffusivity dependence on the CO2 content were considered: in fact, the higher is the CO2 content, the bigger is the mass diffusivity [25] and the shorter is the SLEC time. It is relevant to a fundamental understanding of the PGSS physics to compare the decreasing trend of the SLEC time versus the CO2 content with the experimental measurements regarding the mean particle diameter versus the mixing vessel pressure. In fact, it has been generally found [5–8] that the average particle diameter decreases when the pressure is increased (an exception

123

Fig. 9. Droplet SLEC time vs. initial CO2 mass fraction profile at different initial droplet temperatures (338.15, 348.15, 358.15 K).

is in [8], where the pressure appears not to have any significant effect on the particle diameter). This is a fundamental feature of PGSS, where particle micronization is obtained essentially by means of the CO2 pressure applied to the liquid solution. In the simulations performed, results were obtained in terms of SLEC time and not in terms of mean particle diameter, but these two quantities can be related, at least qualitatively. In fact, during the droplet lifetime after it is formed close to the nozzle exit and before solidification, a droplet can move in the expansion unit entrained by the gas or can collide with other droplets. In this discussion, we do not consider drop break-up, because drop break-up effects are expected to be relevant close to the nozzle exit but not afterwards: in fact, droplet–gas relative velocities, which determine the maximum stable drop size, can be high close to the nozzle exit but afterwards droplets decelerate rapidly to attain the terminal velocity. Therefore, droplet–droplet colli-

Fig. 10. Droplet SLEC time vs. initial droplet temperature profile at different initial CO2 mass fraction (0.02, 0.18, 0.3).

124

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

sions must be primarily considered to understand how SLEC times affect the final particle size. Specifically, the SLEC time should be compared with the droplet–droplet collision time: clearly, the smaller is the SLEC time, the smaller will be the number of droplet collisions (at constant droplet–droplet collision time), which could result in coalescence [18] and in a final particle diameter increase. In this sense, results shown in Fig. 9 provide a fundamental understanding of PGSS: an increase in the mixing vessel pressure means an increase of the initial CO2 content; this determines shorter SLEC times, which reduce droplet coalescence and, therefore, result in smaller particle diameter. The variation of the SLEC time versus the initial droplet temperature is shown in Fig. 10. The SLEC time goes to zero as the temperature approaches the melting point temperature, and increases monotonically with respect to the temperature. Again, this is in agreement with the experimental behaviour [5], that the particle size increases monotonically with respect to the spray temperature: following the same arguments expressed regarding the dependence of the particle size on the CO2 content, a temperature decrease determines shorter SLEC times, which reduce droplet coalescence and result in smaller particle size. The knowledge of the SLEC time can be helpful also in the equipment design, for example, in the expansion unit size computation. In fact, the unit should be large enough to allow the droplet to have the time to solidify. On the one side, given the solution pressure (in the mixing vessel) and temperature, the SLEC time can be computed and the characteristic length of the expansion unit should be equal or bigger than the droplet velocity times the SLEC time. On the other side, if the expansion unit size and the droplet velocity are known, it is possible to compute the minimum pressure to be applied in the mixing vessel and/or the maximum solution temperature such that droplets can solidify before hitting the expansion unit walls. The dependence of the SLEC time on the initial droplet radius can be easily deduced. In fact, governing equations are written in dimensionless form (Section 3), the scale factor for the time being 22 (K1 /ρ1 cP,1 )−1 . Keeping constant in the simulations the ratio R(0)/2 , this results in a quadratic dependence of the SLEC time on the initial droplet radius R(0). In the model developed to describe the single droplet fluid dynamics, no relative motion between the droplet and the surrounding atmosphere has been assumed (Section 1). Such assumption is accepted in the literature [18] for small droplets (droplets whose diameter is below 30 ␮m). For larger droplets, evaluation of the Reynolds number (ρ2 (2R)v∞ )/μ2 (v∞ is the droplet/gas approach velocity) is required to assess the validity of the assumption, but detailed experimental values of v∞ in the PGSS process have not been published yet. However, Reynolds numbers are expected to remain relatively small in most of the PGSS applications, considering that: (1) ρ2 has always a low value in the expansion unit conditions; (2) v∞ is expected to decrease when larger particles are produced and lower CO2 pressures are applied in the mixing vessel. Therefore, convective effects related to the relative motion between droplets and the surrounding atmosphere should not influence significantly the simulation results.

According to the results discussed above, an explanation of the experimental relationship between the particle size and the operative variables has been provided; further, some fundamental features of the PGSS physics have been clarified. However, nozzle hydrodynamics and atomization could play an important role too, particularly to explain the dependence of the final particle size on the process conditions. In spite of the work of Li et al. [12,13], it is still not clear to which extent the formation of particles by crystallization within the nozzle can be considered competitive with respect to the particles formation by atomization and subsequent evaporation/solidification. Further, it would be extremely valuable to clarify the fluid dynamics regime inside the nozzle, particularly which kind of gas–liquid flow pattern should be expected to occur inside the nozzle. In this respect, it could be useful to reconsider the concepts of gas dissolved atomization [33,34]. In fact, although the pressures in PGSS are high compared to those used in some applications of the gas-dissolved atomization [33] and although primary expansion chambers are not used in PGSS, the atomization mechanisms share many analogies. Particularly, in both cases a solution of a low-volatile component is saturated by a high vapour pressure gas and then it is expanded at low pressure. It would be interesting to ascertain whether the mechanisms of bubble formation and bursting/flashing developed in the gas dissolved atomization occur also in the PGSS process. 5. Conclusions The behaviour of an isolated gas saturated solution droplet in a gaseous environment was mathematically modeled. The model describes the transport phenomena within the droplet and in the surrounding atmosphere, considered as incompressible. Inside the droplet the model accounts for species and energy diffusion; in the gas phase it includes convection (due to the evaporation), momentum and energy transfer. At the droplet–gas interface, the evaporation of the dissolved gas causes a remarkable cooling of the droplet and a dissolved gas diffusive flux from the droplet interior. Stefan conditions were developed to determine the boundary conditions at the phase discontinuity and the position of the interface itself. The model was numerically solved until the saturated solution reaches the SLEC at the droplet boundary; the SLEC is given by the melting point temperature curve as a function of the dissolved gas content. It was found that the computed SLEC times decrease monotonically when the initial CO2 content increases and when the initial droplet temperature decreases. This is in fair agreement with experimental data [5–8] showing a particle size decrease at increasing mixing vessel pressures and at decreasing spray temperatures. In the development of the model, some assumptions were made. First, the convective mass flux in the liquid solution and the liquid density variation with respect to the CO2 content were neglected; further, simulations were stopped when the SLEC was achieved, so that the solid phase growth and the crystallization kinetics were not accounted for. In this regard, work is in progress to complete the model, by also including the CO2 bubble nucleation inside the droplet, in order to provide valuable

M. Strumendo et al. / J. of Supercritical Fluids 41 (2007) 115–125

information for the investigation of the final particle morphology and for the particle design. Acknowledgements The authors are thankful to Prof. M. Putti and Prof. L. Bergamaschi of University of Padua (Italy) for the useful discussions regarding the numerical implementation of the model presented. References [1] B. Bungert, G. Sadowski, W. Arlt, Separations and material processing in solutions with dense gases, Ind. Eng. Chem. Res. 37 (1998) 3208–3220. [2] J. Jung, M. Perrut, Particle design using supercritical fluids: literature and patent survey, J. Supercrit. Fluids 20 (2001) 179–219. [3] J. Fages, H. Lochard, J.-J. Letorneau, M. Sauceau, E. Rodier, Particle generation for pharmaceutical applications using supercritical fluid technology, Powder Technol. 141 (2004) 219–226. [4] S.-D. Yeo, E. Kiran, Formation of polymer particles with supercritical fluids: a review, J. Supercrit. Fluids 34 (2005) 287–308. [5] E. Weidner, Z. Knez, Z. Novak, PGSS (particles from gas saturated solutions)—a new process for powder generation, in: G. Brunner, M. Perrut (Eds.), Proceedings for the Third International Symposium for Supercritical Fluids, vol. 3, 1994, p. 229. [6] P. Sencar-Bozic, S. Srcic, Z. Knez, J. Kerc, Improvement of nifedipine dissolution characteristics using supercritical CO2 , Int. J. Pharm. 148 (1997) 123–130. [7] J. Kerc, S. Srcic, Z. Knez, P. Sencar-Bozic, Micronization of drugs using supercritical carbon dioxide, Int. J. Pharm. 182 (1999) 33–39. [8] Z. Knez, High pressure micronisation of polymers, in: G. Brunner, I. Kikic, M. Perrut (Eds.), Proceedings for the Sixth International Symposium for Supercritical Fluids, 2003, p. 1865. [9] P. Kappler, M. Petermann, E. Weidner, Size and morphology of particles generated by spraying polymer melts with carbon dioxide, in: G. Brunner, I. Kikic, M. Perrut (Eds.), Proceedings for the Sixth International Symposium for Supercritical Fluids, 2003, p. 1891. [10] E. Weidner, M. Petermann, Z. Knez, Multifunctional composites by high pressure spray processes, Curr. Opin. Solid State Mater. Sci. 7 (2003) 385–391. [11] N. Elvassore, M. Flaibani, A. Bertucco, P. Caliceti, Thermodynamic analysis of micronization processes from gas-saturated solutions, Ind. Eng. Chem. Res. 42 (2003) 5924–5930. [12] J. Li, H. Matos, E.G. de Azevedo, Modeling of particle formation from a gas-saturated solution process, in: G. Brunner, I. Kikic, M. Perrut (Eds.), Proceedings for the Sixth International Symposium for Supercritical Fluids, 2003, p. 1859. [13] J. Li, H.A. Matos, E.G. de Azevedo, Two-phase homogeneous model for particle formation from gas-saturated solution processes, J. Supercrit. Fluids 32 (2004) 275–286.

125

[14] M. Rantakyla, M. Jantti, O. Aaltonen, M. Hurme, The effect of initial drop size on particle size in the supercritical antisolvent precipitation (SAS) technique, J. Supercrit. Fluids 24 (2002) 251–263. [15] A.K. Jasuja, Airblast atomization of alternative liquid petroleum fuels under high pressure conditions, Trans. Am. Soc. Mech. Eng. 103 (1981) 514– 518. [16] J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, UK, 1984. [17] L.I. Rubinstein, The Stefan Problem, American Mathematical Society, Providence, Rhode Island, 1971. [18] W. Sirignano, Fluid Dynamics and Transport of Droplets and Sprays, Cambridge University Press, Cambridge, UK, 1999. [19] J.O. Werling, P.G. Debenedetti, Numerical modeling of mass transfer in the supercritical antisolvent process, J. Supercrit. Fluids 16 (1999) 167– 181. [20] J.O. Werling, P.G. Debenedetti, Numerical modeling of mass transfer in the supercritical antisolvent process: miscible conditions, J. Supercrit. Fluids 18 (2000) 11–24. [21] M. Henczka, J. Baldyga, B.Y. Shekunov, Modelling of spray-freezing with compressed carbon dioxide, Chem. Eng. Sci. 61 (2006) 2880–2887. [22] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, John Wiley and Sons, New York, 2002. [23] P.L. Chambr`e, On the dynamics of phase growth, Q. J. Mech. Appl. Math. 9 (1956) 224–233. [24] M.W. Formo, Physical properties of fats and fatty acids, in: D. Swern (Ed.), Bailey’s Industrial Oil and Fat Products, John Wiley and Sons, New York, 1979, p. 177. [25] C.A. Lockermann, E.-U. Schlunder, The effect of concentration on the liquid phase diffusivities in the binary systems carbon dioxide–oleic acid, carbon dioxide–methyl myristate and carbon dioxide–methyl palmitate at high pressures, Chem. Eng. Process. 35 (1996) 29–33. [26] W.J. McManamey, J.M. Woollen, The diffusivity of carbon dioxide in some organic liquids at 25◦ and 50 ◦ C, AIChE J. 19 (1973) 667–669. [27] B. Wagner, G.M. Schneider, Investigation of phase transitions in pure and gas-saturated triglycerides up to 20 MPa by thermal analysis according to Smit, Termochim. Acta 274 (1996) 179–186. [28] C. Hirsch, Numerical Computation of Internal and External Flows, John Wiley and Sons, Chichester, 1988. [29] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, New York, 1980. [30] W.E. Schiesser, The Numerical Method of Lines, Academic Press, San Diego, California, 1991. [31] R.M. Furzeland, A comparative study of numerical methods for moving boundary problems, J. Inst. Maths Applics 26 (1980) 411–429. [32] J. Bellan, Supercritical (and subcritical) flow behaviour and modeling: drops, streams, shear and mixing layers, jets and sprays, Prog. Energy Combust. Sci. 26 (2000) 329–366. [33] E. Sher, C. Elata, Spray formation from pressure cans by flashing, Ind. Eng. Chem. Process Des. Dev. 16 (2) (1977) 237–242. [34] S.D. Sovani, P.E. Sojka, A.H. Lefebvre, Effervescent atomization, Prog. Energy Combust. Sci. 27 (2001) 483–521.