9976
Langmuir 2006, 22, 9976-9981
Correct Procedures for the Calculation of Heats of Adsorption for Heterogeneous Adsorbents from Molecular Simulation G. R. Birkett and D. D. Do* Department of Chemical Engineering, UniVersity of Queensland, St. Lucia, Queensland 4072, Australia ReceiVed May 29, 2006. In Final Form: September 5, 2006 Several procedures for calculating the heat of adsorption from Monte Carlo simulations for a heterogeneous adsorbent are presented. Simulations have been performed to generate isotherms for nitrogen at 77 K and methane at 273.15 K in graphitic slit pores of various widths. The procedures were then applied to calculate the heat of adsorption of an activated carbon with an arbitrary pore size distribution. The consistency of the different procedures shows them to be correct in calculating interaction energy contributions to the heat of adsorption. The currently favored procedure for this type of calculation, from the literature, is shown to be incorrect and in serious error when calculating the heat of adsorption of activated carbon.
1. Introduction Equilibrium studies of adsorption are concerned primarily with two main quantities. The first, and most important, is the amount adsorbed, and the other is the heat of adsorption. The heat of adsorption is required for meaningful modeling of many adsorption processes. This can range from interest in the heat generated itself (e.g., adsorption refrigeration cycles) to the impact that the heat generation and dissipation has on adsorption kinetics (e.g., industrial adsorption beds). Commensurate with their relative importance, the amount of experimental and theoretical work dedicated to the adsorption capacities of adsorbents far outweighs that dedicated to heats of adsorption. However, this does not only come about because of the relative importance of each topic but also because the measurement of the heat of adsorption is significantly more difficult than the measurement of an adsorption isotherm. The heat of adsorption must be either approximated using multiple isotherms and Classius-Clapeyrontype equations1 or measured directly using calorimetric methods. The Classius-Clapeyron route requires isotherms at multiple temperatures to calculate the heat of adsorption. It is prone to large calculation errors and has been shown to be inapplicable to strongly adsorbing species or adsorption at high pressure where bulk densities approach those of the adsorbed phase.2,3 Calorimetric methods are more direct but require significantly greater expense than the measurement of adsorption isotherms alone. Simulations provide an appealing way of measuring heats of adsorption because they can be calculated directly at any pressure and temperature. The calculation of the heat of adsorption using simulation is well established for a single simulation.4 To be able to compare with experiment or use in adsorber calculations, these results for a single pore (or arbitrary structure) have to be combined to give the heat of adsorption of the heterogeneous solid of interest, such as activated carbon. The combination of isosteric heats of single pores to calculate the isosteric heat of a solid has been * To whom correspondence should be addressed. E-mail:
[email protected]. Phone: +61 7 3365 4154. Fax: +61 7 3365 2789. (1) Ross, S.; Olivier, J. P. On Physical Adsorption; Interscience Publishers: New York, 1964. (2) Terzyk, A. P.; Rychlicki, G. Adsorpt. Sci. Technol. 1999, 17, 323. (3) Sircar, S.; Mohr, R.; Ristic, C.; Rao, M. B. J. Phys. Chem. B 1999, 103, 6539. (4) Nicholson, D.; Parsonage, N. G. Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: London, 1982.
Table 1. Molecular Parameters molecule or atom
σff (nm)
ff/kB (K)
N2 C H
0.3615 0.331 0.331
101.5 0.01 15.3
presented in the literature.5-7 It is proposed that the method currently used to calculate the solid’s isosteric heat is incorrect. Correct methods of calculating the isosteric heat are presented in this article, and they are applied to the adsorption of nitrogen at 77 K and methane at 273.15 K on a model activated carbon.
2. Details of Simulation Simulations were performed for the adsorption of nitrogen at 77 K and methane at 273.15 K in graphitic slit pores. The molecular model used for nitrogen is a single-site Lennard Jones model from Ravikovitch et al.8 The model for methane is a fivesite model from Chen and Siepmann9 with interaction sites at the carbon and each of the hydrogen atoms. The molecular parameters are given in Table 1 where σff and ff are the Lennard Jones collision diameter and well depth of the interaction energy, respectively, and kB is the Boltzmann constant. The pores of the model activated carbon are treated as infinite slit pores. The potential energy between a pore wall and a fluid interaction site is calculated using the Steele potential.10 The fluid-solid interaction parameters are obtained using the Lorentz-Berthelot mixing rule.11 Simulations were performed in the standard GCMC ensemble.11 The simulation box was bounded in one direction by the finite width, H, of the pore and was given a box length in the remaining two directions of either 3.0 nm or 1.5 times the slit width (whichever was greater). This has been found to be sufficient in avoiding observable finite size effects due to the extent of the box length.12 Periodic boundary conditions were used in the (5) Pan, H.; Ritter, J. A.; Balbuena, P. B. Ind. Eng. Chem. Res. 1998, 37, 1159. (6) He, Y.; Seaton, N. A. Langmuir 2005, 21, 8297. (7) Nicholson, D. Langmuir 1999, 15, 2508. (8) Ravikovitch, P. I.; Vishnyakov, A.; Neimark, A. V. Phys. ReV. E 2001, 64, 1. (9) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 5370. (10) Steele, W. A. Surf. Sci. 1973, 36, 317. (11) Allen, M. P.; Tildesley, T. P. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (12) Jorge, M.; Seaton, N. A. Mol. Phys. 2002, 100, 3803.
10.1021/la0615241 CCC: $33.50 © 2006 American Chemical Society Published on Web 10/28/2006
Procedures for Calculation of Heats of Adsorption
Langmuir, Vol. 22, No. 24, 2006 9977
lateral directions of the pore with fluid-fluid interactions having a lateral cutoff equal to half the simulation box length. For the nitrogen adsorption at 77 K, simulation pressures ranged from 1 × 10-10 to 60 kPa with chemical potentials set to that of an ideal gas at that pressure. For methane at 273.15 K, simulation pressures ranged from 0.05 to 20 000 kPa with the chemical potential set to that of methane according to an accurate equation of state.13 Simulations were started at the lowest pressure from a grid configuration with subsequent isotherm points taking their starting configuration from the final configuration of the previous pressure. Isotherm points were equilibrated until the number of particles was deemed to be at equilibrium. The equilibrium period ranged between 3 × 106 and 150 × 106 attempted Monte Carlo moves depending on the number of particles in the simulation and the stage of adsorption. The collection of ensemble averages was conducted over a range of 3 × 106 to 30 × 106 attempted Monte Carlo moves, depending on the number of particles in the simulation. It was attempted, during the equilibration stage, to set the number of insertions such that a successful insertion was achieved, on average, every five attempted particle moves. This was not always possible, and where acceptance ratios of insertions were very low, the maximum number of insertions was set to 5 times the number of attempted moves. The heat of adsorption can have many definitions, including isosteric heat, differential heat, and equilibrium heat.14 Additionally, there is the difference between absolute and excess heats of adsorption that depend on assumptions of the adsorption system and the adsorbing fluid.15 From a simulation standpoint, all of these quantities require the calculation of the energetic contribution of adsorption to the isosteric heat. There are two major methods of calculating the heat of adsorption using simulation. The first comes from the classical thermodynamic treatment of the adsorption system14 and is given in eq 1
qast(P) ) qst(P) - Hb(P) ) -
( ) ∂Ua(P) ∂Na(P)
(1)
T,Va
where qast is the configurational contribution of the adsorbed phase to the isosteric heat, Hb is the molar enthalpy of the bulk fluid, qst is the absolute isosteric heat, Ua is the total configurational energy of the adsorbed fluid within the pore (simulation), and Na is the total number of particles in the pore with all being a function of pressure, P, at a fixed temperature, T, and pore volume, Va. Now during the course of a GCMC simulation, qast can be calculated directly using fluctuation theory4
qast =
〈NaUa〉 - 〈Na〉〈Ua〉 〈NaNa〉 - 〈Na〉〈Na〉
Figure 1. Pore size distribution used in the simulation study. Pore widths, H, are from carbon center to carbon center of opposing walls, and pore volumes are accessible volumes for nitrogen.
function of pressure, P. We define the absolute amount adsorbed by eq 3
Natot(P) )
(13) Setzmann, U.; Wagner, W. J. Phys. Chem. Ref. Data 1991, 20, 1061. (14) Vuong, T.; Monson, P. A. Langmuir 1996, 12, 5425. (15) Myers, A. L.; Monson, P. A. Langmuir 2002, 18, 10261.
(3)
where F(H, P) is the density of a pore with width H at pressure P and f(H) dH is the volume of pores having a width between H and H + dH. For the practical use of eq 3, the integral is replaced by a sum over discrete intervals using the pore size distribution (PSD) of the solid. The PSD used in this study comes from the fitting of nitrogen adsorption data of Ajax activated carbon at 77 K using density functional theory and is presented in Figure 1. The choice of this pore size distribution is essentially arbitrary because it does not effect the conclusions in this article and is required only for illustration purposes. Although the calculation of the total amounts adsorbed using eq 3 is clear, the calculation of isosteric heat for a heterogeneous solid using individual isosteric heats of individual pores has been calculated incorrectly in many papers. Below we propose the correct methods for calculating the isosteric heat. To do this, we start with a version of eq 1 applicable to the whole solid with the energy and the amount adsorbed by the pore replaced with that of the solid. This is given in eq 4 a (P) qst,tot
(2)
where the brackets, 〈...〉, denote the ensemble average obtained from simulation. Therefore, either eq 1 or 2 can be used to calculate qast for a single pore. Generally, eq 2 is preferable when using GCMC because it can be calculated during the simulation. Equation 1 requires either the fitting of Ua as a function of Na for differentiation or a discrete approximation of the differential using changes in Ua and Na from one isotherm point to the next. All values of qast for single pores are calculated using eq 2 in this study. Because we are interested in the isosteric heat of a solid rather than a single pore, we will now outline the equations used to calculate the adsorption properties of a porous solid. The first quantity of interest is the total amount adsorbed, Natot, as a
∫0∞F(H, P) f(H) dH
)-
( ) ∂Uatot(P) ∂Natot(P)
(4)
T,Vatot
where qast,tot is the total contribution of the adsorbed phase in the solid to the isosteric heat, Vatot is the total volume of the adsorption system, and Uatot is the total energy of the adsorbed phase given by, in analogy to eq 3, the following equation
Uatot(P) )
∫0∞ua(H, P) F(H, P) f(H) dH
(5)
where ua(H, P) is the average configurational energy per mole of particles with a pore of width H at a pressure P. The use of eq 4 to calculate the isosteric heat requires the fitting of Uatot with respect to Natot. This can be fitted with any sort of function as long as it fits the data in a realistic manner (i.e., in a smooth and continuous manner). A discrete form of eq 4 can be used that requires no fitting and can be calculated directly from the simulation data. For this study, we approximate the slope of the energy versus the amount adsorbed at pressure Pi as the average
9978 Langmuir, Vol. 22, No. 24, 2006
Birkett and Do
of the slope between data points at Pi and Pi - 1 and points at Pi+1 and Pi. This discrete approximation of eq 4 is presented in eq 6. a qst,tot (Pi) ≈ -0.5
[(
)
Uatot(Pi) - Uatot(Pi - 1) Natot(Pi) - Natot(Pi - 1)
(
+
Uatot(Pi+1) - Uatot(Pi) Natot(Pi+1) - Natot(Pi)
)]
(6)
As long as there are sufficient isotherm points over the pressure range of interest, eq 6 should be a good approximation of eq 4. Both eqs 4 and 6 are used in this study to calculate the isosteric heat. An alternative to eq 4 is to combine the isosteric heats calculated for the individual pores using eq 2. Now this relies on being able to calculate each pore’s contribution to any incremental increase in the total amount adsorbed. This method is defined by eq 7, whose derivation from eq 4 is given in the Appendix.
a (P) qst,tot
)
∫0∞qast(P, H)(
(
)
∂F(P, H) f(H) dH ∂P T,Va ∂Natot(P) ∂P T,Vatot
)
(7)
The numerator in eq 7 calculates the total heat generated in all of the pores by an incremental change in the pressure. This is divided by the total incremental change in the amount adsorbed to give the isosteric heat of the solid. Now to avoid having to fit and differentiate each pore’s density with respect to pressure, it is more convenient to put eq 7 in a discrete form using the individual results from each pressure, Pi. The discrete form of eq 7 is given by eq 8. a qst,tot (Pi) ≈
∫ q (H, P )0.5{( ∞
0
a st
i
{(
0.5
) (
)}
F(H, Pi) - F(H, Pi - 1) F(H, Pi+1) - F(H, Pi) + Pi - P i - 1 Pi+1 - Pi Pi - Pi - 1
Natot(Pi)
)(
Natot(Pi - 1)
+
Pi+1 - Pi
Natot(Pi+1)
)}
f(H) dH
Natot(Pi)
(8)
Equation 8 is used for the calculation of isosteric heat in this study. This should provide a good approximation of eq 7 as long as there are a sufficient number of isotherm points. Now that we have presented the proposed formulations for the isosteric heat of a solid, we can look at the technique currently used in the literature. The equation commonly used in the literature to calculate the heat of adsorption for a solid5-7 is given by eq 9. a qst,tot (P)
)
∫0∞qast(H, P) F(H, P) f(H) dH Natot(P)
(9)
Now eq 9 gives the isosteric heat of the solid as a weighted average based upon absolute loading at a given pressure, but this is incorrect. The derivation of eq 7 in the Appendix shows that eq 9 cannot be equal to eq 4. This can be demonstrated by a simple example of an idealized solid, pictured in Figure 2, with two pores i and j. The smaller pore, i, adsorbs fluid molecules and releases heat in going from pressure 0 to P1, P1 to P2, and P2 to P3. At pressure P3, the pore is full and adsorbs no additional molecules in going from P3 to P4. Without adsorption, no heat is released by pore i, and it cannot make a contribution to the overall isosteric heat of the solid at pressures above P3. This is
Figure 2. Illustrative diagram of adsorption and heat generation (Q) of an idealized solid with two pores i and j. From top left, pressure increases in a clockwise direction (P1 < P2 < P3 < P4).
reflected by eq 7 because ∂F/∂P ) 0 for pore i at pressures above P3. However, eq 9 would still give significant weight to the isosteric heat of pore i, even though no adsorption is occurring in that pore at pressure above P3. In the next section, we will calculate isosteric heats, using eqs 4, 6, and 8, for the adsorption of the subcritical nitrogen and supercritical methane on an arbitrary activated carbon. These results will be contrasted against the results from using the incorrect eq 9.
3. Results and Discussion Adsorption isotherms and isosteric heats from simulations for nitrogen, at 77 K in selected pore sizes, are presented in Figure 3. This shows typical adsorption behavior for a subcritical nonpolar fluid. The smallest pore shown (0.733 nm) fills at very low pressure whereas the larger pores show layering followed by filling. The isosteric heats obtained are also as expected, with the large pores showing isosteric heat behavior similar to that observed on graphitic surfaces16 up until the discontinuity caused by the condensation in the pore. After condensation, the calculation of the isosteric heat using eq 2 becomes very unreliable because of the very low acceptance of insertions. The calculation is equally unreliable using eq 1 because the change in Na with pressure becomes almost negligible. Adsorption isotherms and isosteric heats from simulation for methane, at 273.15 K in selected pore sizes, are presented in Figure 4. Figure 4 shows typical isotherms for supercritical adsorption in carbon slit pores. For all pores, the isosteric heat initially increases with loading before going through a local maximum. The larger pores then go through a local minimum in the isosteric heat before it increases again as a result of the high fluid density in the pore that gives rise to greater fluid-fluid interaction. Only a maximum is seen in the isosteric heat for the 0.733 nm pore because only a single layer is possible and adsorption at high pressure requires molecules to be close enough to start to repel each other to achieve higher pore densities. Again, as densities increase, the reliability of the isosteric heat calculated using eq 2 decreases but is markedly better than for the case of nitrogen at high density at 77 K. It is important to emphasize that these are absolute isosteric heats and need appropriate adjustments to be able to be matched against experimental data. Now we can look at the isosteric heat calculated for the activated carbon (16) Avgul, N. N.; Kiselev, A. V. Chem. Phys. Carbon 1970, 6, 1.
Procedures for Calculation of Heats of Adsorption
Figure 3. Simulation results for nitrogen at 77 K. (a) Pore density versus simulation pressure and (b) isosteric heat (using eq 2) versus pore density. The different pore widths plotted are 0.733 nm (b), 2.056 nm (0), and 4.039 nm (2). Lines are a guide for the eye only.
Figure 4. Simulation results for methane at 273.15 K. Plots and symbols are the same as in Figure 3.
represented by the pore size distribution in Figure 1. The adsorption isotherm, calculated using eq 3, and isosteric heats, calculated using eqs 4, 6, and 8, for nitrogen at 77 K are presented in Figure 5. Also shown in Figure 5 is the isosteric heat calculated using eq 9.
Langmuir, Vol. 22, No. 24, 2006 9979
Figure 5. Isotherm (a) and isosteric heats (b) for nitrogen at 77 K on the solid with the PSD in Figure 1. The symbols in plot b represent different calculation methods: using eq 4 with the differentiation of the polynomial used to fit the function of Uatot versus Natot (-b-) and eq 6 (0), eq 8 (2) and eq 9 (‚‚)‚‚).
Figure 5 shows the consistency of eqs 4, 6, and 8 in calculating the isosteric heat. It also shows the significant overestimation of the isosteric heat by eq 9. The overestimation by eq 9 is expected because it will overestimate the contribution to the isosteric heat of the smaller pores that have a higher isosteric heat than the larger pores. It can be seen in Figure 3 that the 0.733 nm pore has a significantly higher isosteric heat and pore density, at low pressure, than the larger pores. However, the additional adsorption in the 0.733 nm pore beyond a pressure of 1 × 10-3 is quite small, meaning that its fractional contribution to the incremental total amount adsorbed at higher pressures is negligible. It is a pore’s contribution to incremental, rather than total, loading that determines its observable contribution to the isosteric heat. This is reflected in eqs 4, 6, and 8 but not in eq 9. Figure 6 shows the details of the isosteric heat results for nitrogen at low loading. This shows that eq 9 converges to the correct value of isosteric heat only in the limit of the amount adsorbed approaching zero. The same calculations have been performed with the methane adsorption at 273.15 K, and the results of this are presented in Figure 7. Figure 7 shows similar results to Figure 5. The agreement between the different techniques, using eqs 4, 6, and 8, is very good. The improved agreement and smoothness of all techniques are the results of an increased number of isotherm points and the continuous nature of supercritical adsorption (comparing Figure 4 to Figure 3). The isosteric heat calculated using eq 9 again has a large error. The only place where eq 9 can even approximate the actual isosteric heat is again at very low loading. From Figures 5 and 7, it can be seen that the smoothest function of isosteric heat is from using eq 4 and a fitted function of Uatot with respect to Natot. This is a natural consequence of the fitting process but does not necessarily indicate any greater accuracy. However, it appears to provide a reasonable best fit of the other techniques
9980 Langmuir, Vol. 22, No. 24, 2006
Figure 6. Same plot as in Figure 5b but over smaller ranges of the amount adsorbed and the isosteric heat to show behavior close to zero loading.
Birkett and Do
or by using the discrete changes in energy and the amount adsorbed to approximate the differential. The second technique used the isosteric heats calculated for individual pores. The isosteric heat of the solid was then calculated as the sum of the pores’ isosteric heat weighted by its fractional contribution to the incremental amount adsorbed. All methods were in agreement with each other, and the preference for use is governed by the convenience of calculation. The methods presented are completely transferable to any adsorption system where multiple simulations are required to account for the heterogeneity of the solid. The currently used method in the literature, weighting the isosteric heat with the total amount adsorbed rather than the incremental amount, was shown to be incorrect. Future studies deriving the isosteric heat of a solid through the combination of individual pore isosteric heats should use one of the methods proposed in this article and not the method currently used in the literature. Acknowledgment. This research was made possible by the Australian Research Council, whose support is gratefully acknowledged. Thanks also to the University of Queensland High Performance Computing facility for a generous allocation of computing time.
Appendix Equation 7 is derived starting with eq 4:
a qst,tot (P)
)-
( ) ∂Uatot(P) ∂Natot(P)
(A1)
T,Vtota
Now with the system under constant temperature and at a constant volume, the following derivation is made. (Subscripts to this effect are not used for clarity.) First we apply the product rule to eq A1.
( ) ∂Uatot(P) ∂Natot(P)
Figure 7. Isotherm (a) and isosteric heats (b) for methane at 273.15 K on the solid with the PSD in Figure 1. Symbols are the same as in Figure 5.
( (
) )
∂Uatot(P) ∂P ) ∂Natot(P) ∂P
(A2)
Now the total energy of the system is the sum of energy contributions from each pore. The total energy of the adsorption system is given by
used. From an ease of calculation standpoint, eqs 6 and 8 provide the advantage. These are easily calculated in any spreadsheet and can be updated for any pore size distribution just as easily as the amount adsorbed.
Uatot(P) )
∫0∞ua(P, H) F(P, H) f(H) dH
(A3)
Now we can write
4. Conclusions Two techniques for the calculation of the isosteric heat of a heterogeneous solid from the isosteric heats of individual pores have been presented. The techniques were applied to the adsorption of nitrogen at 77 K and methane at 273.15 K on a model activated carbon using a realistic (but otherwise arbitrary) pore size distribution. The first technique uses the thermodynamic definition of isosteric heat as the differential of the adsorption system energy with respect to the amount adsorbed and applies it to the whole solid. This can be done either by fitting the system energy as a function of the amount adsorbed and differentiating
(
)
∫0∞ua(P, H) F(P, H) f(H) dH
∂ ∂Uatot(P) ) ∂P
∂P
(A4)
Because this is a fixed-volume system, meaning H and f(H) are invariant with respect to P, F, and ua, eq A4 becomes
(
)
∂Uatot(P) ) ∂P
∂{ua(P, H) F(P, H) f(H) dH} ∂P
∫0∞
(A5)
Procedures for Calculation of Heats of Adsorption
Applying the product rule and the definition of Ua to eq A5 gives
∂{u (P, H) F(P, H) f(H) dH} ) ∂P a ∂{u (P, H) F(P, H) f(H) dH} ∂Na(P) ∂Ua(P) ∂Na(P) ) ∂P ∂P ∂Na(P) ∂Na(P) (A6) a
(
)( ) ( )( )
Now the first factor on the RHS of eq A6 is the negative of the isosteric heat of a pore. Substituting this in and replacing the amount adsorbed in the pore with its density and volume, we obtain
Langmuir, Vol. 22, No. 24, 2006 9981
( )( )
(
)
∂Ua(P) ∂Na(P) ∂F(P, H) f(H) dH ) -qast(P, H) ) a ∂P ∂P ∂N (P)
(
-qast(P, H)
Combining eqs A2, A5, A6, and A7 gives
a (P) qst,tot
)
∫0∞qast(P, H)(
)
∂F(P, H) f(H) dH ∂P T,Va ∂Natot(P) ∂P T,Vatot
(
)
such that eq 7 is shown to be equal to eq 4. LA0615241
)
∂F(P, H) f(H) dH (A7) ∂P
(A8)
