Langmuir 2002, 18, 7393-7400
7393
Analysis of Binary Transport Diffusivities and Self-Diffusivities in a Lattice Model for Silicalite David Blanco Maceiras and David S. Sholl* Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213 Received May 21, 2002. In Final Form: July 6, 2002 We have used dynamic Monte Carlo simulations to assess the self-diffusivities and transport diffusivities of binary mixtures in a lattice model for silicalite as functions of hopping rates, mixture composition, and total loading. We compute binary transport diffusivities using the Onsager formalism and subsequently discuss the binary Fickian and Maxwell-Stefan diffusivities. We have used our results to examine the accuracy of several approximations to mixture transport diffusivities that have been suggested previously. Our numerical data strongly support the observation that single-component Maxwell-Stefan diffusivities in this model exactly obey a mean-field expression. Our data also show that a relationship between the binary self-diffusivities and the binary Maxwell-Stefan diffusivities that was claimed by several authors to be exact is in fact only approximate. We have similarly examined several methods that have been proposed for predicting binary self-diffusivities from single-component data.
1. Introduction The diffusion rates of molecules adsorbed inside zeolitic materials are fundamental to the performance of these materials in practical applications. Extensive reviews of this topic are available.1-4 Although essentially all practical applications of zeolites involve the coadsorption of multiple species, many features of multicomponent transport in zeolitic pores are not well understood. For example, the accurate prediction of multicomponent diffusivities in zeolites, even for systems where single-component diffusivities are known, remains challenging. Simulation of molecular transport in zeolites continues to play an important role in advancing our understanding of single-component and multicomponent adsorption and transport in zeolites. Ideally, one could examine the diffusion of molecules and mixtures of interest using fully atomistic models. While there are a number of studies that have examined the diffusions of self-diffusivity5,6 and transport diffusivity7,8 of adsorbed mixtures in this way, the computational expense of fully atomistic methods makes them impractical in many instances. As a result, many groups have examined molecular diffusion in zeolites using lattice models.3,4,9-14 Making rigorous connections * Corresponding author. E-mail:
[email protected]. Fax: 412-268-7139. Telephone: 412-268-4207. (1) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; John Wiley & Sons: New York, 1992. (2) Theodorou, D. N.; Snurr, R. Q.; Bell, A. T. In Molecular Dynamics and Diffusion in Microporous Materials; Alberti, G., Bein, T., Eds.; Comprehensive Supramolecular Chemistry, Vol. 7; Pergamon Press: New York, 1996; p 507. (3) Keil, F. J.; Krishna, R.; Coppens, M. O. Rev. Chem. Eng. 2000, 16, 71. (4) Auerbach, S. M. Int. Rev. Phys. Chem. 2000, 19, 155. (5) Snurr, R. Q.; Ka¨rger, J. J. Phys. Chem. B 1997, 101, 6469. (6) Gergidis, L. N.; Theodorou, D. N.; Jobic, H. J. Phys. Chem. B 2000, 104, 5541. (7) Sanborn, M. J.; Snurr, R. Q. Sep. Purif. Technol. 2000, 20, 1. (8) Sanborn, M. J.; Snurr, R. Q. AIChE J. 2001, 47, 2032. (9) Paschek, D.; Krishna, R. Phys. Chem. Chem. Phys. 2000, 2, 2389. (10) Paschek, D.; Krishna, R. Phys. Chem. Chem. Phys. 2001, 3, 3185. (11) Paschek, D.; Krishna, R. Langmuir 2001, 17, 247. (12) Trout, B. L.; Chakraborty, A. K.; Bell, A. T. Chem. Eng. Sci. 1997, 52, 2265. (13) Coppens, M. O.; Bell, A. T.; Chakraborty, A. K. Chem. Eng. Sci. 1998, 53, 2053.
between these lattice models and more detailed atomistic models remains challenging.15-17 Nevertheless, lattice models have the virtue that they can be used to explore broad ranges of parameters in a computationally efficient way and, as a result, provide insight into possible diffusion behaviors, at least within the assumptions made in these models. In this paper, we present simulation results for binary adsorbed mixtures of components with differing hopping rates in a simple lattice model for silicalite. We have examined both the binary self-diffusivities and the binary transport diffusivities (also known as the Fickian diffusivities). For both properties, several approximations have been previously suggested that predict the mixture properties on the basis of single-component data.11,18 Our simulation data provide a useful way to test the accuracy of these approximations. Our paper is organized as follows. In section 2, we define the lattice model for silicalite and the simulation techniques we have used to measure the relevant diffusivities. The binary transport diffusivities for the lattice model and the accuracy of several approximations to these transport diffusivities are discussed in section 3. We discuss both the matrix of Fickian diffusivities and the formally equivalent matrices of Onsager transport coefficients and Maxwell-Stefan diffusivities. In section 4, the binary self-diffusivities are treated in a similar manner. We summarize our conclusions in section 5. 2. Methods We have simulated binary adsorption and transport using the lattice model first proposed by Trout et al.12 In this model, the three-dimensional channel structure of silicalite is represented by a set of adsorption sites with equal energy lying at channel intersections and in the center of the straight and zigzag channels connecting these intersections. In simulations using two adsorbed species, (14) Coppens, M. O.; Bell, A. T.; Chakraborty, A. K. Chem. Eng. Sci. 1999, 54, 3455. (15) Sholl, D. S. Ind. Eng. Chem. Res. 2000, 39, 3737. (16) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2002, 106, 5058. (17) Tunca, C.; Ford, D. M. J. Chem. Phys. 1999, 111, 2751. (18) Krishna, R. Chem. Phys. Lett. 2000, 326, 477.
10.1021/la025972u CCC: $22.00 © 2002 American Chemical Society Published on Web 08/13/2002
7394
Langmuir, Vol. 18, No. 20, 2002
Maceiras and Sholl
both species were assumed to have the same binding energy in every site. For each species, the hopping rate for molecules in each of the four (two) directions available from intersection (channel) sites was assumed to be equal. The net hopping rate away from intersection sites is double the rate from channel sites. In all our simulations, the hopping rates of Trout et al. were used for species 1. These hopping rates were originally chosen to mimic benzene diffusion at 773 K.12 Species 2 was assumed to have a hopping rate either the same as or slower than that of species 1. The ratio of the hopping rates for species 1 and 2 was characterized by the ratio of the frequencies associated with hopping attempts for the two species, ν1/ ν 2. By construction, binary adsorption in this lattice model satisfies a mixed-species Langmuir isotherm at equilibrium. This lattice model has a saturation loading, Ns, of 12 molecules per unit cell for both adsorbed species. The fractional loading of species i is defined by θi ) Ni/Ns, where Ni is the loading of species i in molecules per unit cell. The mixed-species Langmuir isotherm for this lattice model is
θi )
Kipi 1 + K1p1 + K2p2
(1)
where Ki is the Henry’s law coefficient for species i and pi is the partial pressure of this species in the bulk phase. The characteristics of single-component diffusion in this lattice model are well-known.9,10,12-14,16 The MaxwellStefan diffusivity, ^, decreases linearly with increasing fractional loading. That is, ^(θ)/^(0) ) 1 - θ. The selfdiffusivity, Ds, can be written as Ds(θ) ) (1 - θ)f(θ), where f(θ) e 1 because of correlation effects. We verified that our simulations yielded single-component self-diffusivities in excellent agreement with the values reported by Trout et al.12 The Fickian diffusivity, D, which defines the mass flux generated by a concentration gradient according to
B J ) -D(c)∇c ) -NsFD(θ)∇θ
(2)
is independent of loading. That is, D(θ)/D(0) ) 1. In the expression above, Ns and F are the saturation loading in molecules per unit cell and the adsorbent density, respectively.10 The self-diffusivity, Maxwell-Stefan diffusivity, and Fickian diffusivity are all equal in the limit of dilute loadings, Ds(0) ) ^(0) ) D(0). Molecular diffusion in silicalite is in general anisotropic.5,16,19,20 Because the hopping rates are assumed to be the same in all directions in our lattice model, the only anisotropy in the x and y components of diffusion of this model arises because the unit cell of silicalite has slightly different lengths in these directions. Diffusion in the z direction is considerably slower than diffusion in the x and y directions, as occurs in more detailed models of silicalite.16 For single-component diffusion in this lattice model, the expression due to Ka¨rger16,21
a2 b2 c2 ) + Ds,z Ds,x Ds,y
(3)
is exact for all fractional loadings. Here, a, b, and c are the unit cell dimensions of silicalite. As a result, in the (19) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2001, 105, 3151. (20) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1993, 97, 4173. (21) Ka¨rger, J. J. Phys. Chem. 1991, 95, 5558.
remainder of this paper, we report only orientationally averaged diffusivities. The dynamics of this lattice model were examined using the dynamic Monte Carlo (MC) algorithm described in detail by Trout et al.,12 generalized in the natural way for binary mixtures. Starting from randomly distributed initial conditions, each system was equilibrated by following the dynamics until each particle had successfully hopped 3000 times on average. Subsequently, data for determining the diffusivities were recorded while the trajectory of the system was followed until each particle successfully made 6000 hops on average. Simulations were performed in simulation volumes varying from 4 × 4 × 4 unit cells at high total loadings to 15 × 15 × 15 unit cells at low total loadings. Diffusion coefficients were computed by averaging trajectory data obtained from multiple independent systems at each set of conditions examined. For our studies of loading dependence, the number of independent systems was varied from 600 at low loadings to 100 at high loadings. In simulations where we varied the composition of the adsorbed mixture, 150-200 independent systems were used for each loading. Binary self-diffusivities were computed from our MC simulations using the standard Einstein expression.11 Binary transport diffusivities are most conveniently measured from lattice or atomistic models using the Onsager formulation of binary transport. In this formulation, the flux of each species is related to gradients in the chemical potentials of each species through a symmetric matrix of Onsager transport coefficients,
[ ] [
][ ]
B J1 L11 L12 ∇µ1 B J2 ) - L21 L22 ∇µ2
(4)
with L12 ) L21. The Onsager transport coefficients can be computed from particle trajectories in a system at equilibrium using Einstein expressions analogous to those for self-diffusion.2,7,8,16 In our calculations of the binary selfdiffusivities and the Onsager coefficients, we used a socalled order-N algorithm for data analysis.19,22 Uncertainties in these quantities were estimated by separately analyzing multiple trajectories and calculating the deviation from the mean in the resulting diffusion coefficients. Once the binary Onsager transport coefficients are known, they can be used to generate the binary Fickian diffusivities or Maxwell-Stefan diffusivities, as discussed in section 3. Several groups have examined lattice models of silicalite related to the model of Trout et al. we have used here. Paschek and Krishna have simulated self-diffusion of methane-perfluoromethane mixtures in silicalite using a lattice model with a saturation loading of 24 molecules per unit cell for both species and slightly anistropic hopping rates.11 This model predicts loading-dependent singlecomponent self-diffusivities for both species in excellent agreement with those of atomistic simulations, although it does not correctly reproduce the Fickian diffusivities observed in atomistic simulations.16,19 Paschek and Krishna also examined single-component self-diffusion and transport diffusion in a lattice model for 2-methylhexane in silicalite that allowed adsorption only at channel intersections and in straight channels and again used anisotropic hopping rates.9 This lattice model was subsequently used for performing simulations of diffusion of binary mixtures and comparing these results to those found using a cubic lattice.10 The lattice model of Trout (22) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: London, 1996.
Analysis of a Lattice Model for Silicalite
Langmuir, Vol. 18, No. 20, 2002 7395
Figure 1. Onsager coefficients for a 50:50 adsorbed mixture with ν1/ν2 ) 10. MC results are shown for L11 (squares), L22 (triangles), and L12 (diamonds). Approximations to the Onsager coefficients based on estimating ^12(θ1,θ2) using eqs 9 and 11 are shown as solid curves.
et al. has been extended by Bell and co-workers to examine the influence of strong adsorption sites on singlecomponent diffusion.12-14 3. Binary Transport Diffusivities As described in section 2, we computed the Onsager transport coefficients directly from equilibrium simulations of our lattice model. Example results are shown in Figure 1 for an equimolar adsorbed mixture of species 1 and 2 with ν1/ν2 ) 10 as a function of the total loading. As should be expected, these results are similar to those reported by Paschek and Krishna for a related lattice model.10 The diagonal transport coefficient for species 1, L11, is considerably larger than the diagonal component for species 2, L22, at all loadings because of the higher hopping rate of species 1. It is interesting to note that, in this example, the off-diagonal transport coefficient, L12, is comparable in magnitude to L22 for moderate and high fractional loadings, unlike the case of the silicalite simulations of Paschek and Krishna,10 which used species with more similar hopping rates. L11 and L12 in this example have maximum values at fractional loadings significantly different from 0.5. Once the Onsager transport coefficients are known, it is straightforward to generate the matrix of binary Fickian diffusivities, also known as the transport diffusivities:10
[
] [
][
][
D11 D12 L11 L12 1/θ1 0 Γ11 Γ12 1/θ2 Γ21 Γ22 D21 D22 ) L21 L22 0
]
(5)
Here, Γij are the so-called thermodynamic correction factors, which can be evaluated explicitly for the Langmuir isotherm obeyed by our lattice model;10
Γij ) δij +
θi 1 - θ1 - θ2
(6)
with δij ) 1 if i ) j and 0 otherwise. The binary transport diffusivities derived in this way from the results in Figure 1 are shown in Figure 2. Recall that, for single-component diffusion in this model, the transport diffusivity is independent of loading. The diagonal components of the binary transport diffusivities, D11 and D22, approach the single-component diffusivities for the two species as the fractional loading approaches zero. In the equimolar mixture used for Figures 1 and 2, the diagonal components of the transport diffusivity both decrease as the fractional
Figure 2. Same as Figure 1 but showing the Fickian diffusivities D11 (squares), D22 (diamonds), D12 (circles), and D21 (triangles).
loading increases. In contrast, the off-diagonal transport diffusivities increase with loading. The off-diagonal transport diffusivities become comparable in magnitude to the diagonal components for moderate and high fractional loadings. These results are qualitatively similar to those obtained on a cubic lattice by Paschek and Krishna.10 It is also straightforward to generate the binary Maxwell-Stefan diffusivities once the Onsager coefficients are known. The matrix of Onsager coefficients can be related without approximation to the loading-dependent Maxwell-Stefan diffusivities by10,11
[
] [
L11 L12 B11 B12 L21 L22 ) B21 B22
where
][ ] -1
θ1 0 0 θ2
[
θ2 θ1 1 + B11 B12 ^1 ^12 ^12 θ2 θ1 B21 B22 ) 1 + ^12 ^2 ^12
[
]
(7)
]
(8)
We will first consider the diagonal Maxwell-Stefan diffusivities, ^1 and ^2, as functions of the binary loading, (θ1,θ2). Krishna and co-workers have previously suggested a mean-field expression for these diffusivities, ^i(θ1,θ2) ) Ds,i(0)(1 - θ1 - θ2).11,18 Examination of our numerical results provides very strong evidence that this expression is in fact exact for the lattice model we have simulated. In Figure 3, we plot ^1 and ^2 as determined from our measured Onsager coefficients. Despite the slight scatter in our numerical data, it is clear that the expression ^i(θ1,θ2) ) Ds,i(0)(1 - θ1 - θ2) describes the true result with great precision. For all the data shown in Figure 3, the mean-field expression lies within one standard deviation from the mean of our numerical results. Throughout the remainder of this paper, we will treat this expression as exact. Since we now have an exact expression for the diagonal Maxwell-Stefan diffusivities for binary mixtures in our lattice model, it is interesting to ask if the counterexchange coefficient, ^12(θ1,θ2), can be expressed in a simple form. Two approximations for ^12(θ1,θ2) have previously been suggested. The first is to estimate the counterexchange coefficient as the geometric mean of the single-component Maxwell-Stefan diffusivities,11
^12(θ1,θ2) ) x^1(θ1,θ2) ^2(θ1,θ2)
(9)
7396
Langmuir, Vol. 18, No. 20, 2002
Maceiras and Sholl
Figure 3. Normalized diagonal Maxwell-Stefan diffusivities, ^1/^1(0) (filled symbols) and ^2/^2(0) (open symbols), for binary diffusion in the silicalite lattice model as functions of mixture composition and loading with ν1/ν2 ) 10. Results are shown for mixtures with the mole fraction of species 2 equal to 0.2 (squares), 0.5 (triangles), and 0.8 (diamonds). The solid line indicates the predicted result using ^i(θ1,θ2) ) Ds,i(0)(1 - θ1 θ2).
Since, as shown above, the diagonal Maxwell-Stefan diffusivities both have the same loading dependence for this lattice model, this approximation reduces to
^12(θ1,θ2) ) x^1(0) ^2(0)(1 - θ1 - θ2)
(10)
The second approximation for ^12(θ1,θ2) is the logarithmic interpolation formula of Krishna and Wesselingh,11,23
^12(θ1,θ2) ) ^1(θ1,θ2)θ1/(θ1+θ2)^2(θ1,θ2)θ2/(θ1+θ2)
(11)
Again, substituting the known diagonal Maxwell-Stefan diffusivities leads to a simpler expression:
^12(θ1,θ2) ) ^1(0)θ1/(θ1+θ2)^2(0)θ2/(θ1+θ2) (1 - θ1 - θ2) (12) We first note that when the adsorbed mixture is an equimolar mixture, eqs 9 and 11 are identical. Approximating the counterexchange coefficient using these expressions and subsequently determining the Onsager coefficients and Fickian diffusivities using eqs 7 and 5 lead to the solid lines shown in Figures 1 and 2. The predicted diagonal components of the Onsager and Fickian matrices for the slower species are very accurate. The diagonal components of these matrices for the faster species and the off-diagonal elements of the matrices are somewhat less accurate, particularly for moderate and high fractional loadings. Extending the same predictions to a broader range of mixture compositions gives the results shown in Figure 4 for mixtures where species 2 hops 10 times more slowly than species 1. In this figure, D11 and D12 have been normalized by dividing by D1(0), the single-component Fickian diffusivity of species 1. D21 and D22 have similarly been normalized using D2(0). Using eq 9 to approximate the counterexchange coefficient leads to the predictions shown in these figures as dashed curves, while the solid curves show the results that come from using eq 11. Neither of the two approximations is clearly better than the other over all loadings and compositions. Both (23) Krishna, R.; Wesselingh, J. A. Chem. Eng. Sci. 1997, 52, 861.
approximations are reasonably accurate for D22, the diagonal component of the diffusivity matrix corresponding to the slower diffusing species. The predictions of eqs 9 and 11 are shown as a function of the ratio of the hopping rates of the two species in Figure 5. The data in this figure were normalized in the same way as that in Figure 4. When the slower species (species 2) hops 10 times slower than the faster species (species 1), both D22 and D12 are accurately predicted for all compositions when the total loading is low. Under the same conditions, D11 and D21 are predicted quite accurately for mixtures rich in species 1, but both are underestimated for mixtures rich in species 2. As might be expected, this approximate method is less accurate when the total loading is high. In particular, D11 and D12 (D21 and D22) are overestimated (underestimated) for mixtures of all compositions when the total loading is high. One interesting feature of the results in Figures 4 and 5 is that the accuracy of the predicted values is rather different from that reported by Paschek and Krishna for a binary mixture diffusing on a cubic lattice with a similar ratio of hopping rates.10 In this case, the diagonal diffusivity for the faster diffusing species was underestimated for all mixture compositions when the total loading was 0.48. By contrast, the same quantity is overestimated for simulations of our silicalite lattice model with total loadings of 0.50 (see Figure 1). Also, all elements of the transport diffusion matrix were accurately predicted for results from the cubic lattice with the total loading 0.96. In contrast, at the highest total loading we simulated, 0.8, the only element of the predicted matrix that has an accuracy similar to that seen for the cubic lattice is D12. These quantitative differences between the cubic and silicalite lattices presumably stem from the different correlation effects that are observed on these two lattices.10,13,14 These observations reinforce the conclusions of previous investigators that extending quantitative arguments regarding diffusion from cubic lattices to lattices with reduced connectivity such as silicalite must be done with caution.13,14 4. Binary Self-Diffusivities We denote the self-diffusivity of species i in a lattice with fractional loadings θ1 and θ2 of species 1 and 2 by Ds,i(θ1,θ2). These binary self-diffusivities have been computed for the lattice model of silicalite, as described in section 2. Self-diffusivities for binary mixtures where ν1/ ν2 ) 10 are shown in Figure 6 as a function of total fractional loading and the composition of the adsorbed mixture. The diffusivities in Figure 6 for species i have been normalized by the infinite dilution single-component diffusivity of that species. The main features of Figure 6 are similar to those that have been observed by Paschek and Krishna in their lattice model for methane-perfluoromethane diffusion in silicalite.11 Briefly, at fixed composition, the self-diffusivities of both species decrease as the fractional loading is increased. This also occurs during single-component diffusion and is a result of the steric hindrance of diffusing particles due to other particles in the lattice. At fixed total loadings, the self-diffusivities of both species decrease as the mole fraction of the slower species, species 2, is increased. This so-called acceleration/ deceleration effect has been carefully explained by Paschek and Krishna.11 The fact that this effect acts differently on the two species is made clear by noting that the normalized binary self-diffusivities for the two species are not equal for any total loading or mixture composition. The differences between the two species are particularly marked
Analysis of a Lattice Model for Silicalite
Langmuir, Vol. 18, No. 20, 2002 7397
Figure 4. Normalized binary Fickian diffusivities in the silicalite lattice model as a function of mixture composition and loading with ν1/ν2 ) 10. Each component of the Fickian diffusion matrix is shown: (a) D11, (b) D12, (c) D21, and (d) D22. In each plot, the symbols show results from MC simulations at total loadings of 10% (squares), 40% (triangles), and 80% (diamonds). Approximations based on eqs 9 and 11 are shown as dashed and solid curves, respectively.
for moderate total loadings (see the open and filled triangles in Figure 6a). The result of varying the ratio of hopping rates for species 1 and species 2 is shown in Figure 7, which shows the normalized binary self-diffusivities as a function of ν1/ν2 for simulations with a fractional loading of 0.4. When ν1/ν2 ) 1, the self-diffusivities are of course independent of the composition of the mixture, as shown by the squares in Figure 7. As the ratio of hopping rates is increased, the acceleration/deceleration effects mentioned above become more pronounced. Our main focus here is to examine the accuracy of several approximations that have been suggested for predicting binary self-diffusivities from single-component data. We have examined three approximate methods for estimating binary self-diffusivities. The first approximate method is a logarithmic interpolation formula,11
Ds,i(θ1,θ2) ) Ds,i(0,θ1+θ2)θ1/(θ1+θ2)Ds,i(θ1+θ2,0)θ2/(θ1+θ2) (13) This formula interpolates the binary self-diffusivity between the single-component self-diffusivity of component i at the total loading of the binary mixture and the infinite dilution self-diffusivity of component i in a system where the total loading of the other component is θ1 + θ2. In terms of the mole fraction of species 2 in the adsorbed
mixture, x2 ) θ2/(θ1 + θ2), this expression can be written as
Ds,i(θ1,θ2) ) Ds,i(0,θ1+θ2)
(
)
Ds,i(θ1+θ2,0)
Ds,i(0,θ1+θ2)
x2
(14)
The ability of this interpolation formula to fit the actual binary self-diffusivities for our lattice model is shown in Figures 6a and 7a. As these figures show, eq 14 provides a quantitatively accurate method to correlate the binary self-diffusivities. This has been observed previously by Paschek and Krishna for a lattice model of methaneperfluoromethane adsorbed in silicalite11 and by Snurr and Ka¨rger for atomistic simulations of the same mixture.5 The other two approximations to the binary selfdiffusivities are based on rewriting the self-diffusivities in terms of the single-component Maxwell-Stefan diffusivities, ^i(θ1,θ2), and the counterexchange coefficient, ^12(θ1,θ2). Krishna and co-workers have argued that the binary self-diffusivities can be written as11
Ds,i(θ1,θ2) ) (1 + θ1^2/^12 + θ2^1/^12)-1^i (15) where all the Maxwell-Stefan terms on the right-hand side are evaluated in the binary mixture at loading (θ1,θ2). To use this expression in practice, the single-component Maxwell-Stefan diffusivities and the counterexchange coefficient, ^12(θ1,θ2), must be estimated. Although it has been claimed that this expression is exact for lattice models
7398
Langmuir, Vol. 18, No. 20, 2002
Maceiras and Sholl
Figure 5. Normalized binary Fickian diffusivities in the silicalite lattice model as a function of mixture composition and loading with ν1/ν2 ) 10 (diamonds), 5 (triangles), 2 (circles), and 1 (squares) with the total loading 40%. Each component of the Fickian diffusion matrix is shown: (a) D11, (b) D12, (c) D21, and (d) D22. Approximations based on eqs 9 and 11 are shown as dashed and solid curves, respectively.
that exactly satisfy Langmuir isotherms,11,18 a simple argument shows that this claim is not correct. If eq 15 is exact, then
Ds,1(θ1,θ2)/Ds,2(θ1,θ2) ) ^1(θ1,θ2)/^2(θ1,θ2) (16) Since we showed above that for the lattice model we are considering ^i(θ1,θ2) ) Ds,i(0)(1 - θ1 - θ2), then a direct implication of eq 15 is that
Ds,1(θ1,θ2)/Ds,1(0) ) Ds,2(θ1,θ2)/Ds,2(0)
(17)
That is, when the binary self-diffusivities are plotted in the normalized form used in Figures 6 and 7, the normalized self-diffusivities of the two species are predicted by eq 15 to be equal. It is clear from our numerical data that this is not the case. Another way to state this result is that even if the diagonal Maxwell-Stefan diffusivities, ^i(θ1,θ2), and the Maxwell-Stefan counterexchange coefficient, ^12(θ1,θ2), are known exactly, eq 15 does not correctly predict the binary self-diffusivities. We return to this point below. Despite the fact that eq 15 is not exact, it is interesting to examine the ability of this expression to predict the binary self-diffusivities using the known single-component Maxwell-Stefan diffusivities, ^i(θ1,θ2) ) Ds,i(0)(1 - θ1 θ2), and the approximations described in section 3 for the
counterexchange coefficient, eqs 9 and 11. The binary selfdiffusivities predicted for our lattice model using eq 9 are compared to the actual binary self-diffusivities in Figures 6b and 7b. Since eq 17 applies regardless of the form of ^12(θ1,θ2), the predictions using eq 9 for each species are identical in the normalized form shown in these figures. Examining Figure 6b shows that this approximation systematically underestimates the self-diffusivity of the slower species and overestimates the self-diffusivity of the faster species except for mixtures rich in the slower species at low total loadings. It can be seen from Figure 7b that this approximation is not exact, even when the hopping rates of the two species are identical. In this case, the self-diffusivity is overestimated. Our results are qualitatively similar to those of Paschek and Krishna for their lattice model of methane-perfluoromethane diffusion in silicalite at 200 K11 in that using eq 9 generally underestimates (overestimates) the self-diffusivity of the slower (faster) diffusing species. One exception to this observation occurs in the results of Paschek and Krishna for high fractional loadings, where eq 9 overestimates the self-diffusivities of both species. Using eq 11 to define the counterexchange coefficient leads to the predicted binary self-diffusivities, shown in Figures 6c and 7c. Again, the two species are predicted to have equal normalized self-diffusivities as a consequence of eq 17. Figure 6c shows that using eq 11 gives
Analysis of a Lattice Model for Silicalite
Figure 6. Normalized binary self-diffusivities for the silicalite lattice model with ν1/ν2 ) 10. Filled (open) symbols show the results from MC simulations for species 1 (species 2) with the total fractional loadings 0.1 (squares), 0.4 (triangles), and 0.8 (diamonds). The solid (dashed) curves show the predictions for species 1 (species 2) from (a) eq 14, (b) eq 9, and (c) eq 11.
results that are quite accurate for the self-diffusion of the slower species as the mole fraction of that species approaches zero when the slower species hops 10 times slower than the other species. This agreement does not extend, however, to mixtures where the hopping rates are more similar, as can be seen from Figure 7c. As with the predictions based on eq 9, the self-diffusivities are overestimated for the apparently simple case in which both species have the same hopping rate. When the difference between the hopping rates of the two species is large, as in Figure 6c, the decrease in the self-diffusivities as a function of the mole fraction of the slower species is substantially overestimated. This overestimation is particularly severe for low total loadings. Our results are qualitatively similar to those of Paschek and Krishna for their lattice model of methane-perfluoromethane diffusion in silicalite at 200 K.11
Langmuir, Vol. 18, No. 20, 2002 7399
Figure 7. Same as Figure 6 but for the silicalite lattice model at the total fractional loading 0.4 with ν1/ν2 ) 1 (squares), 2 (circles), 5 (triangles), and 10 (diamonds).
It is interesting to ask whether the discrepancies between the predicted binary self-diffusivities discussed above are due mainly to inaccuracies in the approximated forms of the counterexchange coefficient or to the fact that eq 15 is not exact. To address this issue, we calculated the counterexchange coefficient directly from the binary Onsager transport coefficients (using eqs 7 and 8) and used the resulting values of ^12(θ1,θ2) in eq 15. In this calculation, any discrepancy between the predicted selfdiffusivities and the actual result is due only to the inexactness of eq 15. The results of using eq 15 in this way are shown in Figure 8. The fact that the normalized selfdiffusivities predicted by eq 15 should be equal for the two species (cf. eq 17) provides a good indication of the numerical uncertainties in Figure 8. Predicting the binary self-diffusivities in this way yields results that are arguably more accurate than comparable results determined by approximating ^12(θ1,θ2) using eq 9 or 11.
7400
Langmuir, Vol. 18, No. 20, 2002
Figure 8. Same as Figure 6 except, in this case, the solid (dashed) lines show the binary self-diffusivities for species 1 (2) determined from MC simulations and the filled (open) symbols show the prediction of eq 15 based on computing all components of the Maxwell-Stefan diffusivities directly from the Onsager transport coefficients extracted from the same MC simulations. Results are shown with total fractional loadings 0.1 (squares), 0.4 (triangles), and 0.8 (diamonds).
Nonetheless, it is clear from Figure 8 that predicting binary self-diffusivities on the basis of eq 15 cannot yield results that are quantitatively accurate for both species simultaneously, regardless of how accurately ^12(θ1,θ2) can be predicted. 5. Discussion We have used dynamic Monte Carlo simulations to measure the binary transport diffusivities and selfdiffusivities of adsorbate mixtures in a simple lattice model of silicalite. In this model, the saturation loading of each species is the same and the binding energies of all lattice sites are identical, so adsorption equilibrium for this model is exactly described by a Langmuir isotherm. Our main conclusions can be summarized as follows: 1. The diagonal Maxwell-Stefan diffusivities for binary mixtures in this model appear to exactly follow the meanfield expression ^i(θ1,θ2) ) Ds,i(0)(1 - θ1 - θ2). This expression has previously been assumed to hold,11,18 but
Maceiras and Sholl
we know of no previous numerical data assessing this assumption. 2. The two approximations we have tested for the Maxwell-Stefan counterexchange coefficient, eqs 9 and 11, both give reasonable predictions for some ranges of the model parameters, but neither is quantitatively accurate in all instances. We have presented data that display the variation in the accuracies of these approximate methods as a function of the hopping rates of the two species. 3. As has been observed previously for both lattice11 and atomistic simulations,5 a logarithmic interpolation formula, eq 13, can be used to accurately fit the binary self-diffusivities under all conditions we have examined. Use of this formula requires knowledge of the loadingdependent self-diffusivity of each species and the infinite dilution self-diffusivity of each species as a function of the other species’ loading. 4. Our numerical data shows that eq 15, which relates the binary self-diffusivities to the binary Maxwell-Stefan diffusivities, is not an exact expression, even when the binary Maxwell-Stefan diffusivities are known exactly. Several previous efforts to develop approximate relations between single-component and binary diffusivities in lattice models that obeyed Langmuir isotherms have assumed eq 15 is exact.11,18 Using eq 15 in conjunction with the approximations to the Maxwell-Stefan counterexchange coefficient defined in eq 9 or 11 leads to predicted binary self-diffusivities that are in qualitative but not quantitative agreement with the actual selfdiffusivities. Our analysis shows that the inexactness of eq 15 plays a significant role in these quantitative discrepancies. Acknowledgment. This work was supported by a National Science Foundation CAREER award (CTS9983647). D.B.M was supported by Fundacio´n Pedro Barrie´ de la Maza. D.S.S. gratefully acknowledges an Alfred P. Sloan fellowship and a Camille Dreyfus Teacher Scholar Award. Computing resources were provided by the Pittsburgh Supercomputer Center and through a computer cluster in our department supported by the National Science Foundation (CTS-0094407) and Intel. LA025972U
