G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. of Supercritical Fluids xxx (2014) xxx–xxx
Contents lists available at ScienceDirect
The Journal of Supercritical Fluids journal homepage: www.elsevier.com/locate/supflu
Numerical simulation of particle formation in the rapid expansion of supercritical solution process Jiewei Liu ∗ , Gustav Amberg, Minh Do-Quang Department of Mechanics, The Royal Institute of Technology, 100 44 Stockholm, Sweden
a r t i c l e
i n f o
Article history: Received 18 June 2014 Received in revised form 21 August 2014 Accepted 22 August 2014 Available online xxx Keywords: Supercritical fluid Rapid expansion Particle formation Method of moments Nucleation Condensation Coagulation
a b s t r a c t In this paper, we numerically study particle formation in the rapid expansion of supercritical solution (RESS) process in a two dimensional, axisymmetric geometry, for a benzoic acid + CO2 system. The fluid is described by the classical Navier–Stokes equation, with the thermodynamic pressure being replaced by a generalized pressure tensor. Homogenous particle nucleation, transport, condensation and coagulation are described by a general dynamic equation, which is solved using the method of moments. The results show that the maximal nucleation rate and number density occurs near the nozzle exit, and particle precipitation inside the nozzle might not be ignored. Particles grow mainly across the shocks. Fluid in the shear layer of the jet shows a relatively low temperature, high nucleation rate, and carries particles with small sizes. On the plate, particles within the jet have smaller average size and higher geometric mean, while particles outside the jet shows a larger average size and a lower geometric mean. Increasing the preexpansion temperature will increase both the average particle size and standard deviation. The preexpansion pressure does not show a monotonic dependency with the average particle size. Increasing the distance between the plate and the nozzle exit might decrease the particle size. For all the cases in this paper, the average particle size on the plate is on the order of tens of nanometers. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Fine and monodisperse particles are important in many applications, such as production of pharmaceuticals, dyes and ceramic processing, pigments in paint formulation [1], as well as surface coating [2–4]. There are conventional technologies to reduce the particle size, but these have drawbacks. One is mechanical methods like grinding, milling and crushing, but many solids are difficult to comminute mechanically. Another is the equilibrium-controlled method like liquid phase crystallization, for which a high purity product is hard to obtain. In addition, the sizes of particles obtained from those conventional methods are always nonuniform [1]. The rapid expansion of supercritical solutions (RESS) process, which has received increased attention in the past years, is an innovative and promising method to produce small particles with narrow size distribution [5]. In the RESS process, a supercritical solution at moderate pressure is sprayed from a nozzle into an expansion chamber where the pressure is much lower, and thus the solution expands rapidly. The solvating strength of supercritical fluids is directly related to the fluid density. Generally speaking, the higher
∗ Corresponding author. Tel.: +46 8 790 75 81. E-mail address:
[email protected] (J. Liu).
the density, the stronger the solvating power [6,7]. The rapid expansion in the RESS process leads to a rapid decrease in the solvent’s density and solvating strength, resulting in solute precipitation and the formation of fine particles. Those particles will be collected on a plate which is positioned in front of the spray nozzle. The most common supercritical solvent is carbon dioxide (CO2 ). It has several advantages. First, it has low critical temperature (Tc = 304.1282 K) and mild critical pressure (pc = 7.3773 MPa), which makes its supercritical state easy to reach. Second, it is non-toxic, non-flammable and environmentally friendly, it is also widely available and low-cost. Third, it is in a gaseous state under atmospheric conditions, so the final product is solvent free and has a high purity. Experimental research on the RESS process has been conducted for over 20 years [5], but numerical simulation of the RESS process is still in an early state, and most of the published studies are 1-dimensional (1D). A recent review of numerical investigations regarding supercritical fluid expansion can be found in the work of Moussa and Ksibi [8]. In the following, we are going to briefly review some previous numerical works, first for 1D, and then for 2D. Debenedetti, Kwauk and others investigated particle formation during the partial expansion of a dilute supercritical mixture in a convergent nozzle under subsonic conditions both experimentally and numerically [1,9,10]. They modeled the flow field as a
http://dx.doi.org/10.1016/j.supflu.2014.08.033 0896-8446/© 2014 Elsevier B.V. All rights reserved.
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model SUPFLU-3088; No. of Pages 16 2
ARTICLE IN PRESS J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
one-dimensional steady flow with Peng–Robinson equation of state (EoS), and the particles by the general dynamic equation (GDE) equation. Particle nucleation and condensation is included in this model, while coagulation is neglected since the typical residence time within their simulation domain is much less than the estimated coagulation times. They showed that increasing the pre-expansion temperature would increase particle size while increasing extraction temperature would decrease particle size. They also showed that small changes in surface tension could result in large variations in the predicted particle size. Weber, Russell, and Debenedetti included coagulation mechanics for particle growth in later publications [11], in which an aerosol dynamic equation was solved with a sectional method along an expansion device. But the particles’ sizes found in their experiments were much larger than their numerical results, which might be due to limitations of their model. ¨ In a series of papers [12,5,13,14], Turk, Helfgen and others studied particle formation in the RESS process both experimentally and numerically, for various solutes in different supercritical solvents. The solution they considered is very dilute. They modeled the pure solvent by a one-dimensional, steady state fluid model, and the extended generalized Bender (egB) equation of state was used to close the system. The particle formation process was modeled by the GDE, considering homogenous nucleation, condensation and coagulation effects all together. The GDE was solved by the moment’s method. Their simulation domain included the convergent capillary inlet, the straight nozzle and the large expansion chamber part. They showed [5] that the supersaturation and nucleation rate was very sensitive to the solubility and surface tension of the solute-solvent system, and the classical nucleation theory was not suitable to describe the trend in particle size. The general trends of their numerical results [14] showed good agreement with experiments, but the predicted particle sizes didn’t match the measured particle sizes exactly. Their later work [13] showed that particle formation occurs mainly in the supersonic free jet and coagulation was the main mechanism for particle growth. Other one-dimensional numerical works include Lele et al. [15], Reverchon and Palado [16], Weber and Thies [17], etc. 1D simulations are good attempts to uncover and understand physics in RESS, like the thermodynamic conditions along the centerline. But they might not be able to completely represent the physics of the flow since it cannot capture the processes that happen off-axis, like the large gradient of the thermal conditions in the oblique shock originating from the lip of the nozzle exit. Thus, 1D simulations might be insufficient if we want to explain the experimental results well. Ksibi and Moussa’s group published several papers about numerical simulations on the RESS process in 2D, axisymmetric geometry. Ksibi et al. [18] investigated the rapid expansion of pure supercritical CO2 by modeling it as a Newtonian, viscous compressible fluid, governed by the Navier–Stokes equation, and satisfying the Altumin and Gadetskii equation of state [19]. The system of equations were solved by the Total Variation Diminishing scheme using a finite difference method, with Roe averaged technique [20] for shock capture. Their studies focused on the profiles of thermodynamic variables along the centerline and the normal plate, when the pure supercritical carbon dioxide expanded into the expansion chamber which contained motionless pure carbon dioxide. Following Ksibi’s work, Moussa et al. [21] carried out a parametric study of the nozzle geometry to control the supercritical fluid expansion for carbon dioxide. They showed that the geometry of the nozzle and operating condition would highly influence the jet structure. The aforementioned 2D simulations did not consider particle formation during the RESS process. In 2008, Moussa et al. [22] incorporated particle formation in their two dimensional, axisymmetric supersonic jet during the RESS process. Particle transport and growth were governed by the time dependent aerosol GDE,
without considering nucleation and condensation. A lognormal size distribution function was assumed for the particles and the GDE was solved with a sectional method. They found that particles deposit on the flat plate as a ring since particles located at the jet boundary are subjected to more pronounced Brownian coagulation. They also found that a larger orifice and lower expansion pressure leads to a narrower distribution of particles on the flat plate. Other 2D numerical simulation works on RESS are: Franklin et al. [23] numerically simulated the rapid expansion of the perfluoropolyether–carbon dioxide supercritical solution into air by solving an axisymmetric, Favre-averaged compressible Navier–Stokes equations which was suitable for turbulent flow simulations. The ideal-gas law was applied for the air, and gas-phase mixing rule was applied between CO2 and air. The latticefluid state equation of Sanchez and Lacombe was used for both single-phase CO2 and the solute perfluoropolyether diamide, which was in a liquid state at the conditions of interest. Separate continuity equations were used for air and CO2 . Khali et al. [24] studied the structure of supercritical jet expansion for in viscid, adiabatic pure carbon dioxide flow impacting on a flat plate, but particle formation was not considered. They solved an axisymmetric Euler equation with the Redlich–Kwong equation of state [25], using a two-step Lax–Wendroff finite difference method. They showed that their numerical results compared reasonably well with their experiment if the pressure was not very high. They proposed that a divergentconvergent capillary could be more suitable for the RESS process than a convergent nozzle. Dea et al. [26] investigated the growth of magnetic thin films using carbon dioxide RESS expansion, following the 2D work of Khali et al. [24]. They applied a quasi-one dimensional analysis to thermal properties on the jet centerline, and a kinematic model for the particle formation, which was based on bimolecular collisions. Numerical studies of particle formation in the RESS process are still in an embryonic state [22]. Existing numerical simulations can only roughly interpret the experimental results, sometimes they even contradict experimental results [5]. In order to gain a better understanding of this process and provide some guidance in experimental design to produce small and uniform particles, we numerically study particle formation in the RESS process in a 2D, axisymmetric geometry for CO2 -benzoic acid system. The fluid is described by the classical Navier–Stokes equation, with the usual pressure term being replaced by the generalized pressure tensor [27]. The equation of state is the egB equation of state [5,28]. For the particle model, following previous work of Pratsinis [29] and Helfgen et al. [14], the method of moments is used to solve the general dynamic equation, which accounts for the aerosol dynamics during a simultaneous particle nucleation, condensation and coagulation process. We are going to first study a typical case in detail, and then do a parametric study, to see how things change when we change the preexpansion temperature, the preexpansion pressure and the plate position. The rest of the paper is organized as follows: Section 2 describes the models, including the fluid dynamic equations in Section 2.1 and the particles model in Section 2.2, as well as the moment equations in Section 2.3. In Section 3, we talk about the numerical method, show some numerical results, and discuss some issues. Section 4 is the conclusion. Section 5 is the acknowledgement. Section A is the appendix. 2. Problem formulation The governing equations consist of two parts. One is the fluid dynamic equations, which describe the evolution of density, temperature, pressure and velocity of the flow. The other is the general dynamic equation, which models the evolution of the number density of the particles in the whole domain. The general dynamic
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
3
equation will be solved using the method of moments, using models of the nucleation rate, growth rate, and coagulation coefficients for particles in the free molecular regime and continuum regime, respectively. It will build on an assumption of a lognormal distribution of particle size. As the solution is very dilute, the influence of particles on the fluid flow is negligible. Thus, the fluid dynamic equations are solved for the pure solvent, and this is used as input to the general dynamic equation for the evolution of the particles. In the following, variables with or without ∼ are dimensional and dimensionless, respectively.
In (1)–(6), = / ˜ c is the dimensionless density, T = T˜ /Tc is the dimensionless temperature, v is dimensionless fluid velocity. The coefficients Ai = Ai (1/T, ej ) in (6) are functions of 1/T and ej . The ej are given as correlations and described by 95 general constants and two substance-specific parameters [28]. The total energy eT is the sum of internal energy e and kinetic energy,
2.1. Fluid dynamic equations
e = h − p,
The fluid dynamic equations are based on the conservation of mass, momentum and energy. During the RESS process, the fluid transits from the supercritical state to a two-phase state or a gas state. A good fluid model should be able to capture phase transition if this occurs, and consider the contributions from the phase interface. In 1893, van der Waals introduced a gradient term in the Helmholtz free energy density to describe a gas-liquid interface [30,31] for isothermal, one-component fluids. Such a gradient term began to be widely used since the work of Cahn and Hilliard [32], who studied the interfacial structure for isothermal, binary alloys by adding a gradient term to Helmholtz free energy. Later, Onuki [27] extended van der Waals theory to dynamic processes by considering that phase transitions, for a variety of situations, occur in inhomogeneous temperature or heat flow. He included gradient contributions in both the entropy and the internal energy, leading to a generalized pressure tensor term in the system of equations. We will use this model to calculate the flow field and thermodynamic variables, but instead of van der Waals equation of state, we use the egB equation of state, since previous studies [5,14,33] showed that the cubic equation of state is not suitable for supercritical expansion calculations. With the choice of characteristic variables * = c , p* = pc , T* = Tc , L* = 1m, and v∗ = vs , where (c , pc , Tc ) are taken at the critical point of the fluid, m is the unit of nozzle diameter, vs is the sound speed at the inlet, we obtain the non-dimensional form of the fluid dynamic equations as follows: Mass conservation,
and enthalpy h can be calculated from the equation of state (6). More details can be found in [33]. The dimensionless parameters that appear here are: Reynolds number Re = L∗ v∗ /0 , which measures the relative importance of inertia to viscous stresses, 0 comes from the assumption that ˜ [27]; the dynamic viscosity depends linearly on density = 0 the parameter R = p∗ /∗ v∗2 = a/9mbv∗2 , which is proportional to the ratio of the attractive potential energy to the molecular kinetic energy, the parameter B = 0 ∗ T ∗ /e∗ v∗ L∗ = 8m 0 /3kB v∗ L∗ , the parameter = CT ∗ ∗2 /2m2 p∗ L∗ 2 ; the normalized gravity acceleration G = (mgL* )/a/9b [27]; and the dimensionless parameter C = Rs ∗ T ∗ /p∗ .
∂ + ∇ · (v) = 0, ∂t
(1)
Momentum equation,
∂(v) ˆ + 1 ∇ · ˆ − GRez , + ∇ · (vv) = −R∇ · Re ∂t
(2)
∂eT ˆ − Re−1 R−1 ) + ∇ · (eT v) = −∇ · [( ˆ · v] + B∇ · (∇ T ) − Gvz , ∂t (3) ˆ and viscous stress tensor where the generalized pressure tensor ˆ are
ˆ = (∇ v + ∇ v ) + 1 −
2 3
ˆI ∇ · v.
(4) (5)
Thermodynamic pressure p in (4) is calculated from the egB equation of state, p = CT (1 + A1 + 2 A2 + 3 A3 + 4 A4 + 5 A5 2
2
2
+ (A6 + A7 ) exp(− )).
1 v2 , 2R
(7)
where internal energy e is related to enthalpy h and pressure p of the system, (8)
2.2. General dynamic equation Pratsinis [29] studied simultaneous aerosol nucleation, condensation and coagulation in an isothermal, uniform aerosol reactor, covering the entire particle size range. More than 10 years later, Helfgen et al. [13,14] adopted Pratsinis’s model, to study particle formation during the RESS process in a 1D geometry. They added particle transportation to Pratsinis’s model, and assumed a steady state. We know that nonuniform and nonisothermal conditions will appear in the whole domain during the expansion. Changes that happen off-axis are very important and will have an influence on the final product, but cannot be captured by 1D simulation. In this section we will follow previous studies [34,29,14] to describe particle formation in an axisymmetric, 2D domain in detail. The general dynamic equation for simultaneous aerosol nucleation, condensation, coagulation and convection is [34,29,14] ˜ V˜ , t) ∂n( ˜ V˜ , t)) + · (˜vn( ∂t = ˜J (V˜ c )ı(V˜ − V˜ c ) −
Total energy equation,
ˆ = (p − T (∇ )2 − 2T∇ 2 )ˆI + 2T ∇ ∇ ,
eT = e +
+
1 2
×
V˜
˜ V˜ , t)) ∂(G˜ n( ∂V˜
˜ V˜ − V˜¯ , V˜¯ )n(V˜ − V˜¯ , t)n(V˜¯ , t)dV˜¯ − n(V˜ , t) ˇ(
0 ∞
˜ V˜ , V˜¯ )n(V˜¯ , t)dV˜¯ , ˇ(
(9)
0
˜ ∂t describes the rate of change of particle concenThe term ∂n/ ˜ represents the tration (or the number density). The term · (˜vn) loss or gain of particles due to convection. The term ˜J (V˜ c )ı(V˜ − V˜ c ) accounts for the formation of new particles of critical volume ˜ n) ˜ ∂V˜ V˜ c by homogeneous nucleation at rate ˜J . The term −∂(G/ accounts for the loss or gain of particles by condensation at rate ˜ The last two terms on the RHS account for the gain and loss of G. V˜ ˜ V˜ − V˜¯ , V˜¯ )n(V˜ − particles by Brownian coagulation. The term ˇ( 0
(6)
V˜¯ , t)n(V˜¯ , t)dV˜¯ describes the gain of a new particle of volume V˜ by coagulation of a particle of volume V˜ − V˜¯ with a particle of volume
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
4
V˜¯ . This term has a coefficient 1/2 which is due to the fact that every ∞ ˜ V˜ , V˜¯ )n(V˜¯ , t)dV˜¯ particle is considered twice. The term −n(V˜ , t) ˇ( 0
accounts for the loss of a particle of volume V˜ by coagulation of a particle of volume V˜ with any other particle of volume V˜¯ . In order to solve Eq. (9), three terms on the right hand side: ˜ and coagulation coefficient nucleation rate ˜J , condensation rate G ˜ need to be specified on the whole particle size spectrum. They ˇ will be given in the following subsection, following previous work ¨ [5], Helfgen et al. [14] and Pratsinis [29]. of Turk 2.2.1. Nucleation There are two different nucleation mechanisms: heterogeneous nucleation and homogenous nucleation. The first one occurs when the system is able to overcome its activation free energy barrier with the help of foreign particles or wall surfaces. The latter occurs in very high purity systems where there are no foreign particles or walls, or in very rapid expansions [35]. Homogeneous nucleation is the mechanism considered in this paper. All homogenous nucleation models, although different in detail, start from the general expression ˜J = Ke
− F˜ c /kB T˜
,
(10)
where F˜ c is the free energy change involved in the formation of one critical particle of size r˜ c . The exponential term dominates the behavior of the nucleation rate. kB = 1.3806488e − 23 J/K is the Boltzmann’s constant, T˜ is the absolute temperature. The pre-exponential coefficient K can vary a lot in magnitude between different theories [35]. To form a stable embryo through homogenous nucleation, the system has to overcome a free energy activation barrier F˜ [36,35],
F˜ = 4˜r 2 −
3 3 ˜ d Rs T˜ ln S. ˜r 4
(11)
The positive contribution to the energy barrier 4˜r 2 , represents the surface free energy of an embryo of radius r˜ , with denoting the surface free energy per unit area. The negative contribution to the energy barrier − 43 ˜r 3 R ˜ s T˜ ln S is due to the bulk ˜ d is the density of the condensed phase, T˜ is free energy change. the absolute temperature, S is the supersaturation, Rs is the specific gas constant, Rs = NA kB /Md , NA = 6.02E23 mol−1 is the Avogadro constant, Md is the molar mass of the condensed phase. For S = 1, the second term disappears and F˜ increases monotonically with r˜ . For S < 1, the second term is positive and make F˜ increases even more steeply. For S > 1, the negative contribution from the second term implies the existence of a maximum in F˜ at some critical r˜ c , and this maximum represents the summit of the free energy activation barrier. Any additional growth of the particle beyond the size r˜ c will lead to irreversible growth of the particle at the expense of the supersaturation [36]. ¨ [5] rewrites formula (11) for Using properties of the solute Turk the RESS process as: 4
2
F˜ = 4˜r G −
3V2,molecule
3
r˜ kB T˜ ln S.
2G V2,molecule . kB T˜ ln S
F˜ c =
(13)
16 2/3 G V2,molecule kB T˜ 3
3
1
×
(ln S)2
.
(14)
Given the critical radius r˜ c (13), the height of the free energy activation barrier F˜ c (14), and the preexponential coefficient K, explicitly. Helfgen the nucleation rate ˜J in (10) can be expressed et al. [14] take K = ˛C V2,molecule 22 2G /m2,molecule , which is also adopted in this paper. So,
˜J = ˛C V2,molecule 2 2
2G m2,molecule
exp
−
F˜c kB T˜
,
(15)
Here m2,molecule is the molecular mass of the solute, 2 is the number density of dissolved molecules, which can be calculated as 2 = /M ˜ ˜ is the density of the mixture, Mw is the molar w yE NA [5]. mass of the mixture, yE is the solute mole fraction at extraction conditions. is the non-isothermal factor, which equals 1 in the case of a dilute mixture, and ˛C is the condensation coefficient, which is set to 0.1 [5]. 2.2.2. Solubility and supersaturation The supersaturation S in formula (12) is given by Helfgen et al. [14] as follows S=
ye (T˜e , p˜ e )(T˜ , p˜ , ye ) . y∗ (T˜ , p˜ )(T˜ , p˜ , y∗ )
(16)
ye (T˜e , p˜ e ) is the mole fraction of the solute at extraction temperature and pressure, y∗ (T˜ , p˜ ) is the equilibrium mole fraction of the solute at temperature and pressure of the flow field, (T˜ , p˜ , y) is the solute fugacity coefficient calculated by Eq. (19). Eq. (16) is valid no matter whether the solute is dilute or not [5]. As a measurement of solubility, the mole fraction of the solute y2 can be calculated as follows at a given pressure p and temperature T [37],
vp
y2 =
p˜ 2
p˜ 2 (T˜ , p˜ , y1 , y2 )
exp
v2s p˜ RT˜
,
(17)
where R is the universal gas constant, vs2 is the molar volume of the pure solid solute, vs2 = 0.965E − 4 m3 /mol for benzoic acid at vp 122 ◦ C [37], p2 is the vapor pressure of the pure solid solute at the temperature of interest, which can be calculated with the Antoine equation ln
A−B p˜ vp . = Pa T˜
(18)
Eq. (17) must be solved iteratively because the fugacity coefficient 2 is a function of y2 . In the dilute mixture, the expression for 2 is given by the following formula [37], which is calculated using the PR-EoS. ln 2 (T˜ , p˜ , y1 , y2 ) =
(12)
Here V2,molecule is the molecular volume of the solute, and G is the solid–fluid interfacial energy. Since no experimental data is available on G , an estimated value of 0.02 N/m is used in what follows, the same as used in [5,14]. From formula (12), the critical radius r˜ c , the size that a cluster must achieve to be able to steadily grow, can be expressed as [5] r˜ c =
The corresponding free energy activation barrier is then [5]
b2 (Z − 1) a p˜ (V˜ − b) −√ − ln ˜ b RT 8RTb × ln
a
−
b2 b
V˜ + 2.414b , V˜ − 0.414b
and a1 = 0.4572
R2 Tc2 Pc
b1 = 0.07780 1
2(a y + a y ) 12 1 2 2
1+(0.3746 + 1.5423ω − 0.2699ω2 )
RTc , Pc 1
(19)
1 2 2 ˜ 1−
T Tc
,
(20)
(21)
2
a = (a12 y1 + a22 y2 ) ,
(22)
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
5
where
−1
10
B3 =
(482 )
1/3
a = 7.79232589 Pa*m6/mol2
˜ 2,molecule 2,s (8KB T˜ /m2,molecule )1/2 V . 3
(28)
2
b2 = 1.06589928*10−4 m3/mol
−2
2.2.4. Coagulation The collision frequency function in the free molecule regime is [29,14]
10
y2
˜ fm = B2 ˇ −3
10
1 1 + V˜ V˜¯
1/2
2
1/3 (V˜ 1/3 + V˜¯ ) ,
(29)
where B2 = (3/4) 308.15 K, exp 318.15 K, exp 328.15 K, exp 343.15 K, exp 308.15 K, cal 318.15 K, cal 328.15 K, cal 338.15 K, cal
−4
10
−5
10
100
150
200
250 p (bar)
300
350
400
b = b1 y1 + b2 y2 , 1
a12 = (a1 a2 ) 2 .
dV˜ dt
kB T˜ 2m2,molecule
1/3 (V˜ 1/3 + V˜¯ ),
(31)
where 2KB T˜ , 3 ˜
(32)
˜ B5 , r˜
(33)
2.3. Moment equations Several different methods have been used to solve the GDE, like the sectional method, the “J-Space” method, and the moment’s method. The first two are more complicated to code, and may have difficulties in handling condensation [38]. For simplicity, we use the moment’s method in this paper to solve the general dynamic equation. The kth moment of an aerosol size distribution is
∞
˜ V˜ , t)dV˜ . V˜ k n(
˜ k (t) = M
(34)
0
Since the condensation rate G and the coagulation coefficient ˇ in the free molecular regime and continuum regime are different, the moments equation will be considered in both regimes, in order to analyze the particle size in the entire spectrum [29,14]. 2.3.1. Moment equations in the free molecule regime With the condensation rate Gfm and coagulation coefficient ˇfm specified, we can calculate the first three moments in the free molecular regime [29]. ˜0 dM ˜ −1/2 + 2M ˜ −1/6 + M ˜ 0 ], ˜ 0 ) = ˜J − b0 B2 [M ˜ 2/3 M ˜ 1/3 M ˜ 1/6 M + · (˜vM dt (35)
1/2 .
(26)
2,s = /M ˜ w y2 NA is the number density of dissolved molecules at saturation, which changes along the RESS path. In the continuum regime, the condensation rate is [29,14] ˜ ˜ c = dV = B3 V˜ 1/3 (S − 1), G dt
(30)
(25)
where
C(V˜ ) C(V˜¯ ) + 1/3 V˜ 1/3 V˜¯
C(V˜ ) is the Cunningham correction factor and B5 = 1.257.
fm
B1 = (36)1/3 V2,molecule 2,s
˜ c = B4 ˇ
(24)
= B1 V˜ 2/3 (S − 1),
.
C(V˜ ) = 1 +
2.2.3. Condensation In the free molecular regime, the condensation rate for a particle is [29,14]
1/2
(23)
Z = p˜ V˜ mix /RT˜ represents the compressibility coefficient of the mixture, V˜ is the molar volume of the mixture. a1 and b1 are solvent parameters which can be calculated with its critical data and acentric factor ω. a2 , b2 are solute parameters which can be obtained through nonlinear regression of the experimental data. For a benzoic acid/CO2 solution, we obtain a2 = 7.79 Pa m6 /mol2 , b2 = 1.07E − 4 m3 /mol when using non-linear least squares to fit the experimental data in [37]. Schmitt and Reld obtained a2 = 9.23 Pa m6 /mol2 , b2 = 1.34E − 4 m3 /mol. The sum of squared residuals for 39 experimental data, a residual being the difference between an experimental value and the fitted value, is 0.00113810269 in our model and 0.0034544133 in Schmitt and Reld’s model. Fig. 1 shows the mole fraction of benzoic acid in CO2 at different temperature and pressure when taking a2 = 7.79 Pa m6 /mol2 , b2 = 1.07E − 4 m3 /mol, the solubility calculated with our numerical scheme matches well with the experimental results.
˜ fm = G
(6KB T˜ V2,molecule /m2,molecule )
The collision frequency function in the continuum regime is [29,14]
B4 =
Fig. 1. Solubility of benzoic acid in CO2 .
1/6
(27)
˜1 dM ˜ 1 ) = ˜J V˜ c + B1 (S − 1)M ˜ 2/3 , + · (˜vM (36) dt ˜2 dM ˜ 1/2 ˜ 2 ) = ˜J (V˜ c )2 + 2B1 (S − 1)M ˜ 5/3 + 2b2 B2 [M ˜ 5/3 M + · (˜vM dt ˜ 4/3 M ˜ 7/6 M ˜ 5/6 + M ˜ 1 ]. + 2M
(37)
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
6
2.3.2. Moment equations in the continuum regime With the condensation rate Gc and coagulation coefficient ˇc specified, we could calculate the first three moments in the continuum regime [29]. ˜0 dM ˜ −1/3 ˜ 0 ) = ˜J − B4 [M ˜ 2+M ˜ 1/3 M + · (˜vM 0 dt
4 1/3
+ B5
˜1 dM ˜ 1 ) = ˜J V˜ c + B3 (S − 1)M ˜ 1/3 , + · (˜vM dt
(39)
˜2 dM ˜ 2/3 ˜ 2 ) = ˜J (V˜ c )2 + 2B3 (S − 1)M ˜ 4/3 + 2B4 [M ˜ 2+M ˜ 4/3 M + · (˜vM 1 dt
4 1/3
+ B5
˜ 2/3 + M ˜ 4/3 )]. ˜ 1M ˜ 1/3 M (M
3
1 exp √ ˜ V (t)3 2 ln g (t)
−
density
n(V˜ (t); V˜g (t), g (t), t)
N(t) exp = √ ˜ V (t)3 2 ln g (t)
−
(ln V˜ (t) − ln V˜g (t))
2
distribution
(ln V˜ (t) − ln V˜g (t)) 18(ln g (t))2
.
function
2
(41) n(t) =
,
˜ 0 (t) = N(t) ˜ M represents number concentration, 2 9 ˜ ˜ V˜g (t)e 2 (ln g (t)) N(t)E[ V˜ (t)] = N(t) represents volume
˜ 1 (t) = M
fraction, 2 2 ˜ ˜ 2 (t) = N(t)E[ ˜ and M V˜ (t)2 ] = N(t)( V˜g (t)) e18(ln g (t)) is the sec Var[V˜ (t)] = ond moment. The standard deviation s.d.[V˜ (t)] = 2
E[V˜ (t)2 ] − (E[V˜ (t)]) = E[V˜ (t)] e9(ln g ) − 1. With the lognormal distribution of particle size, we can now rewrite the moments equation in the free molecular regime (35)–(37) as follows (considering the variation of density, from here ˜ k to represent the original M ˜ k , so M ˜ k in the following on we use ˜M represents the quantity per unit mass): ˜ 0) d( ˜M ˜ 0) + · (˜v ˜M dt 2
˜ 0 ) (V˜g ) ˜M = ˜J − b0 B2 ( + 2 exp
5 8
ln2
1/6
exp
+ exp
1 8
25 8 ln2
ln2
,
exp
25 8
ln2
ln2
+ 2 exp
5 8
ln2
.
(46)
˜ 0) d( ˜M ˜ 0) + · (˜v ˜M dt
˜ 0) = ˜J − B4 ( ˜M
1 2
2
2
r˜g
1 + exp(ln2 ) + B5
2
ln (1 + exp(2ln )) ,
(47)
˜ 1) d( ˜M ˜ 1 ) = ˜J V˜ c + B3 (S − 1) ˜ 0 V˜ g1/3 exp + · (˜v ˜M ˜M dt
1 2
ln2 , (48)
2
1/3
˜ 1 V˜ g = ˜J V˜ c + 2B3 (S − 1) ˜M
exp
7 2
4 1/3
× 1 + exp(ln2 ) + B5
3
ln2
(V˜ g )
˜ 1) + 2B4 ( ˜M
−1/3
2
5 exp − ln2 2
.
(49)
(42)
(43)
8
+ exp − ln2 2
9 ˜ k (t) = M ˜ 0 (t)Vg (t)k exp( k2 ln2 g (t)), M 2
2
2
1
ln2
1
˜ where N(t) is the number density at time t. According to formula (77), the kth moment (34) satisfies
3
˜ 1 )2 V˜ g1/6 exp(8ln2 ) + 2b2 B2 ( ˜M
˜ 2) d( ˜M ˜ 2) ˜M + · (˜v dt
18(ln g (t))2
2/3
The moments equations in the continuum regime (38)–(40) are as follows:
× exp
2.3.3. Lognormal distribution and moments equation Assume a log-normal distribution for the particle size V˜ ,
Since the number ˜ we have f (V˜ (t))N(t),
+ exp
(40)
The fractional moments that appear in formula (35)–(40) can be expressed through the zeroth moment, geometric mean and geometric standard deviation if a lognormal distribution function is assumed for the particle sizes.
f (V˜ (t)) =
2
˜ 1 V˜ g ˜M = ˜J V˜ c + 2B1 (S − 1) × exp
˜ 0M ˜ 1/3 M ˜ −1/3 + M ˜ −2/3 )], (38) (M
3
˜ 2) d( ˜M ˜ 2) ˜M + · (˜v dt
(44)
˜ 1) ˜M d( ˜ 1 ) = ˜J V˜ c + B1 (S − 1) ˜ 0 V˜g 2/3 exp(2ln2 ), (45) + · (˜v ˜M ˜M dt
2.3.4. Dimensionless moments equation and harmonic average Helfgen et al. did not discuss about the scaling, and the expressions for the coagulation and condensation coefficients in the dimensionless moment equations were not given in their paper. The characteristic variables that we choose is different from Pratsinis’s choice [29]. The reason is that fluid flow was not considered in Pratsinis’s model, so they could choose the time scale to be the characteristic time for particle growth, which was taken to be 45 s in their paper. The characteristic length and velocity we used to nondimensionalize the fluid equations have already implicitly defined a characteristic time, t ∗ = L∗ /v∗ , and v∗ ∼300 m/s, so t* is on the order of 10−9 s in our case. And also, they chose the characteristic length to be the radius of the monomer, which was usually on the order of 10−9 m. Our smallest simulation domain is almost [2000 m, 2000 m], so the length scale we choose is 10−6 m. Pratsinis studied an isothermal system, so the monomer concentration at supersaturation was constant and was used to scale the number density. In the RESS process, the temperature varies over the domain, so we take the number concentration of dissolved molecules at the extraction condition to scale the ˜ E /Mw yE NA , where E is the density number density, i.e., n∗s = at extraction condition, so M0∗ = n∗s /∗ . Other characteristic variables used in our moments equations are: characteristic volume ∗ V ∗ = 4/3(L∗ )3 , characteristic first moment M1∗ = M0∗ Vol , characteristic second moment M2 = M0∗ (V ∗ )2 , characteristic nucleation rate J ∗ = n∗s /t ∗ . In this paper, the extraction condition is chosen to be TE = 388 K, PE = 15 MPa. Corresponding to this, yE = 1.926e − 3, M0∗ = 1.66e + 22 kg−1 , J* = 2.483e + 33 m3 s−1 .
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
With the above choices, we can nondimensionalize the moments Eqs (44)–(46) in the free molecular regime as follows: zeroth moment equation d(M0 ) + · (vM0 ) = J − f (M0 )2 , dt
where
c = B3
(50) =
where
f = b0 t ∗ ∗ M0∗ (L∗ rg )
× exp
25 8
1/2
ln2
6kB T ∗ TV 2,molecule m2,molecule
5
+ 2 exp
8
ln2
1/2
+ exp
8
ln2
(51)
(52)
= 4L
2/3
exp(2ln2 )
∗ 1/2 ∗ y N 2 A ∗2 kB T T ∗ t
2m1
m1
1
2
∗ y2 NA c = 4L Mw
V2,molecule V∗ (60)
(61)
Mw
7 2
(Vg )2/3
1/2 /3 t ∗ (Vg )1/3
V2,molecule V∗
ln2 ,
V2,molecule exp(2ln2 ). V∗ (53)
(62)
2kB T ∗ T ∗ ∗ ∗ c = M0 t 1 + exp(ln2 ) + B5 ∗ 30 ∗ L rg 1 + exp − ln2 2
d(M2 ) 2 + · (vM2 ) = J(V c ) + 2f (S − 1)M1 + 2f (M1 )2 , dt
(54)
where
k T ∗ T 1/2 B 2m1
8kB T ∗ T m2,molecule
and
Second moment equation
∗ y2 NA 3 V L∗ 2,molecule Mw
t ∗ (Vg )1/3
ln2 .
f =
Mw
d(M2 ) 2 + · (vM2 ) = J(V c ) + 2c (S − 1)M1 + 2c (M1 )2 , dt
× exp Vg
2
ln2
8k T ∗ T 1/2 ∗ y N B 2 A
∗
where
M1∗
4L∗ 3
1
where
d(M1 ) + · (vM1 ) = JV c + f (S − 1)M0 , dt
M0∗ (V ∗ )2/3 t ∗
(V ∗ Vg )1/3 exp
M1∗
Second moment equation .
First moment equation
f = B1
t ∗ M0∗
× exp
1
7
2/3
t ∗ Vg
exp 8ln2 ,
(55)
5 exp − ln2 2
.
(63)
Eqs. (50)–(63) are combined to describe aerosol dynamics over the entire particle size spectrum by using the harmonic average of the condensation and coagulation rates in the free molecule and continuum regimes. d(M0 ) + · (vM0 ) = J − (M0 )2 , dt
(64)
and
f = b2 (L∗ )1/2
× exp
6kB T ∗ TV 2,molecule m2,molecule
25 8
ln2
+ 2 exp
1/2 t ∗ ∗ M0∗ (Vg )1/6 exp
5 8
ln2
+ exp
1 8
ln2
3 2
ln2
1800
F
1600
.
plate
G
1400
(56)
d(M0 ) + · (vM0 ) = J − c (M0 )2 , dt
2kB T ∗ T ∗ ∗ ∗ M0 t 1 + exp(ln2 ) + B5 ∗ 30 ∗ L rg × exp
centerline
1000
open boundary
800 600
where c =
(57)
z (μ m)
1200
Converting the moments Eqs. (44)–(46) in the continuum regime into dimensionless form, we obtain the following: zeroth moment equation
1 2
400 C
ln2 (1 + exp(2ln2 )) .
200
solid wall
A B
E
(58) 0 H
First moment equation d(M1 ) + · (vM1 ) = JV c + c (S − 1)M0 , dt
D
−200
0
inlet
I 500
1000
1500
2000
r (μ m) (59) Fig. 2. RESS geometry.
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
8
d(M1 ) + · (vM1 ) = JV c + (S − 1)M0 , dt d(M2 ) 2 + · (vM2 ) = J(V c ) + 2(S − 1)M1 + 2(M1 )2 , dt
(65)
(66)
where =
f c f + c
,=
f c f + c
,=
f c f c ,= . f + c f + c
(67)
The geometric mean volume Vg and geometric standard deviation g can be calculated through the first three moments,
Vg (t) =
M1 (t)2 [M0 (t)3 M2 (t)]
9ln2 (t) = ln
1/2
,
M0 (t)M2 (t) M1 (t)2
(68)
.
(69)
3. Numerical method and results The numerical method we use to solve the fluid Eqs. (1)–(8) and the moments Eqs. (64)–(67) is exactly the same as what was used in [33]: The characteristic-based splitting scheme plus artificial viscosity to capture the shocks, using a finite element numerical toolbox femLego [39]. At every time step, the fluid equations are solved first, and then the moment equations are solved. The distribution of temperature, density and pressure obtained from the fluid equations are used to calculate the quantities in the moments equation. Fig. 2 shows the geometry considered in this paper. All calculations are made in axisymmetric geometry. HIAB is a convergent capillary inlet with a starting radius of LHI = 600 m and a length of LHA = 330 m. ABCD is a straight nozzle with a diameter of d = 2LAB = 100 m and a length of LAD = 100 m. DCEFG is an outside chamber. A plate GF is positioned at a distance of LDG to the nozzle exit CD, AB is the nozzle inlet, EF is an open boundary, CE is a solid wall, and HG is the axisymmetric axis. For most of cases in our simulation, LDG = 30r, where r is the radius of the nozzle, i.e., r = LAB . LDG = 60r and LDG = 90r are also considered.
Fig. 3. Profiles of (a) temperature T, (b) density (in log scale), and (c) velocity magnitude v on the whole domain at t = 40, 000, under preexpansion conditions T = 388 K, p = 10 MPa.
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
T (r, z, 0) =
To , if zD ≤ z ≤ zG
(r, z, 0) =
Ti , if zH ≤ z ≤ zD
Temperature along the center line
400 350
T (K)
In this section, we will discuss some numerical results. We will first describe in detail the results under preexpansion conditions T = 388 K and p = 10 MPa, where the plate is positioned in front of the nozzle at a distance of Xp = 30r (LDG = 1500 m). We will later show some results for parametric studies, to see how the particle sizes on the plate change if we change the preexpansion temperature, preexpansion pressure and the plate position. The initial distribution of temperature T and density are the following,
9
300 250 t=30000 t=40000 t=50000 t=60000
200 ,
z − zH i + (o − i ), if zH ≤ z ≤ zD zD − zH , o , if zD ≤ z ≤ zG
(70)
150
0
500
(71)
Initial conditions for the moments equations are: on the whole domain, we set geometric standard deviation g = 1.05, geometric mean particle size rg = 10 nm, dimensionless M0 = 10−3 , this leads to ˜ 0 from the inlet on the order of a variation of the number density ˜M 1021 m−3 to the expansion chamber on the order of 1019 m−3 . Inlet condition at the preexpansion chamber is M0 = 1e − 3, g = 1.05, rg = 10 nm. Calculated dimensionless parameters are: Re = 79.50, B = 0.268, R = 0.155, G = 0, = 2.94e − 6. 3.1. Thermal and hydrodynamic properties Fig. 3 shows the profiles of (a) temperature T, (b) density and (c) velocity magnitude ||v|| on the whole computational domain, at dimensionless time t = 40, 000. It shows that thermal conditions within the capillary inlet, which is from z = −130 to z = 200, are almost the same as the inlet condition. Starting from the straight nozzle, the temperature and density begin to decrease, and they continue to decrease very rapidly until the flow reaches the Mach disk, a normal shock after which the jet flow becomes subsonic. Across the Mach disk, temperature and density increase very sharply. Velocity, however, has an opposite trend. It starts to increase inside the preexpansion chamber, and continues to increase and reaches as high as more than 700 m/s before the Mach disk. Across the Mach disk, the velocity becomes very small suddenly and the flow is subsonic. In addition to the Mach disk in the normal direction, we can also see the barrel shock which originates from the nozzle exit and finally interacts with the Mach disk. Across the barrel shock, temperature and density increases sharply, and velocity decreases sharply. We can clearly see the shear layer of the jet in the profiles of temperature in Fig. 3(a) and velocity in Fig. 3(c). Within the shear layer, temperature is relatively low compared to its ambient fluid, and the flow is still supersonic. When the shear layer hits the plate, the fluid is forced to change its direction suddenly. The interaction between the shear layer and the plate makes the flow pattern
1500
Temperature along the plate
340
where Ti , To (i , o ) are temperature (density) at the inlet and outside chamber, respectively. The boundary conditions of the flow are:
320 T (K)
1. At the inlet HI: p = pi , = i , T = Ti . 2. At the plate GF, the solid wall CE, and the nozzle wall IB, BC: v = 0, no mass flux, homogenous Numan boundary condition for T and p. 3. At the open boundary EF: homogenous Numan boundary condition for , T, p. 4. At the centerline AG: v · n = 0, n is the normal vector. Homogenous Numan boundary condition for , T, p.
1000 z (μ m) (a)
300
t=30000 t=40000 t=50000 t=60000
280
260 0
500
1000 r (μ m) (b)
1500
2000
Fig. 4. Profiles of temperature along (a) the centerline and (b) the plate at different time.
after the Mach disk, as well as the thermal conditions on the plate, complicated. As Fig. 4(a) shows, from t = 40, 000 to t = 60, 000, the position of the Mach disk stays at z = 850, and the fluid before the Mach disk has reached a steady state since the temperature (and also other variables) is time independent. However, the fluid after the Mach disk always changes with time. Fig. 4(b) shows the temperature profiles along the plate from t = 40, 000 to t = 60, 000. Due to the impact of the jet boundary on the plate, those profiles are oscillating and change with time.
3.2. Solubility and supersaturation The solubility y2 is a function of the temperature and pressure. The rapid decrease in temperature and pressure before the Mach disk leads to a rapid decrease in solubility y2 (Fig. 5(a)), and thus a rapid increase in supersaturation S (Fig. 5(b)). Profiles of S along the centerline in Fig. 5(c) tells us that, at the inlet of the straight nozzle z = 200, the solution is just saturated with S = 1. Inside the straight nozzle, the solution becomes more saturated. And at the nozzle exit z = 300, supersaturation S achieves 103 , so particle precipitation inside the nozzle might not be ignored. The solubility y2 decreases to as low as 10−15 before the Mach disk, and the supersaturation increases to as high as 1012 . After the Mach disk, the sharp increase in temperature and pressure lead to an increase in solubility y2 and thus, a decrease in supersaturation S. Although the magnitude of S decreases, it is still greater than 1, so the solution is still saturated, and particle precipitation and growth continue after the Mach disk.
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model SUPFLU-3088; No. of Pages 16
ARTICLE IN PRESS J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
10
Fig. 5. Profiles of (a) solubility y2 , (b) supersaturation S on the whole domain, (c) supersaturation along the centerline at t = 40, 000. All graphs are in log scale.
3.3. Nucleation rate and critical nuclei size Classical nucleation theory shows that a large S means a low energy barrier and a small critical radius rc , and thus, might lead to a large nucleation rate J. Fig. 6 shows the profiles of nucleation rate J and the critical nuclei radius rc . To clearly see the profiles within the expansion chamber, values of J and rc near the inlet of the straight nozzle have been cut off in the figure. Fig. 6(a) shows that the nucleation rate J starts to increase very rapidly inside the straight nozzle. The region with maximum nucleation rate is near the nozzle exit. In the expansion jet, J decreases very rapidly until it reaches the Mach disk, and then increases after the Mach disk. The phenomenon that the largest nucleation rate J happens near the nozzle exit rather than near the Mach disk is because the nucleation rate J is not only related to the supersaturation S, but is also proportional to the square of the number concentration of condensable molecules 2 in the solution, where
2 = /Mw y2 NA . Although S is several orders of magnitude higher near the Mach disk than near the nozzle exit, 2 is too small near
the Mach disk to have a large nucleation rate J. After the Mach disk, 2 increases, so the nucleation rate J increases, but it is still several orders of magnitude less than near the nozzle exit. So the main region of particle precipitation should be near the nozzle exit. We also notice that, nucleation rate J in the shear layer cannot be ignored. The variation of rc in Fig. 6(b) and (c) shows that, the critical nucleus size starts to decrease inside the straight nozzle, and continues to decrease during the expansion process. At around z = 270, the critical nucleus radius decreases down to the molecular radius of benzoic acid, rmolecule = 0.33678 nm, meaning that a single solute molecule can serve as a nuclei for particle growth, e.g., homogeneous nucleation. 3.4. Particle properties Fig. 7 shows profiles of the number density M0 , average particle size ra , and geometric standard deviation g along the centerline (a, c, e) as well as along the plate (b, d, f). Fig. 7(a) shows
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model SUPFLU-3088; No. of Pages 16
ARTICLE IN PRESS J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
11
Fig. 6. Profiles of (a) nucleation rate J, (b) critical radius rc on the whole domain, (c) rc along the centerline at t = 40,000. The inset in (c) is a close-up near the Mach disk region.
that the number density starts to increase inside the straight nozzle, it reaches a maximum near the nozzle exit, it then decreases up to the Mach disk at around z = 850 and increases sharply across the Mach disk. After z = 1140, M0 keeps almost a constant. Fig. 7(c) tells us that the average particle size ra before the Mach disk is very small, just a little bit larger than a single solute molecule, which indicates that condensation and coagulation is negligible before the Mach disk. After the Mach disk, ra increases immediately, and especially in the region z > 1000, meaning that condensation and coagulation is important for particle growth after the Mach disk. Fig. 7(e) shows that g within the region where 300 < z < 850 is almost a constant, since ra in this region is also a constant, meaning that the particle size distribution inside this region is very similar. Starting from z = 850, g becomes larger, so particle size diversity increases, especially in the interval where 1100 < z < 1300. Fig. 7(b), (d) and (f) show that the plate can be roughly divided into two parts by the impact point of the shear layer on the plate. The region in the jet (roughly 0 < r < 800) has a relatively higher
number density M0 , lower average particle size ra , and higher geometric standard deviation g . And because of the oscillating thermal properties we saw previously, both ra and g show oscillating behavior in this region. The region outside the jet (r > 800) has a lower number density, larger average particle size and smaller geometric standard deviation. And it seems that the further away from the center, the larger the particle size and smaller the geometric standard deviation. 3.5. The effects of preexpansion temperature Fig. 8(a) and (b) shows profiles of the mean particle size ra and the geometric standard deviation g along the plate with three different preexpansion temperature, T = 318 K, T = 388 K and T = 418 K, and prexpansion pressure p = 10 MPa. It is seen that the profiles of ra and g on the plate are very similar with the results for preexpansion temperature T = 388 K and T = 418 K. The particle radii ra for T = 388 K are between 30 nm and 40 nm, and ra with T = 418 K are between 30 nm and 60 nm, a little bit larger than ra with T = 388 K,
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
12
log10(ρ M0(#/m3))
log10(ρ M0(#/m3))
19.45 19.4
20
3
log10(ρ M0(#/m ))
3 log (ρ M0(#/m )) 10
22
z=850
z=1140
18
16
19.35 19.3 19.25
14
500
1000 z (μ m)
1500
0
500
1000 1500 r (μ m)
(a) ra (nm) along the center line
2500
2000
2500
45
0.4
molecular
ra (nm)
ra (nm)
2000
ra (nm) along the plate
50
30 0.6
radius
500 1000 1500 10
0
2500
(b)
40
20
2000
40 35 30
500
1000 z (μ m)
25 0
1500
500
1000 1500 r (μ m)
(c)
(d)
σ along the center line
σ along the plate
g
g
7
4 3.95
6
g
σ
σ
g
3.9 5
3.85 4
3
3.8
500
1000 z (μ m)
1500
(e)
3.75 0
500
1000 1500 r (μ m)
(f)
Fig. 7. Profiles of the number density M0 , average particle size ra , geometric standard deviation g along the centerline (a), (c), (e) and along the plate (b), (d), (f) at t = 40,000.
while ra and g with preexpansion temperature T = 318 K are much smaller than the other two cases. Moreover, the curve of ra and g are nearly flat with T = 318 K, with ra almost equal to 10 nm everywhere. The reason for the relatively low value of ra and g , as well as of the nearly flat shape of the curve, might be found in Fig. 9. It shows
profiles of (a) the temperature T and (b) the mean particle size ra on the whole domain with preexpansion temperature T = 318 K at t = 40, 000. Fig. 9(a) tells us that, with T = 318 K, the jet expands more both in the radial and axial directions. Furthermore, a very uniform thermal condition is set up in the region that is bounded by the Mach disk, the jet and the plate. The uniform physical condition and
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
13
ra (nm) along the plate
60
r (nm)
40
a
20
0
−20 0
10MPa, 318K 10MPa, 388K 10MPa, 418K 500
1000 r (μ m) (a)
1500
2000
σg along the plate
4
σ
g
3
2
1
0 0
10MPa, 318K 10MPa, 388K 10MPa, 418K 500
1000 r (μ m) (b)
1500
2000
Fig. 8. Profiles of (a) the mean particle size ra , (b) the geometric standard deviation g along the plate at t = 40, 000 with three different preexpansion temperature.
relatively low temperature and density would favor the formation of smaller particles with a narrower size distribution. That is why the mean particle size in Fig. 9(b) is very small in the whole region of the jet core. Those small particles do not have enough time to grow before they deposit on the plate, which leads to smaller values of ra and g in Fig. 8(a). Helfgen et al. [14] experimentally showed that for the benzoic-acid + CO2 system, increasing the preexpansion temperature would increase the particle size, while their 1D numerical simulation shows that the average sizes of particles are almost the same with three different preexpansion temperatures. Our numerical simulation here seems to match the general trend of their experimental results. 3.6. The effects of preexpansion pressure Fig. 10(a) and (b) shows profiles of the mean particle size ra and the geometric standard deviation g along the plate with three different prexpansion pressure, p = 7.5 MPa, p = 10 MPa and p = 12.5 MPa. The preexpansion temperature is kept at T = 388 K. There is no monotonic relationship between the mean particle size ra and the preexpansion pressure. At p = 7.5 MPa, ra increases from 60 nm at the centerline to 100 nm at the distance r = 2000 m. At p = 10 MPa, ra increases from 30 nm at the centerline to 40 nm at the distance r = 2000 m. At p = 12.5 MPa, ra decreases from 100 nm at the centerline to 70 nm at the distance r = 200 m, it then increases to 100 nm at the distance r = 2000 m. The curves of the geometric standard deviation g with p = 7.5 MPa and p = 10 MPa in Fig. 10(b) have similar shapes, especially for the region outside
Fig. 9. Profiles of (a) the temperature T, (b) the mean particle size ra on the whole domain with preexpansion condition T = 318 K, P = 10 MPa.
the jet. However, g with p = 12.5 MPa shows a higher value than the other two cases, indicating a wider size distribution of the particles. In all cases, we see that particles in the jet has a smaller average size and wider size distribution, particles outside the jet have a larger average size and narrow size distribution. Helfgen et al. [14] showed both experimentally and numerically that for the benzoic-acid + CO2 system, decreasing the preexpansion pressure would lead to larger particles, when there is air supplied to the expansion chamber. They set the extraction pressure to be the same as the preexpansion pressure. Decreasing the preexpansion pressure also led to a change of solubility. In our case, we keep the extraction pressure unchanged when the preexpansion pressure is varied. 3.7. The effects of plate position Fig. 11 shows profiles of (a) the mean particle size ra and (b) the geometric standard deviation g along the plate which is positioned in front of the nozzle at a distance of Xp = 30r, Xp = 60r, and Xp = 90r. The profiles of both ra and g with Xp = 60r and Xp = 90r are very much alike, while the profiles of both ra and g with Xp = 30r
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
14
ra (nm) along the plate
ra (nm) along the plate
100
50 40
80 r (nm)
20
a
60
a
r (nm)
30
10
40
X =30r p
0 7.5MPa, 388K 10MPa, 388K 12.5MPa, 388K
20 0
500
1000 r (μ m) (a)
1500
2000
X =60r p
−10
X =90r p
−20 0
500
2000
σg along the plate
4.2
7.5MPa, 388K 10MPa, 388K 12.5MPa, 388K
4.1
1500
(a)
σg along the plate
4.2
1000 r (μ m)
X =30r p
X =60r p
4
X =90r p
σg
σg
4 3.8
3.9 3.6 3.8 3.4 3.7 0
500
1000 r (μ m) (b)
1500
2000
Fig. 10. Profiles of (a) the mean particle size ra , (b) the geometric standard deviation g along the plate at t = 40,000 with three different preexpansion pressure.
0
500
1000 r (μ m)
1500
2000
(b) Fig. 11. Profiles of (a) the mean particle size ra and (b) the geometric standard deviation g along the plate at t = 30, 000. The plate is positioned in front of the nozzle exit at a distance Xp = 30r, Xp = 60r and Xp = 90r, respectively.
Fig. 12. Profiles of (a) the temperature T and (b) the mean particle size ra on the whole domain at t = 30, 000 when Xp = 60r.
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model SUPFLU-3088; No. of Pages 16
ARTICLE IN PRESS J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
shows larger values and has a more oscillating character. With Xp = 60r and Xp = 90r, ra increases from 10 nm near the centerline at r = 0 −−20 nm near r = 2000 m. With Xp = 30r, ra decreases from 50 nm near the centerline at r = 0 −−30 nm near r = 200 m, and then increases to 50 nm near r = 2000 m. The fact that the mean particle size deposited on a plate further away is smaller is unexpected. A larger Xp means a longer time for particles to grow, and would be expected to lead to a larger particle size. However, Fig. 12 might indicate a reason for this behavior. When the plate is positioned further away from the nozzle, the shear layer has a tendency to occupy the central region along the axis. The fluid within the shear layer has relatively lower temperature and a larger nucleation rate. Particles precipitated in the shear layer are very small and will be carried by the shear layer to the regions that the fluid is going to flow through. When the shear layer flows toward the central region, so do those small particles. That is the reason why we see a region with smaller particles in Fig. 12(b) reach as far as z = 2700. When the plate is positioned very close to the nozzle exit, the shear layer is forced to change its direction before it is going to flow toward the central region. So particles within the central region will not be influenced by the shear layer, and will continue to grow to large particles by condensation.
3.8. Discussions In order to obtain good comparison between numerical results and the experimental results for the particle characteristics, accuracy of all models included in this problem and other related issues are important. Several issues that might influence the accuracy of the numerical results are the following: • A good solubility model for the specific solute-solvent system that can describe the large range of thermal conditions in the RESS process will be able to improve the accuracy of the results. A more accurate model than the Antoine equation to estimate the vapor pressure of the solute over a large range of thermal conditions will also help to improve the results. • An accurate estimation of the surface energy of the solid-fluid surface is needed [5]. As already pointed out by Kwauk and Debenedetti [9], small changes in surface energy could cause large changes in the predicted particle size. • The classical nucleation model has been shown to be not accurate in previous studies [5,40,26], De Dea et al. [26] studied particle formation analytically using a kinetic model since they think classical nucleation model is not expected to be valid for small clusters. • Weber and Thies [17] derived three different coagulation growth laws analytically, based on numerous assumptions. A parabolic, exponential or hyperbolic growth law was obtained depending on whether coagulation is driven by Brownian motion, shear stress in laminar or turbulent flow, or the slip between the continuous phase and the dispersed phase. They showed that the accuracy of the coagulation model used to describe particle growth did have a significant impact on the final particle size. They concluded that the effects of coagulation due to slip flow must be accounted for in order to predict more accurately the particle size. • In the real experiment, the expanded jet is surrounded by a mixture of CO2 and air. The entrainment of co-flowing air might also need to be included into the model. • Another problem is the size of the simulation domain and simulation time. The distance between the plate and the nozzle exit in experiments are usually on centimeter scale, much longer than what we have used in the simulation. With our present
15
computational resources, it takes 10 days to run to dimensionless time t = 40, 000 for plate distance Xp = 30r and mesh size lc = 2.5. Taking account of the inaccuracies listed above, we expect that the numerical simulation can predict the main trends in how particle characteristics depend on operating conditions, but we should not expect an exact match with real experiments. We still consider our results to be meaningful. Most of the previous numerical studies on particle formation in the RESS process use one-dimensional models, which cannot capture the changes that happens off-axis and the influence of those off-axis properties. Our systematic 2D simulations might give some hint to explain the difference between experimental results and numerical studies. For example, we find that particle size off-axis is larger than that near the central region. This might explain why the experimentally observed particles are larger than the numerical simulations. We also find that the shear layer of the jet has an important effect on the particle characteristics on the plate.
4. Conclusion In this paper, we investigate particle formation in the RESS process, using 2D, axisymmetric Navier-Stokes equations for the fluid motion, and the general dynamic equation for the particles. The 2D simulation can capture changes off-axis that cannot be described by 1D simulation. The maximal nucleation rate and number density are found to be near the nozzle exit, and particle precipitation inside the nozzle might not be ignored. The secondary major region for particle precipitation is in the shear layer of the jet, which has a relatively low temperature and high velocity. Particles grow mainly outside the shocks where the temperature and density are relatively higher. When the shear layer of the jet impacts on the plate, it will change the local thermal conditions there. And due to the sudden change of the flow direction during the impact, the flow becomes complicated and thermal variables show some oscillatory behavior. Also, on the plate, particles within the jet have smaller average sizes, higher geometric mean and particle properties fluctuate more, while particles outside the jet shows a larger average size, lower geometric mean and less fluctuations. When the preexpansion temperature is low, the jet expands more in the both radial and axial directions, and it seems easier for the fluid to set up a uniform thermodynamic condition in the region that is bounded by the Mach disk, the jet and the plate. This favors the formation of particles with smaller size and narrower distribution. Both the average particle size and standard deviation show a monotonic increase with the preexpansion temperature, but not with the preexpansion pressure. The average particle size is tens of nanometers in all cases. Increasing the distance between the plate and the nozzle exit might decrease the particle size, since the shear layer of the jet has a tendency to flow toward the central region and bring smaller particles to the region that it flows through.
Acknowledgements We are grateful to Prof. Klas Hjort, Prof. Charlotta Turner and Researcher Irene Rodriguez Meizoso for inspiring and interesting discussions. Computer time was provided by the Swedish National Infrastructure for Computing SNIC. Support by the Swedish Research Council VR (2011-5037) is gratefully acknowledged.
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033
G Model
ARTICLE IN PRESS
SUPFLU-3088; No. of Pages 16
J. Liu et al. / J. of Supercritical Fluids xxx (2014) xxx–xxx
16
Appendix A. A.1. Geometric mean and geometric standard deviation The geometric mean of a data set {1 , 2 , . . ., n } is given by: √ g = n 1 2 . . .n , (72) and the geometric standard deviation of the data set {1 , 2 , . . ., n } is
n 2 (ln(i )/(g )) i=1
g = exp
n
.
(73)
A.2. Log-normal distribution In probability theory, if random variable Y is lognormal distributed, it means its logarithm X = log(Y) has a normal distribution. Or, if X is normal distributed, then Y = exp(X) has a log-normal distribution. The probability density function of the random variable X who has a normal distribution is 2
(x−) 1 − e 2 2 , fX (x; , ) = √ 2
(74)
The probability density function of the random variable Y who has a log-normal distribution is 2
(ln y−) 1 − 2 2 e fY (y; , ) = √ , y 2
y > 0,
(75)
the parameters and are expectation (arithmetic mean) and standard deviation of X. Since the parameter relates to the geometric mean g of Y through g = e and the parameter relates to the geometric standard deviation g of Y through g = e . So the probability density function of the random variable Y can be rewritten with g and g as − 1 fY (y; g , g ) = √ e y 2 ln g
(ln y−ln g )2 2ln2 g
,
y>0
(76)
And the kth moment of Y is given by
∞
1 22
yk f (y) dy = ek+ 2 k
E[Y k ] =
1 2 ln2 g
= g k e 2 k
.
(77)
−∞
References [1] J.W. Tom, P.G. Debenedetti, Particle formation with supercritical fluids – a review, Journal of Aerosol Science 22 (1991) 555–584. [2] Y. Chernyak, F. Henon, R.B. Harris, R.D. Gould, R.K. Franklin, J.R. Edwards, J.M. DeSimone, R.G. Carbonell, Formation of perfluoropolyether coatings by the rapid expansion of supercritical solutions (RESS) process: Part 1. Experimental results, Industrial and Engineering Chemistry Research 40 (26) (2001) 6118–6126. [3] N.A. Birkin, S.M. Howdle, U. Geddea, L. Wågberg, C. Tuner, L. Ovaskainen, I. Rodriguez-Meizoso, Towards superhydrophobic coatings made by nonfluorinated polymers sprayed from a supercritical solution, Journal of Supercritical Fluids 77 (2013) 134–141. [4] O. Werner, C. Turner, Investigation of different particle sizes on superhydrophobic surfaces made by rapid expansion of supercritical solution with in situ laser diffraction (RESS-LD), Journal of Supercritical Fluids 67 (2012) 53–59. [5] M. Türk, Influence of thermodynamic behaviour and solute properties on homogeneous nucleation in supercritical solutions, Journal of Aerosol Science 32 (3) (2001) 295–319. [6] J.R. Dean, Applications of Supercritical Fluids in Industrial Analysis, Springer/CRC press Inc., USA/Canda, 1993. [7] M. Jia, Particle formation in supercritical carbon dioxide by the RESS process (M.Sc. thesis), The University of Western Ontario, London, Ontario, 2003. [8] A.B. Moussa, H. Ksibi, A review of numerical investigations regarding the supercritical fluid expansion in the RESS process, International Journal of Emerging Multidisciplinary Fluid Sciences 2 (1) (2010).
[9] P.G. Debenedetti, X. Kwauk, Mathematical modeling of aerosol formation by rapid expansion of supercritical solutions in a converging nozzle, Journal of Aerosol Science 24 (1993) 445–469. [10] X. Kwauk, P.G. Debenedettl, J.W. Tom, S.-D. Yeo, Rapid expansion of supercritical solutions (RESS): fundamentals and applications, Fluid Phase Equilibria 82 (1993) 311–321. [11] P.G. Debenedetti, M. Weber, L.M. Russell, Mathematical modeling of nucleation and growth of particles formed by the rapid expansion of a supercritical solution under subsonic conditions, Journal of Supercritical Fluids 23 (2002) 65–80. [12] M. Türk, Formation of small organic particles by RESS: experimental and theoretical investigations, Journal of Supercritical Fluids 15 (1) (1999) 79–89. [13] B. Helfgen, P. Hils, C. Holzknecht, M. Türk, K. Schaber, Simulation of particle formation during the rapid expansion of supercritical solutions, Journal of Aerosol Science 32 (3) (2001) 295–319. [14] B. Helfgen, M. Türk, K. Schaber, Hydrodynamic and aerosol modelling of the rapid expansion of supercritical solutions (RESS-process), Journal of Supercritical Fluids 26 (3) (2003). [15] A.D. Shine, A.K. Lele, Morphology of polymers precipitated from a supercritical solvent, American Institute of Chemical Engineers Journal 38 (1992) 742. [16] P. Pallado, E. Reverchon, Hydrodynamic modeling of the RESS process, Journal of Supercritical Fluids 8 (4) (1995) 318–328. [17] M.C. Thies, M. Weber, A simplified and generalized model for the rapid expansion of supercritical solutions, Journal of Supercritical Fluids 40 (3) (2007) 402–419. [18] H. Ksibi, C. Tenaud, P. Subra, Y. Garrabos, Numerical simulation of the rapid expansion of supercritical fluid flow, European Journal of Mechanics – B/Fluids 15 (4) (1996) 569–596. [19] O.G. Gadetskii, V.V. Altunin, Method of formulating the fundamental equations of state for pure substances from varied experimental data, Teplofizika Vysokikh Temperatur 9 (3) (1971) 527–534. [20] P.L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics 43 (1981) 357–372. [21] A.B. Moussa, H. Ksibi, C. Tenaud, M. Baccar, Parametric study on the nozzle geometry to control the supercritical fluid expansion, International Journal of Thermal Sciences 44 (2005) 774–786. [22] H. Ksibi1, A. Ben Moussa1, M. Baccar, Simulation of particles transport and coagulation during the RESS process, European Physical Journal – Applied Physics 43 (2008) 253–261. [23] R.K. Franklin, J.R. Edwards, Y. Chernyak, R.D. Gould, F. Henon, R.G. Carbonell, Formation of perfluoropolyether coatings by the rapid expansion of supercritical solutions (RESS) process: Part 2. Numerical modeling, Industrial and Engineering Chemistry Research 40 (26) (2001) 6127–6139. [24] I. Khalil, D.R. Miller, The structure of supercritical fluid free-jet expansions, American Institute of Chemical Engineers Journal 50 (11) (2004) 2697–2704. [25] S.I. Sandler, Chemical, Biochemical, and Engineering Thermodynamics, Wiley, USA, 2006. [26] S. De Dea, D. Graziani, D.R. Miller, R.E. Continetti, Growth of magnetic thin films using CO2 RESS expansions, Journal of Supercritical Fluids 42 (2007) 410–418. [27] A. Onuki, Dynamic van der Waals theory, Physical Review E 75 (2007). [28] B. Platzer, G. Maurer, A generalized equation of state for polar and nonpolar fluids, Fluid Phase Equilibria 51 (1989) 223–236. [29] S.E. Pratsinis, Simultaneous nucleation, condensation, and coagulation in aerosol reactors, Journal of Colloid and Interface Science 124 (2) (1988) 416–427. [30] J.D. van der Waals, Koninklijke Nederlandse Akademie van Wetenschappen, Verhandelingen, Afdeling Natuurkunde. Eerste Reeks 1 (8) (1893) 56. [31] J.S. Rowlinson, Translation of J.D. van der Waals, the thermodynamic theory of capillary under the hypothesis of a continuous variation of density, Journal of Statistical Physics 20 (1979) 197–244. [32] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system: III. Nucleation in a two-component incompressible fluid, Journal of Chemical Physics 31 (1959) 688. [33] G. Amberg, J. Liu, M. Do-Quang, Numerical simulation of rapid expansion of supercritical carbon dioxide, American Institute of Chemical Engineers Journal (2014), http://dx.doi.org/10.1002/aic.14603. [34] S.K. Friedlander, Smoke, Dust and Haze: Fundamentals of Aerosol Behavior, 2nd ed., Oxford University Press, New York, 2000. [35] R.W. Miles, P.J. Erbland, D.P. Rizzetta, Numerical and experimental investigation of CO2 condensate behavior in hypersonic flow, in: AIAA-2000-2379, 21st AIAA Aerodynamic Measurement Technology and Ground Testing Cofnerence, 2000. [36] J.E. McDonald, Homogeneous nucleation of vapor condensation: I. Thermodynamic aspects, American Journal of Physics (1962). [37] R.C. Reid, W.J. Schmitt, Solubility of monofunctional organic solids in chemically diverse supercritical fluids, Journal of Chemical Engineering Data 31 (1986) 204–212. [38] N.A. Webb, J.C. Barrett, A comparison of some approximate methods for solving the aerosol general dynamic equation, Journal of Aerosol Science 29 (1998) 31–39. [39] G. Amberg, R. Tönhardt, C. Winkler, Finite element simulations using symbolic computing, Mathematics and Computers in Simulation 49 (1999) 257–274. [40] B.J. Jurcik, J.R. Brock, Numerical simulation of particle formation and growth in rapidly expanding axisymmetric flows, Journal of Chemical Physics 97 (1993) 323–331.
Please cite this article in press as: J. Liu, et al., Numerical simulation of particle formation in the rapid expansion of supercritical solution process, J. Supercrit. Fluids (2014), http://dx.doi.org/10.1016/j.supflu.2014.08.033