Langmuir 1993,9, 2682-2692
2682
Sorption in Energetically Heterogeneous Model Silica Systems J. M. D.MacElroy Department of Chemical Engineering, University College, Dublin, Belfield, Dublin 4, Ireland Received October 29, 1992. In Final Form: April 6,19934D Grand canonical ensemble Monte Carlo (GCEMC) simulation results for a Lennard-Jones (12-6) vapor adsorbingin model microporoussilicastructures are presented and discussed. The silicamedia are modeled as randomly distributed interconnected nonoverlapping silica microspheres with the latter treated in two different ways: (i) as structureless particles (smearedadmolecule/silica microsphere interactions) and (ii) as fully structured entities (the sorbed particles interact with the individual atoms of the Si02 matrix). Sorption isotherms are determined for fractional monolayer coverages in the range 10-6 < 8 < 0.5 at three different temperatures and porosities and the influence of solid surface structure is investigated using these data in conjunction with the profiles for the admolecule number density as a function of energy within the micropore volume. These density profiles correspond to the kernel of the Fredholm equation for energetically heterogeneous systems and, coupled with complementary studies of the first few local virial coefficientsof the pore fluid, the simulation results reported here suggest that anumber of assumptions currently employed in the characterization of the energetic heterogeneity of microporous materials may require review. 1. Introduction
Since the early work of Langmuir’ the influence of sorbent heterogeneity as well as its characterization has been the subject of a rapidly growing volume of literature, reviews of which may be found in the recent texts by Jaroniec and Madey2 and Rudzinski and Everett.3 Sorbent heterogeneity may physically arise for avariety of reasons; however, regardless of its origin it is now recognized that two fundamental components are needed in order to develop a predictive theory for vapor or liquid sorption in heterogeneous media. The first of these involves the concept of the sorbate/sorbent potential energy distribution function which must be known irrespective of the presence or absence of admoleculeJadmoleculeinteractions while the second relates to the topography of the potential field which requires careful consideration to accurately quantify molecular interactions in the sorbed film. Current theoretical models of adsorption on heterogeneous surfaces are usually based on the assumption that only the distribution of the potential energy minima is needed to predict sorption behavior, and when used in conjunction with an appropriate isotherm equation for a locally homogeneous surface, this approach has enjoyed some measure of success. Although this method of analysis is accurate, in principle, at low temperatures, a serious problem which is generally difficult to overcome is the independent verification of the accuracy of both the local isotherm employed and the energy distribution function obtained from a numerical analysis of sorption equilibrium data. When the intrinsic heterogeneity of the solid surface is further complicated by the structural heterogeneity associated with the configurational properties of the pore network in microporous materials, the difficulties involved in the extraction of reliable information concerning the potential energy distribution from equilibrium sorption data increase dramatically. These difficulties are comAbstract published in Advance ACSAbstracts, August 15,1993. (1) Langmuir, I. J. Am. Chem. SOC.1918,40, 1361. ( 2 ) Jaroniec, M.; Madey, R. Physical Ad8orption on Heterogeneous Solids; Elsevier: Amsterdam, 1988. (3)Rudzinski, W.;Everett, D. H. Adsorption of Gases on Heterogeneous Surfaces;Academic Press: London, 1992. 9
pounded by the lack of a tractable, accurate local sorption isotherm equation for such systems and a frequent reliance on the assumed applicability of Langmuir’s theory1 or on empirical models based on Polanyi’s potential theory? In this paper a nontraditional approach to the problem of sorption in microporous media will be investigated with a view to clarifying some of the underlying assumptions involved in the analysis of sorption equilibrium data for such systems. The method centers on a microscopic analysis using grand canonical ensemble Monte Carlo (GCEMC) computer simulation techniques to determine the equilibrium properties of a Lennard-Jones vapor adsorbing in two model silica media which are structurally heterogeneous over atomic and hyperatomic length scales. One of the most important properties which will be examined here is the energy-dependent local sorption isotherm. This quantity corresponds to the kernel, ( n ( e ) ) , of the Fredholm intergral in eq ICbelow n(z,T) = $Jn(r,z,T) d r = J ( n ( e ) ) d!&) =
J (n(e)) f(4de
(la) (1b) (IC)
where #(e) is the cumulative volume fraction of the medium in which the sorbate/sorbent potential energy is less than or equal to the energy e and f ( e ) is the differential distribution of the sorption energies in the system. The angular brackets in eq 1 denote averaging over all equipotential domains in the medium. In the context of the classical treatment of heterogeneous systems, these domains are either sites or patches on the solid surface and, as noted above, the energy e and the distribution function f ( e ) are usually considered to refer to the potential minima in a given domain. In this paper, however, the distribution of sorptive energies is considered to be characterized by the entire potential energy surface; i.e. the energies e are not restricted to the potential minima within the medium. This approach, which corresponds to the virial description of adsorption, is quite general and it incorporates the classical method cited above as a special (4) Polanyi, M.
Trans. Faraday SOC.1932,28,316.
0743-7463/93/2409-2682$04.00/00 1993 American Chemical Society
Sorption in Model Silica Systems
Langmuir, Vol. 9,No. 10,1993 2683
case. It is also important to note that the quantity (n(e)) appearing in eqs lb,c is not an independent local adsorption isotherm but implicitly includes all possible correlations with neighboring potential energy domains within the system. The influence of the potential energy surface in the neighborhood of r or e is fully accounted for in the local virial coefficients which appear in the z-expansion for the local sorbate number density. The disadvantage of the classical method which assumes e and f ( e ) refers to the potential minima in the system is that, even though it is much easier to obtain comparatively simple mathematical forms for the kernel (n(e)), it will only provide an accurate description of adsorption when the local density n(r)is sharply peaked at the position of the sorbate/sorbent potential minimum. This is strictly true only when le,& kT >> 1.0 where Iemh( is the magnitude of the minimum potential energy. When this criterion is not satisfied and/ or if the occupancy of high energy sites at low coverages in heterogeneous systems complicates the definition of le-(, then the virial approach or an equivalent threedimensional description of the adsorbed fluid should be employed.
2. Models and Simulation Technique The two models of microporous silica investigated in this work are composed of randomly distributed interconnected Si02 microspheres whose size is commensurate with the dimensions of the average elementary silica particle in high surface area silica gel. For silica gel with a specific surface area of 800 m2/gthe estimated diameter of the silica microspheres is 2.76 nm6 and this value was employed in all of the computations reported in this paper. The techniques used to generate the assemblies of interconnected silica particles have been described in refs 5 and 6 and the reader is referred to this earlier work for details. Where the two models of microporous silica differ is in the treatment of the structure of the individual spheres within the assemblies. In the first model it is assumed that the silica particles can be represented as a continuum of interaction sites with each admolecule/site interaction described by the Lennard-Jones (12-6) potential function f$ij(rij)= 4eij(
(
,)12
-
($)
(2)
For a single admolecule A located at a distance rM relative to the center of a sphere S composed of a continuum of sites, it has been shown7 that the smeared interaction potential is given by 21 ~ e A 0 p ~ 3
(rM2- R2)'
\
uA06
(rM2
)
- R2)3
(3)
where UAO and EAO are the admolecule/site interaction parameters, po is the density of sites within a given sphere, and R is the radius of the silica particle. The continuum (6) MacElroy, J. M. D.; Raghavan, K. J. Chem. Phys. 1990,93,2068. (6) MacElroy, J. M. D.; Raghavan, K. J. Chem. SOC.,Faraday Trans. 1981,87,1971. Due to a clerical error, the number of connected silica particles employed in all of the computations reported in this earlier work waa 8 rather than 10, and hence the absolute porosity waa 0.67 instead of 0.69. This does not, however, affect any of the comments or conclusiona expressed therein. (7) Kaminsky, R. D.; Monson, P. A. J. Chem. Phys. 1991,95, 2936.
interaction potential given in eq 3 differs from that employed in earlier work.6 In ref 5 the density of sites within the silica particles was assumed to vary continuously, whereas eq 3 is based on the assumption that the density of sites is uniform at all points within a sphere of radius R beyond which PO drops discontinuously to zero. Both of these continuum models lead to very similarresulta for Henry's law adsorption equilibrium,6J and for ease in computation eq 3 has been employed here. The sites within the Si02 spheres which are of primary concern here are the oxygen atoms since the silicon atoms do not play a significant role in the admolecule/silica dispersion interaction^.^*^ Therefore, based on a bulk vitreous silica density of 2.2 g/cm3,the value of po for the density of oxygen sites to be used in eq 3 is 44 sites/nm3. The smeared silica sphere radius is a fitted parameter determined by comparing the continuum potential function with the potential function for the interaction of an admolecule with a fully structured silica particle. The probe species chosen in ref 7 for this purpose was methane, and a comparison of the positions of the potential minima in both cases provided R = 1.346 nm. For this value of R the continuum interaction potential in eq 3 passes through zero at rM = (us + u ~ ) / 2= 1.56 nm; Le. the effective diameter of the silica particle is us = 2.76 nm and this was verified in this work. Finally, the Lennard-Jones parameters UAO and CAO for the admolecule/oxygen site interactions were estimated using Henry's law sorption data for the methane/silica system10 in conjunction with the Lorentz-Berthelot rules UAO = (UAA + uoo)/2 and CAO = ( e ~ a o o ) l / In ~ . thia paper the admolecules are modeled as structureless particles with Lennard-Jones interaction parameters appropriate to SF6; i.e. U A A = 0.551 nm and e d k = 200.9 K.ll A single oxygen diameter uoo = 0.3 nm is also used here (eq 3 makes no distinction between the bridging and nonbridging oxygens in silica6) and this provides UAO = 0.426 nm. From a Monte Carlo simulation of the methanehilica system, the value of €Ao/k obtained by comparing simulation and experiment was 334 K which is slightly smaller than the value reported in ref 7 due to the interconnectivity of the solid spheres involved in this work. This result provides eoolk = 752 K which, in turn, gives eAo/k = 389 K for the structureless SFdsilicasystem. In the second model investigated in this work the atomic structure of the silica microspheres is explicitly taken into consideration. The modeling techniques involved in this case have been described in detail in refs 5 and 6 and are briefly summarized here. Initially, bulk vitreous silica is simulated via NVT Monte Carlo simulation using the twobody and three-body potential functions for the Si02 system proposed by Feuston and Garofalini.I2 The quenched glass structure generated in this simulation is subsequently used in the creation of model silica spheres which are obtained by randomly cutting spheres of the required dimensions from the bulk solid. The surface of these spheres is then scanned and, if necessary, modified to ensure that all of the silicon atoms in the silica particle are bondedto each of the nearest neighbor oxygens present in the original bulk material. A further condition imposed on the surface structure is the presence only of single and geminal nonbridging oxygen atoms (assumed here to (8)Bezue, A. G.; Kieelev, A. V.; Lopatkin, A. A,; Du, P. Q. J. Chem. SOC.,Faraday Trans. 2 1978, 74,367. (9) Kieelev,A. V.; Lopatkin, A. A.; Shulga, A. A. Zeolites 1985,5,261. (10) Gan~wal,S. K.; Hudpins, R. R.; Silveeton, P. L. Can. J. Chem. Eng. 1979,57,809. (11) Hirechfelder,J. 0.;CurtisS,C. F.;Bird,R.B. TheMokcukar Theory of Gases and Liquids; Wiley: New York, 1964. (12) Feuaton, B. F.; and Garofaliii, 5. H. J . Chem. Phys. 1958,89, 6818.
2684 Langmuir, Vol. 9, No. 10, 1993 represent the isolated and geminal silanol groups known to exist on the surface of real silica materials) and the absence of nonbridging oxygen triplets. The SFdsilica interaction parameters for the structured system were similarly estimated using the LorentzBerthelot mixing rules and the methanehilica Henry's law sorption data reported by Gangwal et al.l0 In this case the bridging and nonbridging oxygens are distinguished by different Lennard-Jones diameters, a00 = 0.27 nm (bridging oxygens) and a00 = 0.3 nm (nonbridging oxygens), while the potential minimum for the methane/ oxygen interaction is assumed to be the same in bothcases. The latter provides eoo/k = 312 K which in turn gives tAO/k = 250.5 K for the SFdoxygen interaction. In the GCEMC simulations, which were carried out using Adams' asymmetric algorithm,13the admolecule/silica and admolecule/admolecule interactions were truncated and ~0 interaction), shifted to zero at R C A=~7 . 5 ~ (structureless RCAO= 2 . 5 ~ ~(structured 0 silica particles), and RCAA= 2 . 5 (SFdSFe ~ ~ interaction). These cut-off radii were also employed in the preliminary methanehilica simulations and therefore the influence of the larger admolecule/silica cut-off radius in the case of the structureless interaction is, at least in part, compensated for in the evaluation of the parameter t ~ o .The fundamental cells employed in the simulations contained 10interconnected silicaparticles with each configuration of microspheres generated by the quenching procedure described in ref 6 (50 configurations were used in the structureless system simulations and 5 (or 1) configurations were employed for the structured system). In the case of the structureless system three different helium porosities $(He) were studied, $(He) = 0.49 (the silica gel used in the experiments reported by Gangwal et al.l0),0.58, and 0.69. The helium porosity is here defined as the fractional volume of the medium into which a hard core helium particle of diameter 0.2193 nm14 can be inserted without overlapping with a silica particle of diameter equal to 2.76 nm and the side lengths of the fundamental cubic cells involved in these simulations were L = 39.67, 42.74, and 47.04, respectively, in units of the silicon-oxygen bond length lsio = 0.162 nm. In the case of the structured silica system only one helium porosity was considered ($(He) = 0.491, where $(He) is now defined as the volume fraction of the medium in which a LennardJones helium particle (UAA = 0.263 nm and t d k = 6.03 Kll) has a potential energy less than or equal to zero (this method was deemed appropriate for the structured system and, as will be shown later, both techniques gave the same SF6 cumulative void fraction for potential energies greater than zero). The fundamental cell length in this case was L = 41.35 in units of lsio = 0.162 nm. For the structureless system, approximately 103GCEMC compound events per SF6 particle were involved in the equilibration stage of a given run followed by 105 events in the equilibrium ensemble (a compound event corresponds to an attempted move and an attempted addition or removal of an admolecule). For 50 silica microsphere configurations the total sampling rate for each p , V (or $(He)), and T under equilibrium conditions was therefore 5 X lo6. In the case of the structured system, the same number of equilibration events was employed and the number of events sampled in the equilibrium ensemble depended on the number of silica particle configurations involved in the simulation. For those simulations employing five silica particle configurations a total of 5 x 106 (13)Adams, D. J. Mol. Phya. 1975, 29, 307. (14)Chapman, S.;Cowling,T. G. The Mathematical Theory of NonUniform Gases, 3rd ed.; Cambridge University Press: London, 1970.
MacElroy
o yr[He] = 0.69
€
-mm
-6m
I I
/ i r
-om
-5m
om
501
mm
E
Figure 1. Cumulative volume fraction as a function of energy for the system composed of structureless microspheres. e is in units of t ~ =0 389k where 12 is Boltzmann's constant. For clarity, only a small fraction of the points computed is shown in this figure. equilibrium events were sampled (one million for each configuration) and for those runs employing only one silica configuration the equilibrium sampling usually involved 2.5 X lo6 events. The properties computed during these simulations included the volume averaged number densities and the internal energies of the adsorbate which were calculated over a range of activities z = exp(p/kT) and temperatures, where p , which is defined here as an "excess" chemical potential, is equal to the true chemical less kT ln(h2/2?rmkT)3/2.Additional quantities which were determined here and which are central to this work are the energy-dependent local density and cumulative volume fraction defined in eq 1. The computation of ( n ( e ) )was carried out as follows. After the execution of each compound event in a given GCEMC run, between ten and one hundred locations within the fundamental cell were randomly selected (the number of locations considered increased with decreasing p ) . The sorbate/sorbent potential energy at each selected position was then calculated and a search was conducted to determine the number of admolecules present within a spherical volume of radius a d 1 6 centered a t the chosen locations. This procedure gave an ensemble of between 5 X 10' and 5 X 108 local number densities which provided reasonable statistics in the low-energy as well as in the high-energy regions of the medium. (A similar procedure was employed in ref 6 in the determination of the radial density profiles of the sorbate, n(rAs), relative to the centers of the individual silica microspheres). The cumulative volume fraction $(e) was similarly determined by a straightforward Monte Carlo insertion technique. 3. Results and Discussion 3.1. Structureless Silica Microspheres. The cumulative volume fraction $(e) is shown in Figure 1 as a function oft for the structureless system at three different porosities. Note that e in this figure is in units of tAo/k = 389 K and that the values t = -14.4, -9.6, and -4.8 correspond to the potential minima for a single admolecule interacting with three, two, and one silica microsphere(s), respectively. The undulations in the function $(e) at these energies are therefore simply related to configurational transitions in the mode of the admolecule/silicainteraction; i.e. in the range -15.0 < t < -10.0 the potential field is located in cavities generated by three silica particles, while potential energies -10.0 < e < -5.0 correspond to locations
Sorption in Model Silica S y s t e m in the vicinity of the minimum potential “ring” at the contact point of two individual silica spheres, and for energies greater than -5.0 the potential field is located at positions near the potential minimum for the isolated spherical surface of the silica particles and extending into the pore volume. It is also of interest to note that experimentally reported “site” energy distribution functions f ( ~ )for silica gel frequently display a bimodal structure and, in at least one case,ls the relative magnitude of the two peak energies is in close agreement with the positions of the steps in $(e) at e = -14.4 and-9.6 observed in Figure 1. Although this comparison can only be considered as qualitative in view of the simplified admolecule/silica smoothed interaction employed here, it is tempting to suggest that the peaks observed in the experimental distribution functions are a result of the particulate (colloidal) structure of silica gel rather than a possible energetic heterogeneity in the oxide and hydroxylated surface of silica itself. Further evidence for this will be provided later in the discussion of the results for the atomistically structured silica system. The influence of the bulk density of the medium as shown in Figure 1 is readily interpreted in terms of the increased openness of the pore structure with increasing helium porosity $(He). In the range e < -14.4 the function + ( E ) at a given energy decreases as $(He) is increased from 0.49 to 0.58 due simply to the widening of the holes generated by triplets of silica particles in the medium. A further increase in the helium porosity to 0.69 however would appear to permit an increase in the effective coordination number of the silica microspheres with the admolecules which is not unreasonable in view of the size of the admolecules. This observation is a simple example of structural heterogeneity a t work and, as will be shown later, what would initially appear to be a comparatively small variation in $(e) with bulk density can have a significant affect on the behavior of the thermodynamic properties of the adsorbate at low coverages. For energies greater than -14.4 the trends in $(e) are as expected with $ ( Edecreasing ) with increasing helium porosity a t a given energy for e < 0 and increasing in parallel with $(He) for e > 0. In Figure 2 the reduced number density (n(e)) is plotted as a function of energy at a number of bulk vapor activities, temperatures, and porosities. For very low and high energies the influence of coarse graining the density within a spherical volume of radius u d 1 6 was found to have a significant effect on the estimated values of (n(e))and for this reason no values are shown in this figure for these energies.18 For intermediate energies however (generally in the range -12.0 < e < -3.0) it was verified that the effects of coarse graining were minor. For each run, regardless of p, T, or $(He) the number density would appear to be independent of e for energies lower than -7.0. As noted above, in the range E < -10.0 the admolecules are primarily located in cavities generated by triplets of silica microspheres and in a “cell” theory of adsorption (e.g. (15) Jaroniec, M.Surf. Sci. 1976,50,553. (16) (n(r)) may also be obtained from a direct computation of d(N(r))/ dr (Bakaev,V. A.; Steele, W. A. Longmuir, 1992,8, 148) which may be readily shown to be equivalent to V(n(r))f(c).Whether this would circumvent the errors associated with the coarse graining procedure employed in the present work at high and low energies, however, remains an open question. The latter group of terma involves the distribution functionf(c)which,inturn,isobtainedbydifferentiating+(c). Theprocess of numericallydifferentiating#(€)atvery low energiesis,however,subject toseriouserror and therefore there is no guaranteethat improvedaccuracy will be obtained using this procedure under these conditions. Similarly, at high energies d(N(c))/+ (or, specifically, d(N(c))/d#(c)) is a small quantity and high accuracies are required to reliably estimate ( n ( r ) )in this energy range particularly for use in plotting eq 4b.
Langmuir, Vol. 9, No. 10,1993 2685 2.50
R
2.25 2.m 1.75
a p = -10.0
1.9
0 /I=
A
h
W
1.25
Ip
V
1.m
-12.0
= -14.0
o u = -16.0
0.75
.. .
-1z.m
-6.m
-9.m
-3.00
E
2.25
(51
o T = 0.95 o T = 1.20 0 T = 1.45
0.25
om
1
-12.m
-6.m
-9.03
-3.m
E
2.50
-
150
A
v
1.25
-
V
1.m 0.75
-
050
-
0.25
-
OM
-12m
I
I
-em
-gm
I -3.m
E
Figure 2. Reduced number density aa a function of energy for the structureless system: (a) $(He)= 0.49, T = 0.95; (b) $(He) = 0.49, p = -10; (c) T = 0.95, p = -10. n(e) is in units of u u a , Tis in units of e d k , p is in units of kT, and t is in units of CAO = 389k.
Langmuir’s classical theory’) one would expect these high energy sites to be completely filled for all of the conditions involved in Figure 2 (for E < -10.0 the lowest value of Ie,i,J/kT involved in this figure is 13.0). Although this is not immediately clear from Figure 2, this is in fact the case. The plateau for each condition shown in this figure arises due to steric exclusion by admolecules occupying the (very few) lowest energy domains in the system coupled with a significant degree of mobility of these particles even at the highest values of l~-l/kT. The latter effect also accounts for the relatively small magnitude of (n(e)) and the comparatively rapid decrease in the plateau density with decreasing p (or z) as shown in Figure 2a.
MacElroy
2686 Langmuir, Vol. 9, No. 10, 1993
0.70
f
-
i o.m
6.m 7.m
6.m
g.m
e
07
------.**
**'
1
t
om rm am um u.m
6m
lem
nm
lem
'AS
Figure 3. Radial density profile for the sorbate relative to the centers of the structureless silica microspheres (#(He)= 0.49, T = 0.95, and p -10). The number density is in units of U M and ~ the radial position r u is in units of lsio = 0.162 nm. All other unita are as in Figure 2.
07
"f u5
An additionalfeature common to each of the plots shown in Figure 2 is the rather slow decrease in ( n ( s ) )with increasingenergy in the range -7.0 < e < -3.0. Of particular interest are the moderately high densities involved in a number of cases for e > -4.8 which, as discussed earlier, is the potential minimum of the least active sorption domain in the system. It is also worth noting that, at ita highest value in these studies, the fractional monolayer coverage is approximately 0.5, and therefore the high densities at high E are not a result of multilayer formation. The evidence for submonolayer conditions is provided in Figure 3 where the local density, n ( r M ) , as a function of position relative to the center of the silica microspheres is illustrated for the lowest temperature and highest activity considered in this work. Within the context of localized or two-dimensional mobile monolayer models of adsorption the rather diffuse character of both (n(e)) and n(rm) is surprising particularly in view of the comparatively high values of the energy barriers for desorption involved here. The general form of the local adsorption isotherm for a variety of models which are currently employed in correlating submonolayer sorption equilibrium data2%can be simply expressed as
where the function g(n) incorporates the effects of admolecule/admoleculeinteractions and whose functional form depends on both the nature of the adsorbed film (localized or mobile) and the topography of the solid surface. The quantity V, is a constant and corresponds to the site (or domain) *volume"in cell theories of localized adsorption or a measure of the molecular volume of the sorbate in theories of mobile sorbed films (e.g. the van der Waals constant b). It is instructive to consider eq 4a in a slightly different form and to replot the results for ( n ( 4) in accordance with
The data shown in Figure 2a,b are replotted in this way in Figure 4, and from this figure it is possible to make the following observations: (i)At high values of z exp(-e/kT) (low E ) the relationship is indeed linear aa implied by eq 4b; however it is quite clear that if this equation is to be accepted, then the
c!.' 0-2
lbl
1 0.2
0.3
n-'
00
0'
a2
z.exp(-a/kT]
03
04
05
08
07
Figure 4. Modified plot of the local sorbate number density for the structureless system as a function of z exp(-JkT), where z is the activity exp(p/kT)in units of UM-? (a) replot of the data shown in Figure 2a; (b) replot of the data shown in Figure 2b. All symbols and units are BB in Figure 2. parameter V, must be considered to be a function of the adsorbate activity and essentially independent of e (V, increases from approximately 1.0 to 10.0, in unita of a h 3 , as the chemical potential drops from -10.0 to -16.0). None of the approximate theories currently in use predicts such a dependence. (ii)At low z exp(-e/kT) (the high-energyregion in Figure 21, the simulation results for the left-hand side of eq 4b also demonstrate a strong dependence on sorbate activity and, to a lesser extent, the temperature of the system. In this regard it is believed that the influence of admolecule/ admolecule interactions as predicted by the term g(n) is at least qualitatively correct. Typically this term is predicted to decrease rapidly with increasing sorbate density and decreasing temperature and in the opposite limit, z 0, g(n) will approach 1.0. Also, since g(n) dominates the right-hand side of eq 4b at high energies, the results shown in Figure 4 suggest that this term approaches a constant value at a given z as d k T increases. Within the context of approximate theories for submonolayer adsorption, this would imply the existence of a random topography of potential domains within the model silica medium. A deeper understanding of the behavior of the systems under consideration here is also possible by examining the leading terms in the virial expansion of the adsorbing vapor. Based on the approach proposed by Kirkwood and
-
Sorption in Model Silica System
Langmuir, Vol. 9, No.10, 1993 2687
Salsburg,l' Steele and Ross18 have shown that the local density of a pure fluid in an external field can be expressed as
dr, dr, where
1
+ ...
(5)
c/kT
h d u,(ri) is synonymous with the energy e. The single integral in the linear z term inside the chain brackets of eq 5 has been computed here via Monte Carlo integration for approximately 70 different positions, r1, within each of the 50 fundamental cells of the silica microspheres employed in the GCEMC simulations. These positions were selected at random from a sequence of equally spaced energy intervals spanning the range of energies from the minimum energy shown in Figure 1up to 0.0. At each position Monte Carlo integration of the expression
0'3
L
bl
C(e)
C(r,) =
svf12 9) exp( -
dr,
=sv(exp(-w)-l)
involved 5 X 104 uniformly selected points within a spherical volume of radius equal to the cut-off radius R, for admolecule/admoleculeinteractions. The results for the lower temperature, T = 0.95eAA/k, and +(He) = 0.49 are provided in Figure 5 and results not shown here for the other temperatures were qualitatively similar. The rather broad distribution in the values obtained for this local third virial coefficient at a given energy are a clear indication of the degree of structural heterogeneity involved in the conceptually simple model of porous silica under consideration here. Also note the change in sign in the average value of C(e) at c = -9.6, the transition energy above which the potential "rings" at the contact points of two individual silica particles become the dominant sorption domains in the system. Above e = -9.6 the positive value of C(c) implies the dominance of admolecule/ admolecule attractive interactions, while at very low energies the repulsion interactions between particles confined within the holes in the medium are most important. The latter observation is consistent with the exclusion effect involved in Langmuir's model of localized adsorption which was rederived in an elegant manner by Steele and Ross18 using eq 5. They have shown that by assuming (a) infinite repulsion between particles located on the same site, (b)zero interaction between particles on different sites, and (c) the sites are, effectively, deep, narrow potential wells with steep repulsive walls, then eq (17)Kirkwood, J. G.;Salsburg, 2.Discuss. Faraday SOC.1953,15,28. (18)Staele, W.A.; Ross, M. J. Chem. Phys. 1960,33,484.
-3zm
-zm
-zim
-mm
-Em
-m
~m
-4m
om
c/kT
Figure 5. Local third vir% coefficient for the Lennard-Jones adsorbate as a function of e/kT for the structureless system with C(c) in units of U A A(T ~ = O . Q S r d k ) : (a) region of negative values of C(e); (b) region of positive values of C(c). The solid line is the average value of C(c) at a given energy. 5 may be written as
(zC(r,,=-O))'
+ ...I (7)
where C(r1,m-O) is the local third virial coefficient for an infinitely repulsive admolecule/admolecule interaction which is assumed to exist within an exclusion volume Vex i.e.
Assumption c above leads to the result that the fourth localvirialcoefficientD(r1,m-O)isequalto P(r1,m-O)with similar relationships for the fifth and higher virial coefficients, and under these conditions eq 7 can be rearranged to give
By considering C(rl,=-O)= -V, exp(-u,(rl)/kT) = (C(s,=0)) where V, is the effective cell "volume" then eq 9 is equivalent to eq 4 with g(n) = 1.0. The third and fourth virial coefficients C(rl,m+) and D(rl,--O) have also been evaluated here by Monte Carlo integration in a manner similar to that described above
2688 Langmuir, Vol. 9,No. 10,1993 1
0'4
:i
ij -3zm
-28m
-z+m
-ma
-Em
-am
am
-rm
om
c/kT
-
MacElroy
-9.6, this virial coefficientis in good agreement with the third virial for the hard core potential C(e,m-O) in both sign and magnitude at low energies, and this supports the possible applicability of eq 9 to low coverage adsorption within the cavities in the medium. To test this, GCEMC simulations were conducted for a range of conditions at very low adsorbate activity. The pVT ensemble results obtained for a number of properties from these simulations and from the simulationsassociated with the density profiles reported above are provided in Table I. The results obtained for ( U M ) and ( u b )suggest that over most of the range of surface coveragethe potential energy of the adsorbing vapor is nearly entirely determined by the admoleculelsilicainteractions. Also note that even at the lowest coverages studied the specific interaction of the admolecules with the solid is much higher than the value associated with the minimum potential shown in Figure 1for the low energy holes (e =-14.4 or, for example, at T = 0.95edk, dkT = -30.0 which should be compared with the lowest value of -21.2 given in Table I for ( u b ) ) . This implies that even at chemical potentials as low as -26.0 the positions of the admolecules in the potential field of the solid are heavily weighted toward higher energies even though one would expect their most probable locations to be the low energy cavities in the medium. Similarly, the abrupt change in the admoleculelsilica potential energy with increasing p in the range -16 C p C -12 is related to filling of the potential rings in the vicinity of the contact points of the individual silica particles. These observations are supported by the trends in the sorption isotherms shown in Figures 7 and 8. The solid lines in both these figures represent the Henry's law sorption isotherms obtained directly from the results for $(e) using the limiting form of eq 5 in eq 1
Figure 6. Local virial coefficients for the hard-core adsorbate as a functionof elkTfor the strudureleeasystem (T = 0.96edk): (a) the local third virial coefficient C(e,--O) in units of uha;(b) the local fourth virial coefficient D(e,-+) in units of ah6. The solid lines in both cases are the average values at a given energy.
for C(rl). In these computations the exclusion volume Vex has been assumed to correspond to the region of repulsive interaction of two of the Lennard-Jonesparticles, i.e. Vex = (4/3)?rau3. Furthermore, as a partial check on the adequacy of assumption c, the fourth local virial coefficient in the form D(r,,--O)
D(e,--O) =S V (u,(rl
J L
+ r') + u&r1+ r' + r"))) dr" dr'
(10)
kT has been explicitly determined. In the evaluation of this double integral, 500 points were sampled in the Monte Carlo integration of each integral and, although this number of points may appear to be quite small, the integration volume in this case is also much smaller than the volume involved in the numerical integration for C(rl). However, to minimize the statistical error in the case of C(rl,m-O), the same sampling rate (5 X 104) as that employed for C(r1) was also used here. The results obtained for C(c,--O) and D(c,m-O) are plotted in Figure 6. The important features displayed here are (i) a distinct stepwise variation in the most probable values of both coefficients with the steps occurring at the transition energies noted above and (ii) to a very good approximation D(e,-4) = C2(e,=,0). Furthermore, although there is a change in sign for C(e) at e
where the lower and upper limits in the second integral are --03 and +m, respectively. It is clear from Figures 7 and 8 that Henry's law in the structureless silica system will only be observed at extremely low activities. Note also the near equality of the Henry's law coefficients for helium porosities of 0.58 and 0.69. This arises due to the apparent increase in the average coordination number of silicamicrospheressurroundinga given admoleculeat very low energies as discussed earlier with reference to the results for $(He) = 0.69 in Figure 1which is counterbalanced by the lower values of $(e) at higher energies for this porosity. The dotted lines in Figure 7 correspond to the predictions of the approximate approach described above using eq 9 with C(r1,m-O) replaced by (C(e)), i.e. the data represented by the solid line shown in Figure 5. Although the trends predicted agree qualitatively with the exact GCEMC simulation results, it is clear that the assumption that the fourth and higher order local virial coefficients can be obtained from a knowledge of C(r1) alone is only reasonable for the structureless silica system at high temperatures and/or at very low coverages where adsorption is dominated by occupancy of the low energy holes. The GCEMC results shown in Figurea 7 and 8for densities in the range 1o-S C n(z)C 1o-S correspond to the uptake in the potential rings at the contact points of the silica particles and for densitieshigher than 109the admolecules are primarily accumulating on the isolated surface of the microspheres.
Sorption in Model Silica Systems
Langmuir, Vol. 9, No. 10,1993 2689 Table I
(a) Equilibrium Properties from GCEMC Simulations for the Structureless System (+(He) = 0.49) !P = 0.95 !P = 1.20 !P = 1.45 n(z)'(lV) 0.0181 0.0562 0.164~ 0.4% 1.74 8.22 30.76 75.09 1341
fib
-26 -24 -22 -20 -18 -16 -14 -12 -10 -8
( U A ~ I ) ~n(z)'(IV)
( U M ) b (lo-')
+.Oh
-21.25 -19.96 -18.75 -17.64 -16.83 -16.23 -15.33 -13.92 -12.42
-0.092
-0.285 -0.677 -2.123 -7.73 -25.28 -572 -962
( U A A ) b (lo-')
(UAnI)'
-0.0028 -0.021 -0).1l& -0.251 -0.763 -2.979 -12.53 -36.85 -71.48
-17.76 -16.54 -15.24 -14.23 -13.32 -12.82 -12.21 -11.21 -9.858
0.00392 0.0171 0.0622 0.2258 0.812 3.929 19.44 62.1e 1281
n(z)'(lV)
-6
( U M ) b (lo-')
(UAm)b
0.00192 0.0104a 0.0472 0.196
-0.022 -0.011 -0.0403 4.161
-14.41 -13.53 -12.43
0.802
-0.542
4.288 22.54 73.71 1481
-2.457 -11.63 -34.67 -661
-10.72 -10.32 -9.92 -8.899 -7.745
-11.53
(b)Equilibrium Properties from GCEMC Simulations for the Structureless System (T= 0.95) +(He) = 0.69
+(He) = 0.58
n(z)(le) 0.00724 0.0221 0.0664 0.1878
Ir
-26 -24 -22 -20 -18 -16 -14 -12 -10 -8
(UM)
( U b )
-0.108 -0.063 -0.164 -0,389 -1.01 -3.12 -11.76 -341 -732
0.562
2.217 10.72 37.07 941
-20.95 -19.24 -18.14
-17.04 -15.93 -15.03 -14.42 -13.22 -11.51
n(z)(1V) 0.00242 0.00794 0.0241 0.0672 0.190~ 0.632 3.008 14.53 50.31
(lo-') -0.11 -0.022 -0.063 -0.298 -0.488 -1.31 -4.75 -181 -471
(UM)
11%
(Uh)
-20.25 -18.95 -17.6, -16.48 -15.13 -13.82 -13.12 -12.32 -10.91 -9.44
-902
temperature Tis in units of t d k = 200.9 K. b The excess chemical potential, p , and the energies per particle, ( U M ) and ( u b ) ,are in units of kT. The volume averaged number density n(z)and the activity z = exp(p/kT) are in units of ~;9. a The
0
0
1 D T = 0.95 (KH = 9.50*107]
i
e
o
0
0
0
D
e
2 o T = 1.20 [KH = 3.53*103) 3 o T = 1.45 [KH = 1.14*104)
*-u
u-t2
-,11
*-a
,-9
a-,
-7,
*.6
,.5
,-4
3.,
e
0-2
2
Figure 7. Sorption isotherms for the structureless system ($-
(He) = 0.49): -, Henry's law; -., eqs 1and 9 with C(rl,m-O) in eq 9 replaced by (C(c)) from Figure 5. Both n(z) and z are in units of U M ~ .
3.2. Structured Silica Microspheres. The cumulative volume fraction as a function of energy (nowin units of ~ A Ofor the structured system) is illustrated in Figure 9. For comparative purposes the corresponding function for the structureless system is also shown in this figure and the differences between both systems are primarily due to the different cut-off radii for sorbatelsorbent interactions involved in both cases. The trends in $(e) are nonetheless quite similar with the exception that the heterogeneity of the atomistically structured surface appears to partially smooth out the undulations associated with the arrangements of the silica microspheres themselves. The two seta of results for the structured system also illustrate that the configurational sampling of the silica particles is only important at very low energies. The open circles in Figure 9 correspond to a single configuration
,-E
-,11
o.a
4 ,
e
o
1 o yr[He] = 0.49 [KH = 9.50*107] 7 2 o y-[He] = 0.58 [K, = 3.05'10 ] 7 3 o yr(He] = 0.69 [KH = 2.40'10 ] ud
*.7
,-E
u-5
u4
)-3
2
Figure 8. Sorption isotherms for the structureless system a t Henry's law. different helium porosities (T= 0.95cJk): -, All units are as in Figure 7.
of ten structured silica particles, whereas the open squares represent the data for five configurations. The local densities obtained at a number of activities and at the lowest temperature investigated here are shown in Figure 10a as a function of energy. For the structured system the effectsof coarse graining the local densitywithin a spherical volume of radius u d 1 6 were found to be much more severe than in the case of the structureless system and reliable estimates of (n(e)) could only be obtained over a significantly narrower range of energies (-14.0 e < -4.0). The reason for this is clearly demonstrated by the very high values of ( d e ) ) shown in Figure 10a at low energies. These high local densities imply the existence of strongly localized monolayer sorption conditions in contrast to the high degree of mobility inferred from the results for the structureless system. Although strong localization of the admoleculeslimits the range of energies in which ( n(e)) may be obtained with reasonable accuracy,
MacElroy
2690 Langmuir, Vol. 9, No.10, 1993
-* Y
I
i
-am
-25m
-2om
-Ism
-am
-5.m
om
5m
om
E
Figure 9. Cumulative volume fraction as a function of energy for the system composed of structured microspheres. e is in units of CAO = 250.5k. The open circles and squares correspond to the results for a single configuration and five configurations of the silica particles, respectively. ,
h
53.m 27.m
-+
I
nm
A
6m
V
12.m 9.m
I
-am
6.m
-m
I
1
-m
-am
I
-6m
I
om
Ek/T
3.m 0.m
-Rm
-um -am
-1100
-om -am .am -7m .am -5m -rm E
I
m3
-nm .am om e/kT Figure 11. Local virial coefficientsaa a function of c/kT for the structured system (2' = 0.95edk): (a) the local third virial coefficientC(e) for the Lennard-Jonesadsorbatein units of UM~; (b) the local third virial coefficient C(e,=-O) for the hard-core adsorbate in units of U M ~ ; (c) the local fourth virial coefficient i?(e,m+) for the hard-core adsorbate in units of UM*. The solid lines represent the average values at a given energy.
-am
0-2 u.3
u.2
0-1
4
z.expy-c/kT]
d
d
o'
d
Figure 10. Reduced number density aa a function of energy for the structured system ($(He) = 0.49, T = 0.95): (a) n(e) aa a function of e; (b) the data in part a replotted in accord with eq 4b. All units are aa in Figures 2 and 4 with the exception that e is in units of CAO = 250.5k. it is still possible to glean important information from these results, in particular when these data are replotted in accord with eq 4b as shown in Figure lob. At high energies (low z exp(elk27) the results shown in Figure 10b display trends similar to those for the structureless system with the exception that a distinct energy dependence of the term g(n) in eq 4 would appear to be involved in the present case. For example, if the data for
-2m
-am
p = -10 alone are considered, then the low-energy results suggest that the parameter V, in eq 4 is approximately equal to O.O~IJM~. If it is assumed that this quantity is independent of T and e (an assumption frequently employed in the literature), then g(n) is the dominant term in eq 4b for z exp(-c/k27 < 10. Within the context of the approximate theories upon which eq 4 is based, g(n) can only depend on c if a patchwise model of adsorption is considered to apply, in which case the density n in g(n) is the local density (n(c)). Unfortunately, it is difficult
Sorption in Model Silica Systems
Langmuir, Vol. 9, No. 10, 1993 2691 Table I1
(a) Equilibrium Properties from GCEMC Simulations for the Structured System0 (T= 0.95, +(He) = 0.49, and Five Configurationsof Silica Microspheres) n(z) (10-9 0.1819 0.734 2.71 9.01 271 682 1402
P
-22 -20 -18 -16 -14 -12 -10
( U h )
(UM.)
-0.022 -0.42 -3.77 -132 -332 -712 -1252
-25.14 -23.68 -21.58 -19.72 -17.82 -16.h -13.91
(b) Equilibrium Properties from GCEMC Simulations for the Structured Systema(+(He)= 0.49 and One Configurationof Silica Microspheres)
-22 -20 -18 -16 -14 -12 -10 -8 -6 4
O.O8€J4 0.492 2.91 9.83 27.6~ 672 1403
-0,078 -0.21 -3.79 -132 -34, -733 -1293
-22.26 -21.98 -21.53 -20.09 -18.02 -16.h -13.92
0.00906 0.062~ 0.401 2.407 11.13 39.96 10Eil
-0,067 -0.021 -0.528 -2.75 -121 -332 -768
-16.42 -16.12 -15.72 -15.43 -14.21 -12.72 -11.01
0.00442 0.0291 0.2269 1.613 9-22 40.5 1232
-0.038 -0.021 -0.122 -1.31 -6.74 -251 -6%
-11.92 -11.4 -11.61 -11.51 -10.91 -9.76 -8.218
All units are as in Table I.
to reconcile the low values of the left-hand side of eq 4b with the low values of ( n ( e ) )illustrated in Figure 10a a t these energies and an explanation other than or in addition to an energy dependence of g(n) needs to be considered. This is provided by the low-energy results for the righthand side of eq 4b shown in Figure lob. As noted above, the parameter V , for p = -10 is approximately O . O ~ U A A and ~ , from the low-energy results in Figure lob, it would initially appear that V, decreases with decreasing p in contrast to the reverse trend observed for the structureless system. However, another and more appropriate explanation for the trend in V , is with regard to a possible dependence on elkT. Such a dependence was originally proposed for localized adsorption by Hill19 who suggested that the parameter V , is approximately proportional to (kTlle1)3/2forIeJlkT>>1. The results shown in Figure 10b qualitatively support this relationship, and therefore at low activities it is possible that the Langmuir form of eq 4 @(n) 1.0) withan energy- and temperaturedependent V,may be suitable as a kernel for the Fredholm integral in eq 1. To explore this possibility further, the local third virial coefficient defiied in eq 6 for the Lennard-Jones adsorbate and the local third and fourth virial coefficients for hardsphere admolecule/admolecule interactions ( Vex = (47d 3 ) u ~ ~defined 3) in eqs 8 and 10were determined via Monte Carlo integration. In this case approximately 180positions, r1, within each of the five fundamental cells were considered in the analysis. In the single integration for the third virial coefficients 2.5 X 104 points were uniformly selectedand in the double integration involved in D(r1,m0) 500 points were employed in the evaluation of each integral. The results obtained for T = 0.95edk are shown in Figure 11. In comparing these results with the virial coefficients for the structureless system, the most notable differences involved are (i) the sign on C(e) (i.e. C(r1)) does not change at elkT = -20.0 as it did for the structureless system although there does appear to be a change in the energy dependence of C(e) at this point and (ii) at low energies both C(4 and C(e,--O) differ by 1to 2 orders of magnitude. Quite clearly admoleculeladmol-
-
(19) Hill, T.L.J. Chem. Phys. 1949, 17, 762.
ecule attractive interactions are much more important in the structured system a t low energies and repulsive interactions play a more significant role at higher energies. The disparity observed in the low energy local third virial coefficients for the Lennard-Jones and hard core interactions raises questions concerning the applicability of the Langmuir model (or, equivalently, eq 9) to the structured system even though the condition D(e,-O) = Q(e,m-O) would again appear to be supported by the resulta shown in parts b and c of Figure 11. As a qualitative test of this, a series of GCEMC simulations at lower activities and a t two additional temperatures were conducted using a single configuration of ten structured silica microspheres (the $(e) relationship for which is provided in Figure 9). The results obtained from these simulations are presented in Table 11. Again, as in the case of the structureless system, the internal energy of the adsorbate is primarily determined by the interactions with the silica particles. ( u b )also reflects a bias in the average locations of the admolecules toward the higher potential energy domains, although the bias is not as severe in the present case (e.g. for the low temperature system involving five configurations of silica particles the minimum value of e recorded in Figure 9 corresponds to e/kT = -29.0 which is not much lower than the value -25.1 reported in Table 11). The sorption isotherms shown in Figure 12 also suggest that the structured system behaves more “ideally” than the structureless system under similar conditions. The GCEMC simulation results in this case approach the Henry’s law limit more closely at low activities and the predictions provided by eq 9 with C(r1,m-O)replaced by C(r1) are quantitatively better over a wider range of conditions. Significant disparities do exist however which point to the inapplicability of eq 9 and hence to a failure of the assumption D(r1) = C2(rl), E(rd = C3(r1), etc. implicit in Langmuir’s theory of adsorption. For this reason the apparent consistency noted earlier between eq 4 and the results shown in Figure 10bat low energies should be considered as fortuitous. In view of the large differences observed between the results for C(e) and C(e,--O) the reason for the failure of eq 9 is in part due to admolecule/ admolecule attractive interactions even at very low cov-
2692 Langmuir, Vol. 9, No. 10, 1993
MacElroy
- n
0-0
9-,,
a,
,.7
,-6
5.,
*-4
,-3
u-2
2
Figure 12. Sorption isotherms for the structuredsystem ($(He) = 0.49): -, Henry's law; -, eqs 1 and 9 with C(rl,=-O) in eq
9 replaced by ( C(c) ) from Figurella The filledcirclescorrespond to the data for five configurationsof silica particles and the open symbols are the data for a single configuration of silica microspheres. Both n ( ~and ) z are in units of U U ~ .
erages. However, it is also clear from a comparison of the isotherms shown in Figure 12 for the simulations involving five configurations and a single configuration of silica microspheres that the structural heterogeneity of the medium itself (which is greater for the data involving five configurationsas illustrated in Figure 9) plays an important role. This observation is in qualitative agreement with an earlier study by Rudzinski20 in which it was demonstrated that the higher the virial considered, the greater is its sensitivity to the degree of heterogeneity of the medium. 4. Summary and Conclusions In this paper adsorption of a simple vapor in two models of microporous silica has been investigated. One of the
models considered was the analogue of systemswhich have been traditionally assumed to exhibit mobile adsorption (structureless surface) and the second model incorporated the effects of localization on surface sites (atomistically structured surface). In bothcases efforts were undertaken to qualitatively compare Monte Carlo simulation results with the predictions of approximate theories for submonolayer adsorption with only minor success. For the structureless system it was observed that the parameter V, (the inverse of which is usually referred to as the Langmuir constant) varied rapidly with sorbate activity. Approximate theories for submonolayer adsorption were found to be at odds with this observation, although they appear to be qualitatively correct in their predictions of the influence of admoleculeladmoleculeinteractions on the behavior of the sorbed vapor. In the case of the structured silica medium the agreement was more satisfactory in that V, was found to vary with sorbatelsorbent potential energy in a manner consistent with the original ideas of Hill.19 However, a deeper analysis based on a study of the first three local virial coefficients of the pore fluid pointed to weaknesses in this approach also. The results of this study demonstrate the need for significant improvements in the theory of inhomogeneous fluids before the characterization of heterogeneous microporous sorbents can be addressed in a truly quantitative manner. An approach based on the virial expansion, while formally exact, is limited due to the large number of terms in the expansion which are required to adequately correlate sorption data even at very low coverages. It is believed, however, that this or an equivalent three-dimensional description is currently the only viable route toward the development of an unambiguous theory of sorption in heterogeneous microporous media. (20) Rudzinski,
W.Phys. Lett. 1972,42A,519.
