Adsorption of Supercritical Fluids on Graphitised Thermal Carbon

Langmuir 2003, 19, 2215-2225

2215

Adsorption of Supercritical Fluids on Graphitised Thermal Carbon Black: Molecular Layer Structure Theory versus Grand Canonical Monte Carlo Simulation D. D. Do,* H. D. Do, and E. Ustinov† Department of Chemical Engineering, University of Queensland, St. Lucia, Queensland 4072, Australia Received September 16, 2002. In Final Form: December 5, 2002 This paper presents a detailed analysis of adsorption of supercritical fluids on nonporous graphitized thermal carbon black. Two methods are employed in the analysis. One is the molecular layer structure theory (MLST), proposed recently by our group, and the other is the grand canonical Monte Carlo (GCMC) simulation. They were applied to describe the adsorption of argon, krypton, methane, ethylene, and sulfur hexafluoride on graphitized thermal carbon black. It was found that the MLST describes all the experimental data at various temperatures well. Results from GCMC simulations describe well the data at low pressure but show some deviations at higher pressures for all the adsorbates tested. The question of negative surface excess is also discussed in this paper.

1. Introduction Supercritical adsorption on open surfaces and porous media has been an interesting problem because of the increasing importance of high-pressure storage of highvalue energy substances such as hydrogen and methane. Supercritical adsorption is different from that of subcritical fluids in a number of ways. First of all, adsorption of supercritical fluids generally forms a monolayer coverage near the solid surface. Second, the adsorption of supercritical fluids when plotted as surface excess (to be defined later) versus pressure will exhibit a maximum at a pressure close to the critical pressure, especially when the temperature is slightly greater than the critical temperature the maximum behaves like a very sharp spike. For example, the adsorption of propane on graphitized thermal carbon black at 371 K1 exhibits such a spike. The critical temperature of propane is 369.8 K. At temperatures lower than 363 K, the adsorption behavior is typical of that for subcritical fluids on nonporous surfaces (type II isotherm according to BDDT classification).2 Similar behavior is also observed for ethylene adsorption on the same material. The sharp spike is observed at a temperature of 283 K, which is slightly greater than the critical temperature of ethylene of 282.4 K. Isotherms at temperatures less than 283 K have type II isotherm behavior, while those for temperatures greater than that exhibit a maximum. This maximum is smaller in magnitude when the temperature increases and the pressure at which this maximum occurs shifts to higher pressure. The maximum in the adsorption isotherm (plotted as surface mass excess versus pressure) has long been recognized; for example, see Menon.3 Subsequently, many * To whom all correspondence should be addressed. Fax: +617-3365-2789. E-mail: [email protected]. † On leave from Saint Petersburg State Technological Institute (Technical University), 26 Moskovsky Prospect, St. Petersburg 198013, Russia. (1) Findenegg, G. H. In Fundamentals of Adsorption; Myers and Belfort, Eds.; Engineering Foundation: New York, 1983; pp 207-218. (2) Gregg, S. J.; Sing, K. S. W. Adsorption, Surface Area and Porosity, 2nd ed.; Academic Press: London, 1982. (3) Menon, P. G. Chem. Rev. 1968, 68, 277.

authors4-14 have also observed this maximum. Various theories have been put forward to explain the maximum behavior. Table 1 summarizes some of these works in the area of supercritical adsorption. As the Table 1 shows, there are a number of theoretical tools proposed to study the problem. It ranges from the classical approach of using standard equation (Langmuir, Freundlich) with the notion of either constant or variable adsorbed phase volume to sophisticated tools such as the DFT theory and the grand canonical Monte Carlo (GCMC) simulation. They all describe the maximum behavior, but it is not clear about the behavior of the surface excess when the pressure is very high. Past work has shown that the surface excess could be negative9 and the GCMC simulations performed by Vermesse and Levesque34 also show that the surface excess once passing through a maximum will go through a minimum and then modestly increase with pressure once again. The negative surface excess was later corrected by Malbrunot et al.35 by using the correct solid density as they pointed out that incorrect use of the density could lead to the negative surface excess. In this paper, we address the same problem by tackling it with two advanced theoretical tools. One is the molecular layer structure theory (MLST) proposed by Ustinov and Do36 and Do et al.37 to study adsorption on surfaces and surface tension as a function of temperature. The other tool is the now well-established grand canonical Monte Carlo simulations. Using these tools, we will study again (4) Ozawa, S.; Kusumi, S.; Ogino, Y. J. Colloid Interface Science 1976, 56, 83. (5) Specovius, J.; Findenegg, G. Ber. Bunsen-Ges. Phys. Chem. 1978, 82, 174. (6) Wakasugi, Y.; Ozawa, S.; Ogino, Y. J. Colloid Interface Sci. 1981, 79, 399. (7) Blumel, S.; Foster, F.; Findenegg, G. J. Chem. Soc., Faraday Trans. 2 1982, 78, 1753. (8) Findenegg, G.; Korner, B.; Fischer, J.; Bohn, M. Ger. Chem. Eng. 1983, 6, 80. (9) Malbrunot, P.; Vidal, D., Vermesse, J.; Chahine, R.; Bose, T. Langmuir 1992, 8, 577. (10) Jiang, S.; Zollweg, J.; Gubbins, K. J. Phys. Chem. 1994, 98, 5709. (11) Benard, P.; Chahine, R. Langmuir 1997, 13, 808. (12) Zhou, L.; Zhou, Y. Ind. Eng. Chem. Res. 1996, 35, 4166. (13) Murata, K.; Kaneko, K. Chem. Phys. Lett. 2000, 321, 342. (14) Zhou, L.; Zhou, Y.; Li, M.; Chen, P.; Wang, Y. Langmuir 2000, 16, 5955.

10.1021/la020786+ CCC: $25.00 © 2003 American Chemical Society Published on Web 02/21/2003

2216

Langmuir, Vol. 19, No. 6, 2003

Do et al.

Table 1. Summary of the Past Work on Supercritical Fluid Adsorption authors

theoretical approaches

Rangarajan et Aranovich and Donohue17 Murata and Kaneko13 Zhou et al.14 Zhou et al.18,19 Kaneko and Murata,20 Amankwah and Schwarz,21 Ozawa et al.4 Tan and Gubbins,22 Kaneko et al.,23,24 Cracknell et al.,25 Gusev et al.,26 Heuchel et al.,27 Suzuki et al.28 Fischer et al.,29 Kim et al.30 Lowry and Olmstead,31 Barrer and Robins,32 Zhou et al.33

simplified local density model (SLD) Ono-Kondo theory Langmuir fitting & VADS is constant Langmuir-Freundlich fitting & VADS is proportional to density transformation of isotherm quasi-saturated vapor pressure

the maximum behavior and explore the existence of minimum as observed by Vermesse and Levesque34 and address the negative surface excess.

jj, σii, σjj using the Lorentz-Berthelot mixing rule

al.15,16

DFT, GCMC potential calculation equation of state

2. The Solid-Fluid Potential Energy Before we present the MLST and GCMC simulation, we would like to briefly describe the potential energy equation used in these methods. Let us ignore the contribution of the electrostatic forces and assume dominant force of adsorption to be purely dispersive. The potential energy of interaction between an adsorbate molecule i with the surface atom j can be calculated by the 12-6 Lennard-Jones potential energy equation

[(

φij ) 4ij

) (

σij |ri - rj|

12

-

σij |ri - rj|

)] 6

(1)

where ri is the position where the center of the molecule i is located. If the interaction is assumed purely additive (which is, of course, questionable for highly dense fluids), the potential energy between an adsorbate molecule and the surface can be obtained by summing all pairwise potential energies (eq 1) between the molecule and all the surface atoms including those lying underneath the outermost surface layer. The parameters ij and σij may be calculated from pure species molecular properties ii, (15) Rangarajan, B.; Lira, C.; Subramanian, R. AIChE J. 1995, 41, 838. (16) Subramanian, R.; Pyada, H.; Lira, C. Ind. Eng. Chem. Res. 1995, 34, 3830. (17) Aranovich, G.; Donohue, M. J. Colloid Interface Sci. 1996, 180, 537. (18) Zhou, L.; Zhou, Y. Chem. Eng. Sci. 1998, 53, 2531. (19) Zhou, L.; Chen, P.; Li, M.; Sun, Y.; Zhou, Y. In Adsorption Science and Technology, Do, D. D., Ed.; 2000; pp 717-721. (20) Kaneko, K.; Murata, K. Adsorption 1997, 3, 197. (21) Amankwah, K.; Schwarz, J. A. Carbon 1995, 33, 1313. (22) Tan, Z.; Gubbins, K. J. Phys. Chem. 1990, 94, 6061. (23) Kaneko, K.; Shimizu, K.; Suzuki, T. J. Chem. Phys. 1992, 97, 8705. (24) Kaneko, K.; Cracknell, R.; Nicholson, D. Langmuir 1994, 10, 4606. (25) Cracknell, R.; Gordon, P.; Gubbins, K. J. Phys. Chem. 1993, 97, 494. (26) Gusev, V.; O’Brien, J.; Seaton, N. Langmuir 1997, 13, 2815. (27) Heuchel, M.; Davies, G.; Buss, E.; Seaton, N. Langmuir 1999, 15, 8695. (28) Suzuki, T.; Kobori, R.; Kaneko, K. Carbon 2000, 38, 623. (29) Fischer, J.; Bohn, M.; Korner, B.; Findenegg, G. Ger. Chem. Eng. 1983, 6, 84. (30) Kim, M.;. Holste, J.; Hall, K.; Slattery, J. J. Colloid Interface Science 1993, 158, 488. (31) Lowry, H.; Olmstead, P. J. Phys. Chem. 1927, 31, 1601. (32) Barrer, R. M.; Robins, A. Trans. Faraday Soc. 1951, 47, 773. (33) Zhou, C.; Hall, F.; Gasem, K.; Robinson, R. Ind. Eng. Chem. Res. 1994, 33, 1280. (34) Vermesse, J.; Lesveque, D. J. Chem. Phys. 1994, 101, 9063. (35) Malbrunot, P.; Vidal, D.; Vermesse, J.; Chahine, R.; Bose, T. Langmuir 1997, 13, 539. (36) Ustinov, E.; Do, D. D. J. Colloid Interface Sci. 2002, 253, 247. (37) Do, D. D.; Ustinov, E.; Do, H. D. Fluid Phase Equilib., in press.

σij ) (σii + σjj)/2 and ij ) (iijj)1/2

(2)

Assuming a pairwise additive potential, the potential energy of a molecule center and a layer of graphite having a two-dimensional density Fss (center/m2) is obtained by summing eq 1 for all combinations of molecule and surface atoms

∑k

φi,L ) 4is

[( ) ( ) ] σis

12

σis

-

ri,k

6

(3)

ri,k

where ri,k is the distance between the molecule i and the surface atom at the position k. Since the two-dimensional density of a graphite surface is dense and if the distance of the adsorbate molecule from the surface is greater than the distance between two adjacent carbon atoms, the process of pairwise summation of 12-6 LJ potential energies (eq 1) can be reasonably approximated by an integration, of which the result is the 10-4 potential.38

φi,L )

( ) [ ( ) ( )] 5 2 σis i,L 3 5 z

10

-

σis z

4

(4)

where i,L is the well depth of the interaction and is given by

i,L )

6 π F sσ 2 5 is s is

(5)

For a single graphite surface, the two-dimensional carbon density is39

Fss ) 3.813 × 1019 m-2

(6)

Carbon black surface is assumed to be comprised of many parallel graphite planes, separated by a distance ∆. For this set of parallel graphite layers, the potential energy of interaction between a molecule and this set of graphite layers is simply obtained by summing the potential of eq 4 between the molecule and all layers. The result is

φi,w )

( ) [ ( ) ( )] 5

3

∞

2

σis

∑ k)0 5 z + k∆

i,L

10

-

σis

z + k∆

4

(7)

One approximation to eq 7 is that proposed by Steele.40 The sum from the second and higher layers is replaced by (38) Everett, D.; Powl, J. J. Chem. Soc., Faraday Trans. 1976, 72, 619. (39) Donnet, J. B.; Voet, A. Carbon Black: Physics and Chemistry and Elastomer Reinforcement; Marcel Dekker: New York, 1976. (40) Steele, W. Surf. Sci. 1973, 36, 317.

Adsorption of Supercritical Fluids

Langmuir, Vol. 19, No. 6, 2003 2217

Table 2. Values of f and g (defined in Eq 9) for Various Values of ∆/σis f -g

0.25

0.30

0.35

0.50

0.65

0.75

0.85

1

1.50

2

0.9336 0.7809

0.9467 0.6619

0.9566 0.5827

0.9750 0.4546

0.9846 0.3958

0.9884 0.3727

0.9912 0.3565

0.9940 0.3402

0.9980 0.3158

0.9991 0.3076

an integration with the lower limit of the integral being z + 0.61∆, and the repulsive terms are neglected for those layers. The result is the well-known Steele 10-4-3 potential for a molecule at a distance z from the center plane of the outermost graphite layer of the carbon black surface

[(

φi,w(z) ) φw

)

1 σis 5 z

10

-

( )

1 σis 2 z

4

-

σis4

]

6∆(z + 0.61∆)3

(8a)

where φw ) 4πFsσis2∆is. Equation 8a is the potential energy of interaction between one molecule with the surface. Thus the energy of interaction between 1 mol and the surface follows the same form with φw defined as

φw ) 4πFsσis2∆isNA

(8b)

where NA is the Avogadro number. The following values are usually adopted in the literature, Fs ) 114 × 1027 m-3, ∆ ) 3.354 × 10-10 m, σss ) 3.4 × 10-10 m, (ss/kB) ) 28 K, and they will be used in this paper. Although the graphite surface comprises of carbon atoms arranged in a hexagonal structure, the structureless 10-4-3 model of Steele (eq 8) works well for many substances including hydrocarbons. This is possibly due to the fact that hydrocarbon molecules spread over the graphite surface covering more than one energetic site on the surface and therefore the periodic surface energy fluctuation is smeared as a result. The potential energy of eq 8 has a single minimum at a position, which is slightly less than one collision diameter (σis) from the plane passing through the centers of carbon atoms residing on the outermost layer of the surface. The minimum position and the potential energy at this point are given by

φ(zmin) zmin ) f(∆/σis) and ) g(∆/σis) σis φw

(9)

The values for f and g for various values of ∆/σis are given in Table 2.41 The solid-fluid molecular parameters, is and σis, are usually obtained by matching the theory with the experimental data of an open surface (usually the Henry constant is used). If the data are not available, they can be estimated from the Lorentz-Berthelot mixing rule of eq 2. Very often we have to introduce a binary interaction parameter (k12) in the calculation of the solid-fluid interaction energy as

is ) (1 - k12)(iiss)1/2

by first imagining that this layer is surrounded by neighboring layers of the same surface density. These layers are equally spaced at a distance of δ at equilibrium. Thus this set of layers behaves as a three-dimensional bulk fluid having a volumetric density of F (mol/m3), which is related to the surface density and layer spacing as

(10)

3. The Molecular Layer Structure Theory (MLST) The MLST is simple in its description.36 It treats the fluids whether it be homogeneous or inhomogeneous, as parallel layers of molecules. The basic crux of this method is the calculation of the intrinsic Helmholtz free energy of an isolated single layer (i.e., there is no presence of surrounding layers). The two-dimensional surface density of this isolated layer is denoted as Fs (mol/m2). We will adopt a quasi-equilibrium approach to obtain this intrinsic Helmholtz free energy of this isolated layer. It is calculated (41) Do, D. D. Adsorption Analysis: Equilibria and Kinetics; Imperial College Press: London, 1998.

F ) Fs/δ

(11)

The equilibrium distance δ is obtained from the mechanical equilibrium condition, details of which will be shown later. The molar Helmholtz free energy of this set of layers (i.e., bulk fluid of density F) can be readily obtained from any suitable equation of state (EOS) of homogeneous bulk fluids, for example, the accurate Bender EOS or the popular cubic Peng-Robinson EOS. Once this is done, the intrinsic Helmholtz free energy of an isolated layer (of surface density Fs) can be obtained by subtracting the interaction energy of that layer with all surrounding layers from the bulk molar Helmholtz free energy. 3.1. Lattice Distance for Bulk Fluid. First let us briefly present the mechanical equilibrium condition, from which we can calculate the equilibrium distance between adjacent layers in the homogeneous bulk fluid. Starting with the 12-6 Lennard-Jones equation for two interacting centers of eq 1, we can calculate the interaction between one molecule with a layer of the same species at the distance z apart.38 Written in terms of potential energy of interaction between 1 mol, instead of one molecule, and a layer, we have

φf,L(z) )

( ) [ ( ) ( )] 5 2 σff f,L 3 5 z

10

-

σff z

4

(12)

where f,L is defined as a function of surface molar density Fs (mol/m2) and molecular properties (collision diameter and energy of interaction) as shown below

f,L )

6 πNAR(ff/kB)Fsσff2 5

(13)

Here, NA is the Avogadro number (6.023 × 1023 molecules/ mol), R is the gas constant (8.314 J/(mol/K)), and kB is the Boltzmann constant (1.3804 × 10-23 J/(molecule/K)). In a homogeneous fluid having layers of molecules separated by a distance δ, the energy of interaction of 1 mol with the jth layer located at a distance of jδ from that 1 mol is

φf,j(jδ) )

( ) [ ( ) ( )] 5 2 σff f,L 3 5 jδ

10

-

σff jδ

4

(14)

Thus, the total energy of interaction between 1 mol with all layers is simply obtained by summing eq 14 with respect to j from j ) 1 to j ) ∞. We have the following result

φf,L‚‚‚L )

( ) [ ( ) ∑ ( )∑ ] 5

3

f,L

2 σff

5 δ

10 ∞

j)1

j-10 -

σff δ

4 ∞

j-4

(15)

j)1

Here the subscript L‚‚‚L denotes for all layers. The two infinite sums in the above equation can be expressed in

2218


Do et al.

energy as follows

closed form as described below ∞

j ∑ j)1 ∞

-10

j ∑ j)1

)

2 -4

)

µs ) ∂(FsAs)/∂Fs

10 1 (2π) B5

10!

Substituting eq 19b into the above equation, we obtain the following expression for µs

4 1 (2π) B2

2

µs(Fs) ) µ(F) + 2CFs

4!

where B2 ) 1/30 and B5 ) 5/66 are Bernoulli numbers. The equilibrium distance between two adjacent layers in the bulk fluid is the distance at which the total energy of interaction is minimum; i.e., it is the distance at which the derivative of the total energy of interaction with respect to δ is zero. Solving this derivative equating to zero, we obtain the following equilibrium distance in the bulk fluid ∞

δ ) σff(

∞

j /∑j ∑ j)1 j)1 -10

-4 1/6

)

( )

) 2π

4!B5

1/6

10!B2

φf,L‚‚‚Leq ) -CFs 3(308) 100

N

σff ≈ 0.987065σff (17)

Fjs[Ajs + φj,L‚‚‚L + φj,w - µ] ∑ j)1

(18)

φj,L‚‚‚L )

1

∞

∑ FksUj,k 2k)1

(23)

k*j

πNAR(ff/kB)σff2

We see that the total potential energy of interaction is proportional to the surface molar density. The unit of the parameter C is (J m2/mol2). 3.3. Intrinsic Helmholtz Free Energy. The intrinsic Helmholtz free energy of an isolated layer of surface density Fs can now be calculated by subtracting the interaction energy of eq 18 from the bulk molar Helmholtz free energy. Let As be the molar intrinsic Helmholtz free energy of an isolated layer of surface density Fs and A be the molar Helmholtz free energy of the bulk fluid of density F ()Fs/δ). We have

(19a)

At a given temperature and a bulk fluid density F, the Helmholtz free energy of the bulk fluid can be accurately calculated from an equation of state (Appendix) or from the p-F-T data. Therefore, the intralayer intrinsic Helmholtz free energy of an isolated layer can be readily obtained from eq 19a since φf,L‚‚‚Leq is given in eq 18, that is

As ) A + CFs

(22)

where N is the number of layer with the first layer being the one closest to the surface and the value for N must be large enough to cover layers in the vapor phase. For supercritical fluids, a choice of N ) 15 is usually adequate. In eq 22, µ is the chemical potential of the bulk gaseous phase. The interaction energy φj,L‚‚‚L is the interaction energy between 1 mol of the layer j with all surrounding layers, and the interaction energy φj,w is the interaction energy between 1 mol of layer j with the surface. The former is defined as

2/3

As ) A - φf,L‚‚‚Leq

(21)

3.5. Adsorption on Nonporous Carbon Black. We have obtained the necessary ingredients in the last section. Here we shall develop the model for adsorption of fluid onto a nonporous surface. The grand potential energy per unit area of surface for a grand canonical (specify T, V, and chemical potential µ) is

Ω)

This equilibrium distance does not depend on the density of the fluid, nor does it depend on the interaction energy. 3.2. Potential Energy of Interaction with Surrounding Layers. Once we know the equilibrium distance δ, the interaction energy between 1 mol and all surrounding layers is simply given in eq 15. By substituting the equilibrium distance of eq 17 into eq 15, we obtain this interaction energy as

C)

(20)

(16)

(19b)

This has completed all the necessary ingredients for the development of the MLST model. 3.4. Intralayer Chemical Potential. Another quantity which is useful for later analysis is the intralayer chemical potential of an isolated molecular layer. It is defined in terms of the intralayer molar Helmholtz free

The factor 1/2 in the above equation is introduced to avoid counting the mutual interaction twice in eq 22. Here Uj,k is the interaction between 1 mol in the layer j with 1 mol of the kth layer, which is assumed to be given by the following 10-4 potential

[ ( ) ( )]

5 2 σff Uj,k ) D 3 5 zj - zk

10

-

σff zj - zk

4

(24)

The variable zj is the coordinate of the layer j, and the parameter D contains only the molecular parameters as shown below

() [ ]

ff 2 6 σ D ) πNAR 5 kB ff

J m2 mol2

(25)

The external potential energy of interaction, φj,w, between 1 mol in the layer j with the surface is assumed to follow the Steele 10-4-3 potential eq 8. The positions of the layers as well as the densities of each layer can be obtained by minimizing the grand potential with respect to the coordinate zj and to the surface molar density Fjs. The partial derivatives of the grand potential with respect to zj equating to zero are

∂Ω ∂zj

(

) Fjs

∂Uj,k

∑ Fks ∂z k)1 k*j

j

+

)

dφj,w dzj

)0

(26)

for j ) 1, 2, ..., N, where ∂Uj,k/∂zj and dφj,w/dzj are given


Langmuir, Vol. 19, No. 6, 2003 2219

4. Grand Canonical Monte Carlo Simulation

by

( )( )[( ) ( ) ] )[ ( ) ( ) σff zk - z j

∂Uj,k 20 D ) ∂zj 3 σff

(

dφj,w 1 φw σsf ) -4 dzj 2 σsf zj

11

+4

11

σsf zj

5

σff zk - zj

5

σsf5

+

(27)

]

∆(0.61∆ + zj)4

(28) The partial derivative of the grand potential with respect to density Fjs equating to zero will provide the second set of equations. We have

∂Ω ) µjs + 2φj,L‚‚‚L + φj,w - µ ) 0 s ∂Fj

(29)

for j ) 1, 2, ..., N. Solving eqs 26 and 29 will yield solution for the positions of all layers and their respective surface densities. Solution is effectively done with the NewtonRalphson method with a special care that the new iterated solutions should not differ from the previous iterated solutions by a few percent. About the starting values for the positions and the surface densities, we use the starting values as follows. For the positions, we choose the position of the first layer (closest to the graphite surface) to be the minimum of the solid-fluid potential energy, and the positions of the second and higher layers are separated by a distance δ, which is the equilibrium distance between layers in the bulk fluid. If the density profile is available for a previous value of bulk density, we use the positions of that profile as the starting values for the next value of the bulk density. For the surface densities, we either choose them to be the same as the bulk density or use the equilibrium density profile obtained from a previous value of the bulk density. Knowing the surface densities of all layers, the molar surface excess is calculated as the sum of surface densities of N layers minus the quantity if the volume of those N layers were to be occupied with molecules of bulk density. N

Γ)

Fjs ∑ j)1

(

zN +

δ

-

2

)

σss 2

Fbulk

(30)

(31a)

where µ is the chemical potential of the gas phase. The density and the Helmholtz free energy at the layer j in the presence of other layers and the surface are

Fj )

Fjs (zj+1 - zj-1)/2

Aj ) Ajs + φj,L‚‚‚L + φj,w

P ) min[1, exp(- ∆U/kT)]

(32)

where ∆U is the difference between the new configuration energy and the old configuration energy. The displacement step is chosen in such a manner that the acceptance rate is between 25 and 50%. 2. Insertion of Particle. The second move is the particle insertion move. A particle is generated at a random position. Its acceptance has the following probability

[

V P ) min 1, 3 exp{[µ - U(N + 1) + Λ (N + 1)

]

U(N)]/kT} (33) where V is the volume of the simulation box and Λ is the thermal de Broglie wavelength. 3. Removal of Particle. The third move in the GCMC simulation is the removal of a particle. A particle is chosen in random and removed. The probability of acceptance of this removal is

[

]

Λ3N exp{- [µ + U(N - 1) - U(N)]/kT} V (34)

P ) min 1,

The absolute amount adsorbed does not make any sense in the case of carbon black because of the fundamental question that where should one draw the boundary of the adsorbed layer. The pressure exerted by the molecules in the layer j can be computed from

Pj ) Fj(µ - Aj)

In the molecular simulation of adsorption in confined space such as pores of adsorbent, the most widely used ensemble is the grand canonical Monte Carlo simulation. In this ensemble, the chemical potential of the fluid µ that the candidate pore is immersed in, the size of the pore, and the temperature are specified. The GCMC simulation is carried out by starting at a very low value of chemical potential. Once this is done, the chemical potential is increased incrementally and the particle configuration of the last run is used as the initial guess for the new chemical potential. In the simulation, the simulation box is constraint by the width of the pore, H, and the lengths Lx and Ly in the x and y directions. The box lengths Lx and Ly are chosen large enough, and usually chosen as twice the cutoff distance. In our work we chose rc ) 5σff, and the usual periodic boundary conditions are applied in the x and y directions. The procedure of the GCMC involves three basic moves.42 1. Displacement of Particle. A particle is selected at random and is given a new conformation. The move is accepted with a probability

In our simulation, we chose the cut off distance of five times the collision diameter. The box lengths in the x and y directions are twice the cut off distance. The number of moves per particle is 30000. The number of moves (usually 3 × 106) is allowed for equilibrium to establish before statistics measurement is carried out to obtain the pore density and its average as a function of chemical potential. The molecular parameters used in the GCMC simulations for the adsorbates tested in this paper are listed in Table 3. 5. Discussion

(31b) (31c)

That concludes our formulation for the adsorption on nonporous carbon surface.

What we will present below are the results of the MLST model and the GCMC simulations. We test them with a (42) Frenkel, D.; Smit, B. Understanding molecular simulation; Academic Press: New York, 2002. (43) Cracknell, R. Mol. Phys. 1993, 80, 885. (44) Reid, R.; Prausnitz, J.; Poling, B. The properties of gases and liquids; McGraw-Hill: New York, 1988.

2220


Do et al.

Table 3. Molecular Parameters of Argon, Krypton, Methane, Ethylene, and SF6 adsorbate

σ (m)

/kB (K)

ref

argon krypton methane ethylene sulfur hexafluoride

3.405 × 10-10 3.685 × 10-10 3.810 × 10-10 4.218 × 10-10 5.128 × 10-10

119.80 164.41 148.10 201.80 222.10

34 34 43 22 44

Figure 2. Prediction of adsorption isotherms of supercritical argon at 253, 273, 298, and 323 K.

Figure 1. Fitting of the 298 K data of argon with the binary interaction as the fitting parameter.

large number of adsorbates, argon, krypton, methane, ethylene, and sulfur hexafluoride. The temperatures and the sources of the experimental data are listed in Table 4. The graphitized thermal carbon black used by Findenegg and co-workers has a specific surface area of 81 m2/g. We first discuss the simplest adsorbate, argon. The molecular parameters used in the MLST are obtained from the matching of the surface tension as predicted from the model with the experimental data,37 while those for GCMC simulations are taken from refs 34 and 45. Argon. We apply our MLST model to predict the adsorption of argon on nonporous carbon black surface at 253, 273, 298, and 323 K. The only parameter required for the prediction of the argon adsorption is the binary interaction parameter, k12, of eq 10 in the calculation of the solid-fluid interaction energy sf. To obtain the optimal value for the binary interaction parameter, we use the adsorption isotherm of 298 K. Figure 1 shows the predictions of the model for 298 K for three values of the binary interaction parameter. The need for this binary interaction parameter is because of the usual overprediction of the Berthelot mixing rule. We see that the model provides an excellent prediction when the binary interaction parameter is 0.2, suggesting that the Lorentz-Berthelot overestimates the solid-fluid interaction energy. The number of layer used in the computation is 15, which is more than sufficient in dealing with supercritical adsorption. Computations with N greater than 15 do not change the result. Using this optimal binary interaction parameter of 0.2, we apply the model to predict the adsorption isotherms of supercritical argon at temperatures 253, 273, and 323 K. Figure 2 shows the excellent prediction by the MLST model. We note the maximum in the 253 K isotherm, and the model also correctly describes not only this maximum at the pressure which the maximum occurs but also the surface excess at this maximum. The surface density profiles for the case of argon at 253 K, the percentage change in the layer spacing, and the (45) Saam, W. F.; Ebner, C. Phys. Rev. A 1978, 17, 1768.

Figure 3. Surface density, pressure, and change in layer spacing profiles for argon at 253 K for three values of bulk pressure: 2.24 (bottom curve), 28.94, and 144.3 bar (top curve).

pressure are shown in Figure 3 for three values of pressure, 2.24, 28.94, and 144.3 bar. Observing the surface density profile, we see that the density distribution gradually increases when the surface is approached. This is what is called the densification in the supercritical fluids. As a result of this densification, the pressure distribution exhibits an increase as the distance from the carbon surface is decreased. The change in layer spacing has a maximum of 16% smaller than the equilibrium distance


Langmuir, Vol. 19, No. 6, 2003 2221

Table 4. Source of Adsorption Data of Argon, Krypton, Methane, Ethylene, and SF6 adsorbate

temperature (K)

references

argon krypton methane ethylene SF6

253, 273, 298, 323 253, 273, 298, 323, and 348 253, 273, 298, and 323 283, 285, 288, 303, 313, and 323 333.29 and 341.96

Specovius and Findenegg5 Blumel et al.7 Specovius and Findenegg5 Findenegg1 Findenegg1

Figure 4. Argon surface excess at 298 K for a wide range of pressure.

in the bulk fluid. We also note from the density profiles that the adsorption of supercritical fluid tends to concentrate near the surface and the adsorbed volume can be seen to vary between about 1 to 4 layers. The same conclusion has been recently derived by Do and Do.46 We also note that the surface density is greatest at the first layer close to the surface, and this density is lower than the liquid density of argon. Taking the bulk liquid density of argon of 1600 kg/m3 and the interlayer spacing of liquid argon as 3.374 × 10-10 m, the surface liquid density is 1.35 × 10-5 mol/m2. The density of argon at the first layer is about 0.8 × 10-5 mol/m2, which is lower than the liquid density. We have also derived this conclusion in Do and Do.46 The power of the MLST model lies in its ability to predict the adsorption isotherms at other temperatures or higher range of pressure. Figure 4 shows the argon adsorption isotherm at 298 K at a much higher range of pressure (50 times greater than the critical pressure), well beyond the range that was covered experimentally by Findenegg and co-workers. Here we see that the model predicts a clear maximum in the isotherm, and this maximum occurs at a pressure (Pmax) just greater than the critical pressure. The bulk density at this maximum pressure is just less than the critical density. This maximum phenomenon occurs because the change in bulk density is very large for a very small change in pressure close to Pmax. What this means is that when the pressure in the bulk phase is increased by a small amount near Pmax, the bulk density is substantially increased, resulting in an increase in the adsorption isotherm. When pressure has passed far away from this Pmax, any large change in pressure results in a small change in density. But this rate of change in the bulk density is greater than the rate of change in the adsorbed density, resulting in a decrease in the surface mass excess. What we see from Figure 4 is that once the surface excess passes the maximum, it decays gradually toward the pressure axis. Even when we increase the pressure much further, there is no evidence of negative surface excess. We will come back to this point later to explore this question of negative surface excess. (46) Do, D. D.; Do, H. D. Carbon, in press.

Next, we would like to compare the MLST model and the GCMC simulation results. To get the correct Henry constant behavior, the binary interaction parameter for the GCMC simulation is found to be 0.2, which is identical to the value obtained in our MLST model. Figure 5a shows the GCMC simulation results for argon adsorption at 253 K. We see that the GCMC simulation results agree well with the data at low pressures (P < 4 MPa), but for pressures above 4 MPa, this method consistently overpredicts the data. The same observation is made when we study the adsorption isotherm of 323 K (Figure 5b). The distinct advantage of the MLST is the very fast computation time. For the GCMC simulation, it takes the order of days to obtain isotherms of four different temperatures, while the MLST takes about few minutes to complete the computation. The obvious question that one would raise is that why the GCMC simulation does not perform as good as (or even better) than the MLST. A possible reason might lie in the fact that the GCMC may not describe well the bulk properties at very high pressures, at which one would expect the multiparticle interactions. Next we would like to investigate to see whether the conclusions made for argon are applicable to all other adsorbates. Krypton. Next we consider another noble gas, krypton. Data for krypton are available from Blumel et al.7 They have obtained isotherms at 253, 273, 298, 323, and 348 K, and the data are available up to 150 bar. The 253 K isotherm is used in the fitting to obtain the binary interaction parameter. It was found to be about 0.2, which is the same as that obtained for argon. The model is then used to predict isotherms at other temperatures. The results are shown in Figure 6. We see that the performance of the MLST model is credible, and even the distinct maximum observed for krypton at 253 K (P ≈ 90 bar) is well described by the model. Next we would like to investigate the behavior of krypton adsorption at extremely high pressure. We performed the calculation for two temperatures: 253 and 273 K. What we observed for argon earlier is repeated here for krypton (Figure 7). We note the distinct maximum and the surface excess decreases gradually after the maximum has been reached, and the surface excess is always positive no matter how high the pressure is. Another feature we note for these two isotherms is that they cross each other after that have passed their respective maxima. The crossing behavior of the isotherm at high pressure is because the change in the bulk density with pressure is different at different temperatures. If we, instead, plot the surface excess in terms of the bulk density, there will be no crossover phenomenon. We next study the comparative performance between the MLST and the GCMC. The isotherm of 253 K is taken for illustrative example (Figure 8). Again, like the case of argon, the GCMC performs well at low pressures but fails to describe the data at higher pressures. Methane. The same analysis is now carried out with methane adsorption data of Specovius and Findenegg.5 Using the 298 K isotherm in the fitting between the MLST model and the data, the binary interaction parameter is

2222


Do et al.

Figure 5. Comparison between MLST and GCMC: (a) T ) 253 K and (b) T ) 323 K.

Figure 6. Plot of the surface excess of krypton versus pressure at 253, 273, 298, 323, and 348 K.

Figure 7. Plot of the surface excess of krypton versus pressure at 253 and 273 K over a wide range of pressure.

found to be 0.1. The agreement between the theory and the data is not as excellent as that observed for argon and krypton, but the agreement still can be regarded as very good (Figure 9a). The negative surface excess is not observed, as was the case for argon and krypton (Figure 9b), and finally the GCMC results, yet again, overpredict the data at very high pressures (Figure 10). Ethylene. We next study the adsorption of ethylene on graphitized thermal carbon black.1 The critical temperature of ethylene is 282.4 K. Our investigation using the MLST model with ethylene adsorption data of 303, 313, and 323 K reveals that what we have observed so far for

Figure 8. Comparison between MLST and GCMC results for krypton at 253 K.

argon, krypton, and methane is equally applied to this case of ethylene. The binary interaction parameter for ethylene was found to be 0.2. What is interesting here is the availability of adsorption data at temperatures just greater than the critical temperature. Findenegg1 reported adsorption data at three temperatures 283, 285, and 288 K, and these isotherms are shown in Figure 11 as symbols. Also presented in those figures are the computed results from the MLST model and the GCMC simulation. First, we note the sharp spike behavior of the adsorption isotherm, and this sharp spike is due to the interplay between the behavior of isotherm “like” subcritical fluids at pressures less than the critical pressure and the behavior of isotherm “like” supercritical fluids at pressures greater than the critical pressure. This sharp spike is well described by the MLST model, the position where the spike occurs and even interesting is the maximum value of this spike (note the different spike heights for the three temperatures 283, 285, and 288 K). The GCMC simulation, on the other hand, fails to describe this spike, especially at 283 K, the temperature very close to the critical temperature of 282.4 K. Thus we have illustrated that the potential of the MLST as a useful tool to study adsorption of supercritical fluids on graphitized thermal carbon black. This potential is finally illustrated with the experimental data of sulfur hexafluoride. This is shown in Figure 12 with data of Findenegg1 as symbols and the MLST results plotted as solid lines (the binary interaction parameter was found to be 0.4). Note the critical temperature of sulfur hexafluoride is 318.96 K.


Langmuir, Vol. 19, No. 6, 2003 2223

Figure 9. Adsorption isotherms of methane at 253, 273, 298, and 323 K: (a) experimental range of pressure (b) wider range of pressure. Table 5. Advantages and Disadvantages of the MLST and GCMC Approaches in the Analysis of Supercritical Adsorption MLST advantages

disadvantages

GCMC

describe well the data, including the sharp spike at temperatures close to critical temperature very short computation time applicable only to rectangular geometry not applicable for fluids where PVT of the bulk fluid is not available

Figure 10. Comparison between the MLST and GCMC for methane at 323 K.

Having performed analysis of argon, krypton, methane, ethylene, and sulfur hexafluoride using the MLST model and the GCMC simulations, we list in Table 5 the advantages and disadvantages of the two approaches.

excellent tools for full predictability power applicable to systems of any geometry and any forms of potential energy of interaction extremely long computation time

The Negative Surface Excess? Having the MLST model as a powerful tool to study the adsorption of supercritical fluids on nonporous carbon black, we would like to use it to study this adsorption system further by changing the system parameters and investigate the sensitivity. First we would like to investigate the effect of the density of the surface carbon atom. The following figure (Figure 13a) shows the surface mass excess of argon adsorption at 298 K when the surface atom density is reduced from its value of 114 × 1027 m-3 by factors of 2 and 5. The curve A is for Fs ) 114 × 1027 m-3, while curves B and C are for Fs ) 57 × 1027 and 22.8 × 1027 m-3, respectively. We see that when the surface carbon atom is less dense enough, the surface mass excess can become negative at some high concentration in the gas phase (curve C in Figure 13a). This is what we would expect physically because if we take one argon molecule in the bulk of high pressure and put it near the surface, this molecule will experience a lesser interaction with the neighboring molecules because of the lesser density of carbon atoms. As a result, the argon density near the surface will be lower than the density in the bulk, leading to a negative surface mass excess. Next we study the effect of the surface

Figure 11. Adsorption isotherms of ethylene on carbon black at temperatures close to the critical temperature: (a) T ) 283 K, (b) T ) 285 K, and (c) T ) 288 K.

2224


Do et al.

Glossary A As C D k1 r

Helmholtz free energy intrinsic layer Helmholtz free energy constant, defined in eq 18 parameter, defined in eq 25 binary interaction parameter position of molecule, or distance between two centers potential, defined in eq 24 distance between adsorbate molecule and the plane passing through centers of carbon atoms of the outermost layer equilibrium distance between adjacent layers in the bulk fluid distance between two adjacent graphite layers interaction energy surface excess, defined in eq 30 chemical potential intralayer chemical potential interaction potential energy volumetric density surface density collision diameter grand potential, defined in eq 22

U z δ

Figure 12. Plots of the surface excess of sulfur hexafluoride versus pressure for various temperatures.

interaction energy. For carbon surface, this interaction energy is ss/kB ) 28 K. We show in Figure 13b the effect of varying this parameter. We use values 56, 28, 14, and 2.8 K. The first value corresponds to a hypothetical surface that is twice stronger in energy than the carbon surface, while the third and fourth values correspond to surfaces that are weaker than the carbon surface by factors of 2 and 10, respectively. We see in this figure that the negative surface excess is also possible when the solid interaction energy is about 10 times less than that of carbon surface. Thus negative surface excess is only possible with surface that has either very low in surface density or very low interaction energy. For typical carbon surface of density Fs ) 114 × 1027 m-3, the negative surface excess is not possible.

∆ Γ µ µs φ F Fs σ Ω

Acknowledgment. This project is supported by the Australian Research Council Appendix: Thermodynamic Properties of Real Fluid To model the thermodynamic properties of a real fluid, the following two equations of state are used: (1) the Peng-Robinson (PR) and (2) Bender (B) equations of state. They have the following form:

6. Conclusions We have presented in this paper a simple and effective new method for describing adsorption isotherms of supercritical fluids onto graphitized thermal carbon black. This method is coined “molecular layer structure theory”. It describes very well the adsorption isotherm and its temperature dependence. We successfully tested our theory with argon, krypton, methane, ethylene, and sulfur hexafluoride.

p)

{

aF2 F RT 1 - bF 1 + 2bF - b2F2

(PR)

(B) FT[R + BF + CF2 + DF3 + EF4 + FF5 + (G + HF2)F2 exp(-a20F2)]

(A1)

In the PR equation of state, the parameters a and b are determined only from the critical properties and the acentric factor. The parameter “a” is also a function of

Figure 13. Effect of (a) surface carbon density and (b) solid interaction energy on the surface excess of argon on graphitized carbon black.


Langmuir, Vol. 19, No. 6, 2003 2225

temperature. In the Bender EOS, the temperaturedependent parameters B, C, ..., H and a20 contain a total of 20 temperature-independent parameters. These 20 parameters of the Bender EOS are found by matching the equation with the P-F-T data, second virial coefficients, etc., and as a result it is more accurate in the calculation of vapor pressure, liquid density, and thermodynamic quantities. Because of this, the Bender EOS is used in our investigation. The two fundamental thermodynamic equations are

dA ) -P dV - S dT

(A2)

dµ ) V dP - S dT

holds, and (the superscript IG denotes ideal gas)

µ(T,P) ) A(T,P) + PV

(A3)

At constant T, the molar Helmholtz and Gibbs free energies of eq A2 are given by

dA ) -P dV ) dµ ) V dP )

P dF F2

(A4)

dP F

where F is the volumetric molar density (F ) 1/V). In the limit of very low density, the ideal gas law of P ) FRT

(A5a)

AIG(T,P) ) µ0(T) + RT ln P - RT

(A5b)

To obtain the molar Helmholtz free energy for the real fluid at a given T, we substitute the Bender EOS (eq A1) in eq A4 and carry out the integration with the boundary condition of eq A5b to get

A ) µ0(T) + RT ln(RTF) - RT +

(

T BF + Combining the two equations of eq A2 and integrating the resulting equation, we obtain the relationship between the molar Helmholtz and Gibbs free energies

µIG(T,P) ) µ0(T) + RT ln P

) ( )[ ) ]

CF2 DF3 EF4 FF5 H T + + + + G+ 2 3 4 5 2a20 a20 H G + HF2 + exp(-a20F2) (A6) a20

(

Knowing the molar Helmholtz free energy of the real fluid as given in eq A6, the molar Gibbs free energy (chemical potential) can be obtained from eq A3, given below

µ ) µ0(T) + RT ln(RTF) +

( ( )[ T 2a20

)

3CF2 4DF3 5EF4 6FF5 + + + + 2 3 4 5 H H G+ - G + HF2 + exp(-a20F2) + a20 a20

T 2BF +

LA020786+

(

)

2

2

2

]

T(G + HF )F exp(-a20F ) (A7)

Adsorption of Supercritical Fluids on Graphitised Thermal Carbon

Recommend Documents