Molecular Dynamics of Carbon Dioxide, Methane and Their Mixtures

J. Phys. Chem. B 2009, 113, 13761–13767

13761

Molecular Dynamics of Carbon Dioxide, Methane and Their Mixtures in a Zeolite Possessing Two Independent Pore Networks as Revealed by Computer Simulations† Marco Sant,‡ Jean-Marc Leyssale,§ George K. Papadopoulos,‡ and Doros N. Theodorou*,‡ School of Chemical Engineering, National Technical UniVersity of Athens, 9 Heroon Polytechniou Street, Zografou Campus, 15780 Athens, Greece and CNRS, Laboratoire des Composites ThermoStructuraux, CNRSsCEAsSnecma Propulsion SolidesUniVersite´ Bordeaux 1, 3 alleé de La Boe¨tie, 33600 Pessac, France

Downloaded by UNIV OF CALIFORNIA SAN DIEGO on September 6, 2015 | http://pubs.acs.org Publication Date (Web): August 4, 2009 | doi: 10.1021/jp902829j

ReceiVed: March 29, 2009; ReVised Manuscript ReceiVed: July 7, 2009

The molecular motion of methane (CH4) and carbon dioxide (CO2) sorbed in the two independent pore networks, being termed hereafter as large cavity (LC) and sinusoidal channel (SC) regions of the siliceous MWWframework-type zeolite ITQ-1, is studied by means of atomistic computer simulation. Equilibrium molecular dynamics predicts different loading dependences for the two gases, for both the self and the collective (Maxwell-Stefan) diffusion coefficients; in particular, the transport coefficients of CH4 go through a maximum as its loading in the zeolite increases, whereas CO2 dynamics exhibits the decreasing trend with loading usually observed in nanoporous materials. The different loading dependence of the self-diffusivity for the two sorbates is attributed to their different probability density distribution within the supercages in the LC system of the ITQ-1 unit cell. The composition and occupancy dependence of the self-diffusivity of each component in their binary mixtures can be explained in terms of the selectivity for CO2 sorption thermodynamics in the zeolite. The collective diffusivity loading dependence of the single and binary sorbate system is explainable on the basis of the strength of intermolecular interactions along the diffusion direction connecting the supercages by invoking the quasichemical mean field theory. 1. Introduction The advent of multiprocessor machines has enriched the field of transport phenomena in microporous media with a tremendous wealth of simulation results, which have appeared in the literature over the last twenty years. In many cases, computational researchers have labored to improve decisively the existing phenomenology, or to devise numerical methodologies for solving complex statistical physics problems pertaining to the modeling of fluid phases, phase transitions, interfaces and thin layers. An insightful analysis of the latter phenomena from the statistical mechanics perspective can be found in the monograph of H. T. Davis.1 The flow of molecular films in micropores of widths commensurate with the size of sorbed molecular size has attracted particular attention. Davis and co-workers adapted for the first time the diffusion theory of Enskog, combined with computer simulation, to an inhomogeneous fluid located in a model slit pore; via both theory and simulations they revealed significant variations of the tracer (self) diffusivity with pore width.2 These predictions were verified by the local average density model developed later by Davis and co-workers.3 Since then, several real nanostructured systems have been studied by means of atomistic simulation; admolecules on surfaces or in macropores and sorbate molecules in meso- or nanopores are typical examples. Sorbed fluids subject to high potential energy fields exerted by the framework atoms, offframework ions or organometallic groups embedded in the structure of a porous medium are encountered in a wide range of materials, ranging from carbon based nanomaterials and †

Part of the “H. Ted Davis Special Section”. * Author to whom correspondence should be addressed. Tel: +30-210772-3157. Fax: +30-210-772-3112. E-mail: [email protected]. ‡ National Technical University of Athens. § CNRSsCEAsSnecma Propulsion SolidesUniversite´ Bordeaux.

zeolite crystals to zeolite imidazolate, metal organic frameworks, and hyperbranched aminosilicas. The molecular modeling of mass transport as a function of the chemical constitution, nanopore geometry and topology, crystal and particle morphology of microporous materials has been the subject of many investigations over the past years. A treatment of mass transport phenomena via atomistic simulation presents the advantage, in comparison with macroscopic phenomenology, of predicting the magnitudes of diffusivities as well as their temperature, loading and composition dependencies. A literature survey on this subject can be found in refs 4-8 and in the rich amount of references therein. In the present article we report results from a simulation study of molecular motion in fluid phases sorbed within a relatively new microporous medium. Our study concerns binary mixtures of carbon dioxide and methane in the recently synthesized ITQ-1 zeolite,9 which is the siliceous analogue of zeolite MCM-22. The purely siliceous nature of the nanoporous medium allowed us to focus on the effects of framework structure and topology and interaction energetics in the absence of complications arising from counterbalancing cations. ITQ-1 belongs to the framework type code MWW; the high void space density of the MWW topology, together with its unique crystal structure, makes it very interesting in the context of gas adsorption and separation. Moreover, it is of particular interest to natural gas separation, as it is a hydrophobic material. The unit cell of ITQ-1 (Figure 1) possesses an unusual framework structure made of two independent large pore systems, offering a variety of possible sorption sites.10 One of the pore systems is a two-dimensional network of sinusoidal ten-member ring channels (hereafter referred to as SC), with a pore diameter of 0.4-0.55 nm, and the other (hereafter referred to as LC) consists of twelve-member ring large cavities (supercages) 1.7 nm in height × 0.8 nm in diameter, connected

10.1021/jp902829j CCC: $40.75  2009 American Chemical Society Published on Web 08/04/2009

13762

J. Phys. Chem. B, Vol. 113, No. 42, 2009

Sant et al. 2

JR ) -

∑ DRβ · ∇Fβ, R ) 1, 2

(2)

β)1

where the transport diffusivity tensor DRβcouples the density gradient of species β to the flux of species R. An analogous expression to eq 2 is given through the phenomenological development due to Onsager which introduces the chemical potential gradient, ∇µ, as the driving force for diffusion; thus, eq 1 takes the form of 2

JR ) -

∑ LRβ · ∇µβ, R ) 1, 2

(3)


β)1

Figure 1. Lowest value isodensity surfaces for CH4 in ITQ-1 at a total loading of 11.4 molecules per unit cell (u.c.) at 300 K and P ) 2.575 atm: 6.6 molecules/u.c. find themselves in the large cavity system (dark gray), while 4.8 molecules/u.c. are in the sinusoidal channel (light gray) system. LC and SC form two independent pore networks in the ITQ-1 unit cell.

by ten-member ring narrow windows. In the following we present an investigation of the loading-dependent diffusivities of CO2, CH4, and their mixtures with reference to their concentration distribution inside these two independent pore systems.

where now L is a symmetric matrix constituting the so-called Onsager transport coefficients. The latter equation is particularly useful in molecular simulations, since the collective diffusion coefficients, L, are readily computed under equilibrium conditions through the Green-Kubo expressions,12 and hence, through the equivalent Einstein form,12,14 from the displacements of the centers of mass of swarms of R or β type molecules with time. In terms of the center of mass position vectors, rRi, of individual molecules i of type R and β, for long correlation times, t, one obtains NR

LRβ )

∑

d 1 [r (0) lim 〈[( tf∞ 2doVkBT dt i)1 Ri Nβ

2. Molecular Dynamics

(4)

j)1

We study the dynamics of molecules sorbed in ITQ-1 by means of equilibrium molecular dynamics (MD) simulations in the NVT ensemble, by utilizing the relationship of transport coefficients and time correlation functions expressed by the Onsager regression hypothesis11,12 as it was later reformulated in the context of linear response theory.12 2.1. Theory. A general expression for the isothermal steadystate single component sorbate flow in the absence of external force fields can be written in Fickian form as follows:

In mixed sorbate systems, we will call the quantity D0,R) LRR kBT/FR the collective diffusivity of species R in the mixture; whereas in the special case of a pure sorbate, we will refer to the above quantity as D0. For a two-component sorbate system the relation between the chemical potential of species R and its fugacity, fR, reads 2

J ) -D · ∇F

∑ [rβj(0) - rβj(t)])]〉

rRi(t)]) · (

∇µR )

(1)

The above dot product provides the column vector J, the overall macroscopic flux, expressed as the mean number of molecules passing per unit cross-sectional area in a directional (anisotropic) medium due to a gradient in the density, F; D is the transport diffusivity tensor. In zeolite crystals or textured materials, transport and, consequently, self-diffusivity and collective diffusivity have a tensorial character; in other words, the flux vector depends upon the spatial orientation of the applied density gradient with respect to the crystal symmetry axes.10 We will use the symbol D to denote the orientationally averaged trace of the diffusivity tensor, namely, (1/3)Tr(D). An alternative to the preceding Fickian phenomenological approach can be found in the Stefan-Maxwell formulation which is based on a balance of forces exerted on the diffusing particles.13 For a two-component sorbate in a 3-D medium, the analogue of eq 1 for the flux of component R is written as

( )

kBT ∂ ln fR Fβ ∂ ln Fβ β)1

∑

∇Fβ, R ) 1, 2

(5)

T,F3-β

Combination of eqs 3 and 5 and comparing with eq 2 reveals the connection between the transport and Onsager diffusivity coefficients as follows:

D ) L·Γ

(6)

where Γ is the square matrix with elements the thermodynamic factors,

ΓRβ )

( )

kBT ∂ ln fR Fβ ∂ ln Fβ

R, β ) 1, 2

(7)

T,F3-β

From eq 6 is seen that tensor D is not necessarily symmetric. The self-diffusivities Ds,R were calculated through the formula

Molecular Dynamics of CO2, CH4 and Their Mixtures

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13763

NR


Ds,R

∑

d 1 ) lim 〈 [r (0) - rRi(t)]2〉 2doNR t f∞ dt i)1 Ri

(8)

The reader is reminded that diffusivities obtained by eqs 4 and 8 are orientationally averaged quantities with respect to dimensionality of the system (do being 1, 2 or 3). 2.2. Computation Details. Zeolite Reconstruction. We digitally reconstructed the ITQ-1 unit cell in accordance to the crystallographic data of Camblor et al.9 For convenience with periodic boundary conditions we used a double unit cell instead of the original hexagonal one. This new unit cell has an orthorhombic shape with axes (N.B. that the following transformation has been misprinted in ref 10 and should be appropriately corrected), a ) 2 sin(60°)a′ ) 2.4609 nm, b ) b′ ) 1.4208 nm, c ) c′ ) 2.4945 nm, where a′, b′, c′ are the original hexagonal unit cell axes.10 In this way, a base-centered orthorhombic unit cell containing 144 silicon and 288 oxygen atoms was created with its apexes located at the centers of large cavities of the zeolite. The whole zeolite was treated as a rigid body in all calculations, as in the preceding work of ref 10. Adopting a rigid zeolite model allowed us to use a pretabulation-interpolation mapping scheme10 for the guest-host interactions, for both the short range dispersive and the long range electrostatic interaction types. Moreover, since the guest-host interactions were tabulated, they had to be computed only once for all potential map grid points and for all sorbate atom-sorbent pairs; this makes the use of Ewald summation very inexpensive in that context. Energetics. The parametrization of the energetics for all interactions in the system under study is presented in Table 1. The dispersion interactions of both sorbates with Si atoms were ignored. Lennard-Jones parameters and partial charges of zeolite atoms were taken from ref 10 (see Tables 1 and 2). The CH4 molecule was described as one united atom. For the majority of molecular dynamics runs, the parameters of ref 15 were used for the methane interactions (first and second rows of Table 1); these parameters resulted from an ad hoc study toward a unified hydrocarbon potential for MFI siliceous structures. An alternative set from ref 16 (third and fourth rows of Table 1) was used in order to examine the sensitivity of simulation results to the parametrization. The CO2 molecule was treated as a rigid dumbbell bearing partial charges with a C-O bond length of 1.149 Å, taken from ref 17 (see Tables 1 and 2). Algorithms. The primary simulation box consisted of a number of unit cells varying from 3 × 4 × 1 at 52 atm up to a maximum of 40 × 70 × 1 at 0.01 atm, so that in all runs a number of sorbed molecules N = 500 was used. Total simulation time varied from 20 to 40 ns moving from low to higher loadings, using a time-step up to 2 fs for the integration of the equations of motion. These parameters ensured statistically satisfactory energy conservation for all runs. The number of cells along the z axis was kept equal to 1, since there is no connectivity in this direction. This is seen in Figure 1 and is also evident from the zero value of the component of the self-diffusivity tensor, Ds,zz, in this direction (see graph in next section). Therefore, since matter is not transferred between the two pore systems, replication of cells along z direction would not contribute to sampling system-sizedependent fluctuations and hence to their effect on system dynamics. Monte Carlo simulations in the grand canonical ensemble (GCMC),10 under constant component fugacities f1, f2, total volume VSC + VLC, and temperature T, were employed for the

TABLE 1: Force Field Parameters for ITQ-1,a CO2,b CH4c and Their Mutuald Interactions pair

ε/kB (K)

σ (Å)

CH4-Oz CH4-CH4 CH4-Oz CH4-CH4 Oz-Oz Si-Si C-C O-O C-O C-Oz O-Oz C-CH4 O-CH4

115.00 158.50 133.30 147.95 89.60 none 28.129 80.507 47.588 50.20 84.93 66.77 112.96

3.47 3.72 3.214 3.73 2.806 none 2.757 3.033 2.892 2.781 2.919 3.2385 3.3765

a Subscript z denotes zeolite oxygen. b C and O denote the CO2 atoms. c Two sets of parameters (see text). d Values in last four rows obtained via the Lorentz-Berthelot combining rules.

TABLE 2: Partial Charges Attributed to Various Atom Typesa

a

atom

q (e)

Oz Si C O

-1 +2 +0.6512 -0.3256

Symbols for atoms as in Table 1.

preparation of the one component (f2 ) 0 or f1 ) 0) as well as the mixed sorbate configurations. Final configurations from the Monte Carlo, after thorough equilibration at fixed f1 and f2, were used as input to the MD runs. At equilibrium, each of the two pore systems, SC and LC, accommodated a certain number of molecules of each species, with both species in both pore systems being at the same fugacity (or chemical potential) with the bulk sorbate phase to which the zeolite was exposed (Figure 2). The integration of the equations of motion was performed via fourth order predictor-corrector algorithm of first order equations of motion, for both translation and rotation of the CO2 molecules altering the original method due to Cheung and Powles.18,19 At the beginning in each MD simulation, translational and angular velocities were assigned random values according to a Gaussian distribution algorithm. These were subsequently readjusted during the initial steps of preliminary microcanonical (NVE) runs, in order to ensure absence of net drift, and also scaled so that the desired kinetic energy (temperature) was attained. However, the lack of molecular population interchange (very low kinetic energy exchange) between the two isolated pore networks in the MWW structure demanded individual thermostatting for each pore network by means of simulations under the NVT ensemble. Therefore, the translational temperature of the molecules in LC and SC was thermostatted independently through the method of a Gaussian thermostat due to Evans and Morriss;20 the rotational temperature was then allowed to reach the desired value indirectly by the kinetic energy equipartition. In molecular dynamics, exceptionally for the fluid-fluid longrange electrostatic interactions, we used conventional Coulomb summation by means of a spherical cutoff in order to avoid an extremely high number of wavevectors in the reciprocal part of Ewald summation.21 This enabled usage of a big simulation box (large number of unit cells).

13764


Sant et al.


Figure 4. Mean square displacement plots of CH4 inside the LC pore system of ITQ-1 along the three orthogonal principal directions, x, y, z; LC loading: 6.6 molecules per unit cell; T ) 300 K.

Figure 2. Indicative equilibrium distribution between the two independent pore systems, LC and SC, of methane (top) and carbon dioxide (bottom) in ITQ-1, at 300 K; circles denote the total loading, triangles the LC loading and squares the SC loading.

Figure 3. Molecular fraction of CO2 in the sorbed phase, x, in ITQ-1 as a function of the molecular fraction in the bulk phase, y, at equilibrium, for various total pressure values.

3. Results and Discussion The sorption results from the GCMC simulations are presented in Figure 2. For methane the potential parameters of ref 15 (Table 1), were used. In the same figure the individual pore system isotherms are shown too. The singlet density distribution for finding a molecule within either of the two networks, SC and LC, F1R(r), was computed by GCMC. This is the average number of R-type sorbate molecules per volume of a given species at a certain position r, averaged over all occupancies and all molecular configurations inside the unit cell, subject to the imposed fugacities and temperature. In Figure 1, low-value isodensity surfaces for methane were selected in order to depict the maximum accessible (void) volume formed by the large cavities and by the sinusoidal channels in the MWW framework type. In Figure 3 the selectivity isotherms for the mixture CH4-CO2, as computed by GCMC, are presented under three total pressures of the bulk phase. As pointed out in ref 10, the zeolite is selective toward CO2, especially at high pressures. An indicative mean squared displacement plot for pure methane in ITQ-1 is depicted in Figure 4. Clearly, a vanishing

component, Ds,zz, of the diffusivity tensor in the z direction is predicted, in agreement with the GCMC contours of Figure 1, which show no communication of the two channel systems in this direction; to an excellent approximation, Ds,xx ) Ds,yy, as expected from the hexagonal symmetry of the ITQ-1 unit cell. Hereinafter the averaged coefficients were taken for do ) 2. Molecular dynamics results for the one component sorbed phase of carbon dioxide and methane are presented in Figure 5. The overall self-diffusivity for CO2 decreases steeply as its concentration in the zeolite increases, controlled mainly by intermolecular collisions. On the contrary, CH4 exhibits a markedly different trend. To analyze this phenomenon better, we studied the dynamics of the two sorbates separately in the LC and SC pore spaces. The results, displayed in the same figure, show that, whereas in the sinusoidal channel of MWW both molecules follow almost exactly the same behavior, in the large cavity the self-diffusivity of methane goes through a maximum before it starts to decrease. Using the potential parameters of ref 16 instead of those of ref 15 (compare Table 1), very similar diffusivities were obtained and the trends pointed out above persisted (Figure 5, open symbols). All the error bars in Figures 5 and 6 are smaller than the symbols’ size. Also, lines are used to guide the eye. The mixture dynamics is presented in Figure 6 for a variety of CH4/CO2 fugacity ratios. Starting from a methane-rich sorbed phase and going up to a carbon dioxide-rich phase we found that the self-diffusivity of CH4 in the presence of CO2 decreases fast with pressure at high pressures, whereas almost no change is observed in the lower pressure region. This can be explained on the basis of the significant increase in selectivity for CO2 for a given fugacity ratio (constant yCO2) as the total pressure rises (see Figure 3). An abundance of sorbed CO2 in the pores slows down CH4 considerably. For low yCO2 values, the variation of selectivity with pressure becomes much more pronounced (cf. Figure 3). The opposite effect appears in the self-diffusivity of CO2 in the presence of CH4. Here (Figure 6, bottom) the self-diffusivity of CO2 increases considerably with CH4 content in the regime of intermediate total pressures but is quite insensitive to composition at high total pressures. The peculiar evolution of the self-diffusivity of CH4 in ITQ-1 with occupancy, whether it is a pure component or part of a mixture with CO2 (Figures 5 and 6), is related to its molecular motion within the LC system of this zeolite (compare maximum in CH4 self-diffusivity in Figure 5, middle). For this reason we analyzed the singlet density distribution of the two sorbates in pure form within this pore system, consisting of supercages connected by narrow necks. Contour plots of this distribution, in xz and xy projections, are shown in the top parts of Figures


Molecular Dynamics of CO2, CH4 and Their Mixtures

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13765

Figure 6. Self-diffusivity vs f ) fCH4 + fCO2, for methane (top), and carbon dioxide (bottom), for various fugacity ratios of the two sorbates (indicated percentages) coexisting in ITQ-1 at 300 K.

Figure 5. Self-diffusivities of pure CH4 (circles) and pure CO2 (triangles) in ITQ-1 at 300 K (top) as functions of loading; contributions from the LC pore system (middle) and from the SC system (bottom) have been separated; open symbols denote results for CH4 by making use of the set of the force field parameters of ref 16 (see Table 1). Symbols from left to right correspond to the following fugacities: 0.128, 0.348, 0.947, 2.575, 6.999, 19.025, and 51.715 atm.

7 and 8. In addition, probability density functions p(z) along the z direction inside the supercage were accumulated. These are shown in the bottom parts of Figures 7 and 8 for three different loadings, corresponding to a very low, an intermediate and a high fugacity value (the maximum in Ds appears around the value of 1 atm). Note that the probability density, p(z), can be converted to a free energy profile, A(z), through the relation

A(z) ) -kBT ln[p(z)]

(9)

The z-dependent histogram, p(z), displayed at the bottom of Figures 7 and 8 is based on an analysis of the contents of the supercage; molecules from the connecting necks between supercages have been excluded from this analysis. To isolate the supercage region, a cylindrical surface of radius R was drawn around the supercage axis. The R value used for each sorbate

Figure 7. Top: xz (left) and xy (right) isodensity surfaces of pure CH4 inside the LC pore system (LC loading: 6.6 molecules/u.c.) under P ) 2.757 atm; dark areas denote higher density. The circled area bounds the supercage in our calculations (R ) 4.8 Å). Bottom: One dimensional probability density distribution along the z-axis inside the supercage of the LC network, and its evolution with increasing pressure.

is indicated in the legends to Figures 7 and 8, and also depicted pictorially in the xy projections of the singlet density distribution plots of the same figures. A larger value of R had to be used for CO2, given its different distribution within the supercage. As shown in Figure 7, at low occupancies methane prefers to reside within culs-de-sac at the top and bottom of the supercage, rather in the middle region from which methane can diffuse along any one of the six narrow channels connecting the supercages. As pressure (occupancy) rises, a smaller fraction of molecules can be accommodated in the energetically favorable culs-de-sac and methane concentration in the middle of the supercage rises spectacularly, facilitating diffusion into neighboring supercages. It is this redistribution, associated with the peculiar geometry and energetics of the LC system, which gives rise to the increasing trend in self-diffusivity displayed in the middle panel of Figure 5 for the large cavity system. At

13766


Sant et al.


Figure 10. Transport diffusivity, D, vs fugacity for pure CH4 (circles) and pure CO2 (triangles), in ITQ-1 at 300 K. Figure 8. Same as Figure 7 for CO2. In the top plot, isodensity surfaces are given for a LC loading of 10.8 molecules/u.c. under P ) 2.757 atm. The radius bounding the supercage is R ) 5.8 Å.

insight of the loading dependence of self-diffusivity, based on a diffusion model previously developed,22 is currently in progress. The collective diffusivities in the single and two component sorbate systems derived from this work are depicted in Figure 9. A qualitative explanation of their loading dependence can be sought by means of the quasichemical mean field theory14,23 through the parameter w, which is directly related to the strength of interactions between neighboring sorbate molecules.12 For pure methane, as pressure increases, the number of methane molecules residing around the windows in the supercages of the LC system also increases (Figure 7); this makes the local intermolecular interactions more “repulsive” in the aforementioned area, resulting in a less negative value of w, which causes a maximum in the collective diffusivity12 of methane. As the sorbed phase becomes leaner in methane, the carbon dioxide energetics dominate, implying much more attractive intermolecular interactions and a more negative value of w; in this case the increase of pressure does not alter the molecular distribution around the windows; consequently, carbon dioxide’s collective diffusivity continues to decrease, as quasichemical theory predicts.13 In Figure 10 is shown the influence of the thermodynamic factor (eq 7, R ) β) in converting D0 to the transport diffusivity, D, through eq 6 for pure methane and pure carbon dioxide, sorbed in ITQ-1. 4. Conclusions

Figure 9. Collective diffusivity (see section 2.1) vs f ) fCH4 + fCO2, for methane (top), and carbon dioxide (bottom), for various fugacity ratios of the two sorbates coexisting in ITQ-1 at 300 K.

even larger methane pressures, crowding around the window region causes the self-diffusivity to drop again. Contrary to methane, CO2 (Figure 8) does not undergo a significant redistribution along the z-direction of the LC system as occupancy rises. Clearly, the p(z) histograms of CO2 in Figure 8 become more structured with increasing pressure. The probability of finding a CO2 molecule in the middle of the supercage actually drops, rather than rising, with increasing pressure. These features lead to the monotonically decreasing self-diffusivity curves of Figure 5. Further work in order to gain

We have studied the single and mixed sorbed phase dynamics of CH4 and CO2 in the zeolite ITQ-1, the siliceous analogue of MCM-22, at room temperature with molecular simulations. The zeolite unit cell was reconstructed atomistically from crystallographic data and GCMC simulations were employed in order to reveal the spatial partitioning of sorbed molecules in its two independent pore networks, the large cavity (LC) and sinusoidal channel (SC) systems. Molecular dynamics verified the independence of the two pore systems. The loading dependence of self-diffusivities was studied for the pure sorbates and their mixtures. While the selfdiffusivity of CO2 exhibited a monotonic decrease with increasing loading, that of CH4 displayed a peculiar loading dependence. Closer examination enabled us to trace this dependence to a maximum in methane’s self-diffusivity along the LC channel system at a pressure around 1 atm. This maximum in methane self-diffusivity can be explained in terms of the redistribution, with increasing loading, of methane within the supercages of the LC system along the z-direction, normal to the direction of diffusion. At low occupancies methane mol-


Molecular Dynamics of CO2, CH4 and Their Mixtures ecules prefer residing at the poles (top and bottom) of the supercage rather than in the equator (middle) from where jumps to other supercages are possible through the narrow necks; the poles of a supercage, then, act as reservoirs of relatively immobile molecules. As occupancy rises, the role of these reservoirs becomes less important because the population in the middle region of the supercage is dramatically enriched facilitating diffusion. This condition lasts up to a pressure where intermolecular collisions become dominant. The loading dependence of the collective (Stefan-Maxwell) diffusion coefficients was also studied. The collective diffusivity of pure methane was found to go through a maximum with occupancy at a pressure around 10 atm, which disappears in the presence of CO2. The collective diffusivity of CO2, either as a pure sorbate or in mixture with CH4, exhibits a decreasing trend with occupancy, behaving very similarly to its selfdiffusivity. The loading dependence of the collective diffusivity was rationalized on the basis of the quasichemical mean field theory of Reed and Ehrlich.12,23 According to that theory, as the number of methane molecules residing around the windows in the supercages of the LC system increases, the local intermolecular interactions become more “repulsive” at high occupancies leading to a maximum; whereas, the different evolution of the distribution of CO2 molecules are much more unfavorable compared to the case of methane (cf. Figures 7 and 8), leading to a decrease in collective diffusivity with loading. Further work toward a quantitative explanation based on a diffusion model previously developed22 is currently in progress. Concluding, it should be noted that experimental verification of the modeling results would be useful. For instance, experiments able to probe the mass transport in a microscopic level (e.g., quasi elastic neutron scattering) or macroscopic experimental techniques (e.g., pulsed field nuclear magnetic resonance) would possibly reveal the contribution of the two morphologically different pore networks to the transport mechanism.

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13767 Acknowledgment. Support by the European Union via the FP6-Marie Curie Research Training Network “INDENS” (MRTN-CT-2004-005503) is gratefully acknowledged. References and Notes (1) Davis, H. T. Statistical Mechanics of Phases, Interfaces and thin Films; Wiley-VCH: New York, 1996. (2) Vanderlick, T. K.; Davis, H. T. J. Chem. Phys. 1987, 87, 1791. (3) Bitsanis, I. A.; Magda, J. J.; Tirrell, M.; Davis, H. T. J. Chem. Phys. 1987, 87, 1733. (4) Snurr, R. Q.; Ka¨rger, J. J. Phys. Chem. B 1997, 101, 6469. (5) Sanborn, M. J.; Snurr, R. Q. AIChE J. 2001, 47, 2032. (6) Dubbeldam, D.; Snurr, R. Q. Mol. Simul. 2007, 33, 305. (7) Bell, A. T.; Maginn, E. J.; Theodorou, D. N. In Handbook of Heterogeneous Catalysis, 1st ed.; Ertl, G., Kno¨zinger, H., Weitkamp, J., Eds.; VCH: Weinheim, 1997; Vol. 3, p 1165. (8) Papadopoulos, G. K.; Theodorou, D. N. Computer simulation of sorption and transport in zeolites. In Handbook of Heterogeneous Catalysis, 2nd ed.; Ertl, G., Kno¨zinger, H., Schu¨th, F., Weitkamp, J., Eds.; WileyVCH: Weinheim, 2008; Vol. 3, p 1662. (9) Camblor, M. A.; Corma, A.; Dıáz-Cabanãs, M. J.; Baerlocher, C. J. Phys. Chem. B 1998, 102, 44. (10) Leyssale, J.-M.; Papadopoulos, G. K.; Theodorou, D. N. J. Phys. Chem. B 2006, 110, 22742. (11) Onsager, L. Phys. ReV. 1931, 37, 405. (12) Papadopoulos, G. K.; Jobic, H.; Theodorou, D. N. J. Phys. Chem. B 2004, 108, 12748. (13) Krishna, R.; Wesselingh, J. A. Chem. Eng. Sci. 1997, 52, 861. (14) Papadopoulos, G. K.; Theodorou, D. N. Mol. Simul. 2009, 35, 79. (15) Dubbeldam, D.; Calero, S.; Vlugt, T. J. H.; Krishna, R.; Maesen, T. L. M.; Smit, B. J. Phys. Chem. B 2004, 108, 12301. (16) Goodbody, S. J.; Watanabe, K.; MacGowan, D.; Walton, J. P. R. B.; Quirke, N. J. Chem. Soc., Faraday Trans. 1991, 87, 1951. (17) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (18) Cheung, P. S. Y.; Powles, J. G. Mol. Phys. 1975, 30, 921. (19) Allen, M. P., Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (20) Evans, D. J.; Morriss, G. P. Comput. Phys. Rep. 1984, 1, 297. (21) Pantatosaki, E.; Papadopoulos, G. K. J. Chem. Phys. 2007, 127, 164723. (22) Sant, M.; Papadopoulos, G. K.; Theodorou, D. N. J. Chem. Phys. 2008, 128, 024504-1. (23) Reed, D. A.; Ehrlich, G. Surf. Sci. 1981, 102, 588.

JP902829J

Molecular Dynamics of Carbon Dioxide, Methane and Their Mixtures

Recommend Documents