Adsorption of Supercritical Fluids and Fluid Mixtures - American

The fluid density profiles are obtained from the inhomogeneous integral equation theory. The theory is shown to be in excellent agreement with Monte C...
0 downloads 0 Views 112KB Size
Download PDF
J. Phys. Chem. B 2001, 105, 6583-6591

6583

Adsorption of Supercritical Fluids and Fluid Mixtures: Inhomogeneous Integral Equation Study† S. A. Egorov* Department of Chemistry, UniVersity of Virginia, CharlottesVille, Virginia 22901 ReceiVed: January 3, 2001

A microscopic statistical mechanical theory of adsorption of supercritical fluids and fluid mixtures on solid surfaces is presented. The fluid density profiles are obtained from the inhomogeneous integral equation theory. The theory is shown to be in excellent agreement with Monte Carlo simulations and to provide substantial improvement over the integral equation theory formulated for homogeneous fluids. A quantitative comparison of theory with the experimental data on the surface excess of supercritical argon and krypton on graphite is presented, and the theory is found to be in good agreement with experiment. In addition, model calculations of adsorption from fluid mixtures are carried out, and the phenomenon of preferential adsorption is discussed.

I. Introduction Supercritical fluids (SCFs) are currently attracting significant industrial and scientific interest due to their unique physical properties.1-3 SCFs are characterized by high dissolving power and enhanced mass-transfer rates, thereby providing an attractive alternative to conventional liquid solvents for a variety of industrial processes.2,4,5 A prominent place among practical applications of SCFs belongs to extraction and separation processes, which generally involve either physical adsorption of SCFs on solid surfaces or selective adsorption/desorption of solutes from/into SCFs. Particular examples of technological applications based on SCF adsorption include transportation and storage of fuel gases,6 separation of lower hydrocarbons,7 remediation of contaminated soils,8 adsorbent regeneration,9 and SCF chromatography.10 Further development of these applications would greatly profit from the ability to predict adsorption isotherms at supercritical conditions. Adsorption of SCFs has been extensively studied experimentally, and certain characteristic patterns in the behavior of adsorption isotherms have emerged.11-13 In particular, the surface excess of the fluid measured as a function of density along a given supercritical isotherm generally passes through a maximum somewhat below the critical density. The height of this maximum depends sensitively on the temperature, and increases rapidly as the critical temperature Tc is approached from above, indicating the multilayer formation, or prewetting of the solid by the fluid. As one gets closer to Tc, the location of the maximum moves toward the critical density. This behavior was observed for numerous adsorbent/adsorbate combinations and, therefore, appears to be fairly universal. In addition to the adsorption of neat SCFs, much attention has been devoted to the adsorption of fluid mixtures, where the phenomenon of preferential adsorption forms the basis of such practical applications as adsorptive separation of hydrocarbons. Experimental studies indicate that under appropriate thermodynamic conditions one can achieve a high degree of selectivity in the adsorption of a given component from the fluid mixture.14 †

Part of the special issue “Bruce Berne Festschrift”

Experimentally observed trends in the density and temperature dependence of adsorption isotherms of neat SCFs and fluid mixtures were successfully reproduced by simulation methods, such as grand canonical ensemble Monte Carlo (MC).14-17 The main advantage of the simulation-based techniques lies in the fact that for a given model these methods produce essentially exact results (albeit with inevitable statistical noise). However, simulation studies are computationally demanding, especially in the near-critical region, where very large system sizes are required in order to obtain converged results. Alternatively, one can employ microscopic statistical mechanical methods, such as density functional and integral equation theories, which are generally less expensive computationally. Unfortunately, these approaches involve unavoidable approximations, and therefore, their accuracy needs to be tested against computer simulations. Some guidance in choosing a reliable theoretical method for treating SCF adsorption can be obtained from the fact that adsorption of a fluid on a solid surface closely resembles in many respects the phenomenon of clustering of a supercritical solvent around a dilute solute, as has been already observed by several authors.18-20 It is now well established that in “attractive” supercritical solutions, where the solvent is attracted to the solute more strongly than to itself, the local density around the solute in the near-critical region is enhanced relative to the bulk value.21-25 Theoretical studies performed for attractive mixtures along near-critical isotherms have shown that the excess local solvent density around the solute increases rapidly with the bulk density in the low-density region, passes through a maximum somewhat below the critical density and then gradually decreases as one moves into the high-density liquidlike region. The height of the maximum in the density decreases for higher temperature isotherms. Thus, density and temperature dependence of the excess local density around an attractive solute exhibits the same general trends as the surface excess of a supercritical adsorbate at a planar wall. The similarity discussed above is hardly accidental because the adsorbing wall can be viewed as an infinitely large attractive solute. Accordingly, several theoretical treatments of clustering in dilute supercritical solutions have adapted the methodology employed in the studies of adsorption. Namely, Kajimoto et

10.1021/jp010007i CCC: $20.00 © 2001 American Chemical Society Published on Web 05/04/2001

6584 J. Phys. Chem. B, Vol. 105, No. 28, 2001 al.18 and Flanagin et al.19 have modeled the local solvent density around the solute in terms of Langmuir adsorption equilibria, while Subramanian et al.20 have described clustering in supercritical solutions on the basis of their simplified local density adsorption model. In the above approaches, the adsorption models are employed to interpret the data on the local density enhancement in supercritical solutions. Conversely, the close analogy between adsorption and clustering can be exploited by going in the opposite direction, i.e., by adapting microscopic statistical mechanical methods of treating local density augmentation around attractive solutes in order to calculate the density profiles of SCFs at adsorbing walls. In our recent work, we have analyzed the accuracy of various integral equation theories in reproducing the local solutesolvent structure of attractive supercritical solutions.26 It was found that standard integral equation approach developed for uniform fluids can be seriously in error in the near-critical region, which is due to its inherent assumption that the pair correlation functions for two solvent particles located near the solute are identical to those in the bulk. This assumption precludes an accurate treatment of the solute-solvent structural properties of attractive near-critical solutions because the latter are characterized by significant solute-induced solvent density inhomogeneities, which can seriously affect the solvent pair correlations in the vicinity of the solute. In view of the above drawback of the uniform fluid approach in treating attractive supercritical solutions, I have subsequently employed an extended version of the integral equation theory,27 which was formulated specifically for inhomogeneous fluids.28-30 This method operates with a nonuniform solvent density profile created by the solute molecule, thereby accounting for the local density inhomogeneities in a consistent way. By performing model calculations for attractive mixtures over a wide range of solvent thermodynamic conditions and by comparing theoretical results with simulation, I have shown that the inhomogeneous fluid approach does provide a very substantial improvement over its homogeneous counterpart in calculating the local solutesolvent structural properties.27 Similarly to an attractive solute, an adsorbing wall creates a strongly nonuniform fluid density profile. Accordingly, the pair correlation functions for two adsorbate particles in the vicinity of the wall can be expected to depend not only on the separation between the two particles (as is the case in the isotropic bulk fluid phase), but also on their positions relative to the wall. Consequently, one could expect the inhomogeneous integral equation theory to be more appropriate for calculating the fluid density profile at the wall in comparison to its homogeneous counterpart. This was indeed confirmed in the study of a supercritical Lennard-Jones (LJ) fluid near a hard wall conducted by Plischke et al.,31 where the results of both homogeneous and inhomogeneous integral equation theories were compared with Monte Carlo simulations, and the latter theory was found to be far more accurate. In view of the successful performance of the inhomogeneous integral equation theory for the model of LJ fluid at a hard wall, it would be of great interest to assess its accuracy for a more realistic model of LJ fluid at LJ wall, which is employed in the vast majority of the simulation studies of SCF adsorption. Some integral equation calculations for this model were carried out,32 but no direct comparison with simulation data was performed. One of the goals of the present work is to perform such a comparison. It will be shown that the inhomogeneous integral equation theory is indeed adequate for calculating the fluid density profile at the wall for the model based on LJ potentials.

Egorov Once the accuracy of the theory is ascertained, it will be employed to analyze the experimental data on SCF adsorption, and to perform model calculations of the adsorbate surface excess and its dependence on the fluid density, temperature, and the strength of the wall-fluid interaction. In addition, the theory will be used to study preferential adsorption from fluid mixtures. The paper is organized as follows. In section II, we briefly outline the inhomogeneous integral equation theory which we use to calculate the fluid density profile at a planar wall. In section III, we compare theoretical results with MC simulations. In section IV, we perform a quantitative analysis of the experimental data on the adsorption of supercritical argon and krypton on graphite. We also carry out model calculations of adsorption isotherms for neat SCFs and fluid mixtures for a wide range of thermodynamic conditions. In section V, we conclude. II. Fluid Density Profiles at the Wall In the present work, we are concerned with a fluid (adsorbate) in contact with a flat wall (adsorbent). We start by considering the case of a neat fluid; generalization of the formalism to the case of a fluid mixture is straightforward and will be discussed at the end of this section. The wall is taken to be located at z ) 0, with the fluid occupying the half-space z > 0. We assume that the particles of the fluid interact with each other via isotropic pair potential φ(r); the wall-particle interaction for an adsorbate particle located at a distance z from the wall is denoted by V(z). As discussed in the Introduction, we use the inhomogeneous integral equation theory to calculate the density profile of the fluid in the vicinity of the wall. To this end, we write the inhomogeneous Ornstein-Zernike (OZ) equation for two fluid r2 ) (x2,y2,z2) particles located at b r ) (x1,y1,z1) and b

h(b r 1, b r 2) ) c(b r 1, b r 2) +

∫dbr 3F(br 3)c(br 1,br 3)h(br 3,br 2)

(1)

where F(r b) is the singlet density of the nonuniform fluid, and h(r b1,r b2) and c(r b1,r b2) are the anisotropic total and direct pair correlation functions of the fluid, respectively. In view of the symmetry of the problem, the fluid density profile varies only along the normal to the wall, i.e., F(r b) ≡ F(z), whereas the pair correlation functions depend on three variables: c(r b1,r b2) ≡ c(z1,z2,s) and h(r b1,r b2) ≡ h(z1,z2,s), where z1 and z2 are the distances of the two particles from the wall, and the coordinate s is defined according to s ) x(x1-x2)2+(y1-y2)2. The inhomogeneous OZ equation can now be simplified by performing a Fourier-Bessel transform31

hˆ (z1,z2,k) ) cˆ (z1,z2,k) +

∫0∞dz3F(z3)cˆ (z1,z3,k)hˆ (z3,z2,k)

(2)

where

∫0∞sdsJ0(ks)f(z1,z2,s)

(3)

∫0∞kdkJ0(ks)fˆ(z1,z2,k)

(4)

ˆf (z1,z2,k) ) 2π and

f(z1,z2,s) )

1 2π

where J0(x) is the zeroth-order Bessel function. The inhomogeneous OZ equation needs to be supplemented with a closure, which relates the total and direct pair correlation functions. In the present work, we use the inhomogeneous

Inhomogeneous Integral Equation Study

J. Phys. Chem. B, Vol. 105, No. 28, 2001 6585

Percus-Yevick (PY) closure

c(z1,z2,s) ) [h(z1,z2,s) + 1][1 - exp(βφ(r12))]

(5)

where β ) 1/kBT and r12 ) |r b1 - b r2|. Since eq 1 contains three unknown functions, two pair correlation functions and the density profile, an additional relation is needed to obtain a closed system of equations. Several exact relations between the density profile and the pair correlation functions are known.33-36 Here, we employ the LovettMou-Buff-Wertheim equation,34,35 which, for our problem, takes the following form

F′(z1)/F(z1) ) -βV′(z1) +

∫0∞dz2F′{z2)cˆ (z1,z2,0)

(6)

where the prime denotes differentiation with respect to argument. By solving eqs 2, 5, and 6 simultaneously, one obtains anisotropic pair correlation functions and the fluid density profile in a self-consistent manner. The numerical procedure for solving the above system of equations involves calculating the correlation functions on a grid for both z and s (or, equivalently, k) variables. By choosing the points si (i ) 1, ‚‚‚, N) and kj (j ) 1, ‚‚‚, N), which satisfy the relations J0(sikN) ) 0 and J0(kjsN) ) 0 for all i and j, the FourierBessel transform can be represented by a discrete sum37 N

ˆf (z1,z2,kj) ) 4π/k2N

f(z1,z2,si)J0(kjsi)/J12(sikN) ∑ i)1

(7)

In our calculations, we have found that taking N ) 100 with a cutoff distance of 6.0 (in units of the fluid particle diameter) was sufficient to obtain converged results. Regarding the z variable, the results were found to be rather sensitive to the choice of the grid spacing and the cutoff distance along the z direction. To achieve convergence while keeping the memory requirements for the computation practical, we have employed a nonuniform grid spacing: we set ∆z ) 0.025 in the vicinity of the wall (up to the distance of two fluid particle diameters) and ∆z ) 0.1 beyond this distance, with the total number of grid points varying between 150 and 180 depending on the thermodynamic conditions. It is instructive to perform a brief comparison of the inhomogeneous integral equation theory represented by eqs 2, 5, and 6 with its homogeneous counterpart. To this end, we note that the homogeneous OZ equation for the fluid pair correlation functions can be obtained from the corresponding inhomogeneous relation by making the following two approximations: (1) the density is taken to be spatially uniform, i.e., F(r b) ≡ F, where F is the bulk fluid density, and (2) the solvent total and direct correlation functions are assumed to be isotropic, e.g., c(r b1,r b2) ≡ c(r12). Under these two assumptions, eq 1 reduces to the OZ equation for an isotropic fluid

∫

r ˜ (r′)c˜ (|b r - b′|) r h˜ (r) ) c˜ (r) + F db′h

(8)

Similarly, eq 5 reduces to the homogeneous PY closure

c˜ (r) ) (h˜ (r) + 1)(1 - exp(βφ(r)))

(9)

Because of the first assumption mentioned above, eqs 8 and 9 are uncoupled from the fluid density profile created by the wall and can be solved independently to yield total and direct isotropic correlation functions of a uniform fluid. The wall-induced fluid density profile, F(z), is treated in the homogeneous integral equation theory as a limiting case of a

solute-solvent pair distribution function in an infinitely dilute solution, where the size of the solute particle becomes infinitely large, i.e., the solute becomes a flat wall. Accordingly, one writes the corresponding wall-particle OZ equation

∫

hw(z) ) cw(z) + F db′h r w(z′)c˜ (|b r - b′) r

(10)

where hw(z) and cw(z) are the wall-particle total and direct correlation functions, respectively. The fluid density profile at the wall can be calculated directly from hw(z) according to F(z) ) F(hw(z) + 1). To solve eq 10, one needs a closure relation between hw(z) and cw(z). We note here (and will return to this point below) that this closure need not be of the same type as the one used in conjunction with the OZ relation for the neat fluid (eq 8). In other words, closure relations other than PY can be employed, such as the hypernetted chain (HNC) closure38

ln[hw(z) + 1] ) -βV(z) + hw(z) - cw(z)

(11)

The singlet (homogeneous) HNC result for the fluid density profile at the wall can be obtained by combining eq 11 with the wall-particle OZ relation eq 10

ln[F(z)/F] ) -βV(z) +

r - F]c˜ (|b r - b′|) r ∫db′[F(z′)

(12)

Incidentally, eq 12 follows immediately from the exact LMBW relation (eq 6) proVided in the latter the anisotropic direct correlation function, c(r b1,r b2), is replaced by its isotropic counterpart, c˜ (r12). It is worth emphasizing that the function c˜ (r) entering eq 12 is obtained by solving eqs 8 and 9 and is, therefore, identical to the direct correlation function of the neat fluid in the absence of the wall. In other words, the homogeneous integral equation approach does not account for the fact that the correlation function for the two fluid particles in the vicinity of the wall can be different from the corresponding pair correlation function in the bulk. By contrast, the inhomogeneous fluid theory operates with anisotropic fluid pair correlations, which depend explicitly on the positions of two particles relatiVe to the wall, thereby accounting for the effect of the inhomogeneous environment induced by the wall on the pair correlation functions. Concomitantly, the relation between the singlet density profile at the wall and the fluid pair correlations utilized in the inhomogeneous approach (such as LMBW equation or other equivalent relations) is exact, whereas the corresponding relation in the uniform fluid theory (such as the singlet HNC result or relations arising from other closures) is approximate. From the above discussion, one would expect the inhomogeneous approach to produce more accurate results for the fluid density profile at the wall in comparison to its homogeneous counterpart. This is confirmed in the next section, where the results from two approaches are compared with MC simulations. Before proceeding to the comparison between theory and simulation, we will briefly outline the generalization of the inhomogeneous integral equation formalism to the case of fluid mixture adsorption. To this end, we consider an n-component mixture in contact with a flat wall. We assume that the fluid particles of species i and j (i,j ) 1, 2, ‚‚‚, n) interact with each other via isotropic pair potential, φij(r); the wall-particle interaction for an adsorbate particle of species i located at a distance z from the wall is taken to be Vi(z). The inhomogeneous OZ equation for two fluid particles of r2 reads species i and j located at b r1 and b

6586 J. Phys. Chem. B, Vol. 105, No. 28, 2001

Egorov

n

hij(b r 1,b r 2) ) cij(b r 1, b r 2) +

∑ ∫dbr 3Fk(br 3)hik(br 1,br 3)ckj(br 3,br 2) k)1

form, Steele has shown that it can be very accurately approximated by the following expression39

(13)

where Fi(r b) is the density profile of the species i, and hij(r b1,r b2) and cij(r b1,r b2) are the anisotropic total and direct pair correlation functions for fluid particles of type i and j, respectively. The generalization of the inhomogeneous PY closure for a multicomponent mixture yields

cij(b r 1,b r 2) ) [hij(b r 1, b r 2) + 1][1 - exp(βφij(r12))] (14) whereas the LMBW equation takes the form n

∑ ∫0 dz2F′k(z2)cˆ ik(z1,z2,0) k)1

F′i(z1)/Fi(z1) ) -βV′i(z1) +

∞

(15)

The system of eqs 13-15 can be solved in exactly the same fashion as outlined above for the case of the neat fluid. These equations will be employed in section IV to perform model calculations of competitive adsorption from a fluid mixture. In the next section, we return to the problem of neat fluid adsorption and compare the results from the inhomogeneous integral equation theory with computer simulations. III. Comparison of Theory with Simulation We consider a neat fluid in contact with a solid substrate. We assume that the particles of the adsorbate interact with each other via the familiar LJ potential

φ(r) ) 4

[(σr ) - (σr ) ] 12

6

(16)

where  is the potential well depth and σ is the effective diameter of the fluid particle. The pair potential for the interaction between a fluid particle and an atom of the solid substrate has the same form but with different  and σ parameters

φfs(r) ) 4fs

[( ) ( ) ] σfs r

12

-

σfs r

6

(17)

To obtain the interaction energy of a given fluid particle with the entire solid substrate, one needs to evaluate the sum ∑iφfs(ri), where ri is the distance between the selected atom of the adsorbate and the ith atom of the substrate, and the summation index i runs over all the atoms of the substrate. This sum can be calculated by dividing the substrate into two-dimensional layers parallel to the surface and summing over the interactions of the adsorbate particle with the atoms in each layer in the layer’s reciprocal space.39 For the model based on the LJ potentials, such calculation has been carried out by Steele.39 The leading term of Steele’s result (corresponding to a flat layer, i.e., no corrugation) gives the following form for the interaction energy of an adsorbate particle with the substrate ∞

V(z) ) 2πFsσ2fsfs

∑

n)0

[ ( ) ( )] 2

σfs

5 z + n∆

10

-

σfs

z + n∆

4

(18)

where z is the distance between the adsorbate atom and the solid surface, ∆ is the distance between the neighboring substrate layers, the summation index n runs over the substrate layers, and Fs is the two-dimensional number density of atoms in each layer. Although the sum in eq 18 cannot be evaluated in closed

V(z) ) 2πFsσ2fsfs

[( ) ( ) 2 σfs 5 z

10

-

σ4fs σfs 4 z 3∆(z + 0.61∆)3

]

(19)

The 10-4-3 wall-fluid potential given by eq 19 has been employed in several simulation studies of supercritical adsorption.15-17 In the present work, we focus on the simulation results obtained by Vermesse et al.16 who have performed grand canonical MC simulations of the density profiles of supercritical argon at a graphite wall for several values of argon density along a particular supercritical isotherm. The following values of LJ parameters for argon were used:  ) 119.8 K and σ ) 3.405 Å. The parameters fs and σfs in the wall-fluid 10-4-3 potential were obtained from the standard Lorentz-Berthelot combining rules: fs ) xs and σfs ) xσσs, where s and σs are the potential well-depth and the effective diameter of the atoms comprising the solid substrate. Using the values s ) 28 K and σs ) 3.400 Å for carbon, gives the following parameters for the argon-graphite 10-4-3 potential: fs ) 57.9 K, σfs ) 3.403 Å.16 The remaining two parameters in this potential correspond to the distance between the neighboring graphite basal planes (∆ ) 3.348 Å) and the carbon atom number density in the graphite basal plane Fs ) 0.382 Å-2). We remark in passing that Vermesse et al. have performed two sets of simulations, one with the argon-graphite potential given by eq 19 with parameters listed above, and another with all interparticle interactions between argon atoms in the fluid and carbon atoms in the substrate (modeled by the LJ pair potentials with the same parameters) explicitly evaluated. The results from the two sets of simulations were essentially equivalent.16 To assess the accuracy of the inhomogeneous integral equation theory in reproducing the wall-induced fluid density profile, we have employed this theory to calculate the density profiles of supercritical argon at the graphite wall at the same thermodynamic conditions for which the simulation results of Vermesse et al.16 have been reported. To present our results, we define dimensionless density and temperature according to: T* ) kBT/ and F* ) Fσ3, and the dimensionless distance variable by z* ) z/σ. We have calculated the argon density profiles along the supercritical isotherm T* ) 2.49 at the following six densities: F* ) 0.051, 0.161, 0.313, 0.519, 0.764, and 0.964. Note that this range of densities spans from a dilute gas to a dense fluid (the LJ triple point density is F/t = 0.85). The density profiles were obtained from the inhomogeneous integral equation theory by solving eqs 2, 5, and 6 as described in the previous section with the potentials given by eqs 16 and 19 with parameters listed above. The corresponding results for F(z) together with the simulation results of Vermesse et al.16 are shown in Figures 1 and 2. For comparison, we also present the homogeneous (singlet) theory results for F(z), as given by eqs 8, 9, and 12. One sees that at all thermodynamic conditions the inhomogeneous theory is in excellent agreement with the simulation, whereas its homogeneous counterpart significantly overestimates the height of the first peak of F(z), except at the lowest density. This situation is very similar to the one observed for the solutesolvent radial distribution functions in the case of dilute attractive supercritical solutions.27 As discussed in the Introduction, this similarity is not accidental because both an adsorbing wall and an attractive solute create strongly inhomogeneous fluid density profiles, and the uniform fluid theory fails to properly

Inhomogeneous Integral Equation Study

Figure 1. Density profiles for argon near graphite wall (basal plane) at three different densities at T* ) 2.49. The open symbols are MC simulation results the solid line is from the inhomogeneous integral equation theory, the dot-dashed line is from the homogeneous integral equation theory.

account for the effect of these local inhomogeneities on the fluid pair correlation functions. To illustrate that this influence can be significant, we show in Figure 3 the anisotropic direct correlation function c(z1,z2,s) calculated from the inhomogeneous integral equation theory at F* ) 0.161 and T* ) 2.49. In plotting this function, we have fixed the values of its first two arguments, i.e., the distances of the two fluid particles from the wall. Specifically, we have set z/1 ) z/2 ) 1.1, which roughly corresponds to the first peak of the fluid density profile at the wall. By varying the variable s, we plot c(z/1 ) 1.1, z/2 ) 1.1, s*) as a function of the separation between the two fluid particles. Also shown in Figure 3 is the isotropic direct solvent correlation function, c˜ (s), obtained from the homogeneous PY theory (eqs 8 and 9) for the same thermodynamic conditions of the fluid. Clearly, the two functions differ substantially at small distances (the magnitude of c˜ (0) is approximately 60% smaller than c(z/1 ) 1.1, z/2 ) 1.1,0)), which demonstrates the pronounced effect of the substrate on the fluid pair correlations. In fact, the anisotropic correlation function is remarkably similar to the isotropic one calculated at the same temperature but at a significantly higher bulk fluid density, F* = 0.5, which is also shown in Figure 3. This is hardly surprising because the local fluid density in the vicinity of an adsorbing wall is expected to be higher than the bulk value. As one would expect, for larger values of the separation between the adsorbate atoms and the wall, the anisotropic correlation function approaches the isotropic one calculated at the actual value of the bulk density. Very similar trends are observed in the behavior of the anisotropic pair distribution function, g(z1,z2,s) ) h(z1,z2,s) + 1, also shown in

J. Phys. Chem. B, Vol. 105, No. 28, 2001 6587

Figure 2. Density profiles for argon near graphite wall (basal plane) at three different densities at T* ) 2.49. The open symbols are MC simulation results the solid line is from the inhomogeneous integral equation theory, the dot-dashed line is from the homogeneous integral equation theory.

Figure 3. One sees that g(z/1 ) 1.1, z/2 ) 1.1, s) differs substantially from the isotropic PY result g˜ (s) ) h˜ (s) + 1 evaluated at the bulk fluid density F* ) 0.161, and at the same time is quite close to g˜ (s) obtained at F* = 0.5. Because the above results indicate that the local fluid density in the vicinity of the wall exceeds the bulk value substantially, it is of interest to obtain some estimate for this local density. It is clear from Figures 1 and 2 that simply taking the value of the actual density profile F(z) by itself is not suitable for determining the local density, since the peak heights of F(z) correspond to density values in excess of those characteristic of close packing. Instead, one can define a coarse-grained local density profile, Fj(z), which is obtained as a weighted average of the actual profile over a microscopic volume determined by the range of the intermolecular potential40

F(z) )

r b r - b′|) r ∫db′F(z′)w(|

(20)

where w(r) is the weight function normalized according to ∫dr bw(r) ) 1. For the model based on LJ potentials, it is convenient to define w(r) in terms of the repulsive part of the fluid pair potential40

exp(-βφrep(r)) - 1

w(r) ) 4π

∫0 drr2(exp(-βφrep(r)) - 1) ∞

(21)

where φrep(r) ) φ(r) +  for r e 21/6σ, and φrep(r) ) 0 otherwise. The above form of the weight function satisfies the condition w(r) g 0, which guarantees that the weighted local density always remains positive.

6588 J. Phys. Chem. B, Vol. 105, No. 28, 2001

Egorov well as for the fluid correlations), which yields the singlet PY result for the density profile

F(z) ) Fexp(-βV(z)) [1 +

r - F]c˜ (|b r - b′|)] r ∫db′[F(z′)

(22)

As expected, no significant improvement over singlet HNC results was observed. In view of this, from now on we will restrict ourselves to the inhomogeneous integral equation theory. IV. Comparison of Theory with Experiment and Model Calculations

Figure 3. (a) Direct pair correlation function for two argon atoms near graphite wall at F* ) 0.161, T* ) 2.49. Solid line is from the inhomogeneous PY theory; dot-dashed line is from the homogeneous PY theory at the actual value of the bulk density, the dashed line is from the homogeneous PY theory at a different (higher) value of the bulk density, as indicated. (b) Radial distribution function function for two argon atoms near graphite wall at F* ) 0.161, T* ) 2.49. Solid line is from the inhomogeneous PY theory; dot-dashed line is from the homogeneous PY theory at the actual value of the bulk density, the dashed line is from the homogeneous PY theory at a different (higher) value of the bulk density, as indicated. (c) Weighted density profile for argon near graphite wall at F* ) 0.161, T* ) 2.49; the dashed line indicates the value of the bulk density.

We have calculated Fj(z) at F* ) 0.161, T* ) 2.49 from eqs 20 and 21, using as input the density profile F(z) obtained from the inhomogeneous integral equation theory; the result is also shown in Figure 3. One sees that the peak of the weighted local density is located around z* ) 1.2, the maximum value attained by Fj*(z) is approximately 0.45, which exceeds the bulk fluid density by nearly a factor of 3. As mentioned above, the anisotropic pair correlation functions for two fluid particles located in the plane z* ) 1.1 are very similar to the isotropic pair correlations evaluated at the bulk density F* = 0.5, which is consistent with the estimate of the local density obtained from the weighted density profile. The results presented in Figures 1 and 2 clearly demonstrate that the accuracy of the inhomogeneous integral equation theory in reproducing the wall-induced fluid density profiles is superior compared to its homogeneous counterpart. As already mentioned, one of the major deficiencies of the uniform fluid theory is its underlying assumption that the pair correlation functions for the fluid particles in the vicinity of the wall are identical to those in the bulk. This assumption is independent of a particular closure relation between the total and direct wall-particle correlation functions and, therefore, choosing a different closure (in place of the HNC relation given by eq 12) is not likely to improve the accuracy of the theory significantly. To confirm this, we have calculated F(z) from the homogeneous theory utilizing the PY closure for the wall-particle correlations (as

SCF adsorption has been studied by several experimental groups.11-14 In the present work, we will analyze the data obtained by Findenegg et al., who have measured adsorption isotherms for supercritical argon and krypton on graphitized carbon black.11,12 These adsorbent-adsorbate combinations are particularly well-suited for theoretical analysis with the model presented in the previous section. First, LJ potentials are appropriate for describing inert-gas adsorbates. Second, graphitized carbon black is a nonporous adsorbent with highly homogeneous, molecularly smooth surface, which justifies the use of Steele’s 10-4-3 form for the wall-fluid interaction potential. The surface excess amount of the adsorbed substance, Γ, is defined as the difference between the actual number of the adsorbate particles per unit surface area of the adsorbent and the number of particles that would be present if the fluid density were uniform up to the substrate-fluid dividing surface (Gibbs surface). For a model adsorbent with perfectly homogeneous flat surface, the surface excess is related to the fluid density profile along the direction normal to the substrate surface by

Γ)

∫0∞[F(z) - F]dz

(23)

where the plane of the surface atoms was chosen as the Gibbs dividing surface.15 To calculate the density profiles F(z) for the experimental systems studied by Findenegg et al., we need to specify the wall-fluid and fluid-fluid interaction potential parameters. We choose the latter by requiring that the reduced experimental critical density and temperature of the adsorbate coincide with the corresponding critical constants of a LJ fluid. The experimental values for the critical parameters of the adsorbates are as follows: Fc ) 13.41 mol/L and Tc ) 150.7 K for argon, and Fc ) 10.92 mol/L and Tc ) 209.3 K for krypton. The best estimates for the critical density and temperature of a LJ fluid are F/c ) 0.316 and T/c ) 1.312.41 Using these numbers, one obtains /kB ) 115 K and σ ) 3.395 Å for argon, and /kB ) 160 K and σ ) 3.636 Å for krypton (note that the parameters for argon obtained through this procedure differ slightly from the Ar parameters used by Vermesse et al. in their simulations). As for the wall-fluid potential parameters fs and σfs, we obtain these from the standard combining rules, using LJ parameters for carbon given in the previous section. This procedure yields fs/kB ) 57 K and σfs ) 3.398 Å for argon, and fs/kB ) 67 K and σfs ) 3.518 Å for krypton. The remaining two parameters in V(z) (∆ and Fs) characterize the substrate; the values of these parameters appropriate for graphite were given in the previous section. With the specified interaction potentials, we calculate the density profiles for argon and krypton at the graphite wall from the inhomogeneous integral equation theory, and then obtain the surface excess for these adsorbates from eq 23. Findenegg et al. have measured surface excess of argon and krypton on graphite as a function of pressure (density) along

Inhomogeneous Integral Equation Study

Figure 4. Surface excess of krypton on graphite as a function of density at three different temperatures. Symbols are the experimental results and lines are from the inhomogeneous integral equation theory.

several supercritical isotherms. Experimental values of Γ for krypton are reproduced in Figure 4 for the following three isotherms: T ) 253, 273, and 298 K, which correspond to the reduced temperatures Tr ) T/Tc ) 1.21, 1.30, and 1.42, respectively. Also shown in Figure 4 are theoretical results for Γ obtained from eq 23. One sees that agreement between theory and experiment is good, except for the high-density portion of the lowest temperature isotherm; it is worth emphasizing that theoretical calculations do not involve any adjustable parameters. In Figure 5 we present a comparison between experimental and theoretical results for the surface excess of argon on graphite for two supercritical isotherms, T ) 253 and 298 K, which correspond to the reduced temperatures Tr ) T/Tc ) 1.68 and 1.97, respectively. Once again, theoretical results for Γ are in good agreement with experiment. The results for the surface excess shown in Figures 4 and 5 illustrate the general trends in the density and temperature dependence of Γ discussed in the Introduction. Namely, in the cases where the experimental measurements cover sufficient density range, Γ as a function of density passes through a maximum somewhat below the adsorbate’s critical density; this maximum becomes more pronounced as the temperature is lowered. We note that the range of densities and temperatures studied experimentally is somewhat limited, the density range does not extend beyond the critical density, and the lowest reduced temperature is Tr ) 1.21. In view of this, we have performed model calculations of Γ for a wider range of thermodynamic parameters. Specifically, we have calculated the surface excess of krypton on graphite for F* extending from 0.01 up to 0.5 for the following three temperatures: T* ) 1.37, 1.59, and 2.50, which correspond to the reduced temperatures Tr ) 1.04, 1.21, and 1.91, respectively. The results of model calculations are displayed in Figure 6, where the dimensionless surface excess Γ* ) Γσ2 is plotted vs fluid density. One sees that as the critical temperature of the adsorbate is approached from the above, the height of the maximum in Γ(F) increases rapidly, and the location of the maximum moves toward the critical density.

J. Phys. Chem. B, Vol. 105, No. 28, 2001 6589

Figure 5. Surface excess of argon on graphite as a function of density at three different temperatures. Symbols are the experimental results and lines are from the inhomogeneous integral equation theory.

Figure 6. Theoretical results for surface excess of krypton and argon on graphite as a function of density at various temperatures.

Turning now to the origins of the density and temperature behavior of the surface excess isotherms, we point out that the general shape of Γ as a function of F is remarkably similar to the density dependence of the excess local density around a dilute attractive solute in a supercritical solvent. The similarity between the SCF adsorption and clustering in supercritical solutions has been already mentioned earlier, it essentially reflects the fact that an adsorbing wall acts as a giant attractive solute. Accordingly, the arguments involved in discussing the density and temperature dependence of excess local density

6590 J. Phys. Chem. B, Vol. 105, No. 28, 2001 around an attractive solute26 can be largely transferred to the case of SCF adsorption (one important difference between the two cases is due to the fact that the excess local density is generally defined for the first coordination shell of the solute and therefore reflects only the short-range part of the solute-solvent radial distribution function, whereas Γ, as defined in eq 23, probes the entire range of the wall-induced fluid density profile). We start by discussing the low-density behavior of the adsorbate surface excess. Because of the strong attractive wallfluid interaction, the local density of the adsorbate near the wall increases with F in the low-density region more rapidly than the bulk density, which results in the initial increase of the surface excess. The plots of Γ(F) given in Figures 4-6 indicate that the rate of this initial increase is quite sensitive to temperature, increasing rapidly when T is lowered. The above behavior can be understood as follows. At very low densities, the fluid density profile reduces to the Boltzmann form: F(z) ) Fexp[-βV(z)]. Substituting this expression into eq 23, one sees that Γ increases linearly with the bulk density, and the slope of this linear dependence is given by the spatial integral over the function exp[-βV(z)] - 1. The Boltzmann factor is peaked at the minimum of V(z), and the height of the peak grows exponentially with the ratio jfs/kBT, where jfs ) 2πFsσ2fsfs. Hence, the lower the temperature, the steeper the initial density dependence of Γ for a given adsorbate-adsorbent combination, which is indeed observed in Figures 4-6. From the above discussion, it follows that for two different adsorbate-adsorbent pairs the low-density behavior of Γ(F) can be expected to be very similar proVided the value of the parameter jfs/kBT (which characterizes the wall-fluid interaction strength relative to the thermal energy kBT) is the same for both systems. To confirm this conjecture, we have calculated the surface excess isotherms for krypton-graphite and argon-graphite systems at jfs/kBT ) 7.8, which corresponds to T* ) kBT/Kr ) 1.59 for Kr, and T* ) kBT/Ar ) 1.75 for Ar. The two isotherms are shown in Figure 6, and one sees that they are indeed very close to each other in the low-density region. Referring to the similarity between SCF adsorption and clustering in supercritical solutions, we note that correlations between the solvent density augmentation around an attractive solute and various measures of the solute-solvent interaction strength (relative to the thermal energy kBT) have been extensively discussed by Maroncelli and co-workers.22 As one proceeds to higher densities, the bulk value of F starts catching up with the local density near the substrate wall, since the repulsive interactions between the adsorbate particles prevent further rapid increase of the local fluid density at the surface. As a result, Γ(F) curve passes through a maximum and gradually decreases. We also note that the adsorption isotherms obtained for krypton-graphite and argon-graphite pairs at the same value of jfs/kBT deviate from each other at intermediate densities, with the maximum in argon surface excess being lower compared to krypton. This behavior can be attributed to the fact that at higher densities the amount of the adsorbed fluid is determined not only by the strength of the wall-fluid interaction but also by the strength of interaction between the particles of the fluid itself, which is determined by the LJ potential well-depth of the fluid relative to the thermal energy. For a given value of jfs/kBT, the ratio Kr/kBT is larger than Ar/kBT, which accounts for the higher maximum in the Γ(F) curve for krypton relative to argon. Up to this point, we have restricted our attention to the adsorption of neat SCFs. However, as already mentioned in the Introduction, it is also very important to study adsorption of SCF mixtures. In particular, it is crucial to understand what

Egorov

Figure 7. Theoretical adsorption isotherms (T ) 298 K) for the argonkrypton equimolar mixture and for pure argon and pure krypton (for pure fluids the surface excess is multiplied by one-half). Solid line gives the adsorption isotherm for krypton from the mixture, and shortdashed line gives the adsorption isotherm for pure krypton. Long-dashed line gives the adsorption isotherm for argon from the mixture, and dotdashed line gives the adsorption isotherm for pure argon.

factors affect the degree of selectivity in the adsorption of a given component from the mixture because preferential adsorption forms the basis for adsorptive separation processes. Once again, one can draw an analogy with supercritical solvation in dilute attractive mixtures. For an attractive solute immersed in a multicomponent supercritical mixture, the local density augmentation around the solute is accompanied by the local composition enhancement, which means that the local mole fraction of a particular solvent species in the vicinity of the solute can be significantly higher than in the bulk.42-45 We have recently employed the inhomogeneous integral equation theory to study the phenomenon of preferential solvation and found that the local composition enhancement for a given species directly correlates with the strength of the solute-solvent interaction.46 As will be seen below, a similar trend holds in the case of preferential adsorption. To study the adsorption of SCF mixtures, we use the corresponding generalization of the inhomogeneous integral equation theory given by eqs 13-15. We apply this theory to analyze the adsorption of an equimolar mixture of argon and krypton at a graphite wall. The interaction between Ar and Kr atoms is modeled via LJ potential with the parameters obtained from the standard Lorentz-Berthelot combining rules. We calculate the surface excess of argon and krypton as a function of the bulk density of the mixture at T ) 298 K. The results of our calculations are shown in Figure 7 together with the surface excess isotherms for neat Ar and Kr at the same temperature (the latter are multiplied by a factor of 0.5 corresponding to the mole fraction of each component in the mixture). One sees that the surface excess of krypton in the mixture generally exceeds the value of Γ(F)/2 for neat krypton (except at very low densities). Conversely, the surface excess of Ar in the mixture is consistently smaller than Γ(F)/2 of neat argon. This indicates that outside the low-density region the density profiles of the two components in the mixture are not the same as the density

Inhomogeneous Integral Equation Study profiles of the respective neat fluids (multiplied by the appropriate mole fractions). It is also clear that there is a pronounced selectivity for the adsorption of krypton from the mixture over the entire density range. This is not surprising, since the parameter fs which governs the strength of the wallfluid interaction is significantly larger for the krypton-graphite jAr pair (jKr fs ) 1986 K) than for the argon-graphite one ( fs ) 1572 K). We close the discussion of preferential adsorption by noting that simulation studies of competitive adsorption available in the literature are primarily concerned with the adsorption of mixtures in slit pores rather than on a single flat surface.14,47,48 In the case of slit pores, the degree of selectivity is determined not only by energetics, but also by geometric factors, such as the relative size of the pores and the adsorbate particles. In this connection, we remark that the inhomogeneous integral equation theory can be easily adapted to calculate the density profile of a fluid in a slit pore or other confined geometry.30 The application of this theory to study the effect of geometric factors on preferential adsorption of SCF mixtures in pores of various form will be the subject of our future research. V. Conclusion In this paper, we have presented a microscopic statistical mechanical treatment of adsorption of supercritical fluids and fluid mixtures on solid surfaces. The inhomogeneous integral equation theory was employed to calculate the wall-induced nonuniform fluid density profile, from which the surface excess of the adsorbate was obtained. The accuracy of the inhomogeneous fluid approach was tested by comparing theoretical results for a particular adsorbate-adsorbent pair (argon on graphite) with simulations. The density profiles obtained from the inhomogeneous integral equation theory were found to be in excellent agreement with the simulation throughout the entire density range spanning from a dilute gas to a dense fluid. It was shown that the nonuniform fluid approach provides a substantial improvement over the theory formulated for uniform fluids due to the fact that the former approach explicitly takes into account the effect of the adsorbing wall on the fluid pair correlation functions. This effect can be very significant, since the local density of the adsorbate in the vicinity of the wall can differ greatly from the bulk. Equipped with an accurate theory of the wall-induced fluid density profile, we have analyzed the experimental data on the surface excess of supercritical argon and krypton on graphite. Theoretical results for the surface excess were obtained without any adjustable parameters and were shown to be in good agreement with experiment. In addition to the adsorption of neat SCFs, we have also considered adsorption of fluid mixtures with particular emphasis on the phenomenon of preferential adsorption. Model calculations were performed for an equimolar argon-krypton mixture in contact with a graphite wall, and it was found that the adsorbed phase is significantly enriched in krypton relative to the bulk. This can be explained by the energetic factors because the interaction strength for krypton-graphite pair is larger compared to argon-graphite pair. We note that the degree of selectivity characterizing the preferential adsorption of a given component from the mixture can also be affected by geometric factors, such as the size of pores of the adsorbent surface. Although the present study was limited to adsorption at a single flat wall, the inhomogeneous integral equation theory can easily be adapted to calculate fluid density profiles in various confined geometries, e.g., in narrow slit pores. The study of preferential

J. Phys. Chem. B, Vol. 105, No. 28, 2001 6591 adsorption under such conditions will be the subject of our future work. Acknowledgment. It is a great pleasure to dedicate this paper to Bruce Berne on the occasion of his 60th birthday. The author is grateful for financial support from the Chemistry Department of the University of Virginia. References and Notes (1) Kiran, E., Brennecke, J. F., Eds. Supercritical Fluid Engineering Science, Number 514 in ACS Symposium Series; American Chemical Society: Washington, DC, 1993. (2) Kiran, E., Sengers, J. M. H. L., Eds. Supercritical Fluids: Fundamentals for Applications; Kluwer: Dordrecht, 1994. (3) Kauffman, J. Anal. Chem. 1996, 68, 248. (4) Johnston, K. P.; Haynes, C. AIChE J. 1987, 33, 2017-2026. (5) Chester, T. L.; Pinkston, J. D.; Raynie, D. E. Anal. Chem. 1998, 70, 301R. (6) Barton, S. S.; Dacey, J. R.; Quinn, D. F. In Fundamentals of Adsorption; Myers, A. L., Belfort, G., Eds.; Engineering Foundation: New York, 1984; p 65. (7) Ju¨ntgen, H. Carbon 1977 15, 273-283. (8) Brady, B. O.; Kao, C. P. C.; Dooley, K. M.; Knopf, F. C.; Gambrell, R. P. Ind. Eng. Chem. Res. 1987, 26, 261-268. (9) Tan, C. S.; Liou, D. C. Ind. Eng. Chem. Res. 1988, 27, 988-991. (10) Strubinger, J. R.; Parcher, J. F. Anal. Chem. 1989, 61, 951. (11) Specovius, J.; Findenegg, G. H. Ber. Bunsen-Ges. Phys. Chem. 1978, 82, 174. (12) Blu¨mel, S.; Ko¨ster, F.; Findenegg, G. H. J. Chem. Soc., Faraday Trans. 2 1982, 78, 1753. (13) Malbrunot, P.; Vidal, D.; Vermesse, J.; Chahine, R.; Bose, T. K. Langmuir 1992, 8, 577-580. (14) Gusev, V. Y.; O’Brien, J. A. Langmuir 1998, 14, 6328-6331. (15) van Megen, W.; Snook, I. K. Mol. Phys. 1982, 45, 629. (16) Vermesse, J.; Levesque, D. J. Chem. Phys. 1994, 101, 9063. (17) Nitta, T.; Shigetta, T. Fluid Phase Equil. 1998, 144, 245. (18) Kajimoto, O.; Futakami, M.; Kobayashi, T.; Yamasaki, K. J. Phys. Chem. 1988, 92, 1347-1352. (19) Flanagin, L. W.; Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1997, 101, 7998-8005. (20) Subramanian, R.; Pyada, H.; Lira, C. T. Ind. Eng. Chem. Res. 1995, 34, 838. (21) Tucker, S. C. Chem. ReV. 1999, 99, 391. (22) Song, W.; Biswas, R.; Maroncelli, M. J. Phys. Chem. A 2000, 104, 6924-6939. (23) Zhang, J.; Lee, L. L.; Brennecke, J. F. J. Phys. Chem. 1995, 99, 9268. (24) Flanagin, L. W.; Balbuena, P. B.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. 1995, 99, 5196. (25) Adams, J. E. J. Phys. Chem. B 1998, 102, 7455. (26) Egorov, S. A.; Yethiraj, A.; Skinner, J. L. Chem. Phys. Lett. 2000, 317, 558. (27) Egorov, S. A. J. Chem. Phys. 2000, 112, 7138. (28) Henderson, D. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; Dekker: New York, 1992; page 177. (29) Attard, P. J. Chem. Phys. 1989, 91, 3072. (30) Kjellander, R.; Sarman, S. Mol. Phys. 1990, 70, 215. (31) Plischke, M.; Henderson, D. J. Chem. Phys. 1986, 84, 2846. (32) Lee, L. L.; Cochran, H. D. J. of Supercr. Fluids 1998, 13, 77-81. (33) Trizenberg, D. G.; Zwanzig, R. Phys. ReV. Lett. 1972, 28, 1183. (34) Lovett, R. A.; Mou, C. Y.; Buff, F. P. J. Chem. Phys. 1976, 58, 1880. (35) Wertheim, M. S. J. Chem. Phys. 1976, 65, 2377. (36) Born, M.; Green, H. S. Proc. R. Soc. London Ser. A 1946, 188, 10. (37) Lado, F. J. Comput. Phys. 1971, 8, 417-433. (38) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic: London, 1986. (39) Steele, W. Surf. Sci. 1973, 36, 317. (40) Zhou, Y.; Stell, G. J. Chem. Phys. 1990, 92, 5544-5550. (41) Potoff, J. J.; Panagiotopoulos, A. Z. J. Chem. Phys. 1998, 109, 10 914. (42) Kim, S.; Johnston, K. P. AIChE J. 1987, 33, 1603-1611. (43) Yonker, C. R.; Smith, R. D. J. Phys. Chem. 1988, 92, 235-238. (44) Zhang, J.; Roek, D. P.; Chateauneuf, J. E.; Brennecke, J. F. J. Am. Chem. Soc. 1997, 119, 9980. (45) Schwarzer, D.; Troe, J.; Zerezke, M. J. Phys. Chem. A 1998, 102, 4207. (46) Egorov, S. A. J. Chem. Phys. 2000, 113, 7502. (47) Cracknell, R. F.; Nicholson, D. J. Chem. Soc., Faraday Trans 1994, 90, 1487. (48) Nicholson, D.; Gubbins, K. E. J. Chem. Phys. 1996, 104, 8126.