Langmuir 1996, 12, 5425-5432
5425
Monte Carlo Simulation Studies of Heats of Adsorption in Heterogeneous Solids T. Vuong and P. A. Monson* Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003 Received April 4, 1996. In Final Form: July 26, 1996X We present calculations of heats of adsorption for molecular models of adsorption in heterogeneous porous solids via grand canonical Monte Carlo simulations. The molecular models treat the adsorbent as a matrix of particles arranged in a predefined microstructure, and the parameters are chosen to represent the adsorption of simple molecules in silica gel. Three different adsorbent microstructures have been considered, allowing us to study the role of translational order and connectivity among the matrix particles in determining the heat of adsorption. We have calculated separately the contributions from the adsorbatesolid and adsorbate-adsorbate interactions to the heat of adsorption. A comparison with experimental data for methane in silica gel is made. During the course of the work we have revisited the derivation of expressions used to determine heats of adsorption and we present a new derivation which is free from some unnecessary assumptions used previously. In addition we test our results for thermodynamic consistency, and the accuracy of different methods for obtaining the isosteric heat from the Monte Carlo simulations is also investigated.
I. Introduction Heats of adsorption are often considered as an indicator of adsorbent heterogeneity in adsorption,1 since they are more sensitive to the adsorbent microstructure than the adsorption isotherm itself. Extensive data for isosteric heats for adsorption in heterogeneous adsorbents are available from experiment,2 but the precise relationship between the observed data and the underlying molecular interactions remains quite poorly understood. Only recently have theoretical calculations for realistic models of the adsorbent microstructure which could help elucidate this issue been made.3-6 The present work focuses on calculations of heats of adsorption for a new class of molecular model of heterogeneous solids which have recently been developed in this group7-10 and by others.6,11-16 These models can provide a realistic three-dimensional treatment of both structural and energetic heterogeneity in the adsorbent. In the model we have developed7 the void space has a complex geometry comparable to that in real materials and the accompanying adsorbate-adsorbent potential energy field is correspondingly nonuniform. Nevertheless the models are so con* Author to whom correspondence should be addressed. X Abstract published in Advance ACS Abstracts, October 1, 1996. (1) Ruthven, D. M. Principles of Adsorption and Adsorption Processes; Wiley Interscience: New York, 1982. (2) Valenzuela, D. P.; Myers, A. L. Adsorption Equilibrium Data Handbook; Prentice-Hall: Englewood Cliffs, NJ, 1989. (3) Bakaev, V.; Steele, W. A. Langmuir 1992, 8, 148. (4) Karavias, F.; Myers, A. L. Langmuir 1991, 7, 3118. Karavias, F.; Myers, A. L. Mol. Simul. 1991, 8, 23. Karavias, F.; Myers, A. L. Mol. Simul. 1991, 8, 51. (5) Woods, G. B.; Panagiotopoulos, A. Z.; Rowlinson, J. S. Mol. Phys. 1988, 63, 49. (6) Segarra, E. A.; Glandt, E. D. Chem. Eng. Sci. 1994, 49, 2953. (7) Kaminsky, R. D.; Monson, P. A. J. Chem. Phys. 1991, 95, 2936. (8) Kaminsky, R. D.; Monson, P. A. In Fundamentals of Adsorption; Suzuki, M., Ed.; Kodansha Publishing Co.: Japan, 1993; p 309. (9) Kaminsky, R. D.; Monson, P. A. Langmuir 1994, 10, 530. (10) Kaminsky, R. D.; Monson, P. A. Chem. Eng. Sci. 1994, 49, 2967. (11) Madden, W. G.; Glandt, E. D. J. Stat. Phys. 1988, 51, 537. (12) Given, J. A.; Stell, G. J. Chem. Phys. 1992, 97, 4573. (13) Fanti, L. A.; Glandt, E. D. J. Colloid Interface Sci. 1990, 135, 385. Fanti, L. A.; Glandt, E. D. J. Colloid Interface Sci. 1990, 135, 396. Fanti, L. A.; Glandt, E. D.; Madden, W. G. J. Chem. Phys. 1990, 93, 5945. (14) Bakaev, V. Surf. Sci. 1988, 198, 571. (15) MacElroy, J. M. D.; Raghavan, K. J. Chem. Phys. 1990, 93, 2068. (16) MacElroy, J. M. D. Langmuir 1993, 9, 2682.
S0743-7463(96)00325-3 CCC: $12.00
structed as to be amenable to theoretical study using methods adapted from liquid state statistical mechanics11,12,17-19 as well as by computer simulation. In previous work we have examined the influence of heterogeneous microstructure upon adsorption isotherms7 and selective adsorption9 within the context of an application to adsorption of simple molecules in silica gel. In this paper we focus on the heats of adsorption.20 We have studied the relative importance of the adsorbatesolid and adsorbate-adsorbate molecular interactions in determining heats of adsorption as well as how these contributions change with the adsorbent microstructure. We have also looked at the accuracy of different methods for obtaining the isosteric heat from the Monte Carlo simulations. As part of our work we have revisited the classical thermodynamic treatment of adsorption in systems like the present ones. We present new derivations of the equations for heats of adsorption from computer simulations, both directly and via Clapeyron relationships. Later we use this analysis to demonstrate the thermodynamic consistency of our computer simulation results. II. Molecular Model and Computer Simulation Techniques The model used in this work is one we have used in a series of studies aimed at elucidating the role of heterogeneous solid microstructure on adsorption equilibrium.7-10 We treat the adsorbing solid as a matrix of particles (initially spheres have been considered, but other shapes could also be treated) in some specified structure.7 To develop the potential, we treat each matrix sphere as a continuum of interaction sites. The potential provides a level of approximation similar to that given by the 9-3 (17) Vega, C.; Kaminsky, R. D.; Monson, P. A. J. Chem. Phys. 1993, 99, 3003. (18) Lomba, E.; Given, J. A.; Stell, G.; Weis, J-J.; Levesque, D. Phys. Rev. E 1993, 48, 233. (19) Rosinberg, M. L.; Tarjus, G.; Stell, G. J. Chem. Phys. 1994, 100, 5172. Kierlik, E.; Monson, P. A.; Rosinberg, M. L.; Tarjus, G. J. Chem. Phys. 1995, 103, 4256. (20) A preliminary account of part of this work has been published in: Vuong, T.; Monson, P. A. In Fundamentals of Adsorption; LeVan, M. D., Ed.; Kluwer: Dordrecht, The Netherlands, 1996; p 1009. (21) Nicholson, D.; Parsonage, N. G. Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: London, 1982.
© 1996 American Chemical Society
5426 Langmuir, Vol. 12, No. 22, 1996
Vuong and Monson
potential22 used for modeling interactions with plane surfaces and reduces to the 9-3 potential in the limit where the matrix to fluid particle size ratio becomes very large. If each matrix particle is modeled as a continuum of Lennard-Jones 12-6 interaction centers, the potential energy between a matrix particle and a fluid molecule is given by
16πgsFsR3 × 3 21 1 d6 + d4R2 + 3d2R4 + R6 σ12 gs σ6gs 5 3 (1) (d2 - R2)9 (d2 - R2)3
ucs(d) )
[(
)
]
where d is the distance from the center of the fluid molecule to the center of the matrix particle, Fs is the density of interaction sites in the matrix particles, R is the matrix particle radius, and σgs and gs are the collision diameter and well depth for the 12-6 potential between the fluid molecule and a matrix particle interaction center. This potential can be regarded as an approximation to a more detailed intermolecular potential developed by MacElroy and Raghavan15 which incorporates the atomic structure of the silica particles making up the matrix. For the interactions between the fluid molecules a Lennard-Jones 12-6 potential truncated at 2.5 collision diameters was used. The matrix-fluid potential was truncated at 4.0σgs beyond the adsorbent surface, as in our previous work. The parameters used correspond to a model of methane adsorbed in a silica xerogel with a surface area of 803.5 m2/g. The parameter values used are those originally used by Kaminsky and Monson:7 R ) 1.3465 nm; σgs ) 0.33 nm; gs ) 339 K; Fs ) 44 nm-3; σgg ) 0.3817 nm; gg ) 148.2 K. The adsorbent microstructure consists of a configuration of 32 matrix particles in periodic boundaries with a volume fraction η ) 0.386. Three types of microstructure have been considered in this work. For most of our studies we have used a configuration from a Monte Carlo simulation of an equilibrium hard sphere (EHS) system. The second type of structure used was an fcc lattice. Comparison between results for these two microstructures reveals the effect of translational order among the matrix particles. The third microstructure considered was generated using configurations from a Monte Carlo simulation of the adhesive hard sphere (AHS) model23 using the method suggested by Seaton and Glandt.24 The surface attraction between the particles in the AHS model creates a connectivity between the matrix particles which is absent in the EHS model at densities below close packing. This connectivity should act to increase the energetic heterogeneity in the system because of the enhanced matrixfluid potential energy field in the neighborhood of the contact points between the matrix particles. We have used the standard grand canonical Monte Carlo method21,25 for our simulation work. To study an adsorption isotherm, we started with an empty matrix and then carried out a series of simulations with successively increasing values of the activity. Each of these simulations was started from the configuration at the end of the previous simulation at a lower activity. Simulations at the lowest adsorbate densities were run for a total of 4 × 106 configurations with 2 × 106 configurations for equili(22) Steele, W. A. The Interaction of gases with solid surfaces; Pergamon: Oxford, 1974. (23) Baxter, R. J. J. Chem. Phys. 1968, 49, 2770. (24) Seaton, N. A.; Glandt, E. D. J. Chem. Phys. 1987, 87, 1785. (25) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987.
bration. Simulations at the highest adsorbate densities were run for a total of 12 × 106 configurations with 6 × 106 configurations for equilibration. A configuration consists of an attempted translation of a molecule followed by an attempted creation or destruction (chosen with equal probability) of a molecule. For the EHS structure, a single configuration of a 32 particle matrix is generally sufficient to obtain a good enough sample of the statistical geometry.7 This is not the case for the AHS structure. In this case all our data are obtained by averaging the results from simulations using eight different configurations from a Monte Carlo simulation of the AHS model using the algorithm of Seaton and Glandt.24 For simplicity in our implementation of the Seaton-Glandt algorithm we restricted the particle moves to include only transitions between particle energy states 0, 1, and 2. Although this means that our configurations are not strictly those of an equilibrium AHS system, the differences are expected to be small, and this implementation is sufficient to generate substantial connectivity between the matrix particles. We have calculated the bulk properties using the equation of state of the 12-6 fluid given by Johnson et al.26 with a correction for the effect of the truncation of the 12-6 potential which we have used. In presenting our results we will primarily use reduced units defined in terms of the collision diameter, σgg, and well depth, gg, of the fluid-fluid Lennard-Jones 12-6 potential. We have for the temperature, density, pressure, activity, configurational energy, and heats of adsorption T* ) kT/gg, F* ) Fσ3gg, P* ) Pσ3gg/gg, ζ* ) ζσ3gg, U* ) U/Ngg, and q* ) q/gg. In making comparisons with experiment we use real units. III. Thermodynamics In this section we present some results from classical thermodynamics which are relevant to our calculations. Recently there has been progress in understanding the thermodynamics for molecular models of fluids in disordered porous materials,17,19 and the present analysis follows a similar point of view. We consider the behavior of a fluid confined in a fixed total volume (solid + void) of an inert adsorbent. We emphasize that the treatment given here is fully consistent and requires no further assumptions beyond the fact that the adsorbent is inert. The treatment is also applicable to molecular models of other adsorbents such as zeolites. The fundamental property relationship for the fluid can be written as
dU(a) ) T dS(a) + φ dV(a) + µ(a) dN(a)
(2)
where the superscript (a) refers to the fluid in the porous material and φ ) Ω(a)/V(a) is the grand potential density. This quantity is -P for a bulk fluid. In our notation an underbar denotes an extensive property (e.g. V(a) ) total volume). Notice that there is no decomposition of the φ dV(a) term into separate terms involving the bulk pressure and surface tension or spreading pressure, as is often done. Such a decomposition is only possible for a bulk fluid in contact with a planar solid surface where the relationship between the excess grand potential and the components of the pressure tensor is well defined. This issue was recently discussed by Cracknell and Nicholson28 in the context of a study of adsorbed solutions in slit pores. The fact that this decomposition is not possible for complex (26) Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. Mol. Phys. 1993, 78, 591. (27) Hill, T. L. Adv. Catal. 1952, 4, 211. (28) Cracknell, R. F.; Nicholson, D. Adsorption 1995, 1, 7.
Heats of Adsorption in Heterogeneous Solids
Langmuir, Vol. 12, No. 22, 1996 5427
porous material geometries should not be of great concern, since it is not necessary for a fully consistent thermodynamic treatment. The above fundamental property relationship implies a Gibbs-Duhem equation of the form
(3)
F(a)
which for constant temperature yields the Gibbs adsorption isotherm
dφ ) -F(a) dµ(a)
(4)
The isosteric heat of adsorption is related to the differential entropy change in a reversible process where fluid is transferred from the porous material at fixed volume and temperature to the bulk phase at a fixed temperature and bulk pressure. This can be written as
[( ) ( ) ] ∂S(b)
qst ) T
∂N
-
(a)
∂N
T,P
dµ(a) )
1 ) (H(b) - µ(b)) T T,P
∂S(a)
(a)
∂N
)
T,V(a)
[( )
1 ∂U(a) T ∂N(a)
- µ(a)
T,V(a)
(6)
( )
]
(7)
[ ] dP dT
[a-b],F(a)
(8) T,V(a)
[( ) ( ) ]
qd ) T
∂S
(b)
∂N
) U(b) ) qst -
∂S
-
T,V(b)
( )
(a)
∂N(a)
T,V(a)
∂U(a)
∂N(a) P
[
) F(b) S(b) +
( )
dT
(13)
F(a)
T,V(a)
F(b)
(9)
The equilibrium heat of adsorption is related to the difference between the molar (molecular) entropy of the fluid in the porous material and that in the bulk phase and is written as
qeq ) T(S(b) - S(a)) ) H(b) - H(a)
(10)
Here we have defined the enthalpy of the fluid in the
dP
(14)
F(b)
( )] ∂µ(a) ∂T
F(a)
[ ( ) ] ∂S(a)
∂N(a)
T,V(a)
)
F(b)qst T
(15)
and
[ ]
This is the expression used in previous studies of heats of adsorption in porous materials via computer simulations. However earlier derivations placed restrictions on the molar volume of the fluid in the porous material and/ or required perfect gas behavior in the bulk. The above analysis shows this expression to be correct without such assumptions. We can similarly obtain the differential heat of adsorption when the bulk fluid is at constant volume via (b)
T
∂µ(a) ∂T
At equilibrium the chemical potential of the bulk fluid and that of the porous material are equal. Using this condition and the above expressions for displacements in the chemical potentials gives
dP dT
∂U(a)
∂N(a)
∂F(a)
dF(a) +
) F(b) S(b) -
Combining these expressions gives
qst ) H(b) -
( ) ∂µ(a)
dµ(b) ) -S(b) dT +
and
( )
(12)
F(a)
We now consider expressions for heats of adsorption of the Clapeyron type using the theory of equilibrium displacements. The chemical potential of the fluid in the porous material can also be expressed as a function of F(a) and T as
T,V(a)
( ) ∂N(b)
φ
qeq ) H(b) - U(a) +
(5)
where the superscript (b) refers to the bulk fluid. This definition of the isosteric heat is similar to that given previously21,22 except that we take special note of the fact that the fluid in the porous material is at constant volume. The partial derivatives of the entropy in the above expression may be expressed as
∂S(b)
(11)
while that in the bulk can be expressed as
∂S(a)
(b)
H(a) ) U(a) - Ω(a) ) U(a) - φV(a) Substituting eq 11 into eq 10, we have
dφ
dµ(a) ) -S(a) dT -
porous material as
[a-b],φ
) F(b)(S(b) - S(a)) )
F(b)qeq T
(16)
where the subscript [a-b] denotes a derivative taken while maintaining equilibrium between the bulk fluid and the porous material. With the aid of the Gibbs-Duhem equation (eq 3), eq 15 can be rewritten as
dP [dT ]
[a-b],F(a)
) F(b)(S(b) - S(a)) -
F(b) ∂φ F(a) ∂T
( )
F(a)
(17)
It follows that
qst ) qeq -
T ∂φ F(a) ∂T
( )
F(a)
(18)
The different expressions for the heats of adsorption serve as a useful thermodynamic consistency check, as will be illustrated later. We should note that the Clapeyron expressions given above agree with those given previously21,22,27 only when the bulk fluid behaves as an ideal gas and the molar volume of the fluid in the porous material is very small relative to that in the bulk. This is a consequence of not decomposing the φdV(a) term in eq 2 into separate terms involving the bulk pressure and surface tension or spreading pressure. In this work we have principally used eq 8 for the calculation of the isosteric heat of adsorption. The partial derivative in eq 8 can be obtained either via numerical differentiation of the configurational internal energy or by using results from the theory of fluctuations in the grand ensemble. From fluctuation theory,21 we have
5428 Langmuir, Vol. 12, No. 22, 1996
Vuong and Monson
〈UN〉 - 〈U〉〈N〉
∂U (∂N )
(19)
〈N2〉 - 〈N〉2
)
T,V
where the angle brackets denote averages in the grand ensemble. The treatment we have presented is based on total properties, since these properties are directly accessible in statistical mechanical calculations. Experimental adsorption measurements, on the other hand, yield the Gibbs surface excess properties.29 It is worthwhile to demonstrate that the isosteric heat determined via eq 8 is the same quantity as that emerging in a treatment based on surface excess properties. To do this we briefly consider our model from the point of view of surface excess properties. Following Sircar,29 with some changes in notation to maintain consistency with that used earlier in this paper, we define a surface excess for a property Z as
Zm ) V(a)[F(a)Z(a) - (1 - η)F(b)Z(b)]
(20)
where 1 - η is the void fraction of the matrix. This factor appears since we have included the volume occupied by the solid in our definitions of V(a) and F(a). The surface excess number of molecules is given by
Nm ) V(a)[F(a) - (1 - η)F(b)]
Figure 1. Isosteric heat of adsorption of methane adsorbed in silica gel at T* ) 1.926 calculated by different methods: (O) eq 8; (4) eq 15; (0) eq 18.
(21)
The isosteric heat is given by
( )
qst ) -T
∂S° ∂Nm
(22) T,P,V(a)
Here S° is the total entropy of the fluid in the porous material and in the bulk, which can be written
S° ) Sm + S(b)(N° - Nm)
(23)
where N° is the total number of molecules in the system. It follows that
( )
qst ) -T
∂Sm ∂Nm
+ TS(b)
(24)
T,P,V(a)
Now using eq 20 we can write
( ) ∂Sm ∂Nm
( ) ( ) ( )
) T,P,V(a)
∂S(a) ∂Nm
∂S(a)
∂N(a)
conditions total density isosteres should be used, but this requires an estimate of the volume of the adsorbed phase, which is not an experimentally accessible quantity.
)
T,P,V(a)
T,V(a)
∂N(a) ∂Nm
( ) ∂S(a)
) T,P,V(a)
∂N(a)
Figure 2. Equilibrium heat of adsorption for methane adsorbed in silica gel at T* ) 1.926 calculated by different methods: (O) eq 12; (4) eq 16.
(25) T,V(a)
so that eq 22 reduces to eq 5. Finally, we note that there is a difficulty in applying a version of eq 15 in which the surface excess rather than the total density is held constant. This is because an isotherm of the surface excess may exhibit a maximum at high pressure for temperatures above the critical temperature of the bulk fluid, and isotherms from different temperatures may cross each other at still higher pressures. This problem was recently discussed by Brauer et al.30 It appears that the determination of isosteric heats from adsorption excess isosteres is only correct when there is a negligible difference between the total adsorption and the adsorption excess (e.g. in regions of low pressure and temperature). At other (29) Sircar, S. J. Chem. Soc., Faraday Trans. 1 1985, 81, 1527. (30) Brauer, P.; Salem, M.; Szombathely, M. V.; Heuchel, M.; Harting, P.; Jaroniec, M. In Fundamentals of Adsorption; LeVan, M. D., Ed.; Kluwer: The Netherlands, 1996; p 101.
IV. Results A. Comparison of the Heats of Adsorption Calculated via Different Methods. We begin by comparing the calculations of the isosteric heats from eqs 8, 15, and 18. Results from the three methods are shown for one temperature in Figure 1. Clearly the agreement is very good. Results for the equilibrium heat of adsorption are shown in Figure 2, and again very good agreement between the direct calculation and the Clapeyron expression is seen. These comparisons lend support to the thermodynamic consistency of our expressions for the heats of adsorption and that of our Monte Carlo simulation results. Next we consider the results for qst obtained using the two methods of evaluating (∂U(a)/∂N(a))T,V(a)sdirect numerical differentiation and the fluctuation method. To evaluate the derivative numerically, we can fit the total internal energy as a function of N(a) via a polynomial and differentiate this polynomial. Alternatively, we can rewrite (∂U(a)/∂N(a))T,V(a) as a function of F(a). That is
Heats of Adsorption in Heterogeneous Solids
Langmuir, Vol. 12, No. 22, 1996 5429
Figure 3. Comparison of the isosteric heat of adsorption via fluctuation theory and numerical differentiation for the EHS structure: (s) differentiation at T* ) 0.8; (- - -) differentiation at T* ) 2.0; (O) fluctuations at T* ) 0.8; (b) fluctuations at T* ) 2.0.
( ) (a)
∂U
(a)
∂N
T,V(a)
) U(a) + F(a)
( ) (a)
∂U
∂F(a)
(26)
T,V(a)
The internal energy per molecule can then be fitted with a polynomial in density, and the polynomial is differentiated. Note that in the grand ensemble simulations U(a) ) 〈U(a)〉/〈N(a)〉. Figure 3 shows the isosteric heat as a function of adsorbed density for two temperatures. The lower temperature is below the Lennard-Jones bulk critical temperature, and the higher temperature is substantially above the bulk critical temperature. The two sets of results appear to be in good agreement although the results from the fluctuation method exhibit more scatter. Although the differentiation method benefits from some data smoothing which accompanies the polynomial fit of the internal energy, this is not sufficient to account for the difference in the scatter. The internal energies obtained from the simulations already show a quite smooth variation with density. The key feature here is the usual difficulty in obtaining precise values of fluctuation properties in Monte Carlo simulations.25 A more limited comparison of this type was previously given by Bakaev and Steele3 for adsorption on a planar heterogeneous surface. They concluded as we have that the numerical differentiation method was better than the fluctuation method, but they further concluded that the use of the Clapeyron expression (eq 15) was the most accurate method. On the basis of our work we would conclude that the Clapeyron expression and the numerical differentiation method yield results of comparable accuracy. However, when using the Clapeyron expression, one needs at least three closely spaced (in temperature) adsorption isotherms, so the calculation requirements are greater. In any event, it is worthwhile to be able to check for thermodynamic consistency, so that using more than one method would seem to be advisable. B. Microstructure Dependence of Heats of Adsorption. We now consider the effect of translational order in the matrix upon the isosteric heat by comparing results for the EHS and fcc matrices. Figure 4 shows a comparison for these two structures at T* ) 0.8, and Figure 5 shows the same comparison for a higher temperature T* ) 2.0. The results for the two matrix structures are markedly different especially at the lower temperature.
Figure 4. Isosteric heat of adsorption of methane adsorbed in silica gel using different microstructures at T* ) 0.8: (O) EHS structure; (b) FCC structure.
At the lower temperature qst decreases as adsorbed density increases for the EHS matrix, while, for the fcc case, qst increases as adsorbed density increases. The results for the EHS matrix at the lower temperature in Figure 4 are closer to what would be typically expected for a heterogeneous adsorbent. The differences are less pronounced at the higher temperature. In this case the results for the EHS matrix show an increase with density at the higher densities. Figure 5 also shows the contribution to the isosteric heat from the bulk enthalpy. Under these conditions there is significant gas imperfection, leading to a pressure dependence in the bulk enthalpy which is sufficient to change the shape of the isosteric heat curves. The behavior of qst shown in Figures 4 and 5 can be partly understood in terms of a decomposition into the contributions from the fluid-fluid (adsorbate-adsorbate) and fluid-matrix (adsorbate-adsorbent) interactions; i.e., we can write
-(∂U(a)/∂N(a))T,V(a) ) (a) (a) (a) -(∂U(a) gs /∂N )T,V(a) - (∂Ugg /∂N )T,V(a) (27)
Figures 6-9 show the two components of the energy derivative as a function of the adsorbed density for the EHS and fcc structures at the two temperatures considered in Figures 4 and 5. From Figures 8 and 9 we notice that the adsorbate-adsorbate contribution increases monotonically and is less sensitive to the difference in microstructure for a given temperature, the microstructure dependence being greatest at the lower temperature. On the other hand, the adsorbate-adsorbent contribution (Figures 6 and 7) decreases with increasing adsorbed density and is much more sensitive to the microstructure especially at low temperature. The decrease seen in Figures 6 and 7 is associated with the energetic heterogeneity in the fluid-solid potential energy field. The most inhomogeneous regions of the matrix are filled by the fluid at low density, and under these conditions the energy derivative is most quickly varying with density. At higher densities the more energetically homogeneous regions of the matrix are filled so that the energy changes less strongly with fluid density. These effects depend strongly on the adsorbate-adsorbent potential energy distribution in the matrix. As we have shown before7 the energy distribution in the EHS matrix is much broader than that
5430 Langmuir, Vol. 12, No. 22, 1996
Vuong and Monson
Figure 7. Adsorbate-adsorbent contribution to the isosteric heat of adsorption of methane adsorbed in silica gel at T* ) 2.0: (O) EHS structure; (b) FCC structure.
Figure 5. Isosteric heat of adsorption and its component contributions for methane adsorbed in silica gel using (a) the EHS structure and (b) the FCC structure at T* ) 2.0: (b) qst; (0) -∂U/∂N; (4) H(b).
Figure 6. Adsorbate-adsorbent contribution to the isosteric heat of adsorption of methane adsorbed in silica gel at T* ) 0.8: (O) EHS structure; (b) FCC structure.
in the fcc structure so that the density dependence of the fluid-matrix contribution to the isosteric heat is greater. Another way in which we have varied the microstructure is to change the porosity of the matrix in the EHS structure.
Figure 8. Adsorbate-adsorbate contribution to the isosteric heat of adsorption of methane adsorbed in silica gel at T* ) 0.8: (O) EHS structure; (b) FCC structure.
In addition to the obvious change in confinement brought about in this way there is also the possibility of microstructural effects because the statistical geometry of the EHS system is density dependent. Figure 10 shows our results for the isosteric heat for three values of the porosity. As might be expected the isosteric heat is lowered with increasing porosity. However, the overall shape of the curves is not changed substantially except perhaps at high density due to the change in the saturation density brought about by the porosity change. A more significant change in the adsorbent microstructure is to change the connectivity of the matrix particles. This we have done by using the AHS model to generate matrix configurations. Figure 11 shows qst as a function of adsorbed density at T* ) 1.926 for both the EHS and the adhesive sphere structures. We see that, for the adhesive sphere structure, there is evidence of a sharper decrease in the isosteric heat with increasing density compared to that of the EHS structure. As can be seen in Figures 12 and 13 the difference between the EHS and AHS results is entirely accounted for in terms of the difference in the adsorbate-adsorbent contribution. C. Comparison with Experiment. There is rather little experimental data for the adsorption of simple
Heats of Adsorption in Heterogeneous Solids
Langmuir, Vol. 12, No. 22, 1996 5431
Figure 9. Adsorbate-adsorbate contribution to the isosteric heat of adsorption of methane adsorbed in silica gel at T* ) 2.0: (O) EHS structure; (b) FCC structure.
Figure 11. Isosteric heat of adsorption of methane adsorbed in silica gel using the EHS and adhesive sphere structures at T* ) 1.926: (O) EHS structure; (4) AHS structure.
Figure 10. Isosteric heat of adsorption of methane adsorbed in silica gel using the EHS structure at T* ) 1.926 for different values of the matrix volume fraction η: (O) η ) 0.386; (0) η ) 0.340; (4) η ) 0.286.
Figure 12. Adsorbate-adsorbate contribution to the isosteric heat of adsorption of methane adsorbed in silica gel at T* ) 1.926. (O) EHS structure; (4) AHS structure.
molecules in silica gel which are suitable for testing the accuracy of the present model. We have compared our simulation results with those of Masukawa31 for methane in silica gel. The adsorbent used in the experiments was a grade 15 silica gel by Davison Chemical. In this case, we used a different adsorbate-adsorbent well depth (gs ) 349 K) than that considered earlier in order to obtain agreement with the experimental Henry’s constant at 298 K. The other parameters in the composite sphere potential were left unchanged. The Lennard-Jones parameters for the fluid-fluid interaction were chosen such that the bulk fluid critical temperature and pressure are the same as that for methane.32 This leads to the values σgg ) 0.39 nm and gg ) 154.8 K. Figure 14 shows a comparison of the adsorption isotherms of methane adsorbed in silica gel at 298 K using (31) Masukawa, S. A Study on Two Phase Equilibria by Use of the Elution Gas Chromatographic Technique: The Methane-Ethane-Silica Gel System and the Methane-Normal Octane System; Rice University: Houston, 1967. Masukawa, S.; Kobayashi, R. J. Chem. Eng. Data 1968, 13, 197. (32) Smith, J. M.; Van Ness, H. C.; Abbott, M. M. Introduction to Chemical Engineering Thermodynamics, 5th ed.; Academic Press: New York, 1996.
the EHS and AHS structures for the matrix in our model. We see that the agreement is quite reasonable given the simplifications in our model, the uncertainty in our knowledge of the actual silica gel microstructure, and the difficulty in estimating the total adsorption from experimental data for the adsorption excess. At low to moderate pressure, the adhesive sphere structure gives some improvement compared to the EHS structure. However, this improvement does not persist to high pressure. Figure 15 shows the comparison of the isosteric heat of adsorption using both the EHS and the adhesive sphere structures and Masukawa’s experimental results. The experimental results were calculated using Masukawa’s estimate of the total adsorption and the Clapeyron expression (eq 15). Perhaps suprisingly the AHS structure does not give substantially better agreement with experiment for the isosteric heat than the EHS structure given that this model has a more inhomogeneous fluid-matrix potential energy field. We have also shown results from our model at a lower temperature which exhibit a decrease with density more in line with that seen in the experimental data at 298 K. If we regard the shape of the experimental isosteric heat curve to be a consequence of the inhomogeneity in the fluid-solid potential field, then Figure 15 suggests that our model does not have sufficient
5432 Langmuir, Vol. 12, No. 22, 1996
Figure 13. Adsorbate-adsorbent contribution to the isosteric heat of adsorption of methane adsorbed in silica gel at T* ) 1.926: (O) EHS structure; (4) AHS structure.
Figure 14. Adsorption isotherm of methane adsorbed in silica gel at T ) 25 °C: (b) experimental data; (O) AHS structure; (4) EHS structure.
inhomogeneity in this field to generate a steep decrease in the isosteric heat with increasing density except at low temperatures. V. Conclusions In this paper we have considered the determination of heats of adsorption from molecular models of fluids in disordered porous materials. We have given a new derivation of the expressions used to calculate the heats of adsorption which is free of assumptions made in previous work. We have used the expressions obtained in this analysis to test our simulation results for thermodynamic consistency. We have presented new results for the heats of adsorption in a molecular model of methane adsorbed
Vuong and Monson
Figure 15. Isosteric heat of adsorption of methane adsorbed in silica gel: (b) experimental data at T ) 298 K; (O) EHS structure at T ) 298 K; (0) AHS structure at T ) 298 K; (4) EHS structure at T ) 124 K.
in a high surface area silica gel. Our experience indicates that numerical differentiation of the internal energy with respect to the adsorbate density provides the most accurate route to the calculation of the isosteric heat for a given investment of computer time. However as a general rule it is probably advisable to perform such calculations by more than one method as a check for thermodynamic consistency, as has been illustrated here. In the present system the isosteric heat is very sensitive to the adsorbent microstructure, as expected. The extent to which this effect is dominated by the adsorbateadsorbent contribution is perhaps surprising, given the inhomogeneity of the fluid density in the system. The adsorbate-adsorbate contribution is relatively insensitive to the microstructure even when a strong coupling between the adsorbate-adsorbate and adsorbate-adsorbent effects would be anticipated. We have also shown that the heat of adsorption can be substantially affected by departures from ideal gas behavior in the bulk. Comparison of our simulation results with experiment showed reasonably good agreement for the adsorption isotherm but not quite so good agreement for the isosteric heat. This suggests to us that the present molecular model may not exhibit sufficient nonuniformity in the adsorbateadsorbent potential energy field. In work in progress we are investigating making the adsorbate-adsorbent potential energy field more nonuniform by roughening the surface of the matrix spheres. Acknowledgment. This work was supported by the National Science Foundation (Grant Nos. CTS-9115297 and CTS-9417649). The authors are grateful to N. Seaton for his assistance with the simulation of the adhesive sphere model. LA960325M
