2212
J. Phys. Chem. C 2007, 111, 2212-2222
Determination of Pore Accessibility in Disordered Nanoporous Materials Thanh X. Nguyen and Suresh K. Bhatia* DiVision of Chemical Engineering, The UniVersity of Queensland, Brisbane, Queensland 4072, Australia ReceiVed: August 29, 2006; In Final Form: NoVember 13, 2006
We propose a new algorithm based on application of cluster analysis to group adsorbate molecules of a highly dense adsorbed phase in the atomistic structural model of a disordered material into connected and disconnected clusters, through which pore network connectivity of the material is identified. Our proposed algorithm is then validated using a synthetic pore structure, as well as the reconstructed structure of a saccharose char obtained in our recent work using hybrid reverse Monte Carlo simulation.1 The algorithm also identifies kinetically closed pores in the latter structural model that are not accessed by adsorbate molecules at low temperature, at which their kinetic energy cannot overcome potential barriers at the mouths of pores that can otherwise accommodate them. The results are validated by transition state theory calculations for N2 and Ar adsorption, showing that N2 can equilibrate in narrow micropores at practical time scales at 300 K, but not at 77 K. Large differences between time scales for micropore entry and exit are predicted at low temperature for N2, the latter being larger by over 3 orders of magnitude, suggesting hysteresis. Similar behavior is predicted for Ar in the same char at 87 K. The results explain several long standing issues such as the observed increase of adsorption of nitrogen with an increase in temperature in coals,2,3 hysteresis phenomena in microporous carbons,4,5 and underprediction of adsorption of supercritical gases using structural parameters extracted from subcritical adsorption of nitrogen.5 Finally, the determination of pore accessibility and connectivity in disordered porous carbons using our proposed model enables one to obtain correct adsorbed quantities as well as selfdiffusivities and transport diffusivities using conventional grand canonical Monte Carlo and molecular dynamics simulations.
1. Introduction The understanding of the microstructure of porous carbon materials in terms of pore morphology and topology plays a key role in the quantitative prediction of adsorption equilibrium and dynamics. Such an understanding is crucial for optimal design of industrial adsorption systems, and for the preparation of advanced porous carbon materials with precisely controlled microstructure for high adsorption capacity and specific mixture separation. Most commonly, it has been based on a slit pore model characterization, obtained through interpretation of the adsorption equilibrium isotherm of simple gases at suitable temperature.6-8 Such a characterization provides a simplified representation of the microstructure of porous carbons as a collection of independent slit pores. The most widely used version of the slit pore model utilizes the infinitely thick wall assumption, but this has been superseded by new development in our laboratory,9,10 whereby both pore size and pore wall thickness heterogeneities are considered. The latter approach has been shown to provide characterization of the pore wall structure of different types of activated carbons consistent with X-ray diffraction and transmission electron microscopy results.11,12 Quantitative prediction of the adsorption equilibrium measurement of simple gases in BPL and Norit R1 Extra carbons at supercritical temperatures has been obtained13 using this approach. Unfortunately, porous carbon materials are not normally crystalline but contain significantly varying contents of amorphous matter. Such amorphous content inevitably leads to * To whom correspondence may be addressed. Phone: +61 7 3365 4263. Fax: +61 7 3365 4199. E-mail:
[email protected].
formation of complex pore networks with imperfect connectivity inside the microstructure. This incomplete connectivity physically or kinetically prevents adsorbate molecules from accessing even pores that are otherwise large enough to accommodate them. A prominent example of this is the difficulty of obtaining useful results from nitrogen adsorption in microporous carbons at 77 K, because of the diffusional limitations at this low temperature and resulting poor accessibility of pore spaces connected by narrow necks or pore mouths.14 Reverse Monte Carlo (RMC)15 and hybrid reverse Monte Carlo16 techniques have recently emerged as promising methods that enable one to generate atomistic models of the microstructure of porous carbons.1,17-19 The reconstructed atomistic model is realistic and permits estimation of equilibrium and transport properties of porous carbons using molecular simulation techniques such as grand canonical ensembles Monte Carlo (GCMC) simulation and equilibrium molecular dynamics (EMD) simulation. Such properties can then be directly compared with corresponding experimental measurements without further assumptions. These appealing features overcome related limitations of the slit pore model. However, similar to real porous carbons, the reconstructed atomistic microstructure also possesses a complex pore network that leads to difficulty in rigorous estimation of the adsorption isotherm using GCMC simulation. This is because adsorption in disconnected regions or closed pores also arises from insertion moves of the technique, which is unphysical. At the same time, the pore connectivity problem also causes difficulty in preparation of initial configurations for EMD simulation using the GCMC simulation or random placement of adsorbate molecules in the reconstructed carbon
10.1021/jp065591f CCC: $37.00 © 2007 American Chemical Society Published on Web 01/12/2007
Pore Accessibility in Disordered Nanoporous Materials structure model, as such methods do not distinguish between open and closed pore spaces. Brennan et al.20 and Biggs et al.21 suggested using geometric criteria to determine pore accessibility in porous materials. While such methods may be suitable for extremely narrow or permanently closed pores, they do not capture the kinetic feature of pore accessibility. In order to overcome the above difficulties, we propose a new algorithm that groups molecules of a close packed adsorbed phase into continuous and disconnected clusters, from which pore network connectivity is determined. Further, we validate the proposed algorithm by the use of a synthetic example of simulated porous carbons with predetermined disconnected and connected regions. Finally, we apply the algorithm to analyze the pore network connectivity of the reconstructed microstructure of a saccharose char, obtained in our recent work using hybrid reverse Monte Carlo simulation.1 Through this technique, kinetically closed pores in the reconstructed microstructure were identified. Subsequently, this kinetic feature of the connectivity was further confirmed in the experimental adsorption time scale by estimation of the crossing time of a single particle through connected neighboring pores using transition state theory (TST), successfully employed by several authors to predict the self-diffusivity of gases in zeolite materials.22,23 From the estimation of the crossing time, the problems arising from apparent pore-blocking effects are elucidated. 2. Mathematical Modeling 2.1. Classification of Closed Pores. We define a closed pore space as space which is inaccessible to the adsorbed molecule when the latter’s kinetic energy is insufficient to overcome the energy barrier at the pore mouth. We distinguish between two kinds of closed pores, physically closed pores and kinetically closed pores, as follows. 2.1.1. Physically Closed Pore. A physically closed pore occurs when the energy barrier at the pore mouth is too high for adsorbate molecules to overcome even at the highest temperature where rigidity of the solid structure is still reasonably valid, or Ea f ∞. In this situation, the physical size of the windows of the physically closed pore is normally smaller than the adsorbate molecular size, and interaction between the adsorbate molecule and solid atoms of the window or pore mouth falls into the steep repulsive region of the fluid-solid potential. 2.1.2. Kinetically Closed Pore. A kinetically closed pore is pore space inaccessible to adsorbate molecules at a particular temperature, due to insufficient adsorbate kinetic energy to overcome the energy barrier at the pore mouth at this temperature. However, the kinetically closed pore becomes open at a higher temperature at which adsorbate molecules gain sufficient kinetic energy to overcome the energy barrier at the pore mouth. The energy barrier, in this case, is contributed not only by the slightly repulsive interaction between solid atoms of the pore mouth and the adsorbate molecule but also by the local potential minimum near the pore mouth. The latter contribution arises from the fact that the carbon surface in porous carbons is not normally smooth, but rough and twisted, as a result of defects. This can lead to the formation of local energetic minima near the pore mouth outside the pore space. Consequently, an adsorbate molecule having low kinetic energy is trapped in these regions before passing through the pore mouth to enter the pore space. However, the adsorbate molecule will jump out of the regions of local energy minimum to access the pore space as it gains sufficient kinetic energy at higher temperature. This leads to a kinetically closed pore, which is explained in more detail in a later section.
J. Phys. Chem. C, Vol. 111, No. 5, 2007 2213 2.2. Determination of Pore Network Connectivity. Our approach to determine pore connectivity analyzes the connectivity of the fluid phase in the model of the solid structure. Since temperature (i.e., kinetic energy) is the main factor responsible for pores being kinetically closed, our analysis is based on the persistence of this effect even at high densities corresponding to that of the liquid phase, that is, when the average intermolecular distance, dha, is approximately equal to the equilibrium separation, 21/6σf, based on the Lennard Jones (LJ) pair potential.24 Here, σf is the fluid LJ size parameter, or collision diameter. This assumption covers most practical applications of gas adsorption using porous carbons, and the high density phase obtained in this condition is considered here as the “close packed” phase. In particular, the close packed adsorbed phase of argon and nitrogen can be obtained using grand canonical Monte Carlo simulation (GCMC) at 77 and 87 K, respectively. For other cases, the position of the first peak of the pair distribution function of the adsorbed phase is used to identify if a close packed adsorbed phase is reached. In our technique, we first fill up all pore spaces in the solid structure, including closed pores, with fluid, to form a close packed phase. Subsequently, the structure of the fluid phase is analyzed to determine if the atomistic solid structure subdivides the guest molecules into discrete clusters. Such discrete clusters signify closed pores, and their absence is indicative of continuous open pore space. Consequently, the determination of the pore network connectivity in the atomistic model of a porous carbon structure is based on analysis of the continuity of the close packed adsorbed phase throughout the solid structure model. This approach is central to our proposed algorithm to determine the pore connectivity in porous carbons. The algorithm is comprised of three principal steps to determine pore connectivity. 2.2.1. Filling All Pore Spaces of the Solid Structure with Close Packed Adsorbate. This step is straightforwardly performed using the conventional grand canonical Monte Carlo (GCMC) simulation technique, commonly used to study the adsorption of confined fluids. In GCMC simulation, the chemical potential, µ, volume, V, and temperature, T, are kept invariant, while the number of particles, N, and associated configurational energy, E, are allowed to fluctuate. In the procedure, microstate configurations are generated through the well-established Metropolis sampling scheme for three trial types: moving, creating, and deleting molecules. The probability for each trial type being accepted is given by Adams’ algorithm.25 In order to obtain close packing of the adsorbed phase, we simulate adsorption equilibrium at the condition of interest such that all the kinetically closed pores are detected at this condition. For distinguishing physically closed pores from kinetically closed pores, the adsorption simulations are performed at a reasonably high temperature at which kinetically closed pores are absent. In the current work, the configuration of the adsorbed phase has been collected after every 10 million GCMC steps to ensure that detection of closed pores is independent of the microstates. The total number of configurations used for every GCMC run is 60 million. 2.2.2. Determination of All Adsorbate Clusters in the Solid Structure Model. In order to determine the clusters, the maximum intermolecular distance between two first nearest neighbor adsorbate molecules in the same adsorbed cluster, RC[i][i], and the minimum intermolecular distance between two adsorbate molecules in two separate adsorbate clusters, RC[i][j], must be assigned. Here, indices i and j inside the square brackets denote clusters i and j. Clearly, RC[i][j] must be greater than
2214 J. Phys. Chem. C, Vol. 111, No. 5, 2007
Nguyen and Bhatia
Figure 1. Pair distribution function of the close packed adsorbed phase of nitrogen in saccharose char1 at 77 K and 1.013 bar.
RC[i][i] in order to obtain complete resolution among the clusters. For a close packed adsorbed phase, RC[i][i] can range from 21/6σf to the mean intermolecular distance for the first nearest neighbor, represented by the first peak of the pair distribution function of the adsorbed phase, as shown in Figure 1. From this figure, it can be seen that the first and second peaks do not resolve completely, indicating some second nearest neighbor molecule can be included if the RC[i][i] is chosen as in Figure 1. At the same time, this may also lead to exclusion of some of the first neighbor molecules. If the first case is true, it does not affect correct determination of clusters, due to the fact that if the first nearest neighbor of a centered adsorbate molecule belongs to a particular cluster its second neighbor must also belong to that cluster. Further, the second case does not occur if closed packing of the adsorbed phase is reached. In all results presented in the current work, it has been found that the following selection for RC[i][i] is satisfactory
RC[i][i] ) 21/6σf + 0.1σf ) 1.2225σf
(1)
For selection of RC[i][j], it can been inferred that RC[i][j] is the shortest intermolecular distance between two adsorbate molecules separated by a graphite sheet perpendicular to the line connecting the centers of these molecules. Accordingly, RC[i][j] is given as
RC[i][j] ) 2σmin sf
(2)
where the second term, σmin sf , in the right-hand side is the closest distance between the adsorbate molecule and the single lattice plane, taken as
σmin sf ) σsf
(3)
Here, the cross Lennard Jones (LJ) collision diameter between a carbon atom and an adsorbate molecule, σsf, is given by the Lorentz-Berthelot mixing rule as
1 σsf ) (σf + σs) 2
(4)
where σs is the carbon LJ collision diameter. Combining eqs 1-4, the largest collision diameter of the adsorbate molecule, σmax f , which provides complete resolution of two clusters is
) 4.495σs σmax f
(5)
Thus, if the carbon LJ collision diameter is chosen to be 0.34 is 1.5 nm. This upper limit of σmax covers nm, the value of σmax f f a wide range of molecules. In practice, large and complex molecules are normally modeled as an array of atom sites linked by chemical bonds. In this case, it can be seen that it is straightforward to apply our proposed algorithm to these large molecules. This is because the chemical bond length between two atom sites in a molecule is normally smaller than the physical bond length between two atom sites in two different molecules. Thus, the atom sites of a large molecule can be treated as part of a cluster. Consequently, the problem of determining whether two large molecules are in the same cluster is now to determine whether any atom site in one molecule belongs to the neighboring cluster and vice versa. This is within the scope of our proposed algorithm. Thus, in this way, the algorithm groups close packed atoms of the adsorbed phase into clusters. However, in complex structures with warped and twisted sheets, situations are possible in which adsorbate molecules do not pack closely in the pore space, because the confining pore space cannot accommodate an additional adsorbate molecule. This may lead to misidentification of adsorbate molecules nominally belonging to part of a different cluster outside this pore space. To avoid this misidentification, more intermediate GCMC configurations of the close packed adsorbed phase are collected to be able to capture the closest distance between two molecules inside and outside the pore space. Finally, we determined all the clusters in the close packed adsorbate phase, using the values of RC[i][i] and RC[i][j] given in eqs 1 and 2, respectively. In particular, the structure of the close packed adsorbed phase is considered as a graph in which adsorbate molecules are vertices. Each vertex connects to its neighboring vertices determined using RC[i][i]. The Breadth-first search algorithm26 was employed to traverse all distinct pathways connected to two arbitrary vertices. During this traversal, visited vertices are marked to guarantee each pathway is not repeatedly visited. Accordingly, two molecules are in the same cluster if there exists at least one pathway to connect them. By this way, all the clusters are determined. 2.2.3. Determination of Continuity of Each Cluster throughout the Structure. After determination of all the adsorbate clusters in the carbon structure model as described above, the final step determines the connectivity of these clusters, and subsequently the pore connectivity of the carbon structure. Due to the periodic boundary conditions (PBC) applied to the primitive unit cell, the carbon structure model is considered as a periodic solid. However, instead of determining the pore connectivity in a periodic unit cell, we perform this task in a periodic supercell constructed from this primitive unit cell. The periodic supercell preserves exactly the pore connectivity feature of the periodic structure. However, the advantage of the periodic supercell over the periodic unit cell is that it explicitly considers structural properties across boundaries. In particular, the relationship between the position of particles in the periodic supercell, rs, and that in the periodic unit cell, ru, is given as
rs ) L‚(ru + n)
(6)
where the matrix L represents the Bravais lattice vectors, defined by the three unit cell vectors Lx ) [Lx, 0, 0], Ly ) [0, Ly, 0], and Lz ) [0, 0, Lz]. Here, Lx, Ly, and Lz are the lengths in the three dimensions of the simulation box. n is a repeating vector, whose three components, nx, ny, and nx, are integers. For the periodic solid, n{} ) 0, 1, without loss of generality. Parts a and b of Figure 2 illustrate a 2D periodic unit cell and the corresponding periodic supercell, respectively. The
Pore Accessibility in Disordered Nanoporous Materials
J. Phys. Chem. C, Vol. 111, No. 5, 2007 2215 time, τ, from a cage A to a neighboring cage B is generally given as
τAfB )
1 kAfB
(8)
where kAfB is a diffusion rate constant, given by TST27-29 as
kAfB ) κ
x
kBT 2πm
∫DSe-βφ (r) d2r ∫V e-βφ (r) d3r sf
(9)
sf
cageA
Here, κ is a transmission coefficient that represents the fraction of particles starting on top of the barrier with velocity toward cage B that successfully reach cage B, kB is the Boltzmann constant, T is the temperature, and m is the mass of the particle. φsf is the interaction potential between the adsorbate particle i at position r and all solid atoms of the adsorbent phase, given as Figure 2. Illustration of the close packed adsorbed phase (black balls) in (a) the periodic unit cell and (b) its corresponding periodic supercell. The black slanted stripes depict the carbon phase.
φsf(r) )
u(|rj - r|) ∑ j)1
(10)
where u is the LJ (12-6) solid-fluid pair potential, given as portions with right slanted black stripes represent the solid phase, and balls depict adsorbate molecules of the close packed adsorbed phase. From these figures, it may be observed that the number of adsorbate molecules in a cluster associated with a closed pore is inVariant from the periodic unit cell to the periodic supercell, while that associated with open space is doubled from the periodic unit cell to the periodic supercell. This observation is the key principle in determining pore network connectivity in the carbon structure model. Accordingly, the degree of connectivity of a pore space, Nf, is given as
u(r) ) 4sf
[( ) ( ) ] σsf r
12
-
σsf r
Nf )
N{i} u
(11)
If we can choose the width, b, of a small cubic box, which contains the dividing surface, to be sufficiently small such that the energy landscape adjacent to the saddle point is similar for surfaces parallel to the dividing surface, the integral term on the right-hand side (R.H.S.) of eq 9 can be rewritten as
∫DSe-βφ (r) d2r 1 ∫V e-βφ (r) d3r ) ∫cageAe-βφ (r) d3r b ∫cageAe-βφ (r) d3r sf
sf
N{i} s
6
box
-1
(7)
{i} where N{i} u and Ns are the number of adsorbate molecules in the ith cluster in the periodic unit cell and in the periodic supercell, respectively. Thus, a cluster is associated with a closed pore if Nf is equal to zero and otherwise to an open pore, for which Nf varies between unity and (2d - 1). Here, d is the number of dimensions in which periodic boundary conditions are imposed. Thus, for a three-dimensional system, Nf can be as high as 7. 2.3. Transition State Theory. 2.3.1. General Theory. As shown in our previous work,1 the atomistic structural model of the saccharose char displays a small fraction of volume inaccessible to nitrogen at 77 K. In the present work, this atomistic model was employed to validate our proposed model for determination of its pore network connectivity. Accordingly, diffusion of nitrogen at 77 K through constricted windows, which connect neighboring pores, may be a rare event whose observation time scale is far beyond that accessible by molecular dynamics (MD), which ranges from 10-8 to 10-7 s. Thus, the estimation of the crossing time of a single molecule from one pore to the neighboring pore at experimental time scales requires other methods. Here, we employed TST23,27-31 to estimate the crossing time of a single adsorbate molecule through pore mouths between cage A and cage B in the structural model of saccharose char, as shown in Figure 9. Such calculations provide insight into the behavior at low densities. According to TST, the crossing
sf
(12)
sf
The two integrals on the R.H.S. of eq 12 can be directly evaluated using the Monte Carlo sampling integration.32 Accordingly, eq 12 is rewritten for the case of one particle in cage A as
Vbox
∫DSe-βφ (r′) d2r′ sf
∫cageAe-βφ (r)/K T d3r sf
)
B
τ′max
∑ e-βφ (r′) sf n
1 τ′max n′)1
(13)
b VcageA τmax τmax
∑ e-βφ (r ) sf
n
n)1
where τ′max and τmax are the number of Monte Carlo (MC) trials. When the particle in cage A falls in the dividing surface, it is expected to vibrate around the saddle point of this surface due to the extremely high energy around this saddle point. Accordingly, for a sufficiently large number of MC trials, eq 13 can be equivalently expressed as
∫DSe-βφ (r) d2r sf
∫cageAe
-βφsf(r)
3
dr
)
1 e-β〈φsf〉N,Vbox,T b
τmax
(14)
e-βφ (r ) ∑ n)1 sf
n
where τmax is the number of grid points in cage A, which is now approximately equal to the ratio of VcageA to Vbox. rn is the
2216 J. Phys. Chem. C, Vol. 111, No. 5, 2007
Nguyen and Bhatia
coordinate of grid points in cage A. Thus, for a single particle, the crossing time, τAfB, from cage A to cage B is finally obtained as
τAfB )
x
b κ
2πm βEa e kBT
(15)
where β ) 1/kBT and Ea is the effective activation energy, given as τmax
Ea ) 〈φsf〉N,Vbox,T + kBT ln(
e-βφ (r )) ∑ n)1 sf
n
(16)
2.3.2. Implementation. From eq 16, it is evident that calculation of the crossing time, τAfB, requires specification of several parameters. These are the transmission coefficient, κ, the width of the box at the barrier, b, and the effective activation energy, Ea. Our choice of value for these is discussed below. 2.3.2.1. Effective Activation Energy (Ea). From eq 16, it can be seen that the effective activation energy, Ea, is determined by the two terms on the right-hand side of this equation. The second term, or the effective free energy of the cage, is evaluated as the summation of Boltzmann factors of the solid-fluid interaction potentials at every grid point in the cage, while the first term or top barrier energy can be evaluated only if the dividing surface is known. Determination of the latter is a nontrivial task for disordered structures. Here, these two terms have been evaluated as follows. a. The Effective Free Energy. In order to evaluate the effective free energy in the cage, we first generate a grid network for the unit cell of the reconstructed saccharose char. In the current work, we employed a grid size of 0.05 nm which is approximately equivalent to the top barrier box dimension. Further, it is noted here that the result of the effective free energy is unaffected by small differences between grid size and the size of the top barrier box for the investigated low-temperature range (
