Langmuir 2000, 16, 8529-8535
8529
Solidification of Lennard-Jones Fluid in Cylindrical Nanopores and Its Geometrical Hindrance Effect: A Monte Carlo Study Hideki Kanda,† Minoru Miyahara,* and Ko Higashitani Department of Chemical Engineering, Kyoto University, Kyoto 606-8501, Japan Received December 20, 1999. In Final Form: July 5, 2000 Grand canonical Monte Carlo (GCMC) simulations have been conducted to investigate freezing phenomena of Lennard-Jones (LJ) methane confined in cylindrical nanopores of reduced pore diameter ranging from 4.5 to 9.5 with saturated vapor as equilibrium bulk condition. The arrangement of the LJ particles in frozen state is found to form a hexagonal array within each concentric circular layer. Nonmonotonic dependence of freezing point against pore diameter, which is observed in strongly attractive pore made of carbon, is interpreted as a result of competition between “geometrical hindrance effect” of a wall’s shape and “compression effect” of a wall’s attractive potential. Geometrical hindrance effect, which originates from excess energy of “cylindrical frozen state” relative to the fcc lattice state, is clarified to become stronger for smaller pore size. The latent heat is also smaller for smaller pore size. We model the shift of freezing point in cylindrical pore considering the “geometrical hindrance effect,” and its reliability is verified successfully through comparison with the GCMC results.
1. Introduction While physical adsorption, capillary condensation, and capillary separation in narrow pores have been much studied theoretically and experimentally, and they are relatively well understood,1-5 freezing and melting in nanopores remain largely unexplored. Freezing and melting in nanopores are basic phenomena of possible application to fabrication of nanomaterials and characterization of porous media. There are a number of experimental studies about freezing/melting phenomena for various combinations of adsorbates and confining materials,6-19 which mostly reported “depression of freez* To whom correspondence should be addressed. Tel: +81-75-753-5582. Fax: +81-75-753-5913. E-mail: miyahara@ cheme.kyoto-u.ac.jp. † Present address: Chemical Energy Engineering Department, Central Research Institute of Electric Power Industry, Nagasaka 2-6-1, Yokosuka City, Kanagawa Prefecture 240-0196, Japan. (1) Evans, R.; Marconi, U. M. B.; Tarazona, P. J. Chem. Phys. 1986, 84, 2376. (2) Miyahara, M.; Yoshioka, T.; Okazaki, M. J. Chem. Phys. 1997, 106, 8124. (3) Yoshioka, T.; Miyahara, M.; Okazaki, M. J. Chem. Eng. Jpn. 1997, 30, 274. (4) Huber, P.; Knorr, K. Phys. Rev. B 1999, 60, 12657. (5) Miyahara, M.; Kanda, H.; Yoshioka, T.; Okazaki, M. Langmuir 2000, 16, 4293. (6) Patrik, W. A.; Kemper,W. A. J. Phys. Chem. 1938, 42, 369. (7) Rennie, G. K.; Clifford, J. J. Chem. Soc., Faraday Trans. 1 1977, 73, 680. (8) Tell, J. L.; Maris, H. J. Phys Rev. B 1983, 28, 5122. (9) Warnock, J.; Awschalom, D. D.; Shafer, M. W. Phys. Rev. Lett. 1986, 57, 1753. (10) Torii, R. H.; Maris, H. J.; Seidel, G. M. Phys Rev. B 1990, 41, 7167. (11) Sokol, P. E.; Ma, W. J.; Herwig, K. W.; M. Snow, W.; Wang, Y.; Koplik, J.; Banavar, J. R. Appl. Phys. Lett. 1992, 61, 777. (12) Strange, J. H.; Rahman, M.; Smith, E. G. Phys. Rev. Lett. 1993, 71, 3589. (13) Unruh, K. M.; Huber, T. E.; Huber, C. A. Phys. Rev. B 1993, 48, 9021. (14) Moltz, E.; Wong, A. P. Y.; Chan, M. H. W.; Beamish, J. R. Phys. Rev. B 1993, 48, 5741. (15) Klein, J.; Kumacheva, E. Science 1995, 269, 816. (16) Duffy, J. A.; Wilkinson, N. J.; Fretwell, H. M.; Alam, M. A.; Evans, R. J. Phys. Cond. Matter 1995, 7, L713. (17) Schaefer, B.; Balszunat, D.; Langel, W.; Asmussen, B. Mol. Phys. 1996, 89, 1057.
ing point.” However most of them treated solidification in silica or similar solids, e.g., silica gels, porous glass, and MCM-41. Recently, Miyahara and Gubbins conducted molecular simulation of freezing/melting of LJ-methane confined in slit-shaped nanopores modeled on graphite,20 and clarified a compression effect, which is caused by strong interaction between pore wall and confined fluid, elevates freezing point. Radhakrishman and Gubbins21 and SliwinskaBartkowiak et al.22 examined freezing point of LJ fluid in slit pores obtained by the Monte Carlo method with the Landau free energy calculation, and demonstrated qualitative agreement with experimental results from the viewpoint of the wall’s potential. Other than slit-shaped pores, Celestini et al. simulated melting of LJ clusters confined in spherical and ellipsoidal geometry and found depression in freezing point, but the effect of external potential field was not examined.23 For cylindrical pores, Maddox and Gubbins pointed out that the geometrical constraint in pore shape hinders freezing.24 The most inner layer confined in cylindrical pore did not freeze when cooled below normal freezing point, even with strong interaction of wall. The freezing point in cylindrical pore was not quantitatively modeled because of its complexity, e.g., geometrical constraint, compression effect and effect of equilibrium pressure in bulk phase. In this paper, we investigate freezing behavior of simple fluid in cylindrical nanopores employing grand canonical Monte Carlo (GCMC) simulations, and modeling of the freezing point in cylindrical pore is attempted taking account of geometrical constraint and compression effect. The reliability of the model is tested with simulation results under various interaction strengths between pore wall and fluid particles. (18) Brown, D. W.; Sokol, P. E.; Ehrlich, S. N. Phys. Rev. Lett. 1998, 81, 1019. (19) Huber, P.; Wallacher, D.; Knorr, K. Phys., Rev. B 1999, 60, 12666. (20) Miyahara, M.; Gubbins, K. E. J. Chem. Phys. 1997, 106, 2865. (21) Radhakrishnan, R.; Gubbins, K. E. Mol. Phys. 1999, 96, 1249. (22) Sliwinska-Bartkowiak, M.; Gras, J.; Sikorski, R.; Radhakrishnan, R.; Gelb, L.; Gubbins, K. E. Langmuir 1999, 15, 6060. (23) Celestini, F.; Pellenq, R. J.-M.; Bordarier, P.; Rousseau, B. Z. Phys. D 1996, 37, 49. (24) Maddox, M. W.; Gubbins, K. E. J. Chem. Phys. 1997, 107, 9659.
10.1021/la991659p CCC: $19.00 © 2000 American Chemical Society Published on Web 10/03/2000
8530
Langmuir, Vol. 16, No. 22, 2000
Kanda et al.
2. GCMC Simulation The GCMC procedure that Miyahara and Gubbins used for solidification in slit pore20 was applied to cylindrical pores: The pore system traced vapor-liquid coexistence line for bulk Lennard-Jones (LJ) fluid by setting corresponding combination of temperature and chemical potential. The method by Adams25-27 was used for the algorithm. A structureless solid wall of cylindrical geometry was employed for the simulation. The employment of a structureless solid is appropriate for a large molecule compared with the scale of potential corrugation on the pore wall,20 such as methane in carbon nanotube. Fluid employed was LJ-methane with ff/k and σff being 148.1 K and 0.381 nm, respectively. The cutoff distance was set as 5σff.
[( ) ( ) ] σff d
u(d) ) 4ff
12
σff d
6
(1)
For the pore potential, we used that for cylindrical pore surrounded by semi-infinite solid, which was derived by Peterson et al.28 with cylindrical coordinate integration. At a radial position r, the interaction between a fluid particle and pore wall with radius R is given as
(r,R) ) πfsFs φcylinder wall
[
]
7σ12 fs K (r,R) - σ6fsK3(r,R) 32 9
(2)
where
∫0πdΘ[-Rr cos Θ + (1 - (Rr )
Kn(r,R) ) R-n
2
1/2 -n
) ]
sin2 Θ
and Fs is the number density of interaction sites in the carbon wall, being 1.14 × 1029 m-3. The parameter ss/k and σss were 28.0 K and 0.340 nm, respectively.29,30 Lorentz-Berthelot combining rules were applied to obtain parameters for the methane-carbon interactions. The reduced pore diameter D* ) D/σff ranged from 4.5 to 9.5. A small difference in pore diameter may cause variation in freezing behavior in connection with the mismatching of pore diameter with particle size, and this is handed over to Appendix A. The reduced pore length Ly* ) Ly/σff was 10, which was confirmed to be sufficient to obtain freezing point correctly: A setting of Ly* ) 30 gave the same freezing point. The periodic boundary condition was imposed in the longitudinal direction. Vapor-liquid coexistence relation for LJ fluid obtained by Kofke and Agrawal31,32 was used to calculate corresponding chemical potential µ at a given temperature T with the following equation.
µ/kT ) 3 ln Λ + ln (P/kT)
(3)
where
Λ)
h
x2πmkT
(25) Adams, D. J. Mol. Phys. 1974, 28, 1241. (26) Adams, D. J. Mol. Phys. 1975, 29, 307. (27) Adams, D. J. Mol. Phys. 1979, 32, 647. (28) Peterson, B. K.; Walton, J. P. R. B.; Gubbins, K. E. J. Chem. Soc., Faraday Trans. 2 1986, 82, 1789. (29) Steele,W. A. Surf. Sci. 1973, 36, 317. (30) Steele, W. A. The Interaction of Gases with Solid Surface; Pergamon: Oxford, 1974. (31) Kofke, D. A. J. Chem. Phys. 1993, 98, 4149. (32) Agrawal, R.; Kofke, D. A. Mol. Phys. 1995, 85, 43.
Figure 1. Change of overall density in cylindrical pore on cooling.
Λ is the thermal de Broglie wavelength and P is bulk phase pressure. The mass of the methane is m ) 2.665 × 10-26 kg. h is Planck’s constant. Preliminary examination on comparison between cooling and heating sequences found some degree of hysteresis that is similar to, or somewhat smaller than, those in slit pores. For slit pores, it was shown that states observed in cooling sequence are closer to the true equilibrium point than those in heating sequence in the previous report.20,33 Here in this study, then, freezing points in cooling sequence starting from sufficiently high temperature were observed. The last configuration obtained in a run was used as the initial configuration for the next run with a lower temperature. Each run consisted of at least 4 × 104 GCMC steps per molecule, each of which contained insertion, deletion and displacement trials. Near the freezing transition calculations up to 3 × 105 GCMC steps per molecule were made. 3. Results of GCMC Simulation for Carbon Pores Results for variation of overall density in various pores of D* ) 4.5-9.5 on cooling from high temperature are shown in Figure 1. The overall density F* of adsorbate in pore is:
F* ) Fσ3ff: F )
〈N) V
and V )
D2 π × 10σff 4
(4)
where 〈N〉 is ensemble average of number of particles in pore. V is the volume of pore space, which is determined as the space between edges of nuclei of surface of carbon atom of wall and includes some dead space in the vicinity of wall where methane cannot penetrate by the repulsive force. For smaller pores, this region forms a larger fraction of the total volume, thus giving a smaller apparent density. Until the temperature drops down to about 120 K, the densities show gradual increase on cooling, and the slope against temperature is almost the same as that for bulk fluid.31,32 With cooling to 110-116 K, overall density in a pore gives a discrete change against temperature. With further cooling, it reaches a plateau and exhibits little change against temperature. The change in density is accompanied by steplike variation in enthalpy at this point. The total internal energy per particle in the pore is written as
UFF + UFW 3 U ) kT + 2 N UFF )
u(rij) ∑ i>j
(5) (6)
Solidification of Lennard-Jones Fluid
Langmuir, Vol. 16, No. 22, 2000 8531
Figure 2. Variation of enthalpy in pores with various pore diameters.
UFW )
(ri,R) ∑i φcylinder wall
(7)
Figure 3. Local density profiles of particles over r-direction. (a) D* ) 7.5, (b) D* ) 9.5
where UFF and UFW are the fluid-fluid (FF) and the fluidwall (FW) contributions. The ensemble average of enthalpy of particles in pore is given as
h ) 〈U + PVP〉
(8)
where P is external pressure and VP is volume per molecule. The second term is so small against the internal energy term that it may be neglected. Enthalpies in the pores are shown in Figure 2 with that for bulk phase obtained by Agrawal and Kofke,31,32 which exhibit discrete changes on cooling. The latent heat is smaller for pore with smaller diameter. The structure of confined fluid at the lowest temperature is analyzed. In the cylindrical pore, frozen methane takes arrangement that differs from the fcc lattice observed in the bulk. The frozen particles have lined up in concentric circular along the wall similarly to those reported by Maddox and Gubbins.24 The local density FL*(r) in the radial direction was calculated as
F*L(r) ) FL(r)σ3ff: FL(r) ) 〈N(r)〉/[Lyπ{(r + ∆r/2)2 - (r - ∆r/2)2}] (9) where N(r) is number of particles in the space confined between r - ∆r/2 and r + ∆r/2. The thickness ∆r was taken to be 0.05σff. Figure 3 demonstrates that, at the temperature with steplike variation in density and enthalpy, the local density profile in r-direction also changes its nature. Every layer completely separates from each other after the change into lower temperature, and the sharpness of the layers increases. Figure 4 shows quasi-two-dimensional snapshots of particles before and after the change, for each concentric layer. The circumferential position zr is given by the following equation:
zr ) riθ
Figure 4. Quasi-two-dimensional snapshots of typical configurations in carbon cylindrical pore of D* ) 7.5: Contacting circular layer (a) T ) 116.6 K, (b) T ) 115.7 K.; second circular layer (c) T ) 116.6 K, (d) T ) 115.7 K.; third circular layer (e) T ) 116.6 K, (f) T ) 115.7 K.
(10)
where ri is the radial position where FL*(r) takes the maximum in each layer, and θ is the angle position in cylindrical coordinate. With slight change in the temperature, all layers of particles in pore show ordered states, which form hexagonal arrays along each concentric circular layer. Each hexagonal arrangement is superimposed in Figure 5 with the abscissa converted from zr into θ. A particle in the second layer is not located at the gravity position of triangle composed of three particles beneath in the contact layer. Similarly, a particle in the third layer is not located at the center of the triangle in the second
Figure 5. Quasi-two-dimensional snapshots of three circular layers superimposed on angular coordinates.
layer. The arrangement shows that interaction between particles within the common concentric layer in θ-direction take precedence over that between the radial direction. In-plane pair distribution function for ith layer gi(dr) and its Fourier transform, or the structure factor Si(k) are analyzed also to examine the freezing behavior, where dr
8532
Langmuir, Vol. 16, No. 22, 2000
Kanda et al.
Figure 6. In-plane pair distribution function in pore of D* ) 7.5.
Figure 7. In-plane structure factor in pore of D* ) 7.5.
denotes distance within a concentric plane. gi(dr) and Si(k) are calculated from quasi-two-dimensional particle configurations within concentric plane, or (y,zr) shown in Figure 4:
gi(dr) ) Si(k) ) 1 + 2πΓi
〈N(dr;ri)〉 2πdr∆drΓi
(11)
∫0∞drJ0(kdr)gi(dr)ddr
(12)
∫(r(r +r+r )/2)/2FL(x)dx
(13)
Γ* i) Γiσ2ff: Γi )
i
i+1
i-1
i
where Γ is in-plane density, N(dr;ri) means the number of particles in the distance between dr - ∆dr/2 and dr + ∆dr/2 from a reference molecule in the ith layer bounded by r ) (ri-1 + ri)/2 and (ri + ri+1)/2. The pair distribution function in Figure 6 shows liquid state at T ) 116.6 K with similar feature as that in bulk phase. With a slight depression in temperature, the function for the contact and second layer indicates solidification. It is judged by the features of the curves: The first minimum goes down to zero, the second peak separates, and the third peak exhibits shoulder. The ordering of the third layer is not clearly recognizable. The two-dimensional structure factor in Figure 7 also shows solidification of the contact and second layer from the sharpness and the height of the first peak: Several studies34-36 proposed, as a criterion for freezing in a purely two-dimensional system, that the first peak in S(k) reaches a value around 4.4-5.0 on solidification. Figure 8 demonstrates difference of inner layers from outer ones in another aspect. The third layer and the innermost layer shows increase in Γi* with further (33) Miyahara, M.; Kanda, H.; Shibao, M.; Higashitani, K. J. Chem. Phys. 2000, 112, 9909. (34) Ranganathan, S.; Pathak,K. N. Phys. Rev. A 1992, 45, 5789. (35) Caillol, J. M.; Levesque, D.; Weis, J. J.; Hansen, J. P. J. Stat. Phys. 1982, 28, 325. (36) Broughton, J. Q.; Gilmer, G. H.; Weeks, J. D. Phys. Rev. B 1982, 25, 4651.
Figure 8. Variation of in-plane density on cooling in pore of D* ) 7.5.
cooling than the freezing point. It means that the fluid in inner layers may not solidify perfectly at the freezing point, which is in accord with results reported by Maddox and Gubbins.24 It is stressed again that the freezing points determined above, from density, enthalpy, arrangement, and structural functions, show nonmonotonic dependence against pore diameter in carbon pores with strongly attractive potential. 4. Results for Various Interaction Strengths between Wall and Fluid To examine dependence of the geometrical difficulty on the pore diameter, we observed solidification phenomena in several pores that have the same interaction strength as that between LJ-methane particles, or “methane pores”. In the methane pores, any compression effect would not arise, and observed shift in freezing point should show directly the contribution of “geometrical hindrance effect.” The procedure was basically the same as that for the carbon pores. Equation 14 was employed as the methane wall’s potential. cylinder φCH (r,R) ) πffFf 4
[
]
7σ12 ff K (r,R) - σ6ffK3(r,R) 32 9
(14)
Solidification of Lennard-Jones Fluid
Langmuir, Vol. 16, No. 22, 2000 8533
Figure 9. Variation of overall density on cooling in cylindrical methane pores.
where Ff is the number density of solid methane, being Ffσ3ff ) 0.962.31,32 With the saturated vapor pressure for the bulk state, methane particles did not stay confined in the methane pore unlike in the carbon pore. Then, in this section, we employed slightly higher value of µ, and set so that µ/kT is constant: The pressure for the bulk was higher than the saturated vapor pressure. It is, strictly speaking, not appropriate for quantitative examination of geometrical hindrance effect, since the contribution of the pressurization is included for the shift of the observed freezing point. However, it is important at least for qualitative understanding of the effect. Figure 9 shows the variations of the overall density. Depression is the case for all pore sizes of methane pores, and the decrease in freezing point is more remarkable for smaller pore diameter. The observed depression in freezing point implies the following: Geometrical hindrance effect, which is thought to originate from excess energy of “cylindrical frozen state” relative to the fcc lattice state, is stronger for smaller pore. Miyahara and Gubbins showed that elevation in freezing point in slit pore is stronger as pore size becomes smaller because of stronger attractive potential.20 Then, nonmonotonic variation of freezing point against pore diameter, which is observed in carbon pores, can be interpreted as a result of competition between “geometrical hindrance effect” of the wall’s shape and “compression effect” of the wall’s potential, both of which being stronger for smaller pores. Similar simulations were carried out with various interaction strengths: ss/k were set to be 9.3, 15.5, 18.6, 21.8, and 24.9 K. Other parameters were the same as those used for the carbon pore. These pores have weaker interaction than the carbon pore and stronger than the methane pore. Results are shown in Figure 10. For some values of ss, freezing points in smaller pores are lower than the bulk even with wall’s interaction stronger than “methane” pore: It corresponds to a condition that geometrical hindrance effect is stronger than compression effect. 5. Simple Model for Solidification in Cylindrical Pores In the modeling of freezing point in cylindrical pores, we consider an excess chemical potential energy of frozen state in cylindrical pores relative to the fcc lattice. The equilibrium for the bulk fluid at T ) Tf is simply described as bulk µbulk S (Tf) ) µL (Tf)
(15)
Figure 10. Fitting of model to observed shifts in freezing point: (a) D* ) 9.5, (b) D* ) 7.5, (c) D* ) 5.5.
where the superscript bulk indicates the chemical potential for the bulk state. The subscripts S and L indicate solid and liquid phase, respectively. Solid phase within the cylindrical pore and liquid phase in the bulk are in equilibrium at T ) Tf + δT, and this condition is expressed by the following equation. bulk µpore C (Tf + δT) ) µL (Tf + δT)
(16)
where the superscript pore indicates particles in pores, and the subscript C indicates frozen state in cylindrical pores. Difference in temperature gives change in chemical potential for the bulk liquid by eq 17, assuming the entropy to be constant over the temperature range δT. bulk µbulk L (Tf + δT) ) µL (Tf) - SLδT
(17)
where S is entropy. Next, we model the chemical potential of the frozen state in cylindrical pores at T ) Tf + δT, taking chemical potential of the bulk solid at T ) Tf to be a reference state. bulk geo - SCδT + ∆φ (18) µpore C (Tf + δT) ) µS (Tf) + µ
The second term in the right-hand side µgeo stands for the excess chemical potential of “cylindrical solid” relative to the fcc lattice at Tf. The third term shows change of the chemical potential by heating of “cylindrical solid”. The last term is the elevating effect by attractive interaction from pore wall. The ∆φ of course can be known for simulation system, but it can be estimated also in an experimental system if so-called the standard isotherm is available. For further detail of ∆φ, Appendix B should be referred.
8534
Langmuir, Vol. 16, No. 22, 2000
Kanda et al.
The following equation is obtained by substituting eqs 17 and 18 into eq 16, and using eq 15.
0 ) µgeo +
∆H′ δT + ∆φ Tf
(19)
where
∆H′ ) (SL - SC)Tf
(20)
The above equation is solved for δT and rearranged to give
∆φ + µgeo δT )Tf ∆H′
(21)
This model equation is a generalization of the model for freezing point in slit pores δT/Tf ) -∆φ/∆Hbulk, proposed earlier:20 For slit pores µgeo ) 0 and ∆H′ is close to ∆Hbulk. 6. Verification of Proposed Model The proposed model is tested by comparison with the GCMC results. Figure 10 shows the relation between ∆φ and δT. The simulation data appear to have a linear trend in accord with the expectation by our model, except those for D* ) 5.5, which is probably because of breakdown of the continuum assumption employed in the model. We attempt a least-squares fitting of the proposed model to the data for the larger pores. The fitting is made for the range of δT > 0, since the model stands for the elevation in freezing point. For pore of D* ) 7.5, slope determined by the fitting is -8.2 × 1022 K/J and the δT-intercept is -17 K, which corresponds to be µgeo ) 0.10ff and ∆H′ ) 0.61ff. Similar fitting gives µgeo ) 0.057ff and ∆H′ ) 0.46ff for D* ) 9.5. It is reasonable that ∆H′ determined for the two systems are smaller than 1.0ff for the bulk state, and roughly in agreement with those observed in simulations. The µgeo indicates larger values for smaller pores, indicating consistency with the expected larger hindrance for smaller pores. Further, it is interesting that the δTintercepts agree with the freezing points in methane pores for both diameters despite its premise of elevated freezing point in the derivation. Summarizing, for cylindrical nanopores, freezing point is successfully explained by the assumption of the change in the latent heat and the excess chemical potential brought by the geometrical hindrance effect. 7. Conclusion Grand canonical Monte Carlo simulations for solidification of LJ-methane fluid in cylindrical carbon pores show that the freezing point in the pores are higher than that for the bulk phase, and the most inner layer of confined phase does not clearly solidify. The frozen phase takes hexagonal arrangement in concentric circular layers along the wall, whose specific structure raise geometrical hindrance effect for solidification. Freezing points are determined by density, enthalpy, arrangement of particles, and structural factor. The change of density is accompanied by steplike variation in enthalpy and ordering of particles at this point. In cylindrical methane pore, where compression effect does not arise, freezing point monotonically decreases as the pore diameter becomes smaller. The tendency proves that the geometrical hindrance effect is larger for smaller pores. Thus, the nonmonotonic dependence of freezing point on pore diameter, which is observed in carbon pores, is the consequence of competition between “geometrical hindrance effect” of wall’s shape
Figure 11. Relationship between freezing points and pore diameter. The freezing points (open circles) approach maxima when the layer separation nears σff.
and “compression effect” of wall’s potential. The latent heat on freezing observed in cylindrical pore is smaller than that for the bulk, and becomes smaller for smaller pore sizes. We model the shift of freezing point in cylindrical pore, considering the “geometrical hindrance effect”. Freezing point predicted by the proposed model shows a linear relation with the potential from pore wall. By comparison between the proposed model and GCMC simulations with several interaction strengths of pore wall, the model successfully shows its validity. Appendix A: Pore-size Hindrance Effect Caused by Mismatching between Pore Size and Particle Size As well as the geometrical difficulty, matching between particle size and pore size is also important for shift of freezing point. Similar simulations were performed for carbon pore of several D* from 6.5 to 7.9 with interval of 0.1, to determine freezing points. The same procedure as described in the text gave overall density, enthalpy, particle arrangement, and structural factor. Figure 11 shows the relationship between freezing points and pore diameter. The average separation between adjacent layers ∆l for a given diameter D* was obtained from the local density distribution as the averaged distance between adjacent peaks. When ∆l nears σff, the freezing point takes a maxima. The freezing point decreased from the maximum with deviation of ∆l from σff, and the solidification itself were not clear in pores where ∆l was most different from σff., or around D* ) 7, even with further cooling from the normal freezing point. This kind of behavior, related to the mismatching, may need additional parameter to be explained precisely by a model. For most cases, however, the freezing point shifts stay around 10-15 K, and general trend can be interpreted as described in the text. Appendix B: Direct Observation of Geometrical Hindrance Effect in Carbon Pores It is naturally supposed that the potential field is an important factor for solidification in cylindrical pore, where geometrical difficulty is also a factor of the shift in freezing point in contrast to slit pore. It is difficult, however, to determine the two effects independently with simple comparison between slit and cylinder, because the difference in geometry brings at the same time difference in potential field. Then, we examined the contribution of the geometrical hindrance effect by GCMC simulations that compared solidification in (1) cylindrical pore and (2) fictitious slit pore with the same potential energy as that in cylindrical geometry. Wall’s potential of cylindrical pore
Solidification of Lennard-Jones Fluid
Langmuir, Vol. 16, No. 22, 2000 8535
have stronger potential profile than that in slit pore because of the overlap of potential energy in a cylindrical pore. Then, a common setting for the two systems in the effective pore wall’s potential ∆φ was used for the comparison. ∆φ, which expresses excess potential energy relative to fluid’s solid state, is given in the following equation,2,5,20,33 and is certainly deducible from experimental isotherm data as demonstrated in ref 37. cylinder ∆φ(r,R) ) φcylinder (r,R) - φCH (r,R) wall 4
(B1)
Then, the fictitious slit pore of width H with common ∆φ as above should have the following potential φi. slit slit (x) + φCH (H - x) + ∆φ(x,R) (B2) φi(x,H) ) φCH 4 4
The first and second term together stand for methane pore of slit geometry, for which we used the LJ(10-4-3) potential of Steele29,30 composed of methane particles. slit (x) ) φCH 4
[( )
2πFfffσ2ff∆ff
2 σff 5 x
10
-
() σff x
4
-
σ4ff
]
3∆ff(x + 0.61∆ff)3
(B3)
where the reduced separation between lattice planes, ∆ff/σff is 0.929. Pore width H was the same as the pore (37) Kanda, H.; Miyahara, M.; Yoshioka, T.; Okazaki, M. Langmuir 2000, 16, 6622.
Figure 12. Shifts of the freezing point δT for cylindrical pores (open circles) and for fictitious slit pores with same ∆φ (squares). Difference between the two shows the geometrical hindrance effect (triangles).
radius of the cylinder in comparison. We took the same procedure as that of GCMC simulation by Miyahara and Gubbins.20 Shifts of the freezing point δT is shown in Figure 12. In the fictitious slit pore, only compression effect by attractive potential contributes to δT. The difference between the two systems implies the extent of the geometrical hindrance effect, which increases as pore diameter becomes smaller, similarly to the results for methane pores. It was also observed that the latent heat in cylindrical pore was smaller by about 25% than that in slit-shaped pore, because solidified particles cannot take the closest packing structure by the geometrical difficulty. LA991659P
