Adsorption of Ethylene on Graphitized Thermal Carbon Black and in

Langmuir 2004, 20, 7103-7116

7103

Adsorption of Ethylene on Graphitized Thermal Carbon Black and in Slit Pores: A Computer Simulation Study D. D. Do* and H. D. Do Department of Chemical Engineering, University of Queensland, St. Lucia, Queensland 4072, Australia Received February 18, 2004. In Final Form: May 18, 2004 In this paper, we studied vapor-liquid equilibria (VLE) and adsorption of ethylene on graphitized thermal carbon black and in slit pores whose walls are composed of graphene layers. Simple models of a one-center Lennard-Jones (LJ) potential and a two-center united atom (UA)-LJ potential are investigated to study the impact of the choice of potential models in the description of VLE and adsorption behavior. Here, we used a Monte Carlo simulation method with grand canonical Monte Carlo (GCMC) and Gibbs ensemble Monte Carlo ensembles. The one-center potential model cannot describe adequately the VLE over the practical range of temperature from the triple point to the critical point. On the other hand, the two-center potential model (Wick et al. J. Phys. Chem. B 2000, 104, 8008-8016) performs well in the description of VLE (saturated vapor and liquid densities and vapor pressure) over the wide range of temperature. This UA-LJ model is then used in the study of adsorption of ethylene on graphitized thermal carbon black and in slit pores. Agreement between the GCMC simulation results and the experimental data on graphitized thermal carbon black for moderate temperatures is excellent, demonstrating that the potential of the GCMC method and the proper choice of potential model are essential to investigate adsorption. For slit pores of various sizes, we have found that the behavior of ethylene exhibits a number of features that are not manifested in the study of spherical LJ particles. In particular, the singlet density distribution versus distance across the pore and the angle between the molecular axis and the z direction provide rich information about the way molecules arrange themselves when the pore width is varied. Such an arrangement has been found to be very sensitive to the pore width.

1. Introduction Vapor-liquid equilibria (VLE) and adsorption behavior on an open surface or in pores have been studied very regularly in the literature by various approaches, ranging from empirical approaches to those that have a basis in classical thermodynamics. Recently, more advanced tools of Monte Carlo simulation,2 molecular dynamics simulation, and density functional theory3,4 have been applied to solve this problem. The availability of the high-speed computer has made the application of these tools more accessible to scientists and engineers. The wider application of these tools is mostly due to the power in predicting VLE and adsorption behavior from a single information of a pairwise potential energy model between two fluid particles or two sites. Usually there are only two parameters in the description of this potential energy model, namely, the collision diameter σ and the well depth of the interaction energy . 1.1. Two-Center Potential Model. In this paper, we will investigate the VLE and adsorption behavior of linear molecules on graphitized thermal carbon black and in slit pores using only the information of pairwise molecular potential. We take ethylene as a model species in our investigation, and the two-center united atom (UA)TraPPE model adopted in this paper is the one developed by Wick et al.1 In this model, the CH2 is grouped together to form one interaction site. Thus, the potential model for one ethylene molecule contains two identical interaction sites. According to their model, the distance between these * Author to whom all correspondence should be addressed. (1) Wick, C. D.; Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8008-8016. (2) Steele, W. A. Appl. Surf. Sci. 2002, 196, 3-12. (3) Olivier, J. In Characterization of Porous Solids III; Rouquerol, J., Ed.; Elsevier: Amsterdam, 1994; p 81. (4) Olivier, J. J. Porous Mater. 1995, 2, 9-17.

two sites is the same as the carbon-carbon double bond length. If we assume that the site-site interaction between two molecules is described by the classical 12-6 LennardJones potential equation, then the interaction potential energy of a site a on a molecule i and a site b on a molecule j is given by

[( ) ( ) ]

(a,b) φi,j ) 4(a,b)

σ(a,b) (a,b) ri,j

12

-

σ(a,b) (a,b) ri,j

6

(1a)

where the subscript denotes for molecule while the superscript is for the interaction site. Each interaction site is characterized by two parameters: collision diameter (σ(a,a)) and well depth of the interaction energy ((a,a)). Because the two groups of ethylene are identical, the ethylene molecule is characterized by one value of the collision diameter, one value of the well depth of the interaction energy, and the distance between the two interaction sites. Knowing the site-site interaction energy, the molecule-molecule interaction energy is then calculated from the following equation: M

φi,j )

M

∑ ∑ φi,j(a,b) a)1 b)1

(1b)

where M is the number of sites on a molecule (M ) 2). The values of these parameters for the UA-TraPPE model1 are σ ) 3.675 Å, /k ) 85 K and ι ) 1.33 Å, where ι is the carbon-carbon double bond length. 1.2. One-Center Potential Model. Before the twocenter potential models became popular, ethylene was modeled with a one-center potential model; that is, it is treated as a “pseudo” spherical particle. Such an approach might be adequate for high-temperature study and when the system is fairly dilute. Here, we take a potential model

10.1021/la0495682 CCC: $27.50 © 2004 American Chemical Society Published on Web 07/17/2004

7104

Langmuir, Vol. 20, No. 17, 2004

Do and Do

with molecular parameters presented in Tan and Gubbins:5 σ ) 4.218 Å, /k ) 201.8 K. This model will be tested for its predictability of VLE against the two-center potential model. 2. Theory The simulation tool used in this paper for the VLE analysis is the Gibbs ensemble Monte Carlo (GEMC) and that for the adsorption study is the grand canonical Monte Carlo (GCMC) simulation. These tools are described in detail in Frenkel and Smit,6 and they are briefly described here. 2.1. GEMC. In this method, the phase coexistence between vapor and liquid can be readily computed without the consideration of an interface between these two phases. Another advantage of the method is that there is no need for the explicit determination of the free energy or chemical potential of the two phases. The Monte Carlo steps are designed so that at equilibrium there will be equality of pressure, temperature, and chemical potential between the two phases. Because there is no need to consider the interface joining the two phases, the two coexisting phases can be simulated in two separate simulation boxes. The liquid box, assigned as box I, and the vapor box, assigned as box II, are specified with volume, number of particles, and temperature. Usually the vapor box is chosen with much greater volume than the liquid box so that we will have a number of particles in that box large enough for acceptable statistical averaging of the Markov chain of configurations. In this GEMC, there are three basic moves. The first move is the displacement of the particle. This can be done by choosing a box at random, and then a particle in that box is chosen randomly. This particle is displaced to a new position, and the acceptance or rejection of this move will follow the usual rule of acceptance commonly used in Monte Carlo simulation:

[

(

)]

∆U kT

pacc ) min 1, exp -

(2)

where ∆U is the difference between the configurational energy after the displacement and that before the displacement. If the energy of the new configuration is lower than the one before the randomly selected particle is displaced, the new configuration is accepted and is counted as one state of the Markov chain. However, if the energy of the new configuration is greater than the old one, this new configuration is either accepted or rejected, depending on if the probability of acceptance of moving from the old configuration to the new configuration is either greater or lower than a random number generated uniformly between 0 and 1. If rejected, the old configuration is counted again in the Markov chain; otherwise, the new configuration is accepted and is counted. The second type of move is the interchange of a particle between the two boxes. In this move, a box is selected at random, say box II, and then a particle is selected at random in this box and moved from this box to a random position in box I. This move has the probability of acceptance as given in eq 2 with ∆U of this particle interchange being given by

∆Uinterchange ) ∆UI + ∆UII + kT ln

[

]

(NI + 1)VII NIIVI

(3)

(5) Tan, Z.; Gubbins, K. J. Phys. Chem. 1990, 94, 6061-6069. (6) Frenkel, D.; Smit, B. Understanding Molecular Simulations; Academic Press: New York, 2002.

where ∆Ui (i ) I, II) is the energy change that occurs in simulation box i and Ni and Vi are the number of particles in (prior to interchange) and the volume of the box i, respectively. The third move in the GEMC is the volume exchange. Let ∆V be the volume change. A box is chosen at random, say box I, its volume is decreased by ∆V, and, consequently, the volume of box II is increased by ∆V so that the total volume of the two boxes remains constant. The probability of acceptance for this move is given as in eq 2 with ∆U being given by

(

) (

VI - ∆V VI VII + ∆V NIIkT ln (4) VII

∆Uvol ) ∆UI + ∆UII - NIkT ln

)

We have described the three basic moves for the GEMC simulation of spherical molecules. For linear molecules, such as ethylene, studied in this paper, we have an additional move to displacement, which is the change in the molecule’s orientation. In this move, a molecule is selected in random and rotated to a new random orientation. The probability of acceptance is the same as for displacement (i.e., eq 2). 2.2. GCMC Simulation. In the GCMC simulation, we specify temperature, volume (pore volume), and the chemical potential. This simulation is ideal to study adsorption where a solid adsorbent (or a single pore) is exposed to a bulk fluid of constant pressure or chemical potential. The common feature shared by all Monte Carlo simulation methods, including the GEMC described in the previous section, is that a Markov chain of molecular configurations is produced. Any properties of interest can be obtained by averaging over this Markov chain. This method was first developed by Norman and Filinov7 and Adams.8 It has been extensively applied to simulate adsorption with great success. In GCMC, there are three different moves used to generate the Markov chain, which is composed of a series of molecular configurations. The first move is the displacement, which is the same as that described earlier for GEMC. For linear molecules, we have an additional move, which is the change of the particle’s orientation. This procedure again is the same as that in GEMC. The two remaining moves in the GCMC are the creation and destruction of a particle. In the creation move, a particle is created at a random position within the simulation box, while in the destruction move a particle is selected at random and removed from the box. These two moves are either accepted or rejected when appropriate probabilities of acceptance are either greater or lower than a random number generated uniformly between 0 and 1. Once the Markov chain (long enough) has been generated, various quantities of interest can be obtained, such as the average number of particles or the configurational energy, and this is simply achieved by averaging their values over the Markov chain. 2.2.1. Solid-Fluid Interaction Energy. In the GCMC for adsorption study, the external interaction exerted by the presence of solid surfaces can be evaluated once a model for the surface or pore is chosen. Pores are generally complex because of the variation in size and shape. Further complexity comes from the corrugation of the surface and the surface chemistry, such as the defect and the functional groups attached to the surface (usually at the edge of the (7) Norman, G. E.; Filinov, V. S. High Temp. (USSR) 1969, 7, 216222. (8) Adams, D. J. Mol. Phys. 1975, 29, 307-311.

Ethylene Adsorption: A Computer Simulation Study

Langmuir, Vol. 20, No. 17, 2004 7105

surface). Attempts have been made in the literature to model the surface or pore as detailed as possible. Such a model increases the effort and the computation time by one or two orders of magnitude. This is only justified when one wishes to study in fine detail the structure of adsorbed molecules close to the surface, and this is particularly important when the temperature is very low. In this paper, we shall take a simple model of a surface and a slit pore, where the surface is assumed to be flat and homogeneous (structureless). It is infinite in the two lateral directions. Because of the infinite nature of the pore, we will work with a finite simulation cell with periodic boundary conditions in the x and y directions. The interaction potential energy between a particle and the homogeneous solid carbon substrate is usually calculated by the wellknown 10-4-3 Steele potential.9 The interaction between a site a on a molecule j and the surface takes the following form: (a) φi,s ) 4πFC(a,s)[σ(a,s)]2∆ ×

{( ) ( ) 1 σ(a,s) 5 z(a) i

10

-

1 σ(a,s) 2 z(a) i

[σ(a,s)]4

4

-

3 6∆(0.61∆ + z(a) i )

}

(5a)

where FC is the volumetric carbon atom density (114 nm-3), ∆ is the spacing between two adjacent graphene layers, and z(a) i is the shortest distance between the site a on the molecule i and the plane passing through all carbon atoms on the outermost layer of the carbon wall. The solid-site molecular parameters, the collision diameter, and the interaction energy are calculated from the LorentzBerthelot mixing rule.

σ(a,s) ) (σ(a,a) + σ(s,s))/2; (a,s) ) (1 - ks)x(a,a)(s,s) (5b) The well depth of the solid-fluid interaction energy (eq 5b) is adjusted with the solid-fluid binary interaction parameter, ks, such that the experimental Henry constant is reproduced by the GCMC simulations. Knowing the site-surface interaction energy as given in eq 5a, the molecule-surface interaction energy is calculated from M

φi,s )

(a) φi,s ∑ a)1

(6)

For a slit pore of width H, defined as the distance between the plane passing through all carbon atoms of the outermost layer of one wall and the corresponding plane of the opposite wall, the solid-fluid potential is

Figure 1. VLE phase diagrams for ethylene: (a) linear scale of density, (b) log scale of density, and (c) vapor pressure versus T.

φi,P ) φi,s(z) + φi,s(H - z)

agreement for saturated vapor density. The good agreement between the simulation results and the data extends to temperatures close to the critical temperature. It is known that in many cases the GEMC simulations are somewhat unstable at temperatures very close to the critical temperature. In such cases, the power law for the critical region can be applied to determine the critical properties. To further check the adequacy of the UA potential model for the description of the VLE of ethylene, we present in Figure 1c the plot of the vapor pressure versus temperature. We see that the description of vapor pressure as a function of temperature is also equally good, with the exception of a small deviation when the temperature is either very close to the critical temperature or close to the triple point. For those temperatures, the molecular parameters need to be slightly adjusted to get the better agreement with the experimental data. This UA potential model will be used later in the description

(7)

3. Results and Discussions 3.1. VLE. We first test the VLE of ethylene with the UA-TraPPE model by performing the GEMC simulations at various temperatures. The parameters that we used in the GEMC are initial number of particles in each of the two simulation boxes (500) and the number of cycles (5000). Figure 1a,b shows the results of saturated densities obtained from the GEMC simulations versus temperature from 120 to the critical temperature of 282.4 K. Figure 1a shows the plot in linear scale to highlight the agreement between the GEMC simulation results for the saturated liquid density and the data, while Figure 1b presents the results in semilog scale to highlight the good (9) Steele, W. A. Surf. Sci. 1973, 36, 317-352.

7106


Do and Do

description of VLE is equally bad. Attempts were made to adjust the molecular parameters so that the VLE description is better with the one-center LJ model; we have found that this adjustment is temperature-dependent, making the potential of one-center LJ models less significant. We will not use the one-center potential model in our subsequent study of adsorption because of its inadequacy in the correct description of vapor pressure and saturated liquid densities because their correct description is essential in the correct description of adsorption isotherms. 3.2. Adsorption on Graphitized Thermal Carbon Black. In the adsorption of linear ethylene molecules on graphitized thermal carbon black, we assume the surface is homogeneous (no defect and no functional group). Although this assumption is ideal, the reproducibility of data obtained by various groups using different sources of graphitized thermal carbon black indeed suggests that the presence of defects or functional groups on this material is minimum. The distance between two interaction sites on the ethylene molecule is longer than the carbon-carbon bond length on the graphite layer, and, therefore, the use of the solid-fluid potential equation of Steele for structureless homogeneous is justified. Thus, in principle the application of the GCMC should describe the adsorption of ethylene on graphitized thermal carbon black. We use the following parameters for the carbon surface: collision diameter σss ) 3.4 Å, ss/k ) 28 K, and ∆ ) 3.35 Å. In the GCMC simulations, we use a slit pore having lengths in the x and y directions equal to 10 times the collision diameter and a cutoff radius of 5 times the collision diameter (half the box length), and the number of cycles is between 50 000 and 100 000. In each cycle, we have N displacements, N rotations and one attempt of either creation or destruction of a particle. Here, N is the number of ethylenes in the simulation box. The pore width used is 80 Å, which is large enough for that slit pore to mimic two independent surfaces. The adsorption isotherm is presented as plot of surface excess versus pressure. The surface excess is defined as

Γ) Figure 2. VLE description of ethylene with the one-center LJ potential model.

of adsorption of ethylene onto graphitized carbon black. The good description of vapor pressure is very essential in the study of adsorption on carbon black and in slit pores over the complete pressure range from 0 to P0 (vapor pressure). This is the case because when pressure is close to the vapor pressure multilayer adsorption occurs; therefore, if the vapor pressure is incorrectly determined from the simulation the multilayer adsorption will be badly described because of the very sharp change in the adsorption isotherm at pressures close to P0. The performance of the one-center LJ potential model is shown in Figure 2, and it is seen from this figure that the description of VLE is inadequate over the range of temperatures investigated. The description of the vapor pressure dependence on temperature is poor (note that we used the semilog plot), and that of the saturated liquid density is even much worse although the saturated vapor density is reasonably described. Even when we use another one-center LJ model suggested in Bird et al.,10 the (10) Bird, R.; Stewart, W.; Lightfoot, E. Transport Phenomena; Wiley: New York, 1960.

[

1 〈Ν〉 - Fb(H - σss) 2 LxLy

]

(8)

The factor 1/2 is because we have two surfaces in the simulation box. Because the solid-fluid interaction does not strictly conform to the Lorentz-Berthelot mixing rule, the binary interaction parameter is introduced as in eq 5b, and this is usually obtained such that the GCMC simulation produces the same Henry constant as the experimental data. We have found that this binary interaction parameter is ks ) -0.05 for the UA-LJ potential model considered in this paper. Figure 3 shows the adsorption of ethylene on graphitized thermal carbon black at 173, 183, and 193 K with the solid lines being the results from the GCMC simulations using the UA-LJ model and symbols being experimental data taken from Avgul and Kiselev.11 As seen from this figure, the description of the adsorption isotherm of ethylene on graphitized thermal carbon black at 173, 183, and 193 K is excellent, suggesting the adequacy of the UA-LJ potential model to describe not just the VLE (as we have shown in the last section) but also adsorption on graphitized thermal carbon black. The important observation is that the temperature dependency of adsorption is described well by the GCMC simulation model with a single set of molecular parameters. To (11) Avgul, N. N.; Kiselev, A. V. Chem. Phys. Carbon 1970, 6, 1-124.


Langmuir, Vol. 20, No. 17, 2004 7107

Figure 3. Plots of surface excess versus pressure for ethylene adsorption on graphitized thermal carbon black at 173, 183, and 193 K with the UA-LJ model of Wick et al.1

investigate the molecules’ orientation, we consider the local density distribution versus distance, the orientationdependent density distribution versus distance, and the angle formed between the molecular axis and the z axis. The local density distribution is defined as

F(z) )

∆N(z) LxLy∆z

(9a)

where ∆N(z) is the number of ethylene molecules whose centers of mass are located in the segment having boundaries at z and z + ∆z and Lx and Ly are the box lengths in the x and y directions, respectively. The orientation-dependent local density is calculated from

F(z, θ) )

∆N(z, θ) LxLy∆z sin θ∆θ

(9b)

where ∆N(z, θ) is the number of ethylene molecules whose centers of mass are located in the segment having boundaries z, z + ∆, θ, and θ + ∆θ. Here, θ is the angle between the axis of the molecule and the z direction. This orientation-dependent density allows us to evaluate the preferential orientation of the molecules at various distances from the wall surface. An angle of 0 means that the ethylene molecule lies vertically against the pore surface, while a value of π/2 means that it lies parallel to the surface. We shall present the density distribution plots as LJ density (nondimensional density) versus distance and angle. The LJ density and LJ orientation-dependent density are defined as

F*(z) ) σ3F(z)

(10a)

F*(z, θ) ) σ3F(z, θ)

(10b)

where σ is the collision diameter of ethylene, used for scaling the density in eq 10. The experimental data are extended only to the pressure range where only two layers are formed on the solid surface. To investigate the layering mechanism of molecules, we show in Figure 4 the local density distributions. Figure 4a shows the local density distribution versus distance from the surface for a specific pressure of 1 × 105 Pa and T ) 193 K, while Figure 4b depicts the threedimensional plot of this local density versus distance and the pressure. The density distribution of ethylene versus distance from the surface of graphitized thermal carbon black in Figure 4 shows that adsorption on the surface is clearly

Figure 4. Plot of the density distribution of ethylene versus distance from the carbon surface for 193 K.

a molecular layering process with the completion of the first layer before the second layer is filled with ethylene molecules. It is interesting to investigate the orientation of ethylene molecules in those layers, and Figure 5 shows the orientation-dependent density distribution plotted versus distance from the surface and the angle between the molecule axis and the z direction for P ) 3 × 105 Pa. We see that the dominant orientation of ethylene in the first layer is parallel to the surface (angle is π/2). On the other hand, the orientation of the second and higher layers does not show any preference between parallel and vertical arrangements (shown as a uniform distribution with respect to the angle for the second layer in Figure 5). This is expected because the molecules in the second and higher layers behave very much like those in the bulk liquid phase, and as a result, they do not show any preferential orientation. The question is whether the GCMC simulation describes correctly the adsorption beyond the second layer. Experimental data are thus needed to further test the potential of the GCMC simulations. However, in the light that, with the same set of molecular parameters, the VLE are wellproduced and the adsorption on graphitized thermal carbon black is produced well up to two layers, one can expect that the GCMC simulations would also describe well the higher layer coverage if the data in that range

7108


Do and Do

Figure 5. Plot of the density distribution of ethylene versus distance from the carbon surface and angle for 193 K.

were available because molecules in the third layer and higher layers behave like those in the bulk liquid and vapor, which we know from the VLE study that the MC simulation with this potential model can describe well. 3.3. Adsorption in Slit Pores. We now study the adsorption of ethylene in a slit pore. Three temperatures are selected in our investigation, 120, 160, and 200 K. The triple and critical temperatures of ethylene are 104 and 282.4 K, respectively. In the GCMC simulations of adsorption in slit pores, we use the following parameters. The box length is 10 times the collision diameter, and the cutoff radius is taken to be 1/2 the box length. The number of cycles in equilibration and in collecting statistics is 100 000, in each cycle we have N displacements (N is the number of particles), and the number of attempts to exchange a particle is chosen such that the acceptance rate is 2%. Maximum Pore Density versus Pore Width. It is known that the way molecules pack in small pores depends on the width of those pores and temperature as well as the shape and size of the adsorbate molecule. This has been a subject of intensive investigation by many workers.12-24 To investigate this, we study the maximum pore density versus pore width. The maximum pore density is the pore (12) Miyawaki, J.; Kaneko, K. Chem. Phys. Lett. 2001, 337, 243247. (13) Ohkubo, T.; Iiyama, T.; Suzuki, T.; Kaneko, K. Tanso 1998, 618619. (14) Ohkubo, T.; Iiyama, T.; Kaneko, K. Chem. Phys. Lett. 1999a, 312, 191-195. (15) Ohkubo, T.; Iiyama, T.; Nishikawa, K.; Suzuki, T.; Kaneko, K. J. Phys. Chem. B 1999, 103, 1859-1863. (16) Ohkubo, T.; Yang, C.; Pinero, E.; Solano, L.; Kaneko, K. Chem. Phys. Lett. 2000, 329, 71-75. (17) Watanabe, A.; Iiyama, T.; Kaneko, K. Tanso 1998, 624-625. (18) Watanabe, A.; Iiyama, T.; Kaneko, K. Chem. Phys. Lett. 1999, 305, 71-74. (19) Suzuki, T.; Iiyama, T.; Gubbins, K.; Kaneko, K. Langmuir 1999, 15, 5870-5875. (20) Kaneko, K. Supramol. Sci. 1998, 5, 267-273. (21) Suzuki, T.; Kaneko, K.; Gubbins, K. Langmuir 1997, 13, 25452549. (22) Iiyama, T.; Suzuki, T.; Kaneko, K. Chem. Phys. Lett. 1996, 259, 37-40. (23) Iiyama, T.; Nishikawa, T.; Suzuki, T.; Otowa, T.; Hijiriyamo, M.; Nojima, Y.; Kaneko, K. J. Phys. Chem. B 1997, 101, 3037-3042. (24) Kaneko, K.; Watanabe, A.; Iiyama, T.; Radhakrishan, R.; Gubbins, K. J. Phys. Chem. B 1999, 103, 7061-7063.

Figure 6. Maximum pore density versus pore width for ethylene with the UA potential model: (a) 120, (b) 160, and (c) 200 K.

density evaluated at the vapor pressure. We define pore density as the number of particles per unit volume available for ethylene molecules, that is,

Fav )

〈N〉 LxLy(H - σss)

(11)

where 〈N〉 is the ensemble average number of particles in the simulation box having lengths Lx and Ly along the x and y directions, respectively. Figure 6 shows the maximum pore density as a function of pore width with temperature as the varying parameter. The oscillation of the pore density versus pore width is an indication of the packing effect in small pores. This


Langmuir, Vol. 20, No. 17, 2004 7109

Figure 7. Local density distribution versus (a) distance and (b) distance and angle for the 6.5-Å pore.

Figure 8. Local density distribution versus (a) distance and (b) distance and angle for the 9-Å pore.

Figure 9. Local density distribution versus (a) distance and (b) distance and angle for the 10-Å pore.

oscillation in pore density has been reported by a number of workers in this area, for example, Olivier25 and Nicholson and Stubos.26 When pores become larger, such oscillation is rapidly diminishing (for pore widths greater than 30 Å). This plot is to show the packing efficiency in each pore. As a rule of thumb, a pore with high packing efficiency is the one that has high pore density, shown as peaks in Figure 6 (hereafter called perfect pores). (25) Olivier, J. Carbon 1998, 36, 1469-1472. (26) Nicholson, D.; Stubos, T. In Recent Advances in Gas Separation by Microporous Ceramic Membranes; Kanellopoulos, N., Ed.; Elsevier: Amsterdam, 2000; p 231.

For comparison, the liquid densities of ethylene at 120, 160, and 200 K are 22.66, 20.75, and 18.62 kmol/m3, respectively (shown as horizontal dashed lines in the same figure). Let us start with the case of 120 K. We see that only small pores that accommodate one layer could have a pore density greater than the liquid density. For example, in the 6.5-Å pore, the maximum pore density is about 25.5 kmol/m3, which is about 13% greater than the liquid density. However, for higher temperatures, 160 and 200 K, the pore density is consistently greater than the corresponding liquid density, particularly for 200 K. This suggests that molecules pack better in the bulk phase at

7110


Do and Do

Figure 10. Evolution of the local density distribution with pore width, from one layer to two layers. Pore widths are 7.5, 8.5, 9.3, and 9.45 Å.


Langmuir, Vol. 20, No. 17, 2004 7111

Figure 11. Evolution of the local density distribution with pore width, from two layers to three layers. Pore widths are 11, 13, 13.1, and 14 Å.

7112


Do and Do

Figure 12. Evolution of the local density distribution with pore width, from three to four layers. Pore widths are 14.5, 16.7, and 18 Å.

low temperatures, while they pack better in the confined space at high temperatures. This is due to the better degree freedom in aligning molecules in the bulk phase at low temperatures than that in the confined space. On the other hand, at high temperatures the high pressures exerted by the pore force the molecules to pack better in the confined space than in the bulk phase. Let us show the local density distribution for the 6.5-Å pore as shown in Figure 7. Here, we see that this very small pore can pack one layer and the packing is very efficient with very high LJ local density as shown in Figure 7a, and the orientation of all molecules is parallel to the wall surface as shown in the dominant peak in Figure 7b occurring at an angle of π/2. This parallel orientation is the only configuration because all other configurations are not possible in this very small pore.

The first imperfect pore is the 9-Å pore, shown as the first minimum in the plot of the maximum pore density versus pore width (Figure 6). The reason for this imperfect packing can be seen in the local density distribution plot as shown in Figure 8. First, we observe the rather broad peak in this pore and the presence of two small shoulders adjacent to the central peak. This is the indication of the onset of two-layer formation in this pore. Because this pore is too large for one layer while too small for two layers, the central peak represents ethylene molecules lying almost vertical to the wall surface, as reflected by the angle of 0 in Figure 8b. The preferential orientation of molecules in the two shoulders is parallel to the wall surface (seen in Figure 8b with an angle of π/2). It is this mix of parallel and vertical orientations that results in imperfect packing in


Langmuir, Vol. 20, No. 17, 2004 7113

Figure 13. Adsorption isotherms of ethylene at 120 K in slit pores having widths 8, 10, 13, 15, 16, and 20 Å.

Figure 14. Adsorption isotherms of ethylene at 200 K in slit pores having widths 10 and 20 Å.

this pore. Let us further highlight the packing effect with the 10-Å pore (perfect pore), which can accommodate exactly two layers (the peak in the plot of maximum pore density versus pore width of Figure 6). Using the argument we have just presented, if this pore is called a perfect pore it must have an integral number of layers and the molecules in those layers lie parallel to the wall surface. Indeed this is the case, as we show in Figure 9 the local density distribution for this 10-Å pore. The local distribution is very sharp and most of the molecules in both layers have an angle of π/2 with the z axis; that is, most ethylene molecules lie parallel to the surface. The maximum pore density for this perfect 10-Å pore is about the same as the liquid density, while that for the “imperfect” 9-Å pore is about 25% lower than the liquid density. Transition of Local Density Distribution from the OneLayer to the Two-Layer Configuration. We have shown the local density distribution for perfect pores of 6.5 Å (one layer) and 10 Å (two layers) and the so-called imperfect pore of 9 Å, falling between the above two perfect pores. It is interesting from the molecular orientation point of view to explore how the local density distribution evolves when pore width varies from 6.5 to 10 Å (i.e., number of layers in the pore varies from one to two). Figure 10 shows this evolution for a number of pore widths ranging from 7.5 to 10 Å.

To investigate this evolution, the ordinates of all plots are set at the same scale so that we can see the relative magnitude of all local densities. Starting with the 7.5-Å pore where we have a single layer as evidenced by very high local density and the dominant parallel orientation in this pore (dominant angle is π/2), the density distribution then shows shoulders in larger pores of 8.5 and 9.3 Å, which is an evidence of molecular rearrangement to start forming two layers. The central peak in these two pores has a vertical orientation with respect to the wall surface while the orientation of molecules of the two shoulders takes the parallel configuration. When the width is large enough (9.3 Å), these two shoulders get larger in magnitude while the central peak representing molecules in the central core (mostly vertical to the wall surface) is decreasing in magnitude. When the pore width is equal to 9.45 Å, the two shoulders now become two dominant peaks and the central peak disappears; this represents the onset of the formation of two parallel layers, and the perfect arrangement of two layers with parallel configuration is finally achieved with the 10-Å pore, as discussed earlier in Figure 9. Transition of the Local Density Distribution from the Two-Layer to the Three-Layer Configuration. The transition from two layers to three layers (H ) 14 Å) is also

7114


Figure 15. Isosteric heat of ethylene in slit pores of various widths versus surface loading.

interesting. We present below in Figure 11 the evolution of the density distribution for pore widths of 11, 13, 13.1, and 14 Å. When the pore width is 11 Å, the two-layer arrangement in the pore remains with some molecules starting to show orientations deviating away from the parallel orientation. When the pore width is increased to 13 Å, two layers are still visible, but now most of the molecules have a vertical orientation with respect to the wall surface. Some molecules are still in the parallel orientation, and they are represented by the two shoulders on the outside of the two major peaks (parallel-oriented molecules are now in the minority). When the pore width is increased by a very small amount, from 13 to 13.1 Å, we see a substantial change in the density distribution profile. The two shoulders in the 13-Å pore now become two major peaks representing molecules in the two contact layers close to the two wall surfaces, while the two major peaks in the 13-Å pore now become two shoulders and one new peak appearing at the center of the pore of 13.1 Å. The orientation of the two peaks close to the surfaces and the central peak is parallel orientation, while that of the two shoulders is vertical orientation. This new pattern persists for pore widths greater than 13.1 Å, with the two major peaks and the central peaks being slightly larger at the expense of the two shoulders; that is, more and more molecules are oriented toward parallel orientation. When the pore is increased to 14 Å, we observe that most molecules in all three layers will adopt the parallel orientation. This 14-Å pore is the perfect pore for accommodating three layers. Transition of the Local Density Distribution from the Three-Layer to the Four-Layer Configuration. The evolution from three layers to four layers is a bit different from the last two evolutions (Figure 12). Starting with the threelayer configuration (H ) 14.5 Å), when the pore width is increased the two peaks adjacent to the wall surfaces remain almost the same, both in magnitude and the parallel orientation, but two small shoulders start to appear on the inside of those two peaks. The orientation of these two small shoulders is vertical to the surface (see pore of 16.7 Å). On the other hand, the central peak starts to deconvolute into two peaks with the orientation of the deconvoluted peaks being parallel while the overlapping portion of these peaks has molecules mostly in a vertical configuration. As the pore is getting larger, these two deconvoluted peaks will become two separated peaks, and the two small shoulders disappear. The perfect pore of four layers is then achieved when the pore width is equal to about 18 Å. Transition of Local Density Distribution from the FourLayer to the Five-Layer Configuration. The evolution from

Do and Do

four to five layers is almost the same as that from three to four layers. When the pore width is increased from H ) 18 Å, the two peaks adjacent to the pore walls remain the same both in magnitude and parallel orientation, and like the previous case there are two very small shoulders appearing on the inside of these two peaks. At the same time, the fifth layer starts to evolve at the center of the pore, and when the pore width is increased further the central peak becomes larger in magnitude and separated from the neighboring peaks. Five perfect layers with dominant parallel orientation are achieved when the pore width is about 22 Å. Adsorption Isotherms at 120 K. Having studied the maximum pore density and the behavior of the local density distribution as a function of pore width, we now turn to the topology of the adsorption isotherms of ethylene at 120 K for a number of pore widths (Figure 13). Generally the smaller is the pore, the sooner is the adsorption, except for pores smaller than 6.5 Å (between 6.1 and 6.5 Å) for which the onset of adsorption is greater for pores having widths smaller than 6.5 Å. This is simply due to the strong repulsion in those small pores. Ethylene molecules cannot enter pores having a width smaller than 6.1 Å, which is defined as the sieving pore width for ethylene. For pores having a width smaller than 9.5 Å (less than perfect two-layer pore of 10 Å), the adsorption isotherm shows a gradual increase in the adsorption isotherm with respect to pressure. The fluid-fluid interaction of mostly single layers in these pores is not strong enough to induce the two-dimensional condensation. However, when the pore width is equal to 10, the pore can accommodate two layers, and for such a pore the fluidfluid interaction is due to those in the same layer as well as those across the layers. This means that the number of neighbors for each molecule is greater than that of smaller pores. It is this “enhanced” fluid-fluid interaction that gives rise to the two-dimensional condensation in this pore (i.e., transition from gaslike adsorbed phase to liquidlike adsorbed phase), and also in this pore we observe the hysteresis, which is consequential of the twodimensional condensation. This two-dimensional condensation is persistent for pores having width greater than 10 A until pore width is equal to 14 A, at which we have three layers residing in the pore (see Figure 11 for the local density distribution). A further increase in the pore width will result in a reversible adsorption-desorption isotherm, and this persists until the pore width has reached 20 Å, where the hysteresis is due to the threedimensional condensation. What we see in this 20-Å pore is the formation of the two contact layers adjacent to the two pore walls, and when these two layers are almost completed the two inner core layers are then completely filled at a pressure of about 150 Pa. Adsorption Isotherms at 200 K. Next, we show the isotherms for two pore widths H ) 10 and 20 Å at a higher temperature of 200 K (Figure 14). At this temperature, the two-dimensional transition in 10 Å pores observed earlier for 120 K no longer exists. This is simply due to the weaker fluid-fluid interaction at this temperature, which is proportional to zw/kT, where z is the number of neighbors and w is the characteristic fluid-fluid energy. If this fluid-fluid interaction energy does not vary with temperature, the strength of the fluid-fluid interaction is inversely proportional to temperature. Thus, at 200 K, this strength is less than that at 120 K, and as a result there is no two-dimensional transition for this temperature. Even for a larger pore of 20 Å, the hysteresis observed for three-dimensional condensation for the case of 120 K no longer exists for this case of a higher


Langmuir, Vol. 20, No. 17, 2004 7115

Figure 16. Plots of isosteric heats contributed by fluid-fluid interaction and fluid-solid interaction versus loadings for pore widths of 7, 9.5, 10, 15, 20, and 50 Å.

temperature of 200 K. Again, this is due to the weaker fluid-fluid interaction. Isosteric Heat. The utility of the GCMC simulations is not only in the generation of the adsorption isotherm and the microscopic configuration of molecules, but also in the calculation of various thermodynamics quantities. One of those is the isosteric heat. Using the fluctuation theory, the isosteric heat is calculated from27

qiso )

〈U〉〈N〉 - 〈UN〉 〈N2〉 - 〈N〉〈N〉

+ kT

where 〈〉 is the ensemble average, N is the number of particles, and U is the configuration energy of the system. This configuration energy can be divided into two contributions; one is due to the fluid-fluid interaction while the other is due to the fluid-solid interaction, that is, U ) Uff + Usf. Therefore, we can write the isosteric heat equation as follows:

qiso )

[

〈Uff〉〈N〉 - 〈UffN〉 〈N 〉 - 〈N〉〈N〉 2

]

+ kT +

above equation is the contribution from the fluid-fluid interaction while the remaining term is the contribution from the interaction between fluid and solid. Figure 15 shows the isosteric heat of ethylene versus surface loading for various pore widths. For small pores, the isosteric heats are greater as a result of the enhancement in the solid-fluid potential energies exerted by the two opposite walls. With respect to loadings, the isosteric heat increases with loading for small pores. This increase is due to the fluid-fluid interaction. For larger pores, the isosteric heat also increases with loading for low values of loadings. When monolayer coverage has been reached, the isosteric heat shows a sharp decrease toward the heat of liquefaction. When the pore is very large, we note that the isosteric heat at zero loading is about 16.8 kJ/mol, and this is in perfect agreement with the experimental value of 16.74 kJ/mol reported by Kiselev and Yashin28 for a graphitized thermal carbon black. Also for high loadings in a large

〈Usf〉〈N〉 - 〈UsfN〉 〈N2〉 - 〈N〉〈N〉

The square bracket term in the right-hand side of the

(27) Nicholson, D.; Parsonage, N. G. Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: London, 1982. (28) Kiselev, A. V.; Yashin, Y. I. Gas Adsorption Chromatography; Plenum Press: New York, 1969.

7116


pore (H ) 50 Å), the isosteric heat approaches the heat of liquefaction (12.1 kJ/mol), shown as an arrow in Figure 15. To show the separate contribution from the fluid-fluid and fluid-solid interactions, we plot in Figure 16 these contributions as a function of loading for pore widths of 7, 9.5, 10, 15, 20, and 50 Å. When the pore width is small (H ) 7 and 9.5 Å), the isosteric heat is dominated by the fluid-solid interaction as we would expect the great enhancement of the solid-fluid interaction due to the overlapping of the potentials exerted by the two opposite walls. When the pore width is greater, the heat contributed by the fluid-fluid interaction can be greater than that contributed by the fluid-solid interaction, as seen in the case of pore widths of 15, 20, and 50 Å. Conclusions We have presented the results of VLE and adsorption of ethylene on graphitized thermal carbon black and in slit pores, using UA-LJ models and one-center LJ models. It was found that the one-center LJ model is inadequate in the description of VLE, while the UA-LJ1 is acceptable models for use in VLE as well as in adsorption studies.

Do and Do

Adsorption behavior in terms of molecular configuration in slit pores is studied in detail, and it was found that the orientation of linear molecules is dominantly parallel in pores that can pack a perfectly integral number of layers. For imperfect pores, the orientation is a mix between parallel and vertical configurations, and it is found to be very sensitive to the pore width. The packing density in the pore depends on the pore width as well as temperature. At low temperatures, the packing density is generally less than the liquid density except for the very fine pore (containing only one layer), while at higher temperatures the packing density is generally greater than the bulk liquid density except for imperfect pores. The isosteric heat of adsorption is a function of loading and pore width. The contribution of the solid-fluid interaction in the isosteric heat of adsorption is dominant in small pores while in large pores it is dominant at low loadings and the fluid-fluid interaction dominates at higher loadings. Acknowledgment. Support from the Australian Research Council is gratefully acknowledged. LA0495682

Adsorption of Ethylene on Graphitized Thermal Carbon Black and in

Recommend Documents