Thermodynamically Consistent Methodology for Estimation of Diffusivities of Mixtures of Guest Molecules in Microporous Materials.

The Maxwell-Stefan (M-S) formulation, that is grounded in the theory of irreversible thermodynamics, is widely used for describing mixture diffusion in microporous crystalline materials such as zeolites and metal-organic frameworks (MOFs). Binary mixture diffusion is characterized by a set of three M-S diffusivities: Đ 1, Đ 2, and Đ 12. The M-S diffusivities Đ 1 and Đ 2 characterize interactions of guest molecules with pore walls. The exchange coefficient Đ 12 quantifies correlation effects that result in slowing-down of the more mobile species due to correlated molecular jumps with tardier partners. The primary objective of this article is to develop a methodology for estimating Đ 1, Đ 2, and Đ 12 using input data for the constituent unary systems. The dependence of the unary diffusivities Đ 1 and Đ 2 on the pore occupancy, θ, is quantified using the quasi-chemical theory that accounts for repulsive, or attractive, forces experienced by a guest molecule with the nearest neighbors. For binary mixtures, the same occupancy dependence of Đ 1 and Đ 2 is assumed to hold; in this case, the occupancy, θ, is calculated using the ideal adsorbed solution theory. The exchange coefficient Đ 12 is estimated from the data on unary self-diffusivities. The developed estimation methodology is validated using a large data set of M-S diffusivities determined from molecular dynamics simulations for a wide variety of binary mixtures (H2/CO2, Ne/CO2, CH4/CO2, CO2/N2, H2/CH4, H2/Ar, CH4/Ar, Ne/Ar, CH4/C2H6, CH4/C3H8, and C2H6/C3H8) in zeolites (MFI, BEA, ISV, FAU, NaY, NaX, LTA, CHA, and DDR) and MOFs (IRMOF-1, CuBTC, and MgMOF-74).

Introduction

Many separation and reaction processes use microporous crystalline materials such as zeolites (alumino-silicates), metal–organic frameworks (MOFs), and zeolitic imidazolate frameworks (ZIFs) as perm-selective membrane layers, adsorbents, or catalysts.[1−8] Separation of nitrogen/oxygen, nitrogen/methane, and propene/propane mixtures in adsorbers packed with LTA-4A, CHA, or ZIF-8 are essentially driven by differences in the pore diffusivities of the guest constituents.[9−16] The conversion and selectivity of several heterogeneous catalytic reactions are influenced by intraparticle transport of reactants and products.[17−19] Pore diffusion characteristics have a significant influence on membrane separation selectivities.[4,20−25]

Mixture diffusion in microporous materials is characterized by the fact that the mobility of any guest constituent is influenced by its partner species;[20,26] the proper modeling of such influences is essential for process development and design.[6,15] It is a common practice to model n-component mixture diffusion by adopting the Maxwell–Stefan (M–S) formulation,[1,6,7,27,28] that has its foundations in the theory of nonequilibrium thermodynamics. The dependence of the intracrystalline molar fluxes, N, on the chemical potential gradients is written in the following form

In eq , R is the gas constant, ρ represents the material framework density, and the component loadings, and q are defined in terms of moles per kg of framework material. The x in eq is the mole fractions of the adsorbed phase components

Two distinct sets of M–S diffusivities are defined by eq , that is phenomenological in nature.[29] The Đ characterize interactions between species i with the pore walls. As established in earlier works,[23,28] the important advantage of the M–S formulation is that the Đ can be identified with the corresponding unary diffusivities, provided the diffusivity data are compared at the same adsorption potential, πA/RT, where A represents the surface area per kg of framework material, and π is the spreading pressure, defined by the Gibbs adsorption equation[3,30,31]

The ideal adsorbed solution theory (IAST) of Myers and Prausnitz[32] enables the calculation of the adsorption potential, πA/RT;[23,28] details are provided in the Supporting Information.

As illustration, Figure presents plots of the transport coefficients, ρĐ/δ, of CO2, CH4, N2, and H2 determined for unary and equimolar binary (CO2/H2, CO2/CH4, CO2/N2, CH4/Ar, CH4/N2, CH4/H2, and N2/H2) mixture permeation across a SAPO-34 membrane of thickness δ.[23,25,33,34] SAPO-34 has the same structural topology as CHA zeolite, consisting of cages of volume 316 Å3, separated by 8-ring windows of 3.8 × 4.2 Å size.[23] Compared at the same value of πA/RT, the magnitude of ρĐ/δ for a binary mixture is comparable to that for the corresponding pure component. Also noteworthy, from the data in Figure is that the dependence of ρĐ/δ on πA/RT is distinctly different for each guest molecule.

Figure 1

Experimental data of Li et al.[25,33,34] for transport coefficients ρĐ/δ of (a) CO2, (b) CH4, (c) N2, and (d) H2 determined for unary and equimolar binary CO2/H2, CO2/CH4, CO2/N2, CH4/Ar, CH4/N2, CH4/H2, and N2/H2 mixture permeation across SAPO-34 membrane at 295 K. The data are plotted as a function of the adsorption potential, πA/RT, calculated for conditions prevailing at the upstream face of the membrane; the calculation details are provided in an earlier work.[23]

The exchange coefficients, Đ, defined by the first right member of eq reflect how the facility for transport of species i correlates with that of species j. The Onsager reciprocal relations impose the symmetry constraint

The magnitude of Đ relative to that of Đ determines the extent to which the flux of species i is influenced by the chemical potential gradient of species j. The larger the degree of correlations, Đ/Đ, the stronger is the influence of diffusional “coupling”. Generally speaking, the more strongly adsorbed tardier partner species will have the effect of slowing down the less strongly adsorbed more mobile partner in the mixture.

For estimation of the exchange coefficient, Đ, the following interpolation formula has been suggested in the literature.[28,35]where the Đ and Đ represent the self-exchange coefficients, that are accessible from molecular dynamics (MD) simulations of self-diffusivities for the constituent unary systems,[26] as will be discussed in a subsequent section. Equation is essentially an adaptation of the interpolation formula for estimation of the M–S diffusivity for binary fluid mixtures.[36]

Specifically, for a binary mixture, that is n = 2, the M–S eq can be rewritten to evaluate the fluxes N explicitly by defining a 2 × 2 dimensional square matrix [Λ]

The elements of [Λ] are directly accessible from MD simulations[28,35,37,38] by monitoring the individual molecular displacements

In 7n and n represent the number of molecules of species i and j, respectively, and r(t) is the position of molecule l of species i at any time t.

Combining eq with eq , the following explicit expression for calculation of the elements of the 2 × 2 dimensional square matrix [Λ] can be derived

The primary objective of this article is to seek validation of the predictive capability of the Maxwell–Stefan formulation by comparing each of the four elements, Λ, determined from MD simulations using eq , with the estimations using eqs and 8; the required data inputs for Đ1, Đ2, Đ11, and Đ22 are determined from MD simulations for the constituent pure components. To meet the stated objective, use is made of the MD simulation data base compiled in an earlier work[28] for a variety of binary mixtures (H2/CO2, Ne/CO2, CH4/CO2, CO2/N2, H2/CH4, H2/Ar, CH4/Ar, Ne/Ar, CH4/C2H6, CH4/C3H8, and C2H6/C3H8) in different host materials. The host materials were chosen to represent a variety of pore sizes, topologies, and connectivities: one-dimensional (1D) channels (e.g., AFI, MgMOF-74), intersecting channels (MFI, BEA, ISV), cages separated by narrow windows (LTA, CHA, DDR), cavities with large windows (FAU, NaY, NaX, IRMOF-1, CuBTC). The Supporting Information provides structural details for all of the host materials considered in this article.

Results and Discussions

Occupancy Dependence of Unary Diffusivities

For any guest molecule, the loading dependence of Đ is strongly influenced also by the pore topology and connectivity and molecule–molecule interactions.[7,8,23,28,39−43] As an illustration, Figure a,b presents data on Đ for the guest species CH4 in a variety of host structures, determined from MD simulations of molecular displacements using the following formula in each of the coordinate direction

Figure 2

MD simulations of (a,b) Maxwell–Stefan diffusivity, Đ, (c) self-diffusivities D, and (d) degree of correlations, Đ/Đ for CH4 in a variety of host structures plotted as a function of the occupancy, θ, determined from eq , where qsat and π are determined from the unary isotherm fits that are specified in the Supporting Information accompanying this article.

The M–S diffusivity, Đ displays a wide variety of dependencies on the fractional occupancy, θ, that serves as a convenient and practical proxy for the adsorption potentialwhere qsat is the saturation capacity.[23,28] The calculations of the adsorption potential in eq use dual-site Langmuir–Freundlich fits of the unary isotherms that are determined from configurational-bias Monte Carlo (CBMC) simulations[44−48] for each guest/host combination; details are provided in the Supporting Information.

For CH4/BEA and CH4/NaX, the Đ appears to decrease almost linearly with occupancy θ till pore saturation conditions, θ = 1, are reached. An appropriate model to describe this occupancy dependence iswhere Đ(0) is the M–S diffusivity at “zero-loading”. Equation is essentially based on a hopping model in which one molecule at any time can jump from one sorption site to an adjacent one, provided it is not already occupied.[23,28,49,50] Using a simple two-dimensional square lattice model, the M–S diffusivity in the limit of vanishingly small occupancies, , where ζ = 4 is the coordination number of the 2D array of lattice sites, λ is the jump distance on the square lattice, and ν(0) is the jump frequency at vanishingly small occupancy.[50]

More generally, molecule–molecule interactions serve to influence the jump frequencies by a factor that depends on the energy of interaction, w. For repulsive interactions, w > 0, whereas for attractive interactions, w < 0. Using the quasi-chemical approach of Reed and Ehrlich[51] to quantify such interactions, the following expression is obtained for the occupancy dependence of the M–S diffusivities[50,52,53]

In 12 the following dimensionless parameters are defined

In the limiting case of negligible molecule–molecule interactions, w = 0, ϕ = 1, β = 1, eqs and 13 degenerate to yield eq . The continuous solid lines in Figure a,b are fits of the MD simulated Đ by fitting the sets of parameters: Đ(0), and ϕ = ϕ0 exp(−aθ). For all of the guest/host combinations, eqs and 13 provide good descriptions of the occupancy dependencies; see Figures S38–S108 of the Supporting Information.

Applying eq to a binary mixture consisting of tagged and untagged species i, that are otherwise identical,[7,50,54] we can derive the following relation between the self-diffusivity, D and the M–S diffusivity, Đ

The self-diffusivities, D may be computed from MD simulations by analyzing the mean square displacement of each species, i for each coordinate direction[7]

By combination of eqs , 14, and 15, we can determine the degrees of correlations due to self-exchange, Đ/Đ. Figure c,d presents MD simulated data on D and Đ/Đ for CH4 in different host materials: MgMOF-74 (1D channels of 11 Å), BEA (intersecting channels of 6.5 Å), ISV (intersecting channels of 6.5 Å), NaX (790 Å3 cages separated by 7.4 Å windows), MFI (intersecting channels of 5.5 Å), and LTA (743 Å3 cages separated by 4.2 Å windows). It is also to be noted that the size of the 1D channels of MgMOF-74 are large enough to preclude single-file diffusion of guest molecules. The degree of correlations is the lowest for the LTA zeolite because the guest molecules jump one-at-a-time across the narrow 4.2 Å windows;[6,55,56] the same characteristics are valid for other cage-window structures with narrow windows, such as CHA, DDR, and ZIF-8.[24,57−59] The variation of Đ/Đ with occupancy is practically linear,[23,28] and the solid lines in Figure d are the linear fits

Equation provides a good description of the occupancy dependence of the degrees of correlations due to self-exchange for all guest/host combinations; see Figures S38–S108 of the Supporting Information.

Occupancy Dependence of [Λ] for Binary Mixture Diffusion

Having established and quantified the occupancy dependence of Đ and Đ/Đ for each guest/host combination, we are in a position to compare the estimations of Λ for binary mixtures using eqs and 8 with the corresponding MD simulated values by monitoring molecular displacements and use of eq ; Figures S38–S108 provide detailed comparisons for each mixture/host combination that was investigated. Figure provides an illustration of the estimation procedure for CH4(1)/C3H8(2) mixture diffusion in NaX zeolite. Figure a show the Reed–Ehrlich model fits for the unary diffusivities Đ1, Đ2 for CH4 and C3H8 in NaX. The linear fits for the degrees of self-exchange Đ1/Đ11 and Đ2/Đ22 are shown in Figure b. In Figure c, the MD simulation data for Λ11, Λ12 = Λ21, and Λ22 for equimolar (q1 = q2) binary CH4(1)/C3H8(2) mixtures are compared with the estimations (shown by the continuous solid lines) using eqs and 8. The Maxwell–Stefan formulation provides very good estimates of dependence of each Λ on the occupancy θ, calculated using eq , wherein the saturation capacity for the mixture is determined fromwhere q1,sat and q2,sat are the saturation capacities of components 1 and 2, respectively. Equation can be derived from the IAST, as detailed in the Supporting Information.

Figure 3

(a,b) MD simulated values of (a) Đ1, Đ2 and (b) Đ1/Đ11, Đ2/Đ22 for guest molecules CH4 and C3H8 in NaX zeolite (86 Al) at 300 K. The continuous solid lines are the fits of the unary diffusivities. (c) MD simulation data for Λ11, Λ12 = Λ21, and Λ22 for equimolar (q1 = q2) binary CH4(1)/C3H8(2) mixtures in NaX zeolite (86 Al) at 300 K, compared with the estimations (shown by continuous solid lines) using eqs and 8.

Similar good estimates of the M–S model are established in Figure for six other equimolar (q1 = q2) binary mixtures: CO2/H2 in MFI, CH4/C3H8 in BEA, CH4/C2H6 in NaY, CH4/CO2 in IRMOF-1, CO2/CH4 in MgMOF-74, and Ne/Ar in CuBTC.

Figure 4

MD simulation data for Λ11, Λ12 = Λ21, and Λ22 for equimolar (q1 = q2) binary (a) CO2(1)/H2(2) in MFI zeolite, (b) CH4(1)/C3H8(2) in BEA zeolite, (c) CH4(1)/C2H6(2) in NaY zeolite (48 Al), (d) CH4(1)/CO2(2) in IRMOF-1, (e) CO2(1)/CH4(2) in MgMOF-74, and (f) Ne(1)/Ar(2) in CuBTC at 300 K, compared with the estimations (shown by continuous solid lines) using eqs , and 8.

A different test of the predictive capability of M–S formulation is to consider diffusion in binary mixtures for which the total loading q1 + q2 is held constant, and the mole fraction of component 1 in the adsorbed mixture, x1, is varied from 0 to 1.[60] One of the earliest investigations of this type were reported by Snurr and Kärger[61] for CH4/CF4 diffusion in MFI zeolite at a total loading of 12 molecules uc–1.

Figure compares the MD simulation data for Λ for binary Ne(1)/Ar(2) mixtures of varying composition x1 in MFI, LTA, CHA, and DDR zeolites. In all four cases, eqs , and 8 provide good predictions of the variation of Λ with composition. It is also to be noted that the off-diagonal elements Λ12 for LTA, CHA, and DDR zeolites are significantly lower, by about an order of magnitude, than the diagonal elements, Λ11 and Λ22. For cage-type zeolites such as LTA, CHA, DDR, ERI with 8-ring windows in the 3.3–4.5 Å size range, the degree of correlations Đ/Đ are negligibly small because the guest molecules jump one-at-a-time across the narrow windows.[20,26−28,55,62]

Figure 5

MD simulation data for Λ11, Λ12, and Λ22 for binary Ne(1)/Ar(2) mixtures of varying composition at constant total loading qt = q1 + q2 in (a) MFI zeolite (qt = 12.5 molecules uc–1), (b) LTA all-silica zeolite (qt = 60 molecules uc–1), (c) CHA all-silica zeolite (qt = 25 molecules uc–1), and (d) DDR zeolite (qt = 40 molecules uc–1) at 300 K, compared with the estimations (shown by continuous solid lines) using eqs , and 8. Note that Λ21 = (x2/x1)Λ12 has not been plotted because it is not independent.

Further evidence of the good predictive capability of the M–S formulation is provided in Figures S38–S108.

Preferential Perching of CO2 in Window Regions of Cage-Type Zeolites

For separation of CO2 from gaseous mixtures containing CH4, H2, N2, Ar, or Ne, cage-type zeolites such as DDR CHA, LTA, and ERI are of practical interest.[8,24,30,31,46,47,55,56] These materials consist of cages separated by narrow windows, in the 3.3–4.5 Å range. CBMC simulations[59] show that the window regions of cage-type zeolites have a significantly higher proportion of CO2 than within the cages. For all four zeolites, CO2 has the highest probability, about 30–40%, of locating at the window regions.[59] The preferential perching of CO2 in the window regions, evidenced by the computational snapshot for CHA (see Figure a), has the effect of hindering the intercage hopping of partner molecules.

Figure 6

(a) Computational snapshot showing the preferential perching of CO2 at the window regions of CHA zeolite.[59] (b–d) MD simulation data for Λ11 and Λ22 for equimolar (q1 = q2) binary CO2(1)/Ne(2) mixtures at 300 K (b) LTA all-silica zeolite, (c) CHA all-silica zeolite, and (d) DDR zeolite, compared with the estimations (shown by continuous solid lines) using eq and assuming that the degrees of correlations are negligible, that is Đ/Đ → 0.

Figure b–d compare the MD simulation data for Λ11 and Λ22 for equimolar (q1 = q2) binary CO2(1)/Ne(2) mixtures in LTA, CHA, and DDR zeolite, with the estimations using eq , assuming that the degrees of correlations are negligible, that is Đ/Đ → 0. For all three zeolites, the MD simulation data for Λ22 are significantly lower than the predictions using eq ; the M–S formulation does not cater for hindering effects caused due to segregated mixture adsorption. Experimental evidence of the importance of hindering effects is provided in published works on CO2/CH4 and CO2/N2 mixture permeation across the DDR membrane.[38,63,64] Analogous hindering effects are also evidenced for CO2/CH4 mixture permeation across ZIF-8 membranes.[24]

Preferential location of branched alkanes and aromatics at the intersections of MFI zeolite often cause intersection-blocking and loss of connectivity;[55,65,66] this leads to failure of the predictions of the M–S model.[67]

Molecular Clustering Due to Hydrogen Bonding

For water/methanol and water/ethanol mixture diffusion in microporous materials, molecular simulations[62,68−71] demonstrate the occurrence of molecular clustering due to hydrogen bonding. As a consequence of cluster formation, the diffusivities of either guest molecule in the mixture is significantly lower than the corresponding unary diffusivities. As illustration of mutual-slowing down effects, Figure presents MD data on the self-diffusivities, D, in water/methanol mixtures in CHA, DDR, and LTA zeolites, plotted as a function of the mole fraction of water in the adsorbed phase, x1. Each of the diffusivities is lowered due to the presence of its partner species. Experimental evidence of mutual-slowing down effects are available for water/alcohol permeation across CHA,[72,73] H-SOD,[74] and DDR[71,75] membranes. Further research is necessary to generalize the M–S formulation in a manner that explicitly allows for cluster formation, by defining a cluster to be a pseudospecies in the mixture.

Figure 7

MD simulations of self-diffusivities, D, of water(1)/methanol(2) mixtures of varying composition in (a) CHA, (b) DDR, and (c) LTA zeolites, plotted as a function of the mole fraction of water in the adsorbed phase, x1. In the MD simulations, the total loading, Θt, expressed as molecules uc–1, is held constant; the values Θt are specified. The MD data are culled from our previous publications.[69−71]

Conclusions

The capability of eqs and 8 for the estimation of the elements of the square matrix of M–S diffusivities, Λ, characterizing mixture diffusion, using input based on unary systems is tested using a large database obtained from MD simulations for a wide variety of guest/host combinations. The key to the estimation methodology is that the estimates are based on comparing the mixture diffusion data with those of the constituent unaries at the same fractional occupancy, θ, that is calculated on the basis of the adsorption potential using eq . For the majority of binary mixtures investigated, 70 in total, summarized in Figures S38–S108, the MD-simulated Λ data are in good agreement with the estimations using the M–S theory. The predictions of the M–S formulation are found to be poor for diffusion of CO2-bearing mixture in cage-type zeolites (LTA, CHA, DDR) wherein the preferential perching of CO2 at the window regions hinders the intercage hopping of partner molecules. The M–S predictions also fail to capture molecular clustering effects in water/alcohol systems that are engendered due to hydrogen bonding.[62]