Langmuir 1997, 13, 4659-4668
4659
Adsorption of Polar Gases on Model Silica Gel Peter A. Gordon and Eduardo D. Glandt* Department of Chemical Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6393 Received April 10, 1997X Two models are developed to study the influence of polar surface groups present in silica gel on adsorption of polar gases. Both are based on the composite sphere model of Kaminsky and Monson (J. Chem. Phys. 1991, 95, 2936-2948). The surface groups are incorporated either as discrete dipoles or as a continuum (smeared) dipolar bilayer. The behavior of each of these descriptions is explored through the calculation of Henry’s constants in various microstructures, water adsorption isotherms, isosteric heats of adsorption, and adsorbate structural correlation functions. Comparison to experimental data shows that while both models produce good qualitative agreement, only the discrete approximation results in a model adsorbent that is sufficiently heterogeneous to predict the correct trends in the isosteric heats of adsorption.
Introduction Materials such as activated carbon and silica gel form the basis for a vast array of separation processes, including oxygen production from air, desiccation, wastewater treatment, and cloud seeding.2,3 Because novel applications often hinge on the properties of the adsorbent, much effort must be invested in characterizing the surface and microstructure of such materials and in assessing their practical merits. Our ability to understand these properties will invariably aid in finding new applications for currently available adsorbents as well as in developing improved materials to fulfill specific needs. There has been much work in the application of molecular theory to the study of fluids adsorbed in micropores. Many studies have employed simple slit or cylindrical geometries to represent pore spaces. More recently, several groups have focused on the behavior of fluids confined within pores of more complex, irregular geometries. Recent approaches to this problem have utilized integral equation theory, density functional theory, other mean field theories, and computer simulations. Most efforts have concentrated on simple fluids, characterized by spherically symmetric Lennard-Jones (LJ) interactions. Relatively little effort has been applied to the fundamental study of hydrogen-bonded fluids in disordered materials. Recently, Segarra and Glandt4 developed a graphitic platelet model of activated carbon that incorporates surface functional groups as dipoles, evenly smeared around the edges of graphitic carbon platelet disks. The arrangement of the model platelets was designed to mimic the tortuous nature of the pore space in activated carbons. In another series of contributions, Muller and co-workers5 have studied the behavior of LJ associating fluids in carbon slit pore models “activated” with specific surface sites capable of forming hydrogenlike bonds with the adsorbate. Both studies saw the characteristic cooperative phenomenon in the adsorption of water vapor on their model carbons. That two qualitatively dissimilar descriptions can produce similar ad* E-mail:
[email protected]. X Abstract published in Advance ACS Abstracts, August 1, 1997. (1) Kaminsky, R. D.; Monson, P. A. J. Chem. Phys. 1991, 95, 29362948. (2) Sircar, S.; Golden, T. C.; Rao, M. B. Carbon 1996, 34, 1-12. (3) Bassett, D. R.; Boucher, E. A.; Zettlemoyer, A. C. J. Colloid Interface Sci. 1970, 34, 436-446. (4) Segarra, E. I.; Glandt, E. D. Chem. Eng. Sci. 1994, 49, 29532966. (5) Muller, E. A.; Rull, L. F.; Vega, L. F.; Gubbins, K. E. J. Phys. Chem. 1996, 100, 1189-1196.
S0743-7463(97)00370-3 CCC: $14.00
sorption behavior suggests that further work needs to be done to clarify the relationship between the intermolecular forces present and the adsorption behavior of such complex systems. Kaminsky and Monson1,6,7 have pursued a promising line of investigation that directly examines the relationship between the structure of the microporous adsorbent and the adsorption behavior. Using their composite sphere (CS) model of amorphous silica gel, they have investigated how adsorbent packing arrangements influence Henry’s constants and adsorption uptake of methane, and selectivities for several simple fluid mixtures. This work, however, has been limited to spherically symmetric nonpolar molecules described by Lennard-Jones potentials. In this paper, we focus on the role of polarity in the behavior of heterogeneous microporous materials. We employ alternative models of silica gel designed to capture the interactions with a hydrogen-bonding fluid such as water. To describe the polar interactions we examine both a discrete-site model and a continuum model, and we compare the resulting energetic environments. We also examine the importance of geometric packing of these model microspheres and are able to identify what factors influence the degree of heterogeneity in the model adsorbent. Finally, we compare the predicted behavior of each adsorbent model to experimental data to evaluate the applicability of the models. Background and Models Employed Silica gel provides a rich framework in which to study the role of polarity in heterogeneous materials. The surface of a typical silica gel particle is characterized by fairly irregular grooves of (predominantly) OH end groups which comprise a hydroxylated surface.8 The density of these groups is highly dependent on adsorbent preparation. At room temperature, a fully hydroxylated surface will readily adsorb moisture and a hydrogen-bonded layer of water will exist on it. Heating to temperatures above approximately 170 °C will drive off most of the adsorbed water, leaving a so-called vicinal anhydrous surface. Further heating above 175 °C will drive off the remaining water and produce dehydrated siloxane patches on the surface. Hydroxylated and pyrogenic (i.e. relatively dehydroxylated) silica surfaces have been shown to have very (6) Kaminsky, R. D.; Monson, P. A. Chem. Eng. Sci. 1994, 49, 29672978. (7) Kaminsky, R. D.; Monson, P. A. Langmuir 1994, 10, 530-537. (8) Iler, R. K. The Chemistry of Silica; John Wiley & Sons: New York, 1979.
© 1997 American Chemical Society
4660 Langmuir, Vol. 13, No. 17, 1997
Gordon and Glandt
different water adsorption characteristics. Bassett and co-workers9 found that the effective water monolayer volume decreases when the distance between OH sites becomes large relative to the size of the water molecule. It has been shown10 that at low adsorbate loadings water tends to orient itself “oxygen down;” i.e., the water dipole is oriented away from the particle surface. As loading increases on the dehydrated surface, clusters of water molecules form before the completion of an adsorbed monolayer. This is due to the energetics of the system: water molecules find more stability by forming clusters away from the surface than by adsorbing separately at isolated silanol sites. The low density of the silanol sites prevents water from forming multiple hydrogen bonds. On the isolated hydroxylated surface, on the other hand, hydroxyl units are packed closely enough to facilitate the formation of hydrogen bonds between water and several sites at once. For this reason, water completes a monolayer before clusters start to form. In our studies, we employ the TIP4P model for water.11 Like other TIPS potential functions, the TIP4P model uses a mixture of partial charges and Lennard-Jones interaction sites at fixed positions to represent water and has been shown to accurately reproduce a wide variety of bulk liquid thermodynamic and structural properties. The potential between two water monomers i and j is given by a Lennard-Jones interaction between the oxygens plus a sum of site-site Coulombic interactions
[( ) ( ) ]
Uij ) 4OO
σOO
12
σOO
-
rOO
rOO
qaqb
6
+
∑a ∑a r
(1)
(
16 π F R3 3 sf s
(
d6 +
σO ) 3.154 Å qM ) -1.04e
O/k ) 78.1 K qH ) 0.52e
rOH ) 0.9572 Å ∠HOH ) 104.52°
rOM ) 0.15 Å
and the center of the composite sphere, σsf and sf are the solid-fluid collision diameter and well depth, respectively, and Fs is the number density of Lennard-Jones sites comprising the adsorbent. The interaction energy diverges as d approaches R, the effective hard core radius of the composite sphere. The physical parameters used here are the same as those determined by Kaminsky and Monson, namely, R ) 1.347 nm and Fs ) 44/nm3. The hard sphere packing fraction is η ) 0.386. The total gas-solid interaction energy is represented as the sum of the symmetric nonpolar term given by eq 2 and an orientationally-dependent polar term,
U(X B ) ) UCS(d) + Uµµ(X B)
(3)
where the contribution Uµµ depends on the specific approximation used for the dipolar interactions. Discrete Dipole Model. A straightforward representation of the presence of polar end groups chemisorbed on an adsorbent particle consists of placing a discrete number of point dipoles on the surface of the composite silica microsphere. The dipoles are oriented radially outward from the center of the particle. The potential between an adsorbate molecule with dipole moment µads at position X B and a surface dipole µs at X B 1 can be expressed in terms of the electric field at X B ,13
ab
The parameters that define the potential are listed in Table 1. We examined two different models of silica microspheres in order to assess how well they represent experimental observations of adsorption of hydrogen-bonding species. Both are based on simulation studies originally performed by MacElroy and Raghavan12 and later Kaminsky and Monson.1 The latter researchers proposed a simplified potential model to investigate the influence of adsorbent microstructure on adsorption characteristics. This model was developed with the goal of qualitatively representing a typical sample of silica xerogel. It consists of a collection of spherical adsorbent particles in an equilibrium hard sphere configuration. Each microsphere is modeled as a continuum of Lennard-Jones 12-6 interaction sites. The interaction of a Lennard-Jones adsorbate particle with one of these composite spheres has been shown to be given by1
UCS(d) )
Table 1. Parameters for TIP4P Potential
)
21 4 2 R6 12 σ d R + 3d2 R4 + 5 3 sf (d2 - R2)9 σ6sf (d2 - R2)3
)
(2)
In this expression, d is the distance between the adsorbate (9) Bassett, D. R.; Boucher, E. A.; Zettlemoyer, A. C. J. Colloid Interface Sci. 1968, 27, 649-658. (10) Klier, K.; Zettlemoyer, A. C. J. Colloid Interface Sci. 1977, 58, 216. (11) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (12) MacElroy, J. M. D.; Raghavan, K. J. Chem. Phys 1990, 93, 20682079.
B) ) E B 1(X
bs‚n b1) - b µs 3n b1(µ |X B-X B 1|3
(4)
where b n1 is the unit vector between X B 1 and X B . The total electric field at X B is simply the vector sum of the fields produced by each surface dipole, Ndip
E B (X B) )
E B i(X B ,X B i) ∑ i)1
(5)
and the dipolar potential experienced by the adsorbate is simply the scalar product of its dipole moment and the resultant electric field,
Uµµ(X B) ) - b µ ads‚E B (X B)
(6)
Smeared Dipole Model. In this approach, the radially-oriented point dipoles are “smeared” over the surface of the composite sphere. We carry out this process in such a way that the total polar strength on the microsphere remains constant. One may imagine increasing the surface density Fd of the dipoles while reducing the magnitude µ of each of them, such that the product Fdµs remains unchanged and expresses the effective polar strength of the dipoles on the microsphere. This approximation is entirely analogous to the smearing of discrete Lennard-Jones interaction sites employed in developing the composite sphere and the well-known 9-3 graphitic carbon model.14 In the limit when the dipoles are smeared over the microsphere, the sum in eq 5 is replaced by a surface integral, (13) Jackson, J. D. Classical Electrodynamics; John Wiley & Sons, Inc.: New York, 1962. (14) Steele, W. Surf. Sci. 1973, 36, 317-352.
Adsorption of Polar Gases on Model Silica Gel
E B (X B) )
∫∫ FµEB(XB) R2 sin θ dθ dφ
Langmuir, Vol. 13, No. 17, 1997 4661
(7)
Integrating over the entire microsphere (0 < θ < π, 0 < φ < 2π) results in a complete cancellation of the electric field. To preserve the polar nature of this model, we have used a “line-of-sight” truncation for the interaction between an adsorbate dipole and the dipoles on the microsphere. In other words, we have only calculated the interaction with the silica gel dipoles that are not blocked by the core of the spherical particle. This implies that the limits on θ in eq 6 vary between 0 and arccos(R/d). The maximum angle corresponds to the ring of dipoles where the vector X B-X B 1 is tangent to the microsphere surface. When these limits are imposed on the integral, the resultant interaction potential becomes
Uµµ(d) ) -
2πFdµsµadsR2 cos φd d2xd2 - R2
(8)
where φd is the angle between the adsorbate dipole vector and d is the vector between the center of the composite sphere and the adsorbate. This potential function is wellbehaved; as d approaches R, the electric field becomes large and the potential becomes infinite in magnitude. In the limit of d . R, the potential behaves as
Uµµ(d.R) ≈ -
2πFDµsµadsR2 cos φd d3
(9)
which is the same decay rate as for a dipole-dipole interaction. This expression gives a simple result that incorporates polarity into the potential description between the adsorbate and adsorbent yet preserves the spherical symmetry of the composite sphere model. Implementation and Properties of the Models. In the case of the discrete dipole model it might seem desirable to achieve a uniform spacing of dipoles on the surface. However, it can be shown that it is not possible to distribute an arbitrary number of points in a regular array over a spherical surface.15,16 Instead, we place them on the surface using a hard particle growth technique. A number of dipoles are arranged on the surface of the composite sphere and assigned an initial hard core radius. The hard core radius of one of the dipoles is then increased by an amount ∆r. If this growth results in overlap with other dipolar cores, the position of the growing core is displaced away from them. If a new position exists for which no overlaps occur, the growth move is accepted. If overlap cannot be avoided, then the local configuration is jammed and the growth move is rejected. After this procedure has been sequentially applied to each dipole, the system is subjected to 5Ndip displacements to keep it from reaching prematurely jammed configurations. To represent the hydroxylated surface of a silica gel microsphere, 122 point dipoles were placed on the surface of the composite sphere. This corresponds to a site density of 5.36 sites/nm2, which lies within 4.6 and 6.2 OH/nm2, the range of hydroxyl site densities reported in the literature.12 The initial positions of the particles were randomly specified, and their hard core radii were set to 1.577 Å. Particles were subjected to growth steps in size increments of ∆r ) 7.9 × 10-6 Å. The growth process was continued over 100 000 cycles, after which the growth and move acceptance rate per cycle had decayed to essentially zero. The average hard core radius of the dipoles on the (15) Mackay, A. L.; Finney, J. L.; Gotoh, K. Acta Crystallogr. A 1977, 33, 98-100. (16) Coxeter, H. S. M. Trans. N. Y. Acad. Sci. 1962, 24, 320-331.
Figure 1. Pair correlation function of dipoles for the discrete model. 122 hard spheres were grown on the surface of a sphere of radius R* ) 4.2691, as described in the text. The data are plotted as a function of the angular distance δ between points and the arc length reduced by the average hard sphere radius 〈σ〉.
surface was 〈σ〉 ) 2.23 ( 0.01 Å, suggesting that the final size distribution of the particles was uniform. Several different initial configurations produced similar structural characteristics of the final dipole arrangement. To examine the positions of the dipoles on the surface, a pair correlation function was calculated as a function of the distance, measured by the arc length on the surface. This function is defined as
g(r) )
no. of pairs separated by a distance r along the surface
N(r) )
no. of randomly distributed pairs separated by r
N(r)ideal
(10)
where
(cos(r/R) - cos((r + dr)/R)) N(r)ideal ) Ndip 2
(11)
In the latter expression, the expected number of pairs separated by a distance r scales with the ratio of the surface area between r and r + ∆r to the surface area of the entire microsphere. Figure 1 shows the calculated pair distribution function for the final configuration of particles on a composite sphere, plotted as a function of both the arc length distance reduced by 〈σ〉 and the angular distance δ. A major feature of interest here is the large initial peak, which reflects both the large number of dipole-dipole contacts and the very slight polydispersity of the cores. If their sizes were monodisperse, this peak would become a Dirac delta function at contact. Mackay and co-workers15 found that in the monodisperse limit the angular distance δ was 18.712 54° and the packing fraction η (fractional coverage of the surface) was 0.811 51. Using the data in Figure 1, the average value of the nearest neighbor distance is 〈δ〉 ) 18.79° ( 0.06. The corresponding fractional coverage is η ) 0.818 67. The slightly higher packing fraction calculated in our system stems from the slight degree of
4662 Langmuir, Vol. 13, No. 17, 1997
Figure 2. Electric field profile around discrete and smeared dipole model spheres. (b) Discrete dipole model, Fs ) 5.36 sites/ nm2, µs ) 1.5. The field is orientationally averaged and presented as a function of d, the reduced distance from the center of the sphere. Standard deviations resulting from the average are also depicted. (2) Smeared dipole model, Fdµs ) 1.67.
polydispersity present, yet it is clear that the general structural characteristics of the two systems are comparable. The arrangement of dipoles shows primarily hexagonal coordination with a small number of heptagonal structures present. The mixture of coordination numbers is reflected in the multiple (twin) peaks representing the second and third nearest neighbors. The structure of the dipoles on the microsphere can be viewed as glasslike, with strong short-range order that decays after just a few diameters. We expect that there will be inhomogeneities in the electric field produced by a single microsphere. The total electric field at a point near the adsorbent particle was calculated via eq 5 at regular intervals in the spherical coordinate system. The average electric field as a function of d, the distance from the center of the particle, is plotted in Figure 2. At each value of d the average and standard deviation were taken over 3600 separate values of φ and θ. For later convenience, we adopt the convention of reducing all quantities with respect to the Lennard-Jones parameters for oxygen in the TIP4P water potential (σ ) 3.154 Å, /k ) 78.1 K). In these units, the composite sphere surface lies at R* ) 4.27. The most evident feature of this plot is that the magnitude of the electric field is extremely inhomogeneous with respect to orientation. The standard deviation is typically 50-60% of the average value of the field at a given distance from the microsphere. This is in sharp contrast to the smeared dipole model, where the field, also shown in Figure 2, is a function only of the radial distance from the microsphere. We also note that the long-range behaviors of the discrete and smeared dipole models are different. A fit of the tail of Figure 2 over the range 6.7 < d < 14.9σff to the functional form AdB yields a power law coefficient of B ) -5.7. In contrast, the electric field of the discrete dipole model decays as B ) -3.3 over the same range. The faster decay rate is a consequence of the cancellation of fields brought about by dipoles on opposite sides of the microsphere.
Gordon and Glandt
We have also calculated the average angle between the resultant electric field vector at X B and the vector connecting the center of the microsphere and X B . For the smeared dipole model, this angle is always zero; there is perfect alignment between the resultant electric field and the center-to-center separation vector. The angle of alignment was computed as a function of d in a manner similar to that for the average electric field calculation. The angle was found to be fairly insensitive to d over the range 4.8 < d < 14.9σff and had an average value of approximately 50°. Thus, we see that even as many as 122 dipoles placed in a relatively regular arrangement on the microsphere can still produce very large variations in the resulting strength and orientation of the electrostatic field. In spite of the spherical symmetry of the smeared dipole model, one must keep in mind that the total electric field at a point in the pore space of the adsorbent is determined by the sum over all microspheres in the matrix. This provides an additional degree of inhomogeneity, since the fields generated by neighboring matrix particles may combine or cancel out at a given position, depending on the local environment. For instance, one would expect cancellation of fields to occur near the centers of the pore spaces. This multiparticle effect may tend to mask the differences observed in the electric field around an isolated smeared and discrete dipole model particle. Henry’s Constants Henry’s constants are a useful tool to characterize the model adsorbents due to their very high sensitivity to the local energetic environment. In order to examine the relationship between the matrix microstructure and the electrostatic environment, we calculate Henry’s constants for TIP4P water in both the smeared and discrete dipole models in a number of different matrix arrangements. In the grand canonical ensemble, the Henry’s constant KH is computed as the average interaction between the adsorbate and adsorbent in the limit of zero activity,
KH ) lim zf0
〈N〉 1 ) zV V
∫ exp(-βφ(rb1)) drb
(12)
where V is the volume of the system, β is 1/kT, φ is the gas-solid interaction energy, z is the gas phase activity of the adsorbate, and 〈N〉 is the average number of molecules in the system. This quantity can also be thought of as an equilibrium partition coefficient between the adsorbed and bulk gas phases, since the activity z approaches the bulk density F in the low-pressure regime where Henry’s law applies. The simulation details are similar to those of the method described by Kaminsky and Monson.1 For the smeared dipole model, we considered the fcc, bcc, and equilibrium hard sphere (ehs) packing arrangements; for the discrete dipole model, we studied the fcc and ehs matrices. The calculated dipole moment of 7.25 × 10-30 C‚m for the TIP4P water model was used in the calculation of the wateradsorbent interaction energy. For the ehs calculations, 32 microspheres were used in the simulation box, and a cycle consisted of 100 insertion and 32 hard sphere displacement attempts. The Henry’s constants for the smeared dipole model were averaged over 107 attempted insertions (100 000 cycles). The Henry’s constants for the discrete dipole model were averaged over a much longer run; between 5 × 106 and 15 × 106 cycles were performed, depending on the system. Uncertainties were calculated by the method of subblock averaging over 10 subsets of the generated Markov chain. Figure 3 shows the computed Henry’s constants for TIP4P water at T* ) 3.854 (301 K) on the smeared dipole
Adsorption of Polar Gases on Model Silica Gel
Langmuir, Vol. 13, No. 17, 1997 4663
Figure 3. Henry’s constants for TIP4P water interacting with smeared dipole model adsorbent for matrix packing fraction η ) 0.386 and reduced temperature T* ) 3.854. The various curves correspond to different matrix packing arrangements: (b) fcc lattice; (~) bcc lattice; (2) ehs packing.
model adsorbent as a function of the polarity of the adsorbent. The packing fraction was 0.386 in all the matrix arrangements considered. At low values of the surface polarity Fdµs, we see little variation between Henry’s constants among the various matrices considered. In this regime, where electrostatic forces make only a small contribution to the total adsorbate-adsorbent interaction energy, the dispersion energy dominates and seems fairly insensitive to the local arrangement of the matrix. Repulsion between the adsorbate and matrix determines the location of the favorable energetic sites. As the polar strength increases, the Henry’s constant for each matrix structure increases at a different rate. The constant for the ehs structure shows the sharpest ascent, while those for the bcc and fcc lattices increase more slowly. As the electrostatic effects become a larger fraction of the total interaction energy, it is clear that the energetic landscape viewed by the adsorbate changes markedly between the ordered and disordered environments. The insertion energy distributions help clarify the changes observed in the Henry’s constants for different structures. A histogram of attractive gas-solid energies was collected over 2 × 107 insertion attempts. Reported in Figure 4 are the gas-solid energy probability distributions for the smeared dipole model in fcc (i.e. ordered) and ehs (i.e. disordered) packing arrangements at several values of the surface polarity. The abscissa has been reduced by the value Umin, the minimum (i.e. most favorable) gas-solid interaction energy between a single model sphere and the adsorbate in its most favorable orientation. For zero polarity, a significant fraction of insertion energies in both the ehs and fcc arrangements result from overlap of potential wells between neighboring matrix particles. This accounts for the tail of the distributions extending to values of E/ Umin > 1. The lower energy insertions correspond to positions near the center of the pore spaces, where the adsorbate position does not coincide with any of the potential wells of the matrix particles.
Figure 4. Insertion energy probability distribution functions for TIP4P water in smeared dipole model micropores. The abscissa is reduced by Umin, the most favorable interaction energy between a water molecule and a single adsorbent microsphere. (a) Packing arrangement is the fcc lattice. The curves correspond to various degrees of adsorbent polarity: (O) Fdµs ) 0; (~) Fdµs ) 0.50; (]) Fdµs ) 2.91; (_) Fdµs ) 7.62. (b) ehs packing arrangement at (O) Fdµs ) 0, (~) Fdµs ) 0.97, (_) Fdµs ) 1.65, and (]) Fdµs ) 2.91.
As the polarity of the microspheres increases, the distribution becomes smoother, decaying more rapidly. There is less attractive cooperativity between neighboring microspheres, as evidenced by the shift of the energy distributions to higher values of E/Umin with increasing polarity. This can be interpreted in terms of a frustration effect. The water molecule cannot simultaneously optimize its orientation with respect to each and every matrix
4664 Langmuir, Vol. 13, No. 17, 1997
Gordon and Glandt
Figure 5. Henry’s constants for TIP4P water interacting with discrete dipole model adsorbent for matrix packing fraction η ) 0.386 and reduced temperature T* ) 3.854. The curves correspond to different matrix packing arrangements: (b) fcc lattice; (2) ehs packing.
particle; it must forgo its most favorable orientation with respect to one microsphere to avoid unfavorable dipolar interactions with other adsorbent particles. As polar terms become the larger contribution to the total potential energy, the frustration effect is magnified both for the fcc and ehs arrangements. Figure 5 shows the Henry’s constants for the discrete dipole model adsorbent in both the fcc and ehs configurations. The conditions are the same as those reported in Figure 4. The controlling parameter for surface polarity is now µdip. We see similar trends for KH as were seen for the smeared dipole model. In the present case, however, there is less difference between the fcc and ehs Henry’s constants at the same value of µdip. Figure 6 displays the corresponding energy distributions for the discrete dipole models. This plot employs an energy scale reduced by 〈Umin〉, the angularly averaged minimum well depth energy between water and the matrix particle. As the magnitude of the surface dipole moment increases, the first low-energy peak increases in size and shifts to energies closer to zero. These values again correspond to the centers of the pore spaces, where the components of the electric field from adjacent matrix particles tend to cancel out just as was seen for the smeared dipole model. One significant difference between the smeared and discrete dipole cases, however, is that the distribution for the latter does not shift significantly to lower values of E/〈Umin〉. In fact, for the fcc arrangement, increasing the surface dipole strength shifts the energy distribution toward more favorable energies. In this case, as the range of electric field overlap increases, regions of energetic “hot spots” also increase. In addition, we note that, at higher values of µdip, the fcc and ehs distributions in Figure 6 are quite similar. Variations in the discrete dipole model packing arrangement (i.e. geometric microstructure) have less impact on the adsorption energetics of the system than those in the case of the smeared dipole model. The probability distributions in Figures 4 and 6 are normalized such that their integrals over the range -∞
Figure 6. Insertion energy probability distribution functions for TIP4P water in discrete dipole model micropores. The abscissa is reduced by 〈Umin〉, the most favorable interaction energy for a water molecule and a single adsorbent microsphere. The energy has been orientationally averaged about the adsorbent. (a) Packing arrangement is the fcc lattice. The curves correspond to various degrees of adsorbent polarity: (O) µs ) 0; (~) µs ) 1.0; (]) µs ) 2.0; (_) µs ) 3.5. (b) ehs packing arrangement at (O) µs ) 0, (~) µs ) 0.88, and (]) µs ) 2.0.
< E < 0 yield the fraction of attractive insertion energies. This quantity is plotted in Figure 7 for each dipole model as a function of polarity. The net effect of increased polarity is clear; this attractive fraction remains almost constant over the polarity values considered for the discrete model, while it drops almost in half for the smeared case. The matrix-matrix interactions are much less cooperative for the former than for the latter.
Adsorption of Polar Gases on Model Silica Gel
Figure 7. Fraction of attractive insertion energies versus adsorbent polarity. Curves correspond to various combinations of dipolar model and packing arrangement: (O) discrete dipole model, ehs packing; (~) discrete dipole model, fcc packing; (4) smeared dipole model, ehs packing; (3) smeared dipole model, fcc packing.
Langmuir, Vol. 13, No. 17, 1997 4665
Figure 8. Henry’s constants of TIP4P water vs temperature on smeared and discrete dipole models. Both models are in ehs arrangements: (O) smeared dipole model, Fdµs ) 10.81; (~) discrete dipole model, µs ) 4.50; (9) data from Pedram et al.
In order to apply these models to realistic systems, we need to determine the surface polarity parameters, Fdµs* and µs*, respectively. Henry’s constants for water vapor on Mobil Sorbead R silica gel were obtained from the lowpressure adsorption data17 at 301, 310, and 326 K. The best-fit values of the surface polarity determined from the data at 301 K were Fdµs* ) 10.81 and µs* ) 4.5. Figure 8 compares Henry’s constants calculated over a range of temperatures with the experimental data of Pedram et al. The agreement is reasonable given the high sensitivity of the Henry’s constant to the adsorbent microstructure and energetics. However, the smeared dipole model slightly overestimates the rate at which KH increases with decreasing temperature. GCMC Simulations To examine the behavior of adsorbed water beyond the infinitely dilute (Henry’s law) regime, we have performed grand canonical Monte Carlo simulations.18 The microstructure employed in both the smeared and discrete dipole models was a representative hard sphere configuration. There will be additional uncertainties in the numbers reported here due to variations in the microstructure. The basic methodology and implementation of these simulations has been described elsewhere.19 An additional insertion biasing scheme was employed to help speed convergence of the simulations. This method was particularly helpful for simulations on the discrete dipole adsorbent. Details are given in the Appendix. The activity was converted to bulk pressure using the Bender equation of state for water.20 Between 4 × 106 and 40 × 106 attempted moves were performed for each value of the (17) Pedram, E. O.; Hines, A. L. J. Chem. Eng. Data 1983, 28, 1114. (18) Adams, D. J. Mol. Phys. 1975, 29, 307-311. (19) Allen, M. P.; Tildesley, D. J. In Computer Simulation of Liquids; Publications, O. S., Ed.; Clarendon Press: Oxford, 1987. (20) Polt, A.; Maurer, G. Fluid Phase Equilibr. 1992, 73, 27-38.
Figure 9. TIP4P water loading isotherms at T ) 326 K. Both models are in ehs arrangements: (O) smeared dipole model; (~) discrete dipole model; (9) data from Pedram et al.
activity before a stable equilibrium was achieved, and an identical number of moves were used to collect thermodynamic averages and structural correlation functions. Figure 9 shows the equilibrium loading isotherms of TIP4P water on the discrete and smeared dipole adsorbents at T ) 326 K. The loading was converted to weight of water per gram of adsorbent using a silica gel density of 1.129 g/cm3, the value employed by Kaminsky and Monson.7 For both adsorbent models, the equilibrium loading curves are qualitatively similar to the experimental values reported by Pedram et al.,17 yet they
4666 Langmuir, Vol. 13, No. 17, 1997
Figure 10. Adsorbate-matrix pair correlation functions for water adsorbed on smeared dipole model adsorbent at 326 K. Curves correspond to different equilibrium loadings. (‚‚‚) 0.28 mg of H2O/g of adsorbent; (- - -) 7.9 mg of H2O/g of adsorbent; (s) 89 mg of H2O/g of adsorbent.
overpredict the loading at intermediate pressures. This disagreement may stem from a variety of physical features of the model adsorbent that were not changed in specifying its parameters. The general trends observed, however, should still be valid when considering the behavior of the real system. As we are primarily concerned with the sub-monolayer regime of loading, an estimate for the monolayer coverage was obtained by calculating the surface area available to particles that reside at the first layer distance from the adsorbent particles.21 The area occupied by a water molecule was estimated from the bulk density of water as 10.5 Å2. This analysis yields a monolayer coverage of approximately 250 mg of H2O/g of adsorbent. Although it has been observed that dense hydroxyl population on the surface of silica gel leads to compression of the physisorbed layer of water,9 we are still far below monolayer coverage even at the highest relative pressures considered. The adsorbate-matrix pair correlation functions for adsorption on the smeared dipole model depicted in Figure 10 reveal that the traditional picture of monolayer adsorption does not hold in this situation. At very low loadings, there is a large peak associated with adsorption at a distance corresponding to the first layer, followed by a depletion region before the correlation approaches unity. As loading increases, we notice the formation of a second peak. Given that we are still far below the predicted monolayer coverage, it is clear that clusters of water molecules are formed and pore filling occurs. This type of behavior is facilitated by the small dimension of the pore spaces relative to the adsorbate particle size. Similar behavior is seen for adsorption on the discrete dipole model as well. The isosteric heats of adsorption shown as a function of loading in Figure 11 mark a departure from the behavior observed between the two adsorbent models. For the (21) Rikvold, P. A.; Stell, G. J. Colloid Interface Sci. 1985, 108, 158173.
Gordon and Glandt
Figure 11. Isosteric heats of adsorption vs adsorbent loading for smeared and discrete dipole model adsorbents at 326 K. Points correspond to loading data reported in Figure 9. (b) smeared dipole model; (~) discrete dipole model; (- - -) heat of vaporization at 326 K.
discrete model adsorbent, the heat of adsorption decreases with increasing loading, approaching the value of the heat of vaporization. The most energetic sites are filled at low pressures, whereas less favorable sites are progressively occupied at higher coverages. The smeared dipole adsorbent, on the other hand, exhibits a relatively uniform heat of adsorption with loading. In this case, cooperative interactions between adsorbate molecules approximately compensate for the energetic heterogeneity of the adsorbent. The relatively uniform heats of adsorption in the smeared dipole model can be explained through Figure 4, where it is seen that increasing the polarity of the smeared dipole model decreases the number of “hot spots” in the adsorbent. The smeared dipole model produces a more homogeneous energetic landscape, behavior not in accord with experimental observations for water adsorption. Conclusions In this work we have presented a discrete and a continuum representation of the functional groups present in silica gel. We have used them in conjunction with the composite sphere model of Kaminsky and Monson to investigate the influence of polarity in water adsorption, and we have compared the model prediction to experimental data. While reasonably good agreement was achieved, it must be pointed out that no attempt was made to match the model to the particular physical characteristics of the gel used by Pedram et al. and that the parameters were chosen independently of the data. Thus, accurate quantitative agreement with the experimental data is not expected. An optimized fit would have included further considerations on adsorbent porosity, microsphere particle size and degree of polydispersity, and pore size distribution, all features that vary from sample to sample in silica gel, depending on the method of preparation. In both models, increasing polarity produced larger Henry’s constants in disordered configurations than in the more ordered bcc and fcc lattice arrangements at the same packing fraction. The electric field produced by a
Adsorption of Polar Gases on Model Silica Gel
single smeared dipole microsphere is longer-ranged and more homogeneous than that of the discrete model. The differences are due to small inhomogeneities in the discrete dipole arrangement and also to the truncation method employed in the smeared dipole model. Correspondingly, the energetic environment within a configuration of smeared dipole microspheres was found to become more homogenous with increasing polarity than that of a bed of discrete-dipole microspheres. These differences were not significant for water adsorption isotherms, where the mechanism was observed to be pore filling rather than simple monolayer adsorption. The effect of cooperativity of the adsorbing water molecules was demonstrated to be significant in this mesoporous model material. The heat of adsorption, however, which is also dependent on the energetic landscape of the adsorbent pore space, showed significantly different behavior at low loadings. The results suggest that the smeared dipole model may be too homogeneous to describe the adsorption characteristics. Given these observations, it appears that the discrete dipole model provides a more complete picture of the hydrophilic adsorbent. The discrete model should also prove to be more robust when extended to partially hydrophilic or hydrophobic materials, where the distribution of sites will influence the clustering effect of polar species. In addition, adsorption from nonpolar and polar mixtures is expected to be more sensitive to the presence of discrete sites and is the subject of ongoing work.
Langmuir, Vol. 13, No. 17, 1997 4667
Figure 12. Evolution of loading during GCMC simulation for insertion biasing schemes. Starting configuration of simulation is an empty pore at T ) 326 K and z* ) 1.0 × 10-5 (from bottom of figure to top): (s) no biasing; (- - -) positional biasing; (- - -) rotational biasing; (‚‚‚) rotational and positional biasing.
∫∫f(rb,ω) drb dω ) 8π2V
Appendix A. Insertion Bias Methods One of the major obstacles toward achieving equilibrium in the grand canonical ensemble is a low insertion and deletion move acceptance probability. Typically, GCMC simulations fail to converge in systems where insertions and deletion moves are difficult to accomplish, such as in dense liquids or long-chain molecules. In our systems, the large orientationally-dependent interaction energies between TIP4P water and the matrix particles make insertions and deletions difficult even at low adsorbent loadings. Insertion attempts fail with greater frequency relative to those for simple fluids; even if a favorable adsorption site is found, a favorable orientation must also be chosen. In addition, molecules residing at favorable adsorption sites will be difficult to dislodge through deletions if the gas-solid interaction energy is high. Thus, the Markov chain of states has a tendency to get locked in small regions of phase space. While the method is still ergodic, the rate at which the Markov chain approaches equilibrium will decrease to the point where it becomes computationally impractical to run such simulations. To achieve greater efficiency in the simulation, we must improve the rate at which the Markov chain moves through phase space. This can be accomplished via insertion biasing methods, a variety of which have been developed by previous workers and tailored for specific applications.22-26 . Following the general language of Cracknell et al., a generalized insertion function f(r b,ω) is defined as directly proportional to the probability of inserting a molecule at a position b r with orientation ω. It is subject only to the normalization constraint (for nonlinear molecules) that (22) Swope, W. C.; Andersen, H. C. J. Chem. Phys. 1995, 102, 28512863. (23) Mezei, M. Mol. Phys. 1980, 40, 901-6. (24) Mezei, M. Mol. Phys. 1987, 61, 565-582. (25) Rao, M.; Pangali, C.; Berne, B. J. Mol. Phys. 1979, 37, 17731798. (26) Cracknell, R. F.; Nicholson, D.; Parsonage, N. G. Mol. Phys. 1990, 71, 931-943.
(A1)
The probability of accepting insertion and deletion moves can be written as
(
Pins ) min 1,
(
))
1 zV exp -β∆Eins + ln N +1 f(r b,ω)
(A2)
and
(
(
b,ω) exp -β∆Edel + ln Pdel ) min 1, f(r
N zV
))
(A3)
The normal unbiased insertion scheme chooses the six degrees of freedom (three positional coordinates and three orientation angles) with equal probability, and f(r b,ω) ) 1 everywhere in the simulation box. We consider here two types of insertion bias functions. The first scheme biases insertion positions only. The simulation box is divided into a large number of subcells of side length δ. A subcell is used for insertion attempts if its center does not overlap with the silica microspheres. If Nins subcells are used for insertions, the insertion function is defined as
f(r b) )
L3 Ninsδ3
(A4)
where L is the simulation box length. In the second scheme, we examine a rotational insertion bias method. We note that we know a priori the electric field at any point b r where we attempt an insertion. Using this information, we construct insertion functions to bias the orientation of the adsorbate dipole in the direction of the electric field. In terms of the Euler angles27 (φ, Θ, ψ), we choose to bias the first two angles according to Gaussian (27) Goldstein, H. Classical Mechanics; Addison-Wesley Publishing Co.: Reading, PA, 1980.
4668 Langmuir, Vol. 13, No. 17, 1997
Gordon and Glandt
distributions centered around the orientation of the electric field vector (φE, ΘE )
f(φ) )
c1
x2πσ1
(
)
-(φ - φE)2
exp
2σ21
,
∫-πf(φ) dφ ) 2π π
(A5)
and
f(cos θ) )
c2
x2πσ2
(
exp
-(cos θ - cos θE) 2σ22
)
2
,
∫-11f(cos θ) d cos θ ) 2
(A6)
The constants c1 and c2 are determined by integrating each (separable) insertion function over the allowable values. The widths of these distributions are specified by the adjustable parameters σ1 and σ2. These distributions may be sampled using a rejection method.28 To compare the efficiency of these biasing schemes, we performed a simulation of water on the discrete dipole model adsorbent at z* ) 1.0 × 10-5 and T ) 326 K starting at zero loading. Figure 12 shows the rates of approach to equilibrium for the unbiased method and each biased method. The required number of moves to achieve a given loading decreases by a factor of nearly 5 with the addition of rotational and positional biasing. We note that it took approximately 6 × 107 moves to reach a stable equilibrium (28) Press, W. H. Numerical Recipes in FORTRAN: The Art of Scientific Computing; Cambridge University Press: Cambridge, 1994.
Table 2. Comparison of Acceptance Rates and Overlaps in Simulations
type of biasing none positional bias rotational bias positional and rotational bias
% gas-solid overlaps
% gas-gas overlaps
% of accepted insertion attempts
% of accepted deletion attempts
54.7 18.2
0.70 1.37
0.62 0.82
0.61 0.81
54.7
1.31
1.21
1.20
18.2
2.25
2.22
2.20
loading. Furthermore, we see that the improvement due to rotational biasing is stronger than that due to positional biasing. Table 2 shows the acceptance rates of insertion and deletion attempts for each of the four runs in Figure 12. Overlaps are defined as having an insertion energy contribution of greater than 100ff (and thus little probability of acceptance). While the differences between the insertion and deletion acceptance rates are the same for both methods, we see that increasing the magnitude of these two quantities increases the rate of approach to equilibrium. Thus, by increasing the acceptance rate of insertions and deletions, we ease the ergodic problems inherent in simulation systems with strong angle-dependent interactions. Acknowledgment. This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences. LA970370+
