Fluid Phase Equilibria 163 Ž1999. 165–173 www.elsevier.nlrlocaterfluid
Monte Carlo simulation for pVT relationship of CO 2 q n-C 4 H 10 and CO 2 q i-C 4 H 10 systems Morio Yamamoto, Yoshio Iwai ) , Yasuhiko Arai Department of Chemical Systems and Engineering, Graduate School of Engineering, Kyushu UniÕersity, 6-10-1 Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan Received 28 August 1998; accepted 10 May 1999
Abstract Monte Carlo simulation has been applied to calculate the pVT relationship of CO 2 q butane Ž n-C 4 H 10 and i-C 4 H 10 . systems at 310.93 K and up to 9.5 MPa. CO 2 is treated as a single-site molecule and butanes are treated as four-site molecules. The Lennard–Jones Ž12–6. potential is used as the site–site potentials and the combining rules proposed by Jorgensen et al. wW.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. 106 Ž1984. 6638.x are adopted for unlike site pairs. The calculated results of the pVT relationship show good agreement with the experimental data wT. Tsuji, S. Honda, T. Hiaki, M. Hongo, J. Supercrit. Fluids 13 Ž1998. 15.x by introducing an intersite interaction parameter between unlike molecules. Furthermore, the radial distribution functions and the number of CO 2 and butanes around butanes are calculated as a fundamental information on the microscopic structures. It is found that the radial distribution functions and the number of CO 2 and n-C 4 H 10 around n-C 4 H 10 are different from those of CO 2 and i-C 4 H 10 around i-C 4 H 10 . q 1999 Elsevier Science B.V. All rights reserved. Keywords: Theory; Monte Carlo simulation; pVT relationship; Carbon dioxide; Butane; Radial distribution function
1. Introduction Supercritical fluids have been given much attention recently as one of the new type of solvents for the separation and the reaction processes. To design such processes, a reliable knowledge of the pVT relationship of binary systems containing supercritical fluids is required. Many researchers have measured the pVT relationship of fluids including their supercritical region and correlated it with equations of state. Tsuji et al. w2x measured the pVT relationship of CO 2 q butane Ž n-C 4 H 10 and )
Corresponding author. Tel.: q81-92-642-3496; fax: q81-92-642-3496; e-mail:
[email protected]
0378-3812r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 Ž 9 9 . 0 0 2 2 8 - 9
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
166
i-C 4 H 10 . systems and correlated it with the Benedict–Webb–Rubin equation of state. These results show that the behavior for the pVT relationship of CO 2 q n-C 4 H 10 system is slightly different from that for CO 2 q i-C 4 H 10 system. It seems interesting to study the difference. Recently, molecular simulation has been applied widely to calculate physical properties owing to the rapid development of high-speed computers. It becomes possible to calculate the pVT relationship and to investigate the phase equilibrium by computer simulation using a personal computer or workstation. Molecular simulation may be helpful in realizing not only the pVT relationship but also the microscopic structures of solvation. In previous works w3–5x, the authors studied chain molecules dissolved in supercritical fluids using Monte Carlo simulation. The static properties of n-paraffins, fatty acids and higher alcohols in supercritical CO 2 and ethane were calculated and discussed. The authors also introduced a group contribution site model to calculate the solubilities of aromatic isomers in supercritical CO 2 w6x. In the model, CO 2 was treated as a single-site molecule, and chain molecules and aromatic compounds were treated as multisite molecules to distinguish isomers. The treatment is the same as common group contribution method such as UNIFAC w7x. In this work, Monte Carlo simulation was applied to calculate the pVT relationship of CO 2 q butane Ž n-C 4 H 10 and i-C 4 H 10 . systems at 310.93 K and up to 9.5 MPa. CO 2 was treated as single-site molecule for simplification. Butanes were treated as four-site molecules to distinguish isomers and the solvation structure around those components was investigated. The Lennard–Jones Ž 12–6. potential was used as the site–site potentials and the combining rules proposed by Jorgensen et al. w1x were adopted for unlike site pairs. Intersite interaction parameters between CO 2 and alkyl group sites in butanes were introduced to fit the calculated results to the experimental pVT data w2x. Furthermore, the radial distribution functions and the number of CO 2 and butanes around butanes are reported as a fundamental information on the microscopic structures.
2. Potential functions In this work, butanes were assumed to be composed of four alkyl group sites with fixed C–C bond length of 0.153 nm w1x. The valence angles between successive pairs Q Ž C–C–C. were kept to be 1128 w1x as shown in Fig. 1. The intermolecular potential energy between molecules a and b, fab , can be calculated by the Lennard–Jones potential between each pair site: na
nb
fab s Ý Ý 4´ i j i
j
si j
12
si j
6
½ž / ž / 5 ri j
y
ri j
Ž1.
where n a and n b are the number of sites in molecules a and b, ri j is the distance between i and j sites, and ´ and s are the energy and size parameters, respectively. The potential parameters for unlike site pairs are expressed by the geometric mean combining rules, according to Jorgensen et al. w1x:
(
si j s sii sj j
Ž2.
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
167
Fig. 1. The geometric structures of n-C 4 H 10 and i-C 4 H 10 .
and
(
´ i j s Ž 1 y k i j . ´ ii ´ j j
Ž3.
The intersite interaction parameter, k i j , was introduced to the combining rule for energy parameter between unlike molecules. The values between the same molecules were fixed to be 0. For Lennard–Jones fluids, Nicolas et al. w8x estimated the reduced critical constant values pUC s Ns 3rVC , TCU s kTCr´ and pUC s s 3 pCr´ to be 0.35, 1.35 and 0.1418, respectively. The parameters for CO 2 were calculated by using 304.2 K for TC and 7.37 MPa for pC w9x, which are more reliable than VC . On the other hand, for alkyl group sites of butanes, the authors adopted the values for those of liquid hydrocarbons proposed by Jorgensen et al. w1x. The potential parameters adopted are listed in Table 1. For n-C 4 H 10 molecule with an internal rotational degree of freedom, it is essential to define the intramolecular potential. According to the potential functions proposed by Jorgensen et al. w1x, the intramolecular potential was taken into account by the following equation: 1 1 1 z Ž F . s z Ž0. q z Ž1. Ž 1 q cosF . q z Ž2. Ž 1 y cos2F . q z Ž3. Ž 1 q cos3F . 2 2 2
Ž4.
where F is the rotational angle in n-C 4 H 10 , and z Ž0. to z Ž3. are the rotational potential parameters of n-C 4 H 10 listed in Table 2.
3. Calculation procedure In calculation, the authors used a system of 269 and eight particles of CO 2 and n-C 4 H 10 Ž2.89 mol%. or 267 and eight particles of CO 2 and i-C 4 H 10 Ž2.91 mol%. enclosed in a central cubic box with periodic boundary conditions. The particle numbers were decided to agree with the experimental concentration of Tsuji et al. w2x. The NVT and NpT ensembles were used in the present simulation.
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
168
Table 1 Potential parameters of CO 2 and butanes
´ r k wKx
Site CO 2 Žcarbon dioxide. CH 3 Ž n-C 4 H 10 . CH 2 Ž n-C 4 H 10 . CH 3 Ž i-C 4 H 10 . CH Ž i-C 4 H 10 . a b
s wnmx
a
0.3910 a 0.3905 b 0.3905 b 0.3910 b 0.3850 b
225.3 88.1b 59.4 b 80.5 b 40.3 b
Nicolas et al. w8x. Jorgensen et al. w1x.
The standard Metropolis method w10x was used to obtain new configurations under the NVT and NpT ensembles. The reptation moves w10x were employed to change the configuration of n-C 4 H 10 molecules. Therefore, n-C 4 H 10 molecules moved in the slithering fashion, that is, the head of n-C 4 H 10 molecule moved to a new position and the rest of n-C 4 H 10 follows like a snake with maintaining bond length, L, and bond angle, Q . Every 40 movements of CO 2 molecules, 10 movements of butane molecules were carried out. In the NpT ensemble, volume was changed after each 1000 movements of butane molecules. On the other hand, i-C 4 H 10 molecules moved and rotated randomly at the position to reduce the influence of molecule geometry. The lengths of the calculation of the pVT relationship were 1.0 = 10 6 configurations and the pVT relationship was obtained by using the virial theorem.
4. Results and discussion The calculated results for the pVT relationship by the NVT and NpT ensembles Monte Carlo simulation are shown in Fig. 2 together with the experimental data w2x. In this case, the intersite interaction parameter, k 12 , between CO 2 and alkyl group sites in butanes was fixed to be 0. As shown in Fig. 2, the calculated results by the NVT ensemble Monte Carlo simulation agree with those by the NpT ensemble. However, they are slightly different from the experimental data w2x. Therefore, the intersite interaction parameters, k 12 , between CO 2 and alkyl group sites in butanes were optimized to fit the calculated results to the experimental data at 6.9 MPa. The values of k 12 for CO 2 q n-C 4 H 10 and CO 2 q i-C 4 H 10 systems were decided to be y0.05 and 0.10, respectively. The calculated results of the pVT relationship by the NVT ensemble MC simulation with the optimized k 12 are shown in
Table 2 a Rotational potential parameters for n-C 4 H 10 Bond
z Ž0. r k wKx
z Ž1. r k wKx
z Ž2. r k wKx
z Ž3. r k wKx
C–C–C–C
0.0
765.9
y158.5
1614
a
Jorgensen et al. w1x.
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
169
Fig. 2. The pVT relationship of CO 2 q n-C 4 H 10 and CO 2 q i-C 4 H 10 systems at 310.93 K: Ž'. n-C 4 H 10 Ž2.89 mol%.; Žv . i-C 4 H 10 Ž2.91 mol%., NpT ensemble Monte Carlo with k 12 s 0; Ž^. n-C 4 H 10 Ž2.89 mol%.; Ž`. i-C 4 H 10 Ž2.91 mol%., NVT ensemble Monte Carlo with k 12 s 0; Ž- - -P. n-C 4 H 10 Ž2.89 mol%.; Ž — . i-C 4 H 10 Ž2.91 mol%., experimental ŽTsuji et al. w2x.; vertical and horizontal lines represent fluctuations in calculated results of compressibility factor and pressure, respectively.
Fig. 3 together with the experimental data w2x. The calculated results showed better agreement with the experimental data by introducing the intersite interaction parameter, k 12 , except at high pressure region.
Fig. 3. The pVT relationship of CO 2 q n-C 4 H 10 and CO 2 q i-C 4 H 10 systems at 310.93 K: Ž^. n-C 4 H 10 Ž2.89 mol%., NVT ensemble Monte Carlo with k 12 sy0.05; Ž`. i-C 4 H 10 Ž2.91 mol%., NVT ensemble Monte Carlo with k 12 s 0.10; Ž- - -P. n-C 4 H 10 Ž2.89 mol%.; Ž — . i-C 4 H 10 Ž2.91 mol%., experimental ŽTsuji et al. w2x.. Horizontal lines represent fluctuations in calculated results of pressure.
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
170
The radial distribution functions, g Ž r ., were also calculated in order to study the difference of the microscopic structures between supercritical CO 2 and butanes around butanes. The NVT ensemble was used to calculate the radial distribution functions. The radial distribution functions for sites in butanes at 310.93 K and 6.9 MPa were averaged from 1.0 = 10 7 to 11.0 = 10 7 configurations. The number of molecules inside sphere of radius r surrounding sites in butanes, nŽ r ., was determined by using the following equation: r
n Ž r . s 4pr
H0 r
2
g Ž r .d r
Ž5.
These results are shown in Figs. 4–7. As shown in Fig. 4Ž a. ,Ž b. , the radial distribution function and the number of CO 2 around CH 2 group in n-C 4 H 10 are similar to those around CH 3 group. On the other hand, the position of the first peaks for CH 3 and CH group are different and the first peak of CH 3 group is lower than that of CH group as shown in Fig. 5Ža.. From Fig. 5Žb., the number of CO 2 , nŽ r ., within a radius of 0.55 nm of CH 3 group is larger than that of CH group though nŽ r . counted from CH 3 group is smaller than that from CH group in the region of 0.55 - r - 0.80 nm. While, the number of CO 2 surrounding n-C 4 H 10 is quite similar to that surrounding i-C 4 H 10 in the region of r ) 0.80 nm. These results show that the microscopic structures of supercritical CO 2 around n-C 4 H 10 are different from those for i-C 4 H 10 . Fig. 6Ža.,Ž b. show the radial distribution functions and the number of sites in n-C 4 H 10 around sites in n-C 4 H 10 , respectively. One can see that all the curves in each figure are almost the same. Fig. 7Ža.,Žb. represent the radial distribution functions and the number of sites in i-C 4 H 10 around sites in i-C 4 H 10 , respectively. As shown in Fig. 7Ža., the first peak of CH group around CH group is the highest and that of CH 3 group around CH 3 group is the lowest. The value of nŽ r . for CH 3 group around CH 3 group within 0.60 nm is the highest as shown in Fig. 7Ž b. . This fact shows that CH 3
Fig. 4. Ža. The radial distribution functions of CO 2 around sites in n-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group; Ž- - -P. CH 2 group. Žb. The number of CO 2 within distance r from sites in n-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group; Ž- - -P. CH 2 group.
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
171
Fig. 5. Ža. The radial distribution functions of CO 2 around sites in i-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group; Ž- - -P. CH group. Žb. The number of CO 2 within distance r from sites in i-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group; Ž- - -P. CH group.
group in i-C 4 H 10 tends to crowd a little more around CH 3 group, because the value of energy parameter of CH 3 group in i-C 4 H 10 is larger than that of CH group. Furthermore, the number density of i-C 4 H 10 near i-C 4 H 10 is slightly higher than that of n-C 4 H 10 near n-C 4 H 10 , because the number of i-C 4 H 10 surrounding i-C 4 H 10 Žabout 0.25 as shown in Fig. 7Žb.. is larger than that surrounding n-C 4 H 10 Žabout 0.22 as shown in Fig. 6Žb.. in the region of r - 0.80 nm. From these results, it is
Fig. 6. Ža. The radial distribution functions of sites in n-C 4 H 10 around sites in n-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group around CH 3 group; Ž- - -P. CH 2 group around CH 3 group; Ž-P-P-P. CH 2 group around CH 2 group. Žb. The number of sites in n-C 4 H 10 within distance r from sites in n-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group around CH 3 group; Ž- - -P. CH 2 group around CH 3 group; Ž-P-P-P. CH 2 group around CH 2 group.
172
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
Fig. 7. Ža. The radial distribution functions of sites in i-C 4 H 10 around sites in i-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group around CH 3 group; Ž- - -P. CH group around CH 3 group; Ž-P-P-P. CH group around CH group. Žb. The number of sites in i-C 4 H 10 within distance r from sites in i-C 4 H 10 at 310.93 K and 6.9 MPa calculated by Monte Carlo simulation: Ž — . CH 3 group around CH 3 group; Ž- - -P. CH group around CH 3 group; Ž-P-P-P. CH group around CH group.
shown that the microscopic structures of n-C 4 H 10 around n-C 4 H 10 are different from those of i-C 4 H 10 around i-C 4 H 10 .
5. Conclusions Monte Carlo simulations have been performed to calculate the pVT relationship of CO 2 q butane Ž n-C 4 H 10 and i-C 4 H 10 . systems at 310.93 K and up to 9.5 MPa. The calculated pVT relationship shows reasonable agreement with the experimental data by introducing the intersite interaction parameter, k 12 , for each pair between unlike molecules. The radial distribution functions of CO 2 and sites in butanes around each sites in butanes are calculated as a fundamental information on their microscopic structures. These results show that the radial distribution functions and the number of supercritical CO 2 and n-C 4 H 10 around n-C 4 H 10 are different from those of supercritical CO 2 and i-C 4 H 10 around i-C 4 H 10 .
6. Nomenclature gi jŽ r . k ki j L N na, n b
radial distribution function Boltzmann constant intersite interaction parameter between sites i and j bond length number of molecules in central cubic box number of sites in molecules a and b
M. Yamamoto et al.r Fluid Phase Equilibria 163 (1999) 165–173
nŽ r . p r ri j T V Z
number of molecules or sites within distance r pressure distance distance between sites i and j absolute temperature volume of central cubic box compressibility factor
Greeks ´i j fab F Q r si j z z Ž0. – z Ž3.
energy parameter between sites i and j intermolecular potential between molecules a and b rotational angle bond angle number density size parameter between sites i and j intramolecular potential intramolecular potential parameters
Subscripts 1 2 C i, j
carbon dioxide butane critical property sites i and j
173
Superscript reduced value
U
Acknowledgements The authors gratefully acknowledge the financial support in part provided by the Japan Society for the Promotion Science ŽJSPS-RFTF96P00401. .
References w1x w2x w3x w4x w5x w6x w7x w8x w9x w10x
W.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. 106 Ž1984. 6638. T. Tsuji, S. Honda, T. Hiaki, M. Hongo, J. Supercrit. Fluids 13 Ž1998. 15. Y. Iwai, Y. Koga, Y. Arai, Fluid Phase Equilib. 116 Ž1996. 267. Y. Koga, Y. Iwai, Y. Arai, J. Chem. Phys. 101 Ž1994. 2283. Y. Koga, Y. Iwai, M. Yamamoto, Y. Arai, Fluid Phase Equilib. 131 Ž1997. 83. Y. Iwai, H. Uchida, Y. Koga, Y. Arai, Y. Mori, Ind. Eng. Chem. Res. 35 Ž1996. 3782. T. Holderbaum, J. Gmehling, Fluid Phase Equilib. 70 Ž1991. 251. J.J. Nicolas, K.E. Gubbins, W.B. Streett, D.J. Tildesley, Mol. Phys. 37 Ž1979. 1429. R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, 4th edn., McGraw-Hill, New York, 1987. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Univ. Press, New York, 1987.