Fluid Phase Equilibria 219 (2004) 33–36
Molecular dynamics studies on clustering process of solute molecules through rapid expansion of supercritical fluids Shin-ichi Furukawa, Satoshi Kato, Tomoshige Nitta∗ Department of Materials Engineering Science, Graduate School of Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan Received 1 September 2003
Abstract The expansion of a supercritical fluid (SCF) solution composed of CO2 and naphthalene from 10 to 1 MPa has been simulated by using both the conventional NVT-MD and one-dimensional pressure-control molecular dynamics (OPC-MD) methods. The characteristic pressure corresponding to a significant decrease in density of the SCF solution (about 6 MPa) is lower than that of pure CO2 (about 8 MPa) due to the strong attractive interactions of NAP molecules with CO2 molecules. The solute clusters dispersed in the SCF start to grow at the characteristic pressure by uniting small clusters to form a larger cluster and gathering dispersed solute molecules to the cluster. In the solute clusters, the number of CO2 molecules gradually decrease as the clusters grow, though the number is still comparable to that of solute molecules at 1 MPa, suggesting that remaining CO2 molecules will escape from the solute clusters at lower pressures. © 2004 Elsevier B.V. All rights reserved. Keywords: Simulation; Supercritical fluid; RESS; Carbon dioxide; Naphthalene
1. Introduction Supercritical fluids (SCFs), the states of which are above their critical points in both pressure and temperature, have received much attention because of their interesting characteristics which are different to gas or liquid; high diffusivity, high solubility and high reactivity. Therefore, studies on engineering applications of SCFs have been intensively carried out as well as on fundamental properties. The rapid expansion of supercritical solution (RESS) process [1,2] is one of such applications for production of small particles with narrow size distribution. The process involves rapid expansion of an SCF solution as it flows through a narrow nozzle toward atmospheric pressure and the rapid expansion brings the SCF solution to be extremely supersaturated, which results in rapid nucleation and growth of the solute particles of small average size. Typical solutes that have been reported are polyaromatic compounds such as naphthalene, phenanthrene and anthracene, while the most popular SCF is CO2 . The particle formation mechanism in the RESS process, however, has not been well understood at the molecu∗ Corresponding author. Tel.: +81-6-6850-6265; fax: +81-6-6850-6265. E-mail address:
[email protected] (T. Nitta).
0378-3812/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2004.01.008
lar level, though some theoretical approaches including the classical nucleation theory, transport equations, and equations of state [3–7] have been utilized to investigate growing behavior of particles. Molecular simulation techniques, consisting of Monte Carlo (MC) and molecular dynamics (MD) methods [8,9], have been recognized to be very powerful in obtaining microscopic and macroscopic information on static and/or dynamic behavior of molecular systems. Some recent simulation studies on SCFs are as follows: solubility [10–12], diffusivity [13–17] and statistical behavior of clustering in SCF solutions [18,19]. However, to our best knowledge, there have been no papers reporting on the dynamical behavior of clustering formation of solute molecules in the RESS process. In the present work, we presume that the solute particles aggregate together and form clusters not outside but inside a nozzle through which an SCF solution flows down with decreasing pressure and density towards the exit. We have developed a one-dimensional pressure-control MD (OPC-MD) method to quickly obtain snapshots of dynamical behavior of clustering with the lowering of pressure. The objective of the present work is, therefore, to provide molecular level information on dynamical aspects of cluster formation of solute particles dissolved in an SCF
34
S.-i. Furukawa et al. / Fluid Phase Equilibria 219 (2004) 33–36
that expands within a narrow nozzle in the process of RESS.
2. Simulation details The simulation cell for an SCF mixture is a rectangular box composed of an expandable side in x-direction and the fixed square in y- and z-directions. The original box lengths are as follows: Lx = 3.42 nm, Ly = Lz = 13.7 nm. We use the conventional periodic boundaries for the simulation cell of the fluid element since the box lengths are small enough compared with the nozzle diameters used in experiments, though they are as small as 0.1–1.0 mm. The simulation method we employed is composed of two types of MD steps: the conventional NVT-MD [8,9] and OPC-MD. Fig. 1 shows a schematic diagram of the simulation cell for the OPC-MD method, in which only Lx is adjusted in such a way that the calculated pressure of x-direction, Px , coincides with a specified pressure P by using the Newton–Raphson numerical method [20]. The specification of pressure P, which is lowering consecutively, is conducted at every 10 MD steps, while Lx is kept constant in other MD steps. The 10 MD step interval was chosen from among five trials of at the 1, 5, 10, 20 and 50 MD steps so as to provide low calculation cost and a minimum of artificial pressure disturbances. The temperature of a fluid in the simulation cell is kept constant by using the velocity scaling technique [8], which corresponds to the assumption that heat flows from nozzle walls during the expansion of fluid. The CO2 and naphthalene (NAP) are chosen as an SCF and a solute, respectively, mainly because the solubility of NAP in CO2 was investigated in both experiments [21–23] and molecular simulations [10–12]. Each molecule was expressed as the spherical Lennard–Jones (LJ) particle with parameters reported by Eya et al. [10]. The LJ cross-parameters are evaluated by using the Lorentz–Berthelot combining rule. In all simulations, the pressure control and the temperature correction are carried out at every 10 MD steps. One MD step time is taken as 10 fs. The temperature is set at 308.2 K, while the pressure decreases from 10 to 1 MPa for 4,500,000 MD steps. The mole fraction of NAP in an SCF solution is set at 0.0523 (276 NAP and 5000 CO2 molecules), which is the experimental solubility of NAP in CO2 at 28 MPa and 308.2 K taken from a paper of Eya et al. [10]. In the preliminary simulations of pressure lowering,
Fig. 1. Schematic diagram of simulation cell for the OPC-MD method.
we found that the densities of the SCF solution between 28 and 10 MPa changed only slightly. Therefore, the starting pressure was chosen as 10 MPa so as to focus on the region of SCF expansion. Before the start of simulations, the initial values for velocities and locations of each molecule are taken from those of the last outcomes of an NPT-MD equilibrium simulation conducted at 308.2 K and 10 MPa for 500,000 MD steps.
3. Results and discussion Fig. 2 shows pressure effects on fluid density of a CO2 –NAP mixture at 308.2 K obtained in this work and that of pure CO2 reported by Eya et al. [10]. The mixture density decreases gradually with decreasing pressure up to 6 MPa, and then it rapidly decreases with decreasing pressure. As compared with the curve of pure CO2 , the pressure dependence of the mixture curve is almost the same, however, the characteristic pressure corresponding to the starting point of rapid decrease in density is appreciably different, i.e., 6 MPa for the mixture while 8 MPa for pure CO2 , which is ascribed to the large attractive potential of NAP molecules that keep the CO2 molecules in high density even at a low pressure of the mixture. In Fig. 2, the box length of x-direction (Lx ) is also displayed against lowering pressure for the CO2 –NAP mixture. The gradient of Lx is proportional to an expansion speed of the box length since a speed of pressure lowering is constant at 2.0 × 108 MPa/s. The expansion speed rapidly increases up to 5.0 m/s at 5.8 MPa, and quickly decreases to 0.8 m/s at 5.7 MPa. From then, the expansion speed increases with decreasing pressure up to 20.0 m/s at 1 MPa. Since the expansion speed is much smaller than average molecular velocities (CO2 : 241 m/s and NAP: 141 m/s at 308.2 K), we can conclude that the speed of pressure lowering we adopted might be appropriate to simulate the growth of clusters.
Fig. 2. Pressure effects on fluid density of a CO2 –NAP mixture, that of pure CO2 , and box length in the x-direction, fluid density of pure CO2 is taken from Eya et al. [10].
S.-i. Furukawa et al. / Fluid Phase Equilibria 219 (2004) 33–36
35
Fig. 3. Local density profiles of component B around a NAP molecule, B: CO2 (a) and NAP (b) and T = 308.2 K, P = 10, 5.5, 1.0 MPa.
Fig. 3a and b show the local density profiles of CO2 and NAP, respectively, around a NAP molecule at pressures of 10, 5.5 and 1.0 MPa and temperature of 308.2 K. At a first glance, we can see that the pressure effects on the local densities of CO2 and NAP are inverse, i.e., the local density of CO2 decreases with a decrease in pressure while that of NAP increases and the increase is significant at 1.0 MPa, which indicates that solute NAP molecules start to gather and form their clusters below 6 MPa. Another note obtained from the two figures is that the local density of CO2 around a NAP molecule at 1.0 MPa is about a half of the density at 10 MPa in the first nearest shell and about a third at distances longer than the first shell, which indicates that many CO2 molecules still exist around each NAP molecule even at 1.0 MPa. In order to pay attention to the dynamical behavior of the NAP cluster formation in more detail, we apply Hill’s cluster definition [24] which judges whether a pair of solute molecules is bound or not, by comparing their kinetic and potential energies, two molecules are defined to be a bound
Fig. 5. Pressure effect on (a) the number of clusters and (b) the number of molecules in clusters at 308.2 K.
pair when the relative kinetic energy (Ekr ) is less than or equal to the potential energy (Ep ) for the two molecules. The Ekr is calculated as Ekr =
1 m i mj (vi − vj )2 2 mi + m j
(1)
where m and v are the molecular mass and the velocity, respectively. The solute cluster is first defined as an assemble of solute (NAP) molecules that are bound directly to other NAP molecules, In addition, solvent molecules (CO2 ) bound with any solute molecules belonging the solute cluster are also included as a member of the cluster. Fig. 4 shows molecular configurations in the simulation cell before and after applying the definition of the solute cluster, from which we can see that molecules of cluster members have been selected properly by the definition of a solute cluster. Fig. 5
Fig. 4. Display of (a) all molecules in the cell and (b) NAP cluster molecules defined by Hill’s cluster criterion at T = 308.2 K, P = 5.5 MPa, full circles stand for naphthalene molecules and open ones for CO2 molecules.
36
S.-i. Furukawa et al. / Fluid Phase Equilibria 219 (2004) 33–36
Fig. 6. Pressure effect on the remaining ratio of NAP molecules in clusters for 10 ps at 308.2 K.
shows pressure effects on the number of clusters and on the number of molecules in clusters. In Fig. 5a, the number of solute clusters is almost constant until the pressure reaches at 6 MPa and then it decreases with decreasing pressure and finally becomes one. In Fig. 5b, at pressures above 6 MPa, the number of CO2 molecules in solute clusters is always much larger than the number of NAP molecules and both numbers are almost constant. However, at pressures lower than 6 MPa, the number of CO2 molecules in the solute clusters decreases with decreasing pressure, while that of NAP molecules increases. Based on these results, we can have a picture for solute clusters to grow, at the characteristic pressure (6 MPa) the solute clusters start to grow by uniting small clusters to form a larger cluster and also gathering solute molecules to the cluster since there have been dispersed many small solute clusters at pressures above 6 MPa. The CO2 molecules in a solute cluster gradually escape from the cluster, but the solute cluster at 1.0 MPa still contained more CO2 molecules than NAP molecules. In order to estimate the molecular mobility in solute clusters, we calculate the remaining ratio of NAP molecules, which is defined as the ratio of the number of NAP molecules that keep remaining in clusters for 10 ps to the total number of NAP molecules in the clusters. Fig. 6 shows the remaining ratio of NAP against pressure. At pressures above 6 MPa, the remaining ratio is about 50%, indicating that almost half NAP molecules have left the cluster and new comers enter in 10 ps, while the remaining ratio sharply increases up to 80% at 6 MPa and then gradually increases with decreasing pressure. 4. Conclusions We have conducted molecular simulations for the expansion of an SCF solution composed of CO2 (solvent) and naphthalene (solute) that decreases in pressure through a
nozzle by using both conventional NVT-MD and (OCP-MD) methods. The characteristic pressure corresponding to the significant decrease in the solution density (about 6 MPa) was found to be lower than that for pure CO2 (about 8 MPa) due to the strong attractive interaction between NAP and CO2 molecules. As a dynamical aspect of clustering, it was found that the NAP solute clusters started to grow at about 6 MPa by uniting small clusters to form a larger cluster and also gathering dispersed solute molecules to the cluster. The solute clusters contain much more CO2 molecules than NAP molecules and the number of CO2 molecules gradually decreases as the clusters grow though the number is still comparable to that of NAP molecules at 1 MPa, which suggests that remaining CO2 molecules will escape from the solute clusters at lower pressures.
References [1] D.W. Matson, J.L. Fulton, R.C. Petersen, R.D. Smith, Ind. Eng. Chem. Res. 26 (1987) 2298–2306. [2] J.W. Tom, P.G. Debenedetti, J. Aerosol Sci. 22 (1991) 555–584. [3] P.G. Debenedetti, J.W. Tom, X. Kwauk, S.D. Yeo, Fluid Phase Equilib. 82 (1993) 311–321. [4] G.R. Shaub, J.F. Brennecke, M.J. McCready, J. Supercrit. Fluids 8 (1995) 318–328. [5] M. Turk, J. Supercrit. Fluids 15 (1999) 79–89. [6] B. Helfgen, P. Hils, C.H. Holzknecht, M. Turk, K. Schaber, J. Aerosol Sci. 31 (2001) 295–319. [7] M. Weber, L.M. Russell, P.G. Debenedetti, J. Supercrit. Fluids 23 (2002) 65–80. [8] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, UK, 1987. [9] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithm to Applications, Academic Press, USA, 1996 (revised ed., 2001). [10] H. Eya, Y. Iwai, T. Fukuda, Y. Arai, Fluid Phase Equilib. 77 (1992) 39–51. [11] Y. Iwai, H. Uchida, Y. Arai, Y. Mori, Fluid Phase Equilib. 144 (1998) 23–244. [12] Y. Iwai, Y. Mori, Y. Arai, Fluid Phase Equilib. 167 (2000) 33–40. [13] H. Higashi, Y. Iwai, H. Uchida, Y. Arai, J. Supercrit. Fluids 13 (1998) 93–97. [14] T. Yamaguchi, Y. Kimura, N. Hirota, Mol. Phys. 94 (1998) 527– 537. [15] J. Zhou, X. Lu, Y. Wang, J. Shi, Fluid Phase Equilib. 172 (2000) 279–291. [16] Y. Zhu, X. Lu, J. Zhou, Y. Wang, J. Shi, Fluid Phase Equilib. 194–197 (2002) 1141–1159. [17] H.D. Ludenmann, L. Chen, J. Phys. Condens. Matter 14 (2002) 11453–11462. [18] I.B. Petsche, P.G. Debenedetti, J. Chem. Phys. 91 (1989) 7075–7084. [19] C.C. Liew, H. Inomata, S. Saito, Fluid Phase Equilib. 104 (1995) 317–329. [20] F.B. Hildebrand, Advanced Calculus for Applications, Prentice-Hall, USA, 1962. [21] Y.L. Tsekhanskaya, M.B. Iomtev, E.V. Mushkina, Russ. J. Chem. 38 (1964) 1173–1176. [22] W.J. Schmitt, R.C. Reid, J. Chem. Eng. Data 31 (1986) 204–212. [23] S. Sako, K. Ogaki, T. Katayama, J. Supercrit. Fluids 1 (1988) 1–6. [24] L.T. Hill, J. Chem. Phys. 23 (1955) 617–622.