7600
Langmuir 2001, 17, 7600-7604
Modeling of Adsorption and Desorption in Pores of Simple Geometry Using Molecular Dynamics L. Sarkisov and P. A. Monson* Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003 Received August 3, 2001. In Final Form: September 27, 2001 We present results from molecular dynamics simulations of adsorption and desorption by diffusive mass transfer into model pores of well-defined geometry. The pore geometries considered were chosen to investigate issues related to classical theories of hysteresis in adsorption. We consider the slit pore with or without a closed end, a wedge or triangular pore, and finally a model with an ink bottle geometry. In each case we show that the molecular dynamics simulations yield the same results for adsorption and desorption as grand canonical Monte Carlo simulations of the same systems. Using computer graphics visualization, we show how hysteresis for the slit and wedge geometries is related to the evolution of vapor-liquid interfaces in a manner largely consistent with the classical picture of these systems. We find that for the ink bottle geometry the large pore can empty during desorption even while the small pore remains filled with fluid. This behavior is different from the classical picture of adsorption and desorption in this geometry based on the concept of pore blocking.
I. Introduction A common feature of adsorption in mesoporous materials is the occurrence of hysteresis between the isotherm on adsorption and that on desorption. A classification of hysteresis loops has been recommended as an aid in the use of adsorption for characterization of porous materials.1,2 The essential phenomenology and the early history of research on hysteresis are described in a review by Everett.3 Perhaps the simplest picture of adsorption hysteresis is that of metastable states in single pores analogous to the supercooled liquid and superheated vapor states in bulk systems. A theoretical treatment of hysteresis for a single infinitely long cylindrical or slit pore was first given by Hill.4 This concept has received a further development to include a size distribution of independent pores in works by, for example, Everett and co-workers,3 who developed the “independent domain” theory, and by Ball and Evans,5 who employed a mean field density functional theory of inhomogeneous fluids in their analysis. Another approach is based on the Kelvin equation and was developed earlier by Cohan.6 In his work Cohan suggested that hysteresis occurs in a single pore because of the difference in geometry of menisci during adsorption and desorption. For example, if a pore is a cylinder open at both ends, capillary condensation is nucleated by the film of fluid on the walls of the cylinder and the relative pressure of this event according to the Kelvin equation corresponds to the radius of the cylinder minus thickness of the film. However, evaporation from the full pore takes place from the hemispherical menisci at each end of the cylinder. The relative pressure of this event will be determined by the radii of curvature of the menisci. Thus, * Author to whom correspondence should be addressed. (1) Rouqerol, F.; Rouqerol, J.; Sing, K. Adsorption by Powders and Porous Solids; Academic Press: New York, 1999; Chapter 7. (2) Sing, K. S. W.; Everett, D. H.; Haul, R. A. W.; Moscou, L.; Pierotti, R. A.; Rouqerol, J.; Siemieniewska, T. Pure Appl. Chem. 1985, 57, 603. (3) Everett, D. H. In The Solid-Gas Interface; Flood, E. A., Ed.; Dekker: New York, 1967; Vol. 2, pp 1055-1113. (4) Hill, T. L. J. Chem. Phys. 1947, 15, 767. (5) Ball, P. C.; Evans, R. Langmuir 1989, 5, 714. (6) Cohan, L. H. J. Am. Chem. Soc. 1938, 60, 433.
condensation and evaporation occur at different relative pressures, which gives rise to the hysteresis. Using similar considerations, one can argue that for a closed end pore, hysteresis should not be observed since both adsorption and desorption are governed by a hemispherical meniscus. This idea has been further developed by Saam and Cole7 who approached the problem from a point of view of stability of adsorbed multilayers. Their theory relates hysteresis in a simple pore to the difference between the critical film thickness tc (on adsorption) and the residual film thickness tr on desorption. Marini Bettolo Marconi and van Swol8 considered the question of the role of meniscus development in hysteresis using a lattice model of a slit pore and a mean field density functional theory. They showed that adsorption behavior of a fluid in a slit pore that is placed in contact with bulk fluid is very different from that of a fluid in a slit pore in periodic boundaries (i.e., an infinite pore). Having a vaporliquid interface introduced on desorption significantly reduces and rounds off the desorption branch of a hysteresis loop. On the other hand, a closed end slit pore has also a reduced adsorption branch of hysteresis, due to a different vapor-liquid interface development path. Later these observations were confirmed using grand canonical Monte Carlo technique9 and density functional theory10 on off-lattice models of open ended and closed ended slit pores. The concepts we have reviewed up to now focus on a single pore or a distribution of independent pores. However, hysteresis in real materials should in principle involve collective phenomena across larger regions of the material. To model this we need to go beyond a single pore length scale and consider interconnected porous structures. One way to approach this problem is to consider a network of pores of a simple geometry. The simplest version of this approach is the so-called “ink bottle” (7) Saam, W. F.; Cole, M. W. Phys. Rev. E 1975, 11, 1086. (8) Bettolo Marini Marconi, U.; van Swol, F. Europhys. Lett. 1989, 8, 531 (9) Papadopoulou, A.; van Swol, F.; Marini Bettolo Marconi, U. J. Chem. Phys. 1992, 97, 6942. (10) Kozak, E.; Chitiel, G.; Patrykiejew, A.; Sokolowski, S. Phys. Lett. A 1994, 189, 94.
10.1021/la015521u CCC: $20.00 © 2001 American Chemical Society Published on Web 10/30/2001
Adsorption in Model Pores
geometry.11 The concept of hysteresis in this kind of a pore is based on the assumption that evaporation of the liquid from a relatively large cavity should be delayed if the only available path for it is through a narrow channel. More sophisticated pore network theories have been developed using percolation concepts.12-14 In all these models hysteresis is thought to occur due to pore blocking effects as explained in detail by Mason12 and Seaton.13 In the last several years novel approaches have been developed to combine the grand canonical ensemble with molecular dynamics.15-18 These approaches have created the possibility to probe directly the dynamics associated with adsorption hysteresis and to address questions that could only be the subject of speculation previously. In this work we use grand canonical molecular dynamics to address the following issues. What is the nature of hysteresis in a simple single pore? How do the processes of adsorption and desorption depend on the geometry of vapor-liquid interfaces in the system? What is the evolution of these interfaces during adsorption and desorption? Can pore blocking effects be observed for a model ink bottle structure? We will also relate our simulations to theoretical predictions and more traditional grand canonical Monte Carlo simulations on the same systems. This will allow us to further investigate the link between the Metropolis Monte Carlo algorithm and diffusive mass transfer mechanisms, which was established in our previous work.19 II. Methodology The approach we use in this paper was originally developed for disordered systems,19 such as silica gels, and is based on grand canonical molecular dynamics with control volume technique.15-18 In this methodology, a system of interest is composed of two regions, the first one containing a model adsorbent and another one being the bulk fluid region. The bulk fluid region serves as a control volume in the terminology of van Swol and co-workers,15 i.e., a source and sink of fluid particles via the conventional grand canonical Monte Carlo scheme. Fluid particles within the whole system are free to move by constant temperature molecular dynamics. Thus, the porous material region is accessible only by diffusion. It is important to notice that for the cases studied here, the bulk region of the system contains a vapor phase at all states of interest, so the events inside the porous material are not skewed by possible transitions in the bulk region. The equilibrium properties of the adsorbed fluid are then sampled from a central part of the porous material region at fixed temperature and chemical potential. This setup essentially imitates a dynamic uptake experiment in that adsorption and desorption take place via diffusive mass transfer mechanisms rather than by addition and deletion of molecules anywhere in the system as in the grand canonical Monte Carlo method. The model pore geometries we consider in this work are illustrated in Figure 1. In each case we use a single twodimensional layer of hexagonally packed spheres for the pore walls. The fluid adsorbing in these structures is a (11) McBain, J. W. J. Am. Chem. Soc. 1935, 57, 699. (12) Mason, G. Proc. R. Soc. London, Ser. A 1988, 415, 453. (13) Seaton, N. A. Chem. Eng. Sci. 1991, 46, 1895. (14) Guyer, R. A.; McCall, K. R. Phys. Rev. B 1996, 54, 18. (15) Papadopoulou, A.; Becker, E. D.; Lupkowski, M.; van Swol, F. J. Chem. Phys. 1993, 98, 4897 (16) Heffelfinger, G. S.; van Swol, F. J. Chem. Phys. 1994, 100, 7548 (17) MacElroy, J. M. D. J. Chem. Phys. 1994, 101, 5274 (18) Cracknell, R. F.; Nicholson, D.; Quirke, N. Phys. Rev. Lett. 1995, 74, 2463 (19) Sarkisov, L.; Monson, P. A. Langmuir 2000, 16, 9857.
Langmuir, Vol. 17, No. 24, 2001 7601
Figure 1. Computer graphics visualization of the pore geometries considered in this work. From top to bottom: a slit pore, a closed end slit pore, a wedge, and an ink bottle pore. The volume on the right outlines the bulk region of the system, serving as a control volume in the grand canonical molecular dynamics simulations. The length of each pore is 40σ (where σ is the diameter of the solid sphere). The pore width (the distance between the planes for the slit or for the closed end slit pore, or the size of the mouth for the wedge) is 10σ. In the direction perpendicular to the picture, the layers of width 10σ are placed in periodic boundaries. For the ink bottle geometry, the length of the whole pore is 40σ and the distance between the layers in the central cavity is 10σ. The width and the length of each of the narrow channels of the ink bottle structure are 5σ and 10σ. respectively.
Lennard-Jones 12-6 fluid. The fluid-fluid and fluidsolid potentials are taken to be identical 12-6 potentials with well depth, , and collision diameter, σ. All potentials were cut and shifted at 2.5σ. The system setup we consider here is similar to that in our previous work,19 but unlike in our previous work, the adsorbent is modeled with a single simple pore. The fluid particles can originate and disappear in the right subvolume (control volume) and then diffuse in the pore. For the closed end slit and wedge geometries, the bulk system is bounded on one side by a repulsive wall. Periodic boundaries were used in the direction normal to the plane of the figure and also for the ink bottle geometry and open slit in the direction parallel to the channels in the plane of the figure. In the molecular dynamics simulations the systems were run up to 85 × 106 dimensionless time steps ∆t/(σ(m/)1/2) ) 4 × 10-3. All the systems were kept at a constant temperature kT/ ) 0.75 (well below the bulk critical temperature of the cut and shifted 12-6 potential) using the Brown and Clark20 velocity rescaling method. At the highest pressures, the systems contained up to 400 particles in the control volume and about 2700, 2500, 1100, and 1200 particles inside the slit, the closed slit, the wedge, and the ink bottle pore, respectively. We have also made grand ensemble Monte Carlo simulations of the systems. For these simulations up to 75 × 106 configurations were used with up to 50 × 106 configurations for equilibration. III. Results A. Adsorption/Desorption in Slit, Closed Slit, and Wedge Pores. In Figure 2 we show an adsorption/ desorption isotherm for a slit geometry calculated via grand canonical Monte Carlo and molecular dynamics. We plot the adsorbed fluid density versus the relative pressure of the bulk vapor, i.e., the pressure divided by (20) Brown, D.; Clark, J. H. R. Mol. Phys. 1984, 51, 1243.
7602
Langmuir, Vol. 17, No. 24, 2001
Figure 2. Adsorption/desorption isotherms of dimensionless density (Fσ3) vs relative pressure (P/P0) at kT/ ) 0.75 for an open ended slit geometry, calculated via grand canonical Monte Carlo (open circles) and molecular dynamics (filled squares, adsorption; filled triangles, desorption).
Figure 3. Computer graphics visualizations of various stages of adsorption and desorption in a slit geometry. The adsorption column corresponds to a sequence of increasing pressures (from top to bottom); the desorption column corresponds to a sequence of decreasing pressures (from top to bottom).
the saturated vapor pressure. There are several important things to notice. There is excellent agreement between the molecular dynamics and Monte Carlo results. Both methods predict the same sequence of states for the system on adsorption and desorption. The nature of this agreement has been discussed in our previous work.19,21 In Figure 3 we show computer graphics visualizations of the key states the system goes through on adsorption and desorption. On adsorption the system goes through several layering stages and then it undergoes a capillary condensation. On desorption, however, we see hemispherical menisci developing on the ends of the pore, followed by a transition as these menisci approach each other. This process was qualitatively predicted by Cohan.6 On the other hand, this picture confirms previous investigations by Marini Bettolo Marconi and co-workers.8,9 Indeed, the development of the menisci is impossible for a slit pore in periodic boundaries. This leads to a smaller and rounded off desorption branch of the isotherm for the open-ended pore in contact with gas phase. Following a similar argument, for a closed end slit pore we should also observe development of the hemispherical meniscus on adsorption, (21) Sarkisov, L.; Monson, P. A. In Fundamentals of Adsorption; Kaneko, K., Ed., in press.
Sarkisov and Monson
Figure 4. Adsorption/desorption isotherms of dimensionless density (Fσ3) vs relative pressure (P/P0) at kT/ ) 0.75 for an open-ended slit geometry (open circles) and for a closed end slit geometry (filled squares, adsorption; filled triangles, desorption), calculated via grand canonical Monte Carlo.
Figure 5. Computer graphics visualization of various stages of adsorption and desorption in a closed end slit geometry. The adsorption column corresponds to a sequence of increasing pressures (from top to bottom); the desorption column corresponds to a sequence of decreasing pressures (from top to bottom).
and as a result the evolution of the states on adsorption and desorption should be very similar; i.e., no significant hysteresis loop should be observed. In Figure 4 we compare adsorption/desorption isotherms from Monte Carlo simulations for a closed end slit pore and for an open end slit pore. Closing one end of the pore dramatically changes the behavior. Now the adsorption branch is also rounded off and almost coincides with the desorption branch. The key configurations that the closed end pore goes through are shown in Figure 5. As one can see, on adsorption the system goes through layering, then it develops a semispherical meniscus similar to that on desorption. On desorption, the process is very similar to that for the openended slit pore system. This picture is in clear accord with the Cohan predictions6 and with more recent work.8,9 In Figure 6 we compare the results from Monte Carlo simulations and from the detailed molecular dynamics approach. The agreement between two methodologies is quite convincing. To take this investigation a step further, we considered adsorption in a generic wedge geometry. One would anticipate the absence of hysteresis or any transitions. This behavior is due to the continuous layering, secured by the geometry of the pore, as was predicted by Everett3 from macroscopic considerations. Indeed, we show in
Adsorption in Model Pores
Langmuir, Vol. 17, No. 24, 2001 7603
Figure 8. Computer graphics visualization of various stages of adsorption and desorption in the wedge geometry. The adsorption column corresponds to a sequence of increasing pressures (from top to bottom); the desorption column corresponds to a sequence of decreasing pressures (from top to bottom).
Figure 6. Adsorption/desorption isotherms of dimensionless density (Fσ3) vs relative pressure (P/P0) at kT/ ) 0.75 for a closed end slit geometry, calculated via grand canonical Monte Carlo (open squares, adsorption; open triangles, desorption) and molecular dynamics (filled circles, adsorption; filled triangles, desorption).
Figure 7. Adsorption/desorption isotherms of dimensionless density (Fσ3) vs relative pressure (P/P0) at kT/ ) 0.75 in the wedge geometry, calculated via grand canonical Monte Carlo (open circles) and molecular dynamics (filled squares).
Figure 7 that this hypothesis is correct. Again we see a good agreement between molecular dynamics and Monte Carlo simulations. On the other hand, the system exhibits no hysteresis and no phase transition, which is very different from slit pore behavior considered earlier. The final point on the adsorption isotherm at high pressure corresponds to a complete filling of the wedge pore with fluid. In Figure 8 we show snapshots of several stages the system goes through on adsorption and desorption. B. Adsorption/Desorption in an Ink Bottle Pore. The ink bottle geometry considered in this work consists of one central rectangular cavity and two narrow slit channels connecting the cavity with the bulk fluid phase. This setup is similar to that recently considered by Maddox et al.,22 but there are significant differences. In their work,
Maddox and co-workers had a system that had no bulk regions and the narrower pores served as control volumes.22 Here, we consider a system placed in a direct contact with bulk gas phase (see Figure 1). The bulk region serves as control volume, which avoids some risk that the location of the control volume in a pore will influence the dynamics. On the other hand, Maddox and co-workers studied a combination of relatively narrow slits and strong solid-fluid field, where micropore filling23 was observed for both larger and smaller slits in the junction rather than a clear capillary condensation. Here we are interested in studying a system where the smallest pore in the junction undergoes a capillary condensation, so that there is a significant change in the adsorbed fluid density between an empty pore and filled pore. According to the traditional view of hysteresis in the ink bottle geometry, the system should go through the following stages. On adsorption, the small cavities become filled with fluid first, then the large central cavity fills. On desorption, fluid in the large cavity cannot evaporate until the narrow channels empty so hysteresis occurs. In Figures 9 and 10 we show the behavior of the system on adsorption and desorption as calculated via molecular dynamics and also by grand canonical Monte Carlo simulations. The system goes through several qualitatively different stages on adsorption and desorption. First, an initial layering of molecules in both the large and small cavities takes place. The second stage, corresponding to a knee in the isotherm at a relative pressure of 0.35, is associated with filling of the small cavities with fluid. At this stage each of the small cavities contains about 300 molecules. This is followed by further diffusion of fluid into the large cavity. The final sharp jump on the adsorption isotherm at a relative pressure of about 0.75 is associated with the complete filling of the large cavity. At this point this cavity contains about 1200 molecules and the whole system about 1800 molecules. On desorption, the large cavity remains filled with liquid until lower pressures relative to those where filling on adsorption occurs. Then evaporation of fluid from the large cavity takes place while the small cavities remain filled with fluid, as can be clearly seen from the isotherm in Figure 9 and from the computer graphics visualization in Figure 10. The next stage is associated with the evaporation from the small cavities, a process that exhibits some hysteresis itself due to some variation in vapor-liquid interfaces geometry on adsorption and desorption. Thus we present a case, where emptying of a small cavity connecting a larger cavity with the gas phase is not a prerequisite for (22) Maddox, M. W.;Quirke, N.; Gubbins, K. E. Mol. Simul. 1997, 19, 267. (23) Gregg, S. J.; Sing, K. S. W. Adsorption, Surface Area and Porosity; Academic Press: New York, 1982; Chapter 3.
7604
Langmuir, Vol. 17, No. 24, 2001
Sarkisov and Monson
The way we have modeled these simple pore geometries using a single layer of solid atoms was chosen for computational simplicity. It might be argued that the phenomenology would be different if we included more layers of atoms in the pore walls, thereby increasing the strength and range of the solid-fluid interactions. This might, for example, act to slow the overall dynamics on desorption, especially from the corners of the large cavity in the ink bottle geometry. Alternatively, it might lead to formation of solid like layers in the small cavities, which would drastically reduce the rate of diffusion through these cavities. Both of these effects were seen in the work of Maddox et al.22 However, the classical picture of the ink bottle geometry is based on there being fluid states in both the small and large cavities, and it is this situation we wanted to address. IV. Conclusions
Figure 9. Adsorption/desorption isotherms of dimensionless density (Fσ3) vs relative pressure (P/P0) at kT/ ) 0.75 for the ink bottle geometry, calculated via grand canonical Monte Carlo (open circles) and molecular dynamics (filled squares, adsorption; filled triangles, desorption).
Figure 10. Computer graphics visualization of various stages of adsorption and desorption for the ink bottle geometry. The adsorption column corresponds to a sequence of increasing pressures (from top to bottom); the desorption column corresponds to a sequence of decreasing pressures (from top to bottom).
the evaporation from the larger cavity. It is important to emphasize here that we have compared the dynamics results with Monte Carlo simulations. The good agreement between these two approaches indicates that the states of the pore fluid along the isotherms in the dynamics calculations correspond to equilibrium with the bulk fluid so that the chemical potentials of those states are thus known and equal to that of the bulk fluid. The hysteresis in our model ink bottle pore is not caused by pore blocking. So what does cause the hysteresis and what is the mechanism by which the liquid can evaporate from the large cavities with the small cavities remaining filled? Our molecular dynamics simulations show that during this process evaporation occurs from the small cavities into the bulk and the molecules lost are replaced by those from the large cavity. This is similar to the process on adsorption where the large cavity fills via mass transfer from the bulk gas through the small cavities. Hysteresis occurs because on desorption there is essentially nothing in the system to nucleate a large scale evaporation process in the large cavity until the pore liquid reaches a sufficiently expanded state (analogous to a thermodynamic stability limit in the bulk) that spontaneous local density fluctuations can lead to cavitation.
We have shown that many classical concepts of adsorption and desorption in simple pore geometries can be carefully tested using a molecular dynamics method in which adsorption and desorption occur by diffusive mass transfer. Using this technique, we have shown the importance of the evolution of vapor-liquid interfaces in the development of hysteresis in a single slit pore consistent with classical ideas. We have further strengthened the link between the Metropolis Monte Carlo algorithm24 in the grand ensemble and grand canonical molecular dynamics. It is clear that the Metropolis algorithm is capable of correct reproduction of the vaporliquid interface development for the geometries considered in this work. Arguably the most significant feature of this work has been the use of the molecular dynamics technique to test the pore blocking concept of hysteresis for the case of an ink bottle geometry. The traditional picture of adsorption and desorption in an ink bottle cavity is not supported by our computer simulations. Mass transfer of evaporating liquid from the large cavity can occur even though the neck of the ink bottle remains filled with liquid. The wider significance of this observation is that if pore blocking does not occur in the ink bottle geometry, its importance in other situations is now subject to question. Theories of hysteresis that are built on the concept of pore blocking have until now provided the most plausible explanation for adsorption hysteresis phenomena for complex porous media.5,12-14 Interestingly, recent work has shown that a very plausible alternative theory of hysteresis, which does not invoke the concept of pore blocking, can be developed by using more realistic models for the porous material microstructure and statistical mechanics in the grand canonical ensemble.25,26 The fact that results from the molecular dynamics and Monte Carlo simulations are essentially identical to each other in both the present case and for other more sophisticated models of fluids in porous materials18 lends further support to this approach. Acknowledgment. This work was supported by the National Science Foundation (Grant No. CTS-9906794). LA015521U (24) Metropolis, N.; Rosenbluth, A. W.;Rosenbluth, M. N.; Teller, A. N.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (25) Kierlik, E.; Monson P.; Rosinberg M. L.; Sarkisov, L.; Tarjus, G. Phys. Rev. Lett. 2001, 87, 055701. (26) Sarkisov, L.; Monson, P. A. Phys. Rev. E, in press.
