Prediction of Permeation Properties of CO2 and N2 through Silicalite

Jan 5, 2001 - Location of CO 2 on silicalite-1 zeolite using a single-crystal X-ray method. Shinjiro Fujiyama , Natsumi Kamiya , Koji Nishi , Yoshinob...
2 downloads 7 Views 138KB Size
Download PDF
J. Phys. Chem. B 2001, 105, 777-788

777

Prediction of Permeation Properties of CO2 and N2 through Silicalite via Molecular Simulations Konstantinos Makrodimitris,†,‡ George K. Papadopoulos,† and Doros N. Theodorou*,†,§ Institute of Physical Chemistry, Demokritos National Research Center, GR 15310 Ag.ParaskeVi, Athens, Greece, Department of Chemical Engineering, National Technical UniVersity of Athens, GR 15773 Zografos, Athens, Greece, and Department of Chemical Engineering, UniVersity of Patras, GR 26500 Patras, Greece ReceiVed: August 8, 2000; In Final Form: October 27, 2000

The sorption isotherms and self-diffusivities of CO2 and N2 in silicalite have been calculated via grand canonical Monte Carlo and equilibrium molecular dynamics simulations over a wide range of occupancies, using various force fields proposed in the literature. Predictions for the sorption thermodynamics are in very favorable agreement with the experiment, especially when detailed point-charge models are used to represent the interaction of the quadrupole moments of the sorbate molecules with the lattice field and with each other. They indicate that the zeolite cannot be in its para (P212121) form under the conditions of the measurements. Permeabilities corresponding to a perfectly crystalline membrane have been estimated for CO2 and N2, as well as for methane, examined in past simulation work, from the predicted sorption isotherms and lowoccupancy self-diffusivities by invoking the Darken equation. The ratios of pure component permeabilities obtained in this way agree very well with actual macroscopic values obtained from carrying out permeation measurements for the different pure sorbates in the same silicalite membrane. Absolute magnitudes of the permeabilities, however, exceed by more than 2 orders of magnitude the reported macroscopic values, which themselves vary widely among different experimental investigations. The large, morphology-dependent nonuniformity in membrane thickness of actual supported silicalite membranes is proposed as a plausible reason for this disparity.

1. Introduction There are strong incentives for simulating, at the molecular level, the thermodynamics and dynamics of sorbed phases of pure gases or mixtures in zeolites, because oftentimes complicated experimental measurements can be guided or even replaced by a rigorously designed computer experiment. An additional advantage of simulations is that one can save time and money by measuring properties on a computer under conditions which are difficult or impossible to attain in the laboratory. Molecular modeling, starting from the atomic scale, is an efficient tool for understanding quantitatively and qualitatively structure-property relations, elucidating the mechanisms of microscopic phenomena, and even retrieving quantitative information about zeolite-sorbate systems. Atomistic simulations not only complement experiments but also constitute an appropriate tool for checking the assumptions of approximate theories. The work presented in this paper aims at the molecular simulation of the sorption thermodynamics and steady-state permeation kinetics of nitrogen and carbon dioxide in silicalite. Isotherms have been predicted by grand canonical Monte Carlo (GCMC) simulations. Self-diffusivities have been extracted from equilibrium molecular dynamics (EMD) simulations. Several pairwise additive effective potentials, proposed in the literature, have been invoked to represent the sorbate molecules and * To whom correspondence should be addressed. Telephone: +3061 997 398. Fax: +3061 993 255. E-mail: [email protected]. † Demokritos National Research Center. ‡ National Technical University of Athens. § University of Patras.

silicalite structure. Three different crystal structures of silicalite have been considered in the calculations. Simulated isotherms are compared directly with gravimetric measurements of sorption carried out in both the low-occupancy (Henry’s law) and high-occupancy regions. Simulation results for the permeability are discussed in light of microscopic [EMD,1-3 pulsed field gradient NMR (PFG-NMR),4,5 and quasielastic neutron scattering (QENS)]6,7 and macroscopic methods [supported membrane permeation,8-17 gravimetric,18 chromatographic,19 Wicke-Kallenbach (WK), single-crystal method (SCM),20 zero-length column (ZLC),21 transient analysis of products (TAP),22,23 and frequency response (FR)]24 that have been used in the literature to measure diffusion in zeolites and permeation through zeolite membranes. Eventually, we have attempted to estimate the permeability of silicalite membranes under steady-state conditions to pure CO2, N2, and CH4 on the basis of the atomistic-simulation results we are reporting herein for CO2 and N2 and analogous results obtained in the past for CH4. Pure component permeability ratios, which are relevant to the selectivity exhibited by silicalite membranes in separations, are compared with those of the experiment. Absolute values of predicted permeabilities are discussed in light of existing microscopic and macroscopic measurements of self-diffusivities and membrane permeation experiments in order to elucidate the role of the structural parameters of the sorbent on observed transport rates. In general, the transport of sorbate through a zeolite layer is described17,25 in terms of a five-step model: 1. Adsorption from the gas phase to the external surface. 2. Mass transfer from the external surface into the zeolite pores.

10.1021/jp002866x CCC: $20.00 © 2001 American Chemical Society Published on Web 01/05/2001

778 J. Phys. Chem. B, Vol. 105, No. 4, 2001

Makrodimitris et al.

TABLE 1: Crystallographic Data on the Three Forms of Silicalite Considered in This Work properties

ortho (ZSM-5)30

mono (H-ZSM-5)31

para (H-ZSM-5)32

a, Å b, Å c, Å R, deg β, deg γ, deg space group asymm unit/unit cell max pore size, straight channel, Å min pore size, straight channel, Å max pore size, sinusoidal channel, Å min pore size, sinusoidal channel, Å

20.07 19.92 13.42 90 90 90 Pnma 1/8 5.6 (circular) 5.4 (circular) 5.5 (circular) 5.1 (circular)

20.107 19.879 13.369 90.67 90 90 P21/n11 1/4 5.78, 5.83 (ellipt) 5.18, 5.27 (ellipt) 5.89, 5.35 (ellipt) 5.78, 5.01 (ellipt)

20.121 19.82 13.438 90 90 90 P212121 1/4 6.06, 6.18 (ellipt) 5.07, 4.80 (ellipt) 6.37, 6.15 (ellipt) 4.76, 4.58 (ellipt)

3. Intracrystalline diffusion. 4. Mass transport out of the zeolite pores to the external surface. 5. Desorption from the external surface to the gas phase. Intracrystalline diffusion is the only step that a microscopic method would usually measure. It entails molecules moving within the pores of bulk zeolite. Macroscopic methods, on the other hand, may be subject to resistances associated with all of the five steps mentioned above. To avoid the diffusion limitations on the outside of the sorbent particles, a macroscopic experimental technique termed Multitrack,22,23 inspired by past TAP measurements, is performed under conditions such that extracrystalline diffusion is in the Knudsen regime (low pressure, high temperature, and low occupancy). On the other hand, there have been efforts to simulate with microscopic computations (molecular dynamics) the outcome of a macroscopic permeation experiment. As pointed out by Sun et al.,20 a concentration gradient exists in the macroscopic-supported membrane permeation and in FR techniques, whereas the PFG-NMR technique measures the mobility of molecules in the pores under equilibrium conditions, in the absence of any gradients. Simulations in the presence of concentration gradients can be performed with nonequilibrium molecular dynamics.1 Thus, Pohl et al.26 applied a dual control volume grand canonical molecular dynamics technique to simulate pressure-driven gas permeation through a porous membrane. This computer experiment resembles actual experiments carried out in testing molecular-sieving membranes. Predicted values of permeabilities for CH4 in silicalite at 300 K were higher than the experimental values of Yan et al.8 by a factor of 40 and lower than those estimated by equilibrium MD1,2 and other equilibrium microscopic methods by a factor of about 5.4,6,27 A simulation strategy which would provide direct molecularlevel insight into permeation through zeolite membranes would be to analyze all transport processes occurring in different regions of the membrane (e.g., adsorption, surface diffusion, intracrystalline diffusion, and desorption). Unfortunately, such a study has to overcome the difficulty of describing the membrane surface and the interfacial regions between crystallites in the membrane; surface and interfacial structures are currently not known to sufficient detail. In atomistic-simulation work, silicalite is usually modeled as a rigid crystal; in this work, too, lattice vibrations have been neglected. In addition, sorbate molecules have been represented as rigid. Vibrational motion could be introduced for a more realistic description.28 On the basis of past studies, however, we expect that the introduction of vibrations into the crystal of silicalite and into the sorbate molecules would leave the sorption thermodynamics practically unaffected and would also not make a big difference in the estimated diffusivities for small molecules

such as N2 and CO2, which do not experience a close fit in the zeolite channels. 2. Model Description Construction of our atomistic model requires the geometry of the sorbate molecules and silicalite lattice and the potentials of interactions between atoms or functional groups. 2.1. Geometry. 2.1.1. Sorbent. Silicalite is a crystalline solid with orthorhombic or monoclinic symmetry depending on the temperature and the amount and the kind of adsorbate present. In highly pure, crystalline silicalite (Si/Al > l000), a reversible monoclinic (P21/n11)-orthorhombic (Pnma) phase transition occurs between 350 and 363 K.29 Silicalite adopts the three crystalline forms shown in Table 1(i.e., the monoclinic form with P21/n11 symmetry (mono), the orthorhombic form with Pnma symmetry (ortho), and the orthorhombic form with P212121 symmetry (para)). There are three intracrystalline environments: interiors of straight channels, interiors of zigzag channels, and channel intersections. The interactions of sorbate molecules with the framework in each of these environments influence their diffusion and adsorption behavior. The calculations reported here have shown that N2 and CO2 reside preferentially between the channel intersections. The distinctive difference between silicalite on the one hand and ZSM-5 or more specifically H-ZSM-5 on the other is the Si/Al ratio in the framework. All three zeolites belong to the structure-type code MFI. ZSM-5 has a Si/Al ratio ranging from 10 to infinity; when this ratio exceeds 1000, the zeolite is known as silicalite. Typical values for Si/Al ratios in an H-ZSM-5 material range from 30 to 300.27,29 The coordinates of Si and O in the unit cell are usually assumed to be approximately the same for the three zeolites. In this work, the coordinates of the atoms in the asymmetric unit of silicalite, required for the sorbent potential map calculation, have been taken from ZSM-5 (ortho)30 and H-ZSM-5 (mono and para);31,32 Al was assumed to be completely absent in the calculations. For the calculation of the Coulombic part of the potential energy (see below), partial charges of -1 and +2, in units of the elementary electronic charge, were placed on oxygen and silicon atoms of the zeolite framework, respectively. These charges are based on semiempirical and local density-functional theory electronic embedded-cluster calculations33 and have given satisfactory results in past computations with aromatic molecules in silicalite.34 The presence of the code “CB” in the abbreviations we use here for denoting the silicalite model indicates that partial charges in the framework atoms are taken into account. The presence of the code “Si” indicates inclusion of silicon atoms in the calculation of the non-Coulombic potential energy of

Permeation Properties of CO2 and N2

J. Phys. Chem. B, Vol. 105, No. 4, 2001 779

TABLE 2: Models and Interaction Parameters Used To Represent the Zeolite Framework LJ and BUCK parameters name LJ.JBT36 LJCB.JBTLC33,36

LJCBSi.WSAGM37,38

BUCKCB.JBTLCRSO33,36,39

authors

O

Si

June Bell Theodorou June Bell Theodorou Lonsinger Chakraborty Wantanabe Stapleton Austin Genechten Mortier June Bell Theodorou Lonsinger Chakraborty Reischman Schmitt Olson

/kB ) 89.6 K σ ) 0.2806 nm

partial charges O

Si

0e

0e

-1 e

+2 e

-1 e

+2 e

-1 e

+2 e

/kB ) 89.6 K σ ) 0.2806 nm

/kB ) 101.65 K σ ) 0.2708 nm

/kB ) 18.62 K σ ) 0.0677 nm

/kB ) 89.6 K σ ) 0.2806 nm

A ) 287 521 kJ b ) 3.96 Å-1

interaction with the sorbates; otherwise, only oxygen atoms35,36 are considered for that part of the energy. The symmetry of the silicalite crystal, namely, para (P212121), mono (P21/n11), or ortho (Pnma), is also denoted explicitly. Previous simulation works on silicalite have considered the ortho as its prevalent form, except for a few calculations with para.34 In the present work, we have examined the sorption and diffusion of nitrogen by modeling the ortho, para, and mono space groups as well. The simulation box used in this work, after testing for system size effects, consisted of 8 unit cells for GCMC simulations and 27 cells for EMD simulations. The larger box size in the case of EMD was chosen to ensure the collection of good statistics on the motion and to avoid artificial correlations in motion because of periodic boundary conditions. 2.1.2. Sorbates. The diatomic N2 molecule was modeled as a dumbbell with a rigid interatomic bond. The triatomic linear molecule of CO2 was modeled as two consecutive dumbbells sharing the central C atom, arranged on a straight line. To describe the interactions, a simplified representation was used, cast in terms of Lennard-Jones (LJ) sites on the atoms and partial charges on the molecular axis; electronic clouds were not considered explicitly. Partial charges were distributed around each molecule so as to reproduce experimental quadrupole moments (QN2 ) -4.67 × 10-40 C m2 and QCO2 ) -13.67 × 10-40 C m2). 2.2. Interatomic Potentials. Table 2 presents the four models with the corresponding parameters used for the silicalite-sorbate interactions. Figures 1 and 2 present in detail the sorbate potential parameters. The models shown in the above table and figures were classified according to (a) the way we handle short and longrange forces and (b) the types of the atoms involved in the potential map calculation. Figures 1 and 2 refer to carbon dioxide and nitrogen, respectively. In the abbreviations, the first part declares the type of potential employed, e.g., LJ or Buckingham (BUCK), and the presence of partial charges (CB). The authors’ initials or the model code name, if it exists, is included in the last part of the abbreviation (e.g., the notation 3LJ3CB.EPM2 is derived from the EPM2 model developed by Harris and Yung).40

Figure 1. Two CO2 pairwise additive effective potentials chosen from the literature.

2.2.1. Short-Range Forces. In this category belong the nonCoulombic forces exerted by the atoms of the zeolite on the sorbate molecules and among the sorbate molecules. The potential energy V(r) here is either a 12-6 LJ potential

V(r) ) 4[(σ/r)12 - (σ/r)6]

(1)

or an exp-6 “modified” Buckingham potential (N2 case), in the sense that V(r) is the sum of a London dispersion-type r-6 attractive part taken from LJ, whereas the repulsive part is an exponential Born-Mayer term:

V(r) ) A exp(-br) - 4(σ/r)6

(2)

780 J. Phys. Chem. B, Vol. 105, No. 4, 2001

Makrodimitris et al. by a test charge at the nodes of a fine cubic mesh constructed throughout the asymmetric unit of silicalite. The partial charges attributed to oxygen and silicon atoms of the zeolite framework are shown in Table 2 for the various zeolite models used in this work. The designation “CB” signifies the inclusion of Coulombic interactions in the model. A map of the electrostatic potential was created via the simple and rigorous direct summation of the Coulombic interactions between the test charge and all charges on the atoms of an infinite zeolite lattice. In this calculation, the surrounding unit cells were added progressively in order of decreasing proximity to the central unit cell containing the test charge. In this way, about 10 000 cells were placed around the central cell (2.5-3 million Si and O atoms of the silicalite crystals participated in the calculations); resulting potentials were fully converged and did not change further upon the addition of more unit cells. 2.2.3. Atoms Considered in Representing the Zeolite Framework. In most calculations reported here, Si atoms were neglected in calculating the short-range (non-Coulombic) part of the sorbate-sorbent interaction energy. In these calculations, we have accepted the values used by June et al.36 for the LJ parameters of zeolite oxygens (see Table 2; models LJ.JBT,36 LJCB.JBTLC,33,36 and BUCKCB.JBTLCRSO33,36,39). The optimal parameters for the zeolite oxygen were obtained by June et al.36 by comparing computed and experimental data for the Henry’s constant of methane at 300 K. Alternatively, we adopted the model of Watanabe et al.,37 in which Si atoms are taken explicitly into account in the calculation of the short-range terms. These authors obtained optimal parameters for both O and Si after fitting to the experimental Henry’s constants for argon, nitrogen, and oxygen (model LJCBSi.WSAGM37,38). 3. Simulation Methods

Figure 2. Four N2 pairwise additive effective potentials chosen from the literature.

The cross-interaction LJ potential parameters were calculated through the Lorentz-Berthelot combining rules. For the exponential part in the BUCK relation, we have used the same combining rules; Aij ) (AiAj)1/2 and bij ) (bi + bj)/2. The way the zeolite-sorbate interactions were calculated is described in detail elsewhere.36 2.2.2. Long-Range Forces. This category includes forces of an electrostatic nature. The sorbate-sorbate electrostatic potential energy and, hence, the forces were calculated through the Ewald summation method. We have used the potential energy and the self-energy due to partial charges within a molecule from ref 46 (eqs 5.20 and 5.24, respectively). The parameter κ determines the width of the Gaussian distribution, which screens interactions between neighboring charges in the Ewald method, making them short-ranged. The number of reciprocal space vectors used has to be sufficient for the chosen value of κ and as small as possible to keep central processing unit (CPU) time low. With our choice of parameters, the direct space term converges within the cutoff distance (minimum image convention). The choices for CO2, N2 (3 partial charges), and N2 (4 partial charges) were respectively κ ) 8/Lmax, 8.2/ Lmax, and 9/Lmax, where Lmax is the length of the longest side of our rectangular simulation box. The sorbate-zeolite electrostatic interactions were calculated on the basis of a pretabulation of the electrostatic potential felt

3.1. GCMC. We have used the Adams GCMC algorithm,47,48 where three types of moves are attempted with equal probability: random displacement and reorientation (a compound move), random insertion, and random deletion.49 The Adams algorithm makes use of a quantity B related to the excess chemical potential by the equation

B)

1 ex µ + ln 〈N〉µVT kBT

(3)

B is related to the fugacity f of the bulk sorbate phase according to

βVf ) exp(B)

(4)

To sample the phase space as satisfactorily as possible, we adjusted the maximum allowed displacement of a particle so that the acceptance ratio was about 20%, and not 50% as usually practiced.50 The efficiency of the sampling process is inversely proportional to the total number of attempted moves needed to obtain statistically significant ensemble averages of macroscopic quantities. For the rotation of the linear molecule, we needed to generate a random unit vector starting at the origin, whose end is uniformly distributed on the surface of a sphere centered at the origin; we used Marsaglia’s51 algorithm for this purpose, with a maximum allowed reorientation angle for a better sampling of the phase space. The CPU time was over 3 weeks on a Pentium III (450 MHz) for the lower temperatures and higher occupancies, merely because of the Ewald summations (N2, 2LJ3CB.MSKM model in LJCB.JBTLC zeolite, 1 atm gas-phase pressure, 29 molecules/

Permeation Properties of CO2 and N2

J. Phys. Chem. B, Vol. 105, No. 4, 2001 781

cell mean occupancy, T ) 77 K, and 10 million cycles of which 6 million were used for equilibration before initiating the averaging). For the same molecular model and fugacity but at 298 K, only a few hours were required to obtain a wellequilibrated average loading in the isotherm. 3.2. EMD. The dynamical behavior of the system was studied in the NVE ensemble, in the absence of any thermodynamic forces that could lead to fluxes within the box. The starting configuration for the EMD simulation was taken from configuration files stored in the previous GCMC simulations for the respective temperatures and species. In the EMD, the particles followed classical trajectories, whose time evolution is governed by Newton’s equation of motion. (The neglect of quantum effects is fully justified at the temperatures of interest here.) For continuous potentials, we need a finite difference approximation to solve the differential equations of motion for the particles. The algorithm chosen for the solution of the differential equations of motion belongs to the category of extrapolation methods and, more specifically, to Sto¨rmer leapfrog integrators. Because of its simplicity and exceptional stability, this algorithm affords a long integration time step and low computational cost. Its stability has been attributed to the correspondence between the order of the finite difference equation and its analytical analogue, to their time reversibility (the algorithm preserves the volume in the phase space, i.e., is symplectic), and to the oscillating nature of the fundamental solutions (periodicity).52 Sto¨rmer leapfrog integrators give fixed time averages and a low drift of the total energy, even with a long integration time step, though they conserve instantaneous total energy rather poorly. Thus, it appears that the trajectory constantly deviates from the initial hypersurface but always finds a way back. Mazur52 supports the view that the leapfrog scheme is distinguished among the symplectic operators in terms of both the order of truncation errors it produces and the conservation of the total energy it affords and that the trajectory really manages to sample from a correct hypersurface in the phase space with larger step sizes than commonly recommended. The center of mass velocities v and the positions x for each molecule in Cartesian space are calculated with the leapfrog algorithm53

vn+1/2 ) vn-1/2 + an δt

(5)

and

T)

∑R rR × fR ) eˆ × ∑R dRfR ) eˆ × g

with g being the “turning force”, which can be determined from the nonbonded forces fR on each atom, the position vectors rR of each interaction site, and the (algebraic) distances dR of each atom R from the center of mass of the linear molecule. In a linear molecule, g can be replaced by its component perpendicular to the molecular axis, gp, without affecting the torque. gp is defined as

gp ) g - eˆ (eˆ ‚g)

(6)

where δt is the integration time step and a is the acceleration, computed from the total force on each molecule. For a linear molecule with fixed internal geometry, we would need holonomic constraints introduced into our equations of motion if we were to perform the integration in Cartesian coordinates. However, one of the principal moments of inertia vanishes, and the other two are equal. Thus, we can track rotational motion in generalized coordinates using a simple quaternion algorithm with one of the body-fixed angular velocity components (rotation about the molecular axis) being always kept to zero. The total angular velocity and the torque must always be perpendicular to the molecular axis.54 Fincham has proposed four algorithms as a solution to this problem using a leapfrog scheme.54 We chose the LEN algorithm, which is recommended by the author as superior to the other three (ORT, EXP, and IMP), because LEN has the smallest drift energy, even for long time steps. If eˆ is the unit bond vector fixed along the molecule axis, the torque T on the molecule can be written as

(8)

There is an alternative to the use of the angular velocity, and that is used by the LEN algorithm, which works in terms of the velocity of the axis vector:

u ) deˆ /dt

(9)

In this algorithm, to apply the constraint that the length of the eˆ vector remain unity, we need to calculate an undetermined multiplier λ.54 This is achieved exactly with the following steps in the subroutine or function which implements the LEN algorithm in our Fortran90 or C++ programming code, respectively: (1) Find an auxiliary vector, e´ , which would result from eˆ at time t + δt if one ignored temporarily the unit length constraint:

e´ ) eˆ n + un-1/2 δt + (gp/I) δt2

(10)

where I is the moment of inertia of the molecule around any axis drawn through its center of mass perpendicular to the molecular axis. (2) The Lagrange multiplier λ is associated with the requirement that the length of the eˆ vector remain equal to unity

λ ) -e´ ‚eˆ n/2 + [(e´ ‚eˆ n/2)2 - e´ 2 + 1]1/2

(11)

(3) The new unit bond vector orientation is given by

eˆ n+1 ) e´ + λeˆ n

(12)

(4) The new midstep rate of change of the axis vector is

un+1/2 ) (eˆ n+1 - eˆ n)/δt xn+1 ) xn + vn+1/2 δt

(7)

(13)

(5) Finally, if we want to check for energy conservation (the part of rotational kinetic energy), un is found by

un ) (un-1/2 + un+1/2)/2

(14)

Of course, to be completed, the LEN algorithm needs the implementation of the linear velocities and positions (leapfrog algorithm; eqs 5 and 6) and the calculation of the translational kinetic energy. The LEN algorithm showed remarkable stability for N2 and CO2 linear molecular models, with no appreciable energy drift compared to the ORT54 algorithm, which was also tested. The time steps we used in the EMD simulations ranged between 1 and 5 fs. The LEN algorithm behaved stably even with a time step of 10-15 fs, but in this case, it led to slightly different results for the diffusivities. The energy drift for the two gases and all models was lower than 0.2-0.4% during all EMD runs (Table 3). Because of the Ewald summations of sorbate-sorbate interactions, the EMD runs were very timeconsuming. A promising approach to save CPU time would be

782 J. Phys. Chem. B, Vol. 105, No. 4, 2001

Makrodimitris et al.

Figure 3. High-pressure adsorption isotherms of N2 in silicalite at 305 K.

TABLE 3: Self-Diffusivities Obtained from EMD Simulations for CO2 at 303 K and N2 at 300 K in Silicalite over a Range of Occupancies molecules/cell

time step (fs)

simulation time (ps)

diffusivity (10-8 m2 s-1)

1.075 2.25 3.185

CO2 (3LJ3CB.EPM2) 1.0 30.0 5.0 35.0 5.0 40.0

0.66 0.6 0.55

0.33 1.89 3.89

CO2 (3LJ3CB.EP) 2.0 180.0 2.0 80.0 2.0 45.0

0.92 0.73 0.5

1.37 3.11 6.444

1.0 1.0 1.0

N2 (2LJ.BP) 100.0 180.0 51.0

1.35 1.05 0.61

1.07 2.11 3.22

N2 (2LJ3CB.MSKM) 1.0 55.0 5.0 45.0 5.0 45.0

2.05 1.67 1.31

to handle the long-range electrostatic forces computed via the Ewald technique through a special multiple time step scheme.55 The test runs indicated that reliable self-diffusivities could be obtained from production runs of 30-180 ps (Table 3), so we chose not to implement a more sophisticated multiple time step scheme. In some cases where the model did not incorporate electrostatic forces, we generated production runs of a longer duration to confirm the convergence of our diffusivity results. EMD runs were carried out on a R10000 SGI Origin 200 workstation. Computational requirements depended on the system, model, and loading examined. A 100 ps run with the 3LJ3CB.EPM2 model of CO2 at an occupancy of 2 molecules/ cell required 30 CPU hours, whereas a 100 ps run with the 2LJ3CB.MSKM model of N2 at an occupancy of 1 molecule/ cell required half of a CPU day. 4. Results and Discussion 4.1. Sorption Equilibria for CO2 from GCMC. We have compared our isotherms, obtained from GCMC simulation with models 3LJ3CB.EPM2 and 3LJ3CB.EP in the zeolite model LJCB.JBTLC, against experimental adsorption measurements from Yamazaki et al.56 in silicalite. As we can see in Figures 3 and 4, our results for both models are in excellent agreement

Figure 4. Low-pressure adsorption isotherms of CO2 in silicalite at 303 K.

with those of the experiment. Moreover, the rigid model 3LJ3CB.EPM2 (proposed by Harris and Yung40 also in a flexible version) exhibits better behavior, both in the Henry’s law region and at high gas-phase pressures (saturated silicalite), in comparison with that of the experiment. 4.2. Sorption Equilibria for N2 at 298 K from GCMC. We have predicted room-temperature sorption isotherms for N2 in silicalite in the ortho, para, and mono phases using various effective potentials for N2, with or without taking into account the Si atoms in the calculation of short-range interactions. Our results have been compared with the experimental measurements of Chaffee et al., presented in ref 37, and with the GCMC simulations of Watanabe et al.37 in ortho silicalite, who have included the short-range interactions between silicon atoms and N2. Our results are shown in Figure 5. Figure 5a shows excellent agreement of the experimental isotherms with the simulated isotherms for both the mono and ortho phases of silicalite, as well as a big discrepancy with the simulation values obtained in the para phase. This leads us to the conclusion that silicalite exposed to nitrogen at room temperature is in the ortho or mono phase or in a biphasic state consisting of coexisting ortho and mono phases. The model we use here for nitrogen is that suggested by Murthy et al.42 and also used by Watanabe et al.37 For the zeolite, for all symmetries under consideration, the results shown in Figure 7a are based on the LJCB.JBTLC model.33,36 Figure 5b presents results obtained with different effective potentials for nitrogen in ortho silicalite, using the LJCB.JBTLC and LJ.JBT models33,36 to represent the zeolite. One can discern a systematic deviation of the 2LJ.BP model, which does not include any Coulombic interactions. This model produces systematically lower values relative to those of the experimental isotherm. Thus, neglect of charge separation and the resulting Coulombic interactions for N2 in silicalite seems unjustified. On the other hand, the expensive 2LJ4CB.BAMOM potential (consideration of 4 partial charges/N2 molecule requires large CPU time because of the Ewald summation) exhibits excellent agreement with the experimental isotherm. In Figure 5c is presented the sorption of N2 in the mono phase of silicalite, using a potential parameterization identical to that employed by Watanabe et al.,37 i.e., the 2LJ3CB.MSKM model for the sorbate molecules and the LJCBSi.WSAGM model for the zeolite. In this case, reasonably good agreement is obtained with the experiment and also with the past simulation results of Watanabe et al.37 in the ortho phase, which are also shown.

Permeation Properties of CO2 and N2

J. Phys. Chem. B, Vol. 105, No. 4, 2001 783

Figure 5. N2 adsorption isotherms in silicalite at 298 K in the three phases of the silicalite crystal (a), for various sorbate models in ortho silicalite (b), and in silicalite where Si interactions are included (c).

Figure 6. N2 isotherms at 77 K in two phases of silicalite with (a) and without (b) partial charges.

A comparison of parts a and c of Figure 5 shows that Theodorou and co-workers’33,36 parameterization of short-range sorbatezeolite interactions (LJCB.JBTLC), which neglects the Si atoms in the LJ summations, and Watanabe et al.’s37 parameterization (LJCBSi.WSAGM), which includes the Si, produce very similar results; the Theodorou and co-workers representation gives predictions somewhat closer to those of the experiment. 4.3. Sorption Equilibria for N2 at 77 K from GCMC. Douguet et al.57 simulated the adsorption of nitrogen in silicalite at 77 K with the GCMC method and compared their predictions with the experimental isotherm of Pellenq. At this temperature, the experimental isotherm displays two steps: the first step occurs at a pressure of 0.002 atm and an occupancy of 20.521.8 molecules/unit cell (“fluid to lattice-fluidlike transition”) and the second at pressures above 0.2 atm and occupancies of 22-28 molecules/unit cell, using a high silica H-ZSM-5 (Si/Al > 1000) and large monocrystals. Mu¨ller et al.58 have measured an adsorption isotherm which shows only one transition step from 24 to 31 molecules/cell at P > 0.15 atm. Although the work of Douguet et al.57 involved a detailed potential function including Ewald sums, induction terms, Born-Mayer repulsion terms, and even three-body terms, it reproduced only the first small step not at the correct pressure but at 0.02 atm. In this work, we tried a different approach to this challenging problem of calculating the low-temperature sorption isotherm of nitrogen, which includes consideration of

different crystal symmetries and different effective potentials for N2. The N2 potentials we considered ranged from a simple one without partial charges (2LJ.BP) to a very detailed one (2LJ4CB.BAMOM), which includes a modified Buckingham potential, 4 partial charges/molecule, and Ewald sums but no three-body terms. Figure 6a shows results obtained with the 2LJ3CB.MSKM model in the ortho form of silicalite using the LJCB.JBTLC model33,36 for the zeolite. The model produces a transition step at 0.2 atm from 25.5 to 29.5 molecules/unit cell; this prediction is in good agreement with Mu¨ller et al.’s experimental measurements.58 It should be noted that at 77 K the zeolite adsorbs nitrogen very strongly and a reliable calculation of the sorption isotherm demands a large amount of CPU time. At pressures below 10-4 atm, this model overpredicts the sorbed amount at this temperature. GCMC results with the same model but in the mono symmetry of silicalite have been obtained only in the lowoccupancy region and are also shown in Figure 6a. The results are practically indistinguishable from those obtained in simulations with the ortho form; no significant differences are seen if we change the symmetry of the zeolite. For the para symmetry, similar calculations (not shown) have given sorbed amounts much higher than those of the experiment; this means that we have to exclude this symmetry at 77 K, as we did at 298 K.

784 J. Phys. Chem. B, Vol. 105, No. 4, 2001

Makrodimitris et al.

Figure 8. N2 isotherms at 77 K in the ortho phase of silicalite.

Figure 7. N2 isotherms at 77 K in two phases of silicalite with and without Si interactions (a) and for various models in ortho silicalite (b).

In Figure 6b, we adopt a coarse approach to modeling the adsorption of nitrogen at 77 K by using the effective potential proposed by Baranello and Panagiotopoulos.45 This potential (2LJ.BP) describes N2 as a LJ dumbbell, without any partial charges. Its parameters were determined so as to give a good representation of the low-temperature vapor-liquid coexistence properties of N2. The isotherms obtained with this model, using the LJ.JBT model for the zeolite, are in very good agreement with those of the experiment for both the ortho and mono phases of silicalite, especially in the Henry’s law region. The agreement is worse at high occupancies. There exists a slight step at 0.03 atm in the simulated isotherm, whose shape deviates considerably from the one obtained experimentally. We have also performed low-temperature N2 sorption calculations, using the potential proposed by Watanabe et al. (2LJ3CB.MSKM for the sorbate and LJCBSi.WSAGM for the zeolite);37 as already mentioned, this includes contributions from silicon atoms in the short-range zeolite-sorbate energy. However, for both the ortho and mono phases, no significant differences are observed with this model relative to the zeolite representation of Theodorou and co-workers,33,36 which takes into account only oxygen in the short-range zeolite-sorbate interactions. This is seen characteristically in Figure 7a. We also tried the more detailed and expensive 2LJ4CB.BAMOM model of N2 with 4 partial charges/molecule within a LJCB.JBTLC zeolite in the ortho phase, in the Henry’s law region, and found no significant differences from the simpler 2LJ3CB.MSKM model (Figure 7b).

Finally, a Buckingham potential with 4 partial charges, suggested by Berns and van der Avoird44 was used. In this case, silicalite oxygens were represented in terms of the potential parameters of June et al.36 for the attractive part and a BornMayer term suggested by Reischman et al.39 for the repulsive part (BUCKCB.JBTLCRSO representation in Table 2). As we can see in Figure 7b, this expensive potential gives excellent agreement with the experiment in the Henry’s law region. Figure 8 shows experimental results found elsewhere57,58 together with our simulated 2LJ3CB.MSKM model in the ortho phase of silicalite represented by the LJCB.JBTLC model and the simulated isotherm of Douguet et al.57 The maximum occupancies obtained are 23 molecules/unit cell in the model of Douguet et al.,57 28-31 molecules/unit cell in the experiment, and 28-29.5 molecules/unit cell in our model. 4.4. Diffusivities from EMD. Table 3 presents the selfdiffusivities for N2 and CO2 (two models used for each in conjunction with the LJCB.JBTLC or LJ.JBT representation of the zeolite) at room temperature for a variety of finite, relatively low occupancies. Self-diffusivities at infinite dilution have been obtained by extrapolation of these values to zero occupancies and are given in Table 4. Table 4 also includes the self-diffusion coefficients at infinite dilution (after an approximate extrapolation) for methane from the work of Maginn et al.,1 Goodbody et al.,68 and Schuring et al.2 and for n-butane and n-hexane from the work of June et al.3 The intracrystalline diffusivities for methane are in excellent agreement with those obtained from the microscopic method of PFG-NMR.4,27 The value from PFG-NMR for methane (10-8 m2 s-1) shown in the table corresponds to an occupancy of 4 molecules/unit cell; in the limit of infinite dilution, it would be approximately 1.6 × 10-8 m2 s-1, very close to the value obtained by June et al.,59 with whose calculations PFG-NMR data agree very well at finite occupancies. The latter value from PFG-NMR is the intracrystalline diffusivity we have used in order to produce an estimated intracrystalline permeability and then to compare it to the measured permeabilities of the supported membranes8 in Table 5. 4.5. Permeabilities from GCMC and EMD. The methodology and assumptions we have followed in this estimation of permeability from the corresponding diffusivities and isotherms computed through our simulations are described in the Appendix. The basic hypotheses being made are the validity of

Permeation Properties of CO2 and N2

J. Phys. Chem. B, Vol. 105, No. 4, 2001 785

TABLE 4: Predicted Self-Diffusivities Ds at Infinite Dilution in Silicalite Compared with Experimental Diffusivitiesa sorbate

Ds

technique

T (K)

N2 N2 CO2 CO2 CO2 CO2 CH4 CH4 CH4 CH4 CH4 CH4 n-butane n-butane n-butane n-butane n-butane n-butane hexane hexane hexane hexane

250 160 70 95 28 24 ≈140 100 10-20 ≈1 41 9 32 1 0.037 10 10 10 22 7 7.5 × 10-6 6.9

EMD (model 2LJ3CB.MSKM) EMD (model 2LJ.BP) EMD (model 3LJ3CB.EPM2) EMD (model 3LJ3CB.EP) FR66 PFG NMR67 EMD1,2,68 PFG-NMR4,5 (4 mol/cell), QENS6,7 FR24 supported membrane methods9 SCM (steady state)69 SCM (transient state)20 EMD3 FR24 supported membrane methods9 QENS6,7 TAP22,23 PFG-NMR4,5 EMD3 gravimetric18 chromatography19 QENS6,7

298 298 303 303 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300

a

All Ds values are in 10-10 m2 s-1.

TABLE 5: Pure Gas Permeabilities (Pe) (10-8 mol/(m s atm)) as Measured in Membrane M3 (ZSM-5) and as Predicted for EMD Simulations in Silicalite sorbate

Pe

technique

sorbent

T (K)

N2 N2 N2 CO2 CO2 CO2 CH4 CH4 CH4 CH4 n-C4H10 n-C4H10

3.2 382 572 8.9 2002 3010 3.85 250 925 1050 0.75 29.5

supported membrane8 EMD (model 2LJ3CB.MSKM) EMD (model 2LJ.BP) supported membrane8 EMD (model 3LJ3CB.EPM2) EMD (model 3LJ3CB.EP) supported membrane8 supported membrane17 EMD1,2 PFG-NMR4, QENS6 supported membrane8 supported membrane17

ZSM-5 (M3) silicalite silicalite ZSM-5 (M3) silicalite silicalite ZSM-5 (M3) silicalite-1 silicalite silicalite ZSM-5 (M3) silicalite-1

303 298 298 303 303 303 303 300 300 300 303 300

the Darken equation, as suggested by many authors,1,27 and the assumption that the corrected diffusivity is independent of concentration.1 Because intracrystalline diffusivities of CO2 and N2 in silicalite from past equilibrium microscopic simulation or experimental methods (EMD, PFG-NMR, and QENS) are unavailable, our results can only be compared to the macroscopic permeabilities measured by Yan et al.8 for supported (porous R-Al2O3) ZSM-5 membranes. For the purpose of this comparison, it is more reasonable to use permeability ratios for different sorbates than absolute permeabilities, so that the influence of the morphological factors of the membrane on the adsorbate mobility can be weakened. Predicted single-component permeabilities for N2, CO2, and CH4 are given in Table 5, and various ratios formed from these permeabilities are shown in Table 6. As seen in the latter table, the calculated permeability ratios agree well with corresponding experimental values; they follow the order CO2 > CH4 > N2, as has been measured by van de Broeke et al.11 and Kusakabe et al.12,13 Also, they are reasonably close to the permselectivity reported by Kusakabe et al.12 for the binary mixture CO2/N2 (≈5-7). It should be noted that permeability ratios measured on a given membrane practically coincide with the corresponding permeance ratios, i.e., the effect of the membrane thickness is largely factored out in taking the ratio.

Although the predicted permeability ratios for the three pure gases through silicalite under low-occupancy conditions are in reasonable agreement with those of the experiment, the absolute values of the predicted permeability are higher by a factor of 120-340 than the macroscopically measured experimental values (Table 5). These discrepancies in absolute values can be attributed to the following reasons: (1) In the present simulation analysis, only the intracrystalline diffusion resistance has been taken into account in calculating the permeabilities. In a macroscopic-supported membrane experiment or in an adsorption/desorption measurement, external heat- and mass-transfer limitations can limit transport rates. Vroon et al.70 have studied the effect of silicalite particle size for CH4 and found that an increase of the interparticle interfaces seems to be beneficial for the flux. Experimental studies have shown a large resistance to mass transfer at the surface of microporous hosts.71 Hargreaves et al.72 simulated an uptake experiment of ethane in MFI zeolites and found that the contribution of surface diffusion to the rate of adsorbate uptake is significant but not dominant; on the other hand, benzene, which possesses a larger kinetic diameter, exhibits a slow uptake, because sorption is governed mainly by intracrystalline diffusion.73 Heat transfer is still a possibility and may reduce the overall permeation rate very significantly.27 Thus, transport rates from EMD simulations, which are sensitive only to the resistance presented by intracrystalline diffusion in a perfect zeolite crystal under equilibrium conditions, tend to be considerably higher than macroscopically measured effective diffusivities but quite close to diffusivities measured by the “microscopic” PFG-NMR and QENS methods. This is seen very characteristically in the results presented for the three hydrocarbons (methane, n-butane, and n-hexane) in Table 4. (2) NMR techniques and an EMD simulation measure the mobility of the molecules inside pores under equilibrium conditions. On the contrary, permeation experiments involve the presence of concentration gradients. Most of the transient methods introduce a disturbance to the system, so the fastest diffusivity that can be measured is restricted by the response time of measurement devices, which is limited below to 10-6-10-7 cm2 s-1.69 Thus, transient experiments are subject to the detrimental role of heat and mass transfer. More recent techniques, such as FR, ZLC, and membrane methods, can reduce heat effects. On the other hand, microscopic methods are suitable for fast-diffusing processes, and the lowest diffusivity that the methods can measure is 10-5-10-7 cm2 s-1, which is the upper limit for macroscopic techniques. However, the steady-state SCM covers a wide range of micropore diffusivities and eliminates heat transfer.69 As we can see in Table 4, the diffusivities from this method agree with those of microscopic methods. (3) In the literature, there is a tendency for reported diffusivities obtained from macroscopic methods to increase with time. Using a ZLC macroscopic experiment (a transient uptake flow method which eliminates considerably the resistances to external heat and mass transfer),62 for example, Loos et al.21 measured larger diffusivities (up to an order of magnitude) for benzene and ethylbenzene in silicalite than for those of previous macroscopic experiments.63-65 Also, Sun et al.,20 using an improved experimental apparatus, procedure, and method of data analysis, obtained diffusivities in silicalite crystals that are only 1 order of magnitude lower than PFG-NMR and EMD values for alkanes, are 1 order of magnitude higher than previous membrane measurements, and agree well with FR methods.24

786 J. Phys. Chem. B, Vol. 105, No. 4, 2001

Makrodimitris et al.

TABLE 6: Predicted and Experimental Ratios in Silicalite sorbates

Pe ratios

technique

sorbent

T (K)

CH4/N2 CO2/CH4 CO2/N2 CH4/N2 CO2/CH4 CO2/N2 CH4/N2 CO2/CH4 CO2/N2 CH4/N2 CO2/CH4 CO2/N2

1.2 2.3 2.8 1.6 2.2 3.4 2.5 2.2 5.3 1.6 3 4.6

supported membrane permeation8 supported membrane permeation8 supported membrane permeation8 EMD (sphere/2LJ3CB.MSKM) EMD (3LJ3CB.EPM2/sphere) EMD (3LJ3CB.EPM2/2LJ3CB.MSKM) EMD (sphere/2LJ.BP) EMD (3LJ3CB.EPM2/sphere) EMD (3LJ3CB.EPM2/2LJ.BP) EMD (sphere/2LJ3CB.MSKM) EMD (3LJ3CB.EP/sphere) EMD (3LJ3CB.EP/2LJ3CB.MSKM)

ZSM-5 (M3) ZSM-5 (M3) ZSM-5 (M3) silicalite silicalite silicalite silicalite silicalite silicalite silicalite silicalite silicalite

303 303 303 300/298 303/300 303/298 300/298 303/300 303/298 300/298 303/300 303/298

We can see also in Table 4 the large differences between macroscopic results for n-butane (Van-den-Begin et al.24 (1989) against Hayhurst and Paravar9 (1988)) and for hexane (Beschmann et al.18 (1990) against Wu and Ma19 (1984)). Recently, Millot et al.10 reported macroscopic-supported membrane permeation measurements operating at steady state with no concentration gradient along the surface of membrane. They compared their measurements of permeation of isobutane against microscopic QENS values and found good agreement, within experimental error. (This has been made possible thanks to recent developments of tubular membranes of high quality.) Another interesting development is the macroscopic experimental method proposed by Nijhuis et al.,22,23 which operates at very low pore occupancies. This method uses an experimental apparatus referred to as Multitrack, based on the TAP technique; we denote it here as TAP. A result from its application is given in Table 4. Diffusivity values for n-butane based on TAP22,23 are lower by only a factor of 3 than those computed from EMD simulations and very close to PFG-NMR and QENS results; there is still some discussion in the literature concerning the interpretation of TAP measurements, however.60 (4) We have simulated silicalite as a perfect, fully siliceous crystal. In macroscopic measurements, defects and aluminum (acid) sites are inevitably present, leading to a decrease in the mobility of sorbates. Intercrystalline boundaries, gel remaining between crystals, or nonidealities in the crystals have been proposed as sources of additional resistance.74 Such defects are probably responsible for the fact that PFG-NMR and QENS diffusivities obtained in real silicalite crystals for n-butane and n-hexane are lower than EMD values by a factor of 3 or so (see Table 4). (5) Most importantly, macroscopic measurements of permeability through zeolite membranes exhibit large variations, depending on the type of zeolite layer and supporting material used.14-16 The large differences seen in Table 5 between CH4 permeabilities measured in different macroscopic-supported membranes (Yan et al.8 and Bakker et al.17) illustrate this point characteristically. Clearly, the transport resistance in the case of the metal support (used by Bakker et al.17) is much lower than that measured in the case of ceramic supports (used by Yan et al.8). Pochusta et al.75 presented permeances for CO2 and CH4 that are 2 orders of magnitude larger than those of Yan et al.8 and concluded that the differences were attributable to different relative amounts of zeolitic and nonzeolitic pores. Coronas et al.76 have also found permeances for N2 that are more than 1 order of magnitude higher than those of Yan et al.;8 they provided an explanation in terms of transport through pathways formed by the grain boundaries. In zeolite membranes grown on ceramic supports, there is intimate and intricate intergrowth of the zeolite with the porous support, with the result being that the effective thickness of the zeolite phase may be much larger than the nominal thickness L, typically attributed

to it.8 The simulation gives the permeability of a hypothetical membrane composed of perfect silicalite crystals, which is an upper bound for nonleaky real membranes. It is noteworthy that the permeability estimated for CH4 based on the experimental isotherm and the PFG-NMR or QENS diffusivity is actually in excellent agreement with the simulated value (Table 5). 5. Conclusions GCMC predicts the sorption equilibria of CO2 and N2 successfully. In agreement with the experiment, CO2 was found to exhibit greater thermodynamic affinity (more steeply rising isotherms) for the zeolite than N2. A variety of models have been used here for the representation of sorbate-zeolite and sorbate-sorbate interactions. In the case of CO2, the 2LJ3CB.EPM2 model has slightly better behavior than that of 2LJ3CB.EP in the Henry’s law region. In the case of N2, all models cast in terms of partial charges and LJ short-range interactions agree reasonably well with each other. They reproduce very well experimental sorption equilibria near 300 K, as well as the high-pressure part of the 77 K isotherm, including the steplike feature it exhibits between 22 and 28 molecules/unit cell. At 77 K and low pressures, the LJ plus partial-charge models of N2 overpredict the sorbed amount. This deficiency can be corrected by resorting to a computationally more demanding model cast in terms of modified BUCK shortrange interactions plus four partial charges (2BUCK4CB.BA). A very simple and computationally inexpensive model of N2 (2LJ.BP) performs surprisingly well in predicting the experimental 77 K isotherm, especially at low occupancies; this is probably related to the fact that the 2LJ.BP model has been optimized with respect to vapor-liquid equilibria at low temperatures. The same charge-free model underpredicts sorption at 300 K. Inclusion of short-range interactions with the framework silicon atoms does not generally improve the quality of isotherm predictions. On the basis of these results, we would recommend the use of rigid polyatomic models with partial charges, such as 3LJ3CB.EPM2 for CO2 and 2LJ3CB.MSKM or 2BUCK4CB.BA for N2 in conjunction with a rigid framework representation such as LJCB.JBTLC or BUCKCB.JBTLCRSO for the prediction of sorption equilibria in silicalite. Calculations of N2 sorption in the three crystal forms (ortho, mono, and para) of silicalite unambiguously exclude the possibility that the zeolite may be in the para form under the conditions considered; the sorption behavior obtained in the ortho and mono forms, on the other hand, is very similar. EMD was employed for the calculation of the self-diffusivities of CO2 and N2 in ortho silicalite at 300 K for a variety of pairwise additive effective potentials. At room temperature and low occupancy, N2 was found to diffuse faster (by a factor of 2-3) than CO2. Thus, the thermodynamic and diffusivity factors of selectivity for low-pressure separation of N2 and CO2 in

Permeation Properties of CO2 and N2

J. Phys. Chem. B, Vol. 105, No. 4, 2001 787

silicalite go in opposite directions, with the thermodynamic factor being the dominant one. Using the Darken equation and assuming a concentration-independent corrected diffusivity D0, we have calculated intracrystalline permeabilities and their ratios for the three gases CO2, N2, and CH4 from our simulation results; in the case of CH4, this calculation was based on past simulation investigations using the LJ.JBT model for the zeolite. Simulated permeability ratios turn out to be in very good agreement with experimental ones, following the order CO2 > CH4 > N2. Again, models of CO2 and N2 incorporating partial charges prove superior to charge-free models in the comparison of permeability ratios against those of the experiment. Simulation-based absolute values of the permeabilities, however, are too high by a factor of 120-340 relative to reported experimental values from silicalite membrane permeation measurements, which themselves may vary by 2 orders of magnitude between different membranes. The main reason for this discrepancy is thought to be the complex morphology of supported silicalite membranes. Intergrowths of the zeolite with the mesoporous support material may result in only a small fraction of the membrane having the nominal thickness that is attributed to it and used in reporting experimental permeability values. Furthermore, structural defects that may be present within the crystallites or in the grain boundaries between the crystallites of real membranes but are certainly absent from the computer models may contribute additionally to the disparity between the experimental and simulated permeability values.

To determine the thermodynamic behavior (sorption) of the fluid in a membrane, we must measure or predict the solubility coefficient S, defined61 as

S ) c/f

We define the local permeability coefficient Pe, which gives us information on the relationship between flux and driving force (fugacity gradient) locally at the considered point as a product of fundamental thermodynamic and dynamic properties:

Pe ) D0S

1 Ds ) lim 〈|r(t) - r(0)|2〉 6 tf∞

(A1)

J ) -D0S∇f ) -Pe∇f

(A7)

On the other hand, macroscopic methods, such as the measurement of the permeation rate, report the transport diffusivity Dt, defined through Fick’s law

J ) -Dt∇c

(A8)

A comparison of eqs A5-A8 gives

Dt ) D0S

c df ∂ ln f df ) D0 ) D0 (Darken equation) dc f dc ∂ ln c (A9) Dt ) P e

df dc

(A10)

The thermodynamic correction ∂ ln f/∂ ln c approaches 1 as c f 0 (linearity between fugacity and concentration in the Henry’s law region); thus, transport and corrected diffusivity become equal in this limit. The corrected diffusivity can also be expressed microscopically in terms of velocity time-correlation functions for the same or different molecules as1

D0 )

or equivalently by the Green-Kubo equation

(A6)

Now the flux density can be written as

Appendix From Self-Diffusion Coefficient and Sorption Isotherm to Permeability. From an EMD simulation we can predict the selfdiffusion coefficient Ds, which characterizes the displacement of an individual molecule in a “sea” of other molecules and is defined by the Einstein equation:

(A5)

1

N

N

∫∞dt ∑∑〈vi(0)‚vj(t)〉 3N 0

(A11)

i)1 j)1

Ds )

1

N

∞ dt ∑〈vi(0)‚vi(t)〉 ∫ 0 3N

(A2)

i)1

where N is the number of molecules considered, r(t) is the position vector of a molecule at time t, vi(t) is the velocity of sorbate molecule i, and the angular brackets denote averaging over all equilibrium dynamical trajectories. In the membrane, the molecules are considered to move under isothermal conditions with average velocity v subject to a driving force ∇µ (chemical potential gradient of the penetrant in the membrane; µ ) µ0 + RT ln f in terms of the fugacity f). So the flux (molecular current density) of fluid will be59

c c J ) cv ) - ∇µ ) -D0 ∇f ζ f

(A3)

where c is the concentration of fluid (molecules/unit volume) at fugacity f and temperature T. D0 is the themodynamic diffusion coefficient or “corrected” diffusivity61

D0 ) RT/ζ

(A4)

and ζ a friction coefficient which governs the dynamic behavior of the system.

From equations A11 and A2, we can easily see that the corrected diffusivity D0 is a collective property that depends on the total flux of all sorbate molecules, while the self-diffusion coefficient Ds depends on the autocorrelation of individual molecule velocities. At low occupancies, we can assume that the velocities of different molecules are uncorrelated; from eqs A2 and A11, we can conclude that, under these conditions, the two coefficients are equal, Ds ) D0. To estimate this low-occupancy value of D0, we compute Ds by EMD simulations at a number of low, finite occupancies and then extrapolate the results to zero occupancy in order to determine D0 ) Ds at zero coverage. On the basis of past simulations1 and experimental findings, we assert that D0 is not strongly occupancy independent; we therefore use the value of D0 estimated by extrapolation of the simulation values of Ds to zero occupancy within eq A9 to obtain Dt at every occupancy. In an actual steady-state permeation experiment, we measure not the local permeability but an average permeability Pe through the entire membrane. This is because we know the overall fugacity drop across the membrane of thickness L, ∆f ) f1 - f2, but not the local values of fugacity inside the membrane, as we do in a simulation. To estimate the average permeability Pe, we integrate eq A8 across the membrane,

788 J. Phys. Chem. B, Vol. 105, No. 4, 2001

Makrodimitris et al.

using the Darken equation (eq A9) and the assumption of a concentration-independent D0. So we have

Pe ≡

D0 J )∆f/L ∆f

D

∫cc ∂∂ lnln cf dc ) - ∆f0∫f f 2

1

c(f) df f

2

1

(A12)

where c(f) is the function expressing the adsorption isotherm of our gas (CO2 or N2), as obtained from GCMC simulations. Average permeabilities computed in this way are compared to experimental values, with f1 and f2 being the fugacities on the high- and low-pressure sides of the membrane, respectively, used under the conditions of the experiment. When f1 is low enough to lie within the Henry’s law region of the isotherm, the ratio c(f)/f is constant over the domain of integration, and the permeability turns out to be equal to D0 times that constant, i.e., the low-occupancy slope of the isotherm. References and Notes (1) Maginn, Ed. J.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1993, 97, 4173. (2) Schuring, D.; Jansen, A. P. J.; van Santen, R. A. J. Phys. Chem. B 2000, 104, 941. (3) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1992, 96, 1051. (4) Caro, J.; Bu¨llow, M.; Schirmer, W.; Ka¨rger, J.; Heink, W.; Pfeifer, H.; Zdanov, S. J. Chem. Soc., Faraday Trans. 1985, 81, 2541. (5) Heink, W.; Ka¨rger, J.; Pfeifer, H. J. Chem. Soc., Faraday Trans. 1992, 88, 3505. (6) Jobic, H.; Be´e, M.; Caro, J.; Ka¨rger, J. In Zeolites for the Nineties, Recent Research reports; Jansen, J. C., Moscou, L., Post, M. F. M., Eds.; Elsevier: Amsterdam, The Netherlands, 1989. (7) Jobic, H.; Be´e, M.; Caro, J. In Procceedings from the Ninth International Zeolite Conference; von Ballmoos, R., Higgins, J. B., Treacy, M. M. J., Eds.; Butterworth-Heineman: Boston, MA, 1993. (8) Yan, Y.; Davis, M. E.; Gavalas, G. R. Ind. Eng. Chem. Res. 1995, 34, 1652. (9) Hayhurst, D. T.; Paravar, A. R. Zeolites 1988, 8, 27. (10) Millot, B.; Methivier, A.; Jobic, H.; Moueddeb, H.; Be´e, M. J. Phys. Chem. B 1999, 103, 1096. (11) van den Broeke, L. J. P.; Bakker, W. J. W.; Kaptein, F.; Moulijn, J. A. Chem. Eng. Sci. 1999, 54, 245. (12) Kusakabe, K.; Yoneshige, S.; Murata, A.; Morooka, S. J. Membr. Sci. 1996, 116, 39. (13) Kusakabe, K.; Sakamoto, S.; Toshiyuki, S.; Morooka, S. Sep. Purif. Technol. 1999, 16, 139. (14) Chau, J. L. H.; Tellez, C.; Yeung, K. L.; Ho, K. J. Membr. Sci. 2000, 164, 257. (15) Bakker, W. J. W.; van den Broeke, L. J. P.; Kaptein, F.; Moulijn, J. A. AIChE J. 1997, 43, 2203. (16) Lovallo, M.; Tsapatsis, M. AIChE J. 1996, 42, 3020. (17) Bakker, W. J. W.; Kaptein, F.; Poppe, J.; Moulijn, J. A. J. Membr. Sci. 1996, 117, 57. (18) Beschmann, K.; Fuchs, J.; Riekert, L. Zeolites 1990, 10, 798. (19) Wu, P.; Ma, Y. H. In New DeVelopments in Zeolite Science and Technology, Procceedings of the 6th International Zeolite Conference; Olsen, D., Bisio, A., Eds.; Butterworth: Guildford, U.K., 1984. (20) Sun, M. S.; Talu, O.; Shah, D. B. AIChE J. 1996, 42, 3001. (21) Loos, J.-B W. P.; Verheijen, J. T.; Moulijn, J. A. Chem. Eng. Sci. 2000, 55, 51. (22) Nijhuis, T. A.; van den Broeke, L. J. P.; van de Graaf, J. M.; Kapteijn, F.; Makkee, M.; Moulijn, J. A. Chem. Eng. Sci. 1997, 52, 3401. (23) Nijhuis, T. A.; Linders, M. J. G.; Makkee, M.; Kapteijn, F.; Moulijn, J. A. Chem. Eng. Sci. 2000, 55, 1939. (24) Van-den-Begin, N.; Rees, L. V. C.; Caro, J.; Bu¨llow, M. Zeolites 1989, 9, 287. (25) Barrer, R. M. J. Chem. Soc., Faraday Trans. 1990, 86, 1123. (26) Pohl, P. I.; Heffelfinger, G. S.; Smith, D. M. Mol. Phys. 1996, 89, 1725. (27) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; Wiley-Interscience: New York, 1992. (28) Demontis, P.; Suffritti, G. B.; Fois, E. S.; Quartieri, S. J. Phys. Chem. 1992, 96, 1482. (29) Mentzen, B. F.; Sacerdotte-Perronnet, M. Mater. Res. Bull. 1993, 28, 1017. (30) Olson, D. H.; Kokotailo, G. T.; Lawton, S. L. J. Phys. Chem. 1981, 85, 2238. (31) van Koningsveld, H.; Tuinstra, F.; van Bekkum, H.; Jansen, J. C. Acta Crystallogr. 1989, B45, 423.

(32) van Koningsveld, H.; van Bekkum, H.; Jansen, J. C. Zeolites 1990, 10, 235. (33) Lonsinger, S. R.; Chakraborty, S. R.; Theodorou, D. N. Catal. Lett. 1991, 11, 209. Cook, S. J.; Chakraborty, A. K.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1993, 97, 6679. (34) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1993, 97, 13742. (35) Kiselev, A. V.; Lopatkin, A. A.; Meier, W. M. Zeolites 1985, 5, 261. (36) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 1508. (37) Watanabe, K.; Austin, N.; Stapleton, M. R. Mol. Simul. 1995, 15, 197. (38) van Genechten, K. A.; Mortier, W. J. Zeolites 1988, 8, 273. (39) Reischman, P. T.; Schmitt, K. D.; Olson, D. H. J. Phys. Chem. 1988, 92, 5165. (40) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (41) Errington, J.; Panagiotopoulos, A. Z. Unpublished data, http:// thera.umd.edu/jerring/gcmc/examples/co2/index.htm. (42) Murthy, C. S.; Singer, K.; Klein, M. L.; McDonald, I. R. Mol. Phys. 1980, 41, 1387. (43) Murthy, C. S.; O’Shea, S. F.; McDonald, I. R. Mol. Phys. 1983, 50, 531. (44) Berns, R. M.; van der Avoird, A. J. Chem. Phys. 1980, 72, 6107. (45) Baranello, P.; Panagiotopoulos, A. Z. Unpublished data, http:// thera.umd.edu/ppe/old/pbaranel/paper.htm, 1997. (46) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987. (47) Adams, D. J. Mol. Phys. 1974, 28, 1241. (48) Adams, D. J. Mol. Phys. 1975, 29, 307. (49) Norman, G. E.; Filinov, V. S. High Temp. 1969, 7, 216. (50) Mountain, R. D.; Thirumalai, D. Physica A 1994, 210, 453. (51) Marsaglia, G. Ann. Math. Stat. 1972, 43, 645. (52) Mazur, Al. J. Comput. Phys. 1997, 93, 171. (53) Frenkel, D.; Smit, B. Understanding Molecular simulations; Academic Press: London, 1996. (54) Fincham, D. Mol. Simul. 1993, 11, 79. (55) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 11990. (56) Yamazaki, T.; Masahiro, K.; Sentaro, O.; Yoshisada, O. Mol. Phys. 1993, 80, 313. (57) Douguet, D.; Pellenq, R. J.-M.; Boutin, A.; Fuchs, A.; Nicholson, D. Mol. Simul. 1996, 17, 255. (58) Mu¨ller, U.; Reichert, H.; Robens, E.; Unger, K. K.; Grillet, Y.; Rouquerol, F.; Rouquerol, J.; Donfgeng, P.; Mersmann, A. Fresenius’ Z. Anal. Chem. 1989, 333, 433. (59) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 8232. (60) Brandani, S.; Ruthven, D. Chem. Eng. Sci. 2000, 55, 1935. (61) Petropoulos, J. H. AdV. Polym. Sci. 1985, 64, 93. (62) Ruthven, D.; Brandani, S. In Physical Adsorption: Experiment, Theory and Application; Fraissard, J., Ed.; NATO Advances Study Institute Series 491, p 261; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1997. (63) Shen, D.; Rees, L. Zeolites 1991, 11, 666. (64) Zikanova, A.; Bu¨llow, M.; Schlodder, H. Zeolites 1987, 7, 115. (65) Ruthven, D.; Eic, M.; Richard, E. Zeolites 1991, 11, 647. (66) Shen, D.; Rees, L. J. Chem. Soc., Faraday Trans. 1994, 90, 3011. (67) Ka¨rger, J.; Pfeifer, H.; Stallmach, F.; Feoktistova, N. N.; Zhdanov, S. P. Zeolites 1993, 13, 50. (68) Goodbody, S. J.; Watanabe, K.; Gowan, D. M.; Walton, J. P. R. B.; Quirke, N. J. Chem. Soc., Faraday Trans. 1991, 87, 1951. (69) Talu, O.; Sun, M. S.; Shah, D. B. AIChE J. 1998, 44, 681. (70) Vroon, Z. A. E. P.; Keizer, K.; Burgraaf, A. J.; Verweij, H. J. Membr. Sci. 1996, 116, 39. (71) Ka¨rger, J.; Pfeifer, H. J. Chem. Soc., Faraday Trans. 1991, 87, 1989. (72) Hargreaves, M.; McLeod, A. S.; Gladden, L. F. In Fundamentals of Adsorption 6; Meunier, F., Ed.; Elsevier: Amsterdam, The Netherlands, 1998. (73) Shah, D. B.; Hayhurst, D. T.; Evanina, G.; Guo, C. J. AIChE J. 1988, 34, 1713. (74) Graaf, J. M.; Kapteijn, F.; Moulijn, J. Microporous Mesoporous Mater. 2000, 35-36, 267. (75) Pochusta, J. C.; Noble, R. D.; Falconer, J. L. J. Membr. Sci. 1999, 160, 115. (76) Coronas, J.; Noble, R. D.; Falconer, J. L. AIChE J. 1997, 43, 1797.