Testing Predictions of Macroscopic Binary Diffusion

Langmuir 2006, 22, 3707-3714

3707

Testing Predictions of Macroscopic Binary Diffusion Coefficients Using Lattice Models with Site Heterogeneity David S. Sholl* Department of Chemical Engineering, Carnegie Mellon UniVersity, Pittsburgh, PennsylVania 15213, and DiVision of Chemical Engineering, UniVersity of Queensland, St. Lucia, QLD, 4072, Australia ReceiVed December 15, 2005. In Final Form: February 16, 2006

Quantitatively predicting mass transport rates for chemical mixtures in porous materials is important in applications of materials such as adsorbents, membranes, and catalysts. Because directly assessing mixture transport experimentally is challenging, theoretical models that can predict mixture diffusion coefficients using only single-component information would have many uses. One such model was proposed by Skoulidas, Sholl, and Krishna (Langmuir, 2003, 19, 7977), and applications of this model to a variety of chemical mixtures in nanoporous materials have yielded promising results. In this paper, the accuracy of this model for predicting mixture diffusion coefficients in materials that exhibit a heterogeneous distribution of local binding energies is examined. To examine this issue, single-component and binary mixture diffusion coefficients are computed using kinetic Monte Carlo for a two-dimensional lattice model over a wide range of lattice occupancies and compositions. The approach suggested by Skoulidas, Sholl, and Krishna is found to be accurate in situations where the spatial distribution of binding site energies is relatively homogeneous, but is considerably less accurate for strongly heterogeneous energy distributions.

Introduction The diffusion of molecules in porous materials plays a crucial role in the performance of adsorbers, membranes, and other applications of these materials.1-5 In essentially all practical applications, adsorbed molecules are present in chemical mixtures, so the core task of describing diffusion in these applications is one of multicomponent diffusion. Unfortunately, accurate experimental measurement of multicomponent diffusion in microporous materials such as zeolites, activated carbons, and so on is difficult. Models that can accurately predict the multicomponent diffusivities associated with net mass transport in these materials from single-component data would be of considerable utility. In the past several years, atomistic simulations of molecular diffusion in nanoporous materials have begun to produce data that can complement experimental approaches to describing multicomponent diffusion in ordered microporous materials. Although molecular dynamics (MD) simulations have been used for some time to examine self-diffusion in adsorbed chemical mixtures,2,6-9 the use of MD to obtain information related to macroscopic diffusion is more recent. MD simulations of macroscopic mixture diffusion have now been performed for a variety of chemical mixtures in various all-silica zeolites4,10-17 and in single-walled carbon nanotubes.18-20 * E-mail: [email protected]. (1) Ka¨rger, J.; Ruthven, D. Diffusion in Zeolites and Other Microporous Materials; John Wiley & Sons: New York, 1992. (2) Keil, F. J.; Krishna, R.; Coppens, M. O. ReV. Chem. Eng. 2000, 16, 71. (3) Yang, R. T. Gas Separation by Adsorption Processes; Butterworth: Boston, MA, 1987. (4) Skoulidas, A. I.; Sholl, D. S. AIChE J. 2005, 51, 867. (5) Auerbach, S. M. Int. ReV. Phys. Chem. 2000, 19, 155. (6) Snurr, R. Q.; Ka¨rger, J. J. Phys. Chem. B 1997, 101, 6469. (7) Jost, S.; Ba¨r, N.; Fritzsche, S.; Haberlandt, R.; Ka¨rger, J. J. Phys. Chem. B 1998, 102, 6375. (8) Gergidis, L. N.; Theodorou, D. N. J. Phys. Chem. B 1999, 103, 3380. (9) Sholl, D. S.; Fichthorn, K. A. J. Chem. Phys. 1997, 107, 4384. (10) Sanborn, M. J.; Snurr, R. Q. Sep. Purif. Technol. 2000, 20, 1. (11) Sanborn, M. J.; Snurr, R. Q. AIChE J. 2001, 47, 2032. (12) Skoulidas, A. I.; Sholl, D. S.; Bowen, T. C.; Doelling, C.; Falconer, J. L.; Noble, R. D. J. Membr. Sci. 2003, 227, 123.

The availability of MD data for mixture diffusion has been useful in exploring models for predicting multicomponent diffusivities from single-component information. For example, Skoulidas, Sholl, and Krishna used MD data for the diffusion of CH4/CF4 mixtures in silicalite to introduce a mixing theory that accurately predicted the macroscopic diffusivities of this mixture using the single-component diffusivities for the same species, which was also measured using MD.13 This mixing theory has subsequently been applied to a range of adsorbed mixtures, including hexane and butane isomers in all-silica mordenite,14 linear alkanes in silicalite15 and all-silica faujasite,16,17 and a variety of gases adsorbed in single-walled carbon nanotubes.20 In all cases, including the small number of ternary and quaternary mixtures that have been examined,16 the predictions of the mixing theory are quantitatively accurate. The results from the mixing theory mentioned above are very promising. The characteristics of the adsorbed mixtures that have been examined with this theory, however, do not span the full range of properties that can occur in practical materials. All of the examples listed above involve crystalline, nanoporous adsorbents that are homogeneous in the sense that the chemical structure of the pore walls is uniform. Unlike these examples, many practical adsorbents have chemically inhomogeneous pore walls. In zeolites, for example, this inhomogeneity arises when framework substitutions and charge-balancing cations are present.1,2,5 In porous protein crystals, this inhomogeneity may arise from the varied functional groups exposed by the protein forming the pore wall.21 Complications can also arise in real (13) Skoulidas, A. I.; Sholl, D. S.; Krishna, R. Langmuir 2003, 19, 7977. (14) van Baten, J. M.; Krishna, R. Microporous Mesoporous Mater. 2005, 84, 179. (15) Krishna, R.; van Baten, J. M. Chem. Phys. Lett. 2005, 407, 159. (16) Krishna, R.; van Baten, J. M. J. Phys. Chem. B 2005, 109, 6386. (17) Chempath, S.; Krishna, R.; Snurr, R. Q. J. Phys. Chem. B 2004, 108, 13481. (18) Chen, H.; Sholl, D. S. J. Am. Chem. Soc. 2004, 126, 7778. (19) Chen, H.; Sholl, D. S. J. Membr. Sci. 2006, 269, 152. (20) Krishna, R.; van Baten, J. M. Ind. Eng. Chem. Res., submitted for publication, 2006. (21) Malek, K.; Odijk, T.; Coppens, M. O. Ind. Eng. Chem. Res. 2006, 45, 2084.

10.1021/la053405b CCC: $33.50 © 2006 American Chemical Society Published on Web 03/15/2006

3708 Langmuir, Vol. 22, No. 8, 2006

materials because of spatial inhomogeneity in the shape of the pores within a material,22 as in disordered materials such as activated carbons, or from nonisotropic interactions between adsorbed molecules arising from, for example, hydrogen bonding. The aim of this paper is to explore the ability of the mixing theory introduced by Skoulidas, Sholl, and Krishna to describe diffusion in materials with energetically inhomogeneous binding sites for adsorbed molecules. Atomically detailed simulations of molecular adsorption and diffusion in cationic zeolites are possible,2,5,23-25 but they are more time-consuming than analogous simulations of all-silica materials. Since my aim here is to test the general applicability of a mixing theory for multicomponent diffusion, it is useful to adopt an approach that can rapidly generate well-defined data for mixture diffusion, namely, lattice models. Lattice models have a long history in the phenomenological study of diffusion.26-35 Of special relevance to this work, lattice models have been used to show that the characteristics of molecular diffusion in porous materials that include localized binding sites that are strongly favored for adsorption are considerably different than those in models where the binding energies of all the lattice sites are homogeneous.36-39 In some instances, it is possible to rigorously connect lattice models with diffusion in real materials.40-44 My intent here, however, is to use lattice models merely as a means for examining the influence of site energy heterogeneity on mixture diffusion in a phenomenological manner. It is useful to state the main conclusion of this work before launching into a detailed description of my results: the mixing theory of Skoulidas, Sholl, and Krishna performs less accurately for models with strong energetic heterogeneity than it does for models that are relatively homogeneous. A judgment about how accurate a theory like this needs to be is of course dependent on the specific context in which it will be used. The results below, however, indicate that situations may exist in real materials where improvement in the accuracy of the current mixing theory would be highly desirable. In addition to providing a test for this specific mixing theory, the lattice model results that are described in the remainder of this paper should also be useful in assisting the development of improved mixing theories in the future. (22) Malek, K.; Coppens, M.-O. Phys. ReV. Lett. 2001, 87, 125505. (23) Goj, A.; Sholl, D. S.; Akten, E. D.; Kohen, D. J. Phys. Chem. B 2002, 106, 8367. (24) Beerdsen, E.; Smit, B.; Calero, S. J. Phys. Chem. B 2002, 106, 10659. (25) Beerdsen, E.; Dubbeldam, D.; Smit, B.; Vlugt, T. J. H.; Calero, S. J. Phys. Chem. B 2003, 107, 12088. (26) Mak, C. H.; Andersen, H. C.; George, S. M. J. Chem. Phys. 1988, 88, 4052. (27) Uebing, C.; Gomer, R. J. Chem. Phys. 1994, 100, 7759. (28) Paschek, D.; Krishna, R. Phys. Chem. Chem. Phys. 2000, 2, 2389. (29) Paschek, D.; Krishna, R. Phys. Chem. Chem. Phys. 2001, 3, 3185. (30) Paschek, D.; Krishna, R. Langmuir 2001, 17, 247. (31) Sholl, D. S.; Fichthorn, K. A. Phys. ReV. E 1997, 55, 7753. (32) Sholl, D. S. Chem. Eng. J. 1999, 74, 25. (33) Ka¨rger, J.; Demontis, P.; Suffritti, G. B.; Tilocca, A. J. Chem. Phys. 1999, 110, 1163. (34) Demontis, P.; Suffritti, G. B.; Tilocca, A. J. Chem. Phys. 2000, 113, 7588. (35) Maceiras, D. B.; Sholl, D. S. Langmuir 2002, 18, 7393. (36) Trout, B. L.; Chakraborty, A. K.; Bell, A. T. Chem. Eng. Sci. 1997, 52, 2265. (37) Coppens, M. O.; Bell, A. T.; Chakraborty, A. K. Chem. Eng. Sci. 1998, 53, 2053. (38) Coppens, M. O.; Bell, A. T.; Chakraborty, A. K. Chem. Eng. Sci. 1999, 54, 3455. (39) Coppens, M.-O.; Iyengar, V. Nanotechnology 2005, 16, S442. (40) Kamakoti, P.; Morreale, B. D.; Ciocco, M. V.; Howard, B. H.; Killmeyer, R. P.; Cugini, A. V.; Sholl, D. S. Science 2005, 307, 569. (41) Tunca, C.; Ford, D. M. J. Chem. Phys. 2004, 120, 10763. (42) Sholl, D. S.; Skodje, R. T. Phys. ReV. Lett. 1995, 75, 3158. (43) Uberuaga, B. P.; Voter, A. F.; Sieber, K. K.; Sholl, D. S. Phys. ReV. Lett. 2003, 91, 106901. (44) Krishna, R.; van Baten, J. M. Chem. Eng. Technol. 2005, 28, 160.

Sholl

Multicomponent Diffusion The results below are for lattice models that phenomenologically describe the adsorption and diffusion of binary chemical mixtures on a structured adsorbent. The adsorbent is defined by a series of discrete sites, with molecules of both chemicals adsorbing on the same set of sites. The occupancy of any site is defined by the normalized probability of finding a particle in that site. The net occupancy of an adsorbed mixture is denoted as B θ ) (θ1,θ2), and the total occupancy of the same mixture is θ ) θ1 + θ2. The adsorption isotherm for these models is defined by θi ) θi(p1,p2), where pi is the partial pressure of species i in the gas phase in equilibrium with the adsorbent. The various mathematical formalisms available for describing multicomponent diffusion have been described in detail elsewhere.1,2,5,30,45-47 Here, I summarize the definitions and relationships between the Onsager, Fick, and Maxwell-Stefan (MS) diffusivities for the diffusion of a binary mixture in a form that is correct when both species have identical saturation occupancies. Several minor modifications of the expressions below are necessary in cases where the saturation occupancies of the two species are not equal.13 Perhaps the most concise definition of binary diffusive transport is to use the Onsager formalism, which relates the species fluxes, N B i, to gradients in chemical potential:

(N) ) -[L](∇µ)

(1)

Here, the Onsager transport coefficients, Lij, define a symmetric matrix. An alternative representation that uses gradients in concentration is the Fickian formalism:

(N) ) -[D](∇θ)

(2)

The binary Fickian (or transport) diffusion coefficients, Dij, are not symmetric. [L] and [D] are related without approximation by

Dii )

[

]

∂ ln pj kT ∂ ln pi L + Lij ; θi ii∂ ln θi ∂ ln θi Dij )

[

]

∂ ln pj kT ∂ ln pi L + Lij (3) θj ii∂ ln θj ∂ ln θj

where j * i. An equally useful approach to defining binary diffusive transport is to use the MS description. With this formalism,

(N) ) -[B]-1[Γ](∇θ)

(4)

where

Bii )

θj θi 1 + corr; Bij ) - corr, j * i }i } }

(5)

θi ∂ ln pi θj ∂ ln θj

(6)

ij

ij

and

Γij )

Comparing eqs 2 and 4 shows that the Fickian and MS diffusivities (45) Krishna, R.; van den Broeke, L. J. P. Chem. Eng. J. 1995, 57, 155. (46) Wesselingh, J. A.; Krishna, R. Mass Transfer in Multicomponent Mixtures; Delft University Press: Delft, The Netherlands, 2000. (47) Cussler, E. L. Multicomponent Diffusion; Elsevier Scientific Publishing Co.: Amsterdam/New York, 1976.

Testing Binary Diffusion Coefficient Predictions

Langmuir, Vol. 22, No. 8, 2006 3709

are related without approximation by

[D] ) [B]-1[Γ]

(7)

A Lattice Model With Site Heterogeneity I have used a lattice model originally suggested by Mak, George, and Andersen26 as a simple example of site heterogeneity. Particles of two species are allowed to diffuse via hops between neighboring sites on a two-dimensional square lattice with periodic boundary conditions with the distance between adjacent sites defined as unity. The fraction of sites filled with species i is denoted θi. A random fraction, ft, of the lattice sites is defined to be trap sites at which particles of one (or both) species bind more favorably than at all other lattice sites. The energy of a particle of species i in a trap (nontrap) site is defined as i,t (i,o) with the convention that these energies are negative. Hops are only allowed into neighboring sites that are unoccupied by other particles. Hops occur for species i at a rate defined by exp(i/kT), where i is the binding energy in the currently occupied site. This hopping rate can be thought of as occurring in a situation where the hopping rate is defined by transition state theory, and the transition states between all pairs of sites have the same energy for both species. Although this definition is somewhat artificial, it captures the idea that more tightly bound species hop away from their binding sites more slowly than weakly bound species. An advantage of this definition is that the single-component diffusion coefficient of species i in the limit θi f 0 can be found in closed form:26

D(θi f 0) ) [ft exp(i,t/kT) + (1 - ft) exp(i,o/kT)]-1

(8)

Mak, George, and Andersen used MC simulations to explore macroscopic single-component diffusion in this lattice model for various trap densities and energies.26 They also examined situations in which the binding site energies formed a continuous distribution. The general features of diffusion in this case can be described in a satisfactory way by considering models based on the simple two-site energy distribution. This model does not include the possibility that the hopping rates for individual particles can be influenced by other nearby particles, although it is certainly possible to efficiently incorporate these effects into kinetic Monte Carlo (KMC) simulations of diffusion in lattice models.44 Single-component and binary diffusion coefficients were obtained for the lattice model defined above using methods analogous to those used in MD simulations of diffusion in nanoporous materials.10,11,13,48,49 A simple KMC scheme was used to allow realizations of each lattice model to evolve in time while defining the correct local hopping rates.26 This method defines the position of each molecule of species i as a function of time, b ri(t). The Onsager transport coefficients for binary mixtures diffusing on an L × L lattice were obtained using10,11,13

Lij ) lim tf∞

1 2

4L kTt

〈[

∑i ∆br i(t)]‚[∑j ∆br j(t)]〉

(9)

where t is time, ∆r bi(t) ) b ri(t) - b ri(0), the sums are over all particles of species i and j, and the angled brackets indicate an average over a large number of independent realizations. For most of the results discussed below, data for this process was averaged over 5000 independent simulations using a 30 × 30 (48) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2001, 105, 3151. (49) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2002, 106, 5058.

lattice of length t ) 4000 after equilibrating an initially random distribution of particles for t ) 400. A small number of simulations for low values of θ1 + θ2 were performed with 100 × 100 lattices. Uncertainties in the computed transport coefficients were estimated by analyzing the N independent trajectories in 10 groups of N/10 trajectories and computing the standard deviation of the resulting values for Lij. Numerical tests were performed to ensure that the computed values of the transport coefficients were not significantly affected by the specific choices listed above for simulation size and length. For each example that is described below, binary diffusion was simulated at all coverages where θ1 and θ2 were multiples of 0.1, and 0 < θ1 + θ2 e 0.9. For single-component diffusion, the expression above for the Onsager transport coefficient is related in a simple way to the single-component corrected diffusivity:10,11,13

D0,i )

kT L θ ii

(10)

During the KMC simulations of single-component diffusion, the self-diffusion coefficients were also computed, using

Ds,i ) lim tf∞

1

∆b r i2(t)〉 ∑ 4t i 〈

(11)

This expression was also used to calculate the self-diffusivities of each species during KMC simulations of mixture diffusion. Once the matrix of binary Onsager transport coefficients, [L], was determined for a lattice model, the matrix of Fickian diffusivities, [D], was defined using eq 3. The binary adsorption isotherm for the lattice model defined above can be expressed in closed form:

θi ) (1 - ft) [1 +

K exp(-i,o/kT)pi

+

∑ K exp(-j,o/kT)pj]

j)1,2

K exp(-i,t/kT)pi ft [1 +

∑

(12) K exp(-j,t/kT)pj]

j)1,2

Here, K is a constant that involves the impingement rate of molecules on each lattice site. Without loss of generality, K is set to be unity for all calculations below. From this expression, derivatives of the form ∂ ln θi/∂ ln pj can be evaluated analytically. To evaluate the derivatives of the form ∂ ln pi/∂ ln θj needed in eq 3, I solved eq 12 numerically to determine the pressures p1 and p2 associated with a particular B θ ) (θ1,θ2) of interest and then evaluated the desired derivatives using the set of derivatives ∂ ln θi/∂ ln pj at those conditions.50

The SSK Method For Predicting Binary Diffusion The method that I will test below for predicting binary diffusivities from single-component information is the approach described by Skoulidas, Sholl, and Krishna;13 for brevity I will denote this approach the SSK method. This method assumes that information is available on the single-component selfand corrected diffusivities of both species of interest. These two single-component quantities are used to define an (50) Chen, H.; Sholl, D. S. Langmuir 2006, 22, 709.

3710 Langmuir, Vol. 22, No. 8, 2006

Sholl

Table 1. Parameters Defining the Four Binary Mixtures Discussed in the Text example #

ft

1,o/kT

1 2 3 4

0 0.15 0.40 0.15

0 0 0 0

1,t/kT

2,o/kT

2,t/kT

-1.0 -1.0 -5.0

-ln(5) -ln(5) -ln(5) -ln(5)

-ln(5) - 3.0 -ln(5) - 3.0 -ln(5)

auxiliary quantity, }corr ii (θ), via

1 1 θ ) + Ds,i(θ) Do,i(θ) }corr ii (θ)

(13)

To define the mixture MS diffusivities (see eq 5), it is necessary to estimate }i(θ1,θ2) and }corr ij (θ1,θ2) for j * i. The first quantity is estimated by assuming that

}i(θ1,θ2) = D0,i(θ1 + θ2)

(14)

That is, the MS diffusivity in the mixture is replaced by the single-component corrected diffusivity at the same total occupancy as in the mixture. The mixture correlation diffusivities are estimated using the interpolation formula θ1/θ [}jj(θ1 + θ2)]θ2/θ (15) }corr ij (θ1,θ2) = [}ii(θ1 + θ2)]

where }ii(θ) and }jj(θ) are defined using eq 13. Several factors involving the saturation occupancies appear in this expression if the saturation occupancies of the individual species are different.13 Once the binary diffusivities have been estimated using eqs 14 and 15, the Onsager and Fickian descriptions of binary diffusion can be recovered using eqs 7 and 3. To make these conversions in a way that uses only single-component information, the mixture adsorption isotherm must be estimated from the single-component isotherms. This task can be accomplished, for example, with ideal adsorbed solution theory (IAST).13 To focus attention on the approximations made with the SSK method defined above, I chose to use the exact mixture adsorption isotherm, eq 12, for the lattice models of interest here rather than introduce another level of approximation by using IAST.

Results and Discussion To examine the effect of site heterogeneity on the accuracy of predictions of the SSK method, I have used four examples that have progressively stronger site heterogeneity. The parameters defining these four examples are listed in Table 1. Example 1 has no heterogeneity; it defines a uniform lattice on which the single-component diffusion of species 2 is 5 times slower than that of species 1 at any fixed concentration. The second and third examples are defined so that trap sites exist that favorably bind particles of both species, although species 2 is trapped more strongly than species 1. The only difference between these examples is that, in Example 2 (Example 3), 15% (40%) of the lattice sites are defined as trap sites. In the final example, a more heterogeneous situation is considered in which species 1 is very strongly bound at the trap sites that form 15% of the lattice, but species 2 is bound with equal strength on all sites. For each example, the SSK approximation was used to predict the matrix of binary Fickian diffusivities, [DSSK], for all occupancies where θ1 and θ2 were multiples of 0.1, and 0 < θ1 + θ2 e 0.9. This set of 36 separate occupancies spans the full range of binary compositions that can exist for these lattice models. One challenge with this comprehensive approach is to find numerical quantities that succinctly capture the accuracy of the

Figure 1. Single-component self-diffusivities (filled symbols) and corrected diffusivities (open symbols) for species 1 and 2 of Example 1. The solid and dashed curves show the fitting functions used in describing these diffusivities in predicting mixture diffusion coefficients.

approximate method. For each diffusion matrix element, the mean and standard deviation of DSSK ij /Dij were calculated with respect to the complete set of occupancies. These values give one estimate of the global accuracy of the SSK approach. The maximum and minimum values of DSSK ij /Dij are also reported to provide another way to describe this accuracy. Example 1: Homogeneous Lattice. The single-component self- and corrected diffusivities for both species in this initial example are shown in Figure 1. The corrected diffusion coefficient is known analytically for this example.26 For species 1, D0(θ) ) 1 - θ. The self-diffusion coefficient of species 1 is well described by Ds(θ) ) 1 - 1.572θ + 0.572θ2. This function, as with the other functions used below to fit single-component diffusion coefficients, obeys the known limits of the diffusion coefficient exactly for θ ) 0 and θ ) 1. By construction, the diffusion coefficients of species 2 in this model are one-fifth the values of species 1. Figure 1 shows error bars on both the selfand corrected diffusion coefficients obtained from MC simulations, although the error bars on Ds are smaller than the symbol sizes. For clarity, error bars on the matrix elements of [D] are not shown below; they are similar in size to those on the singlecomponent correct diffusivity. The numerical values of each element of [D] and the uncertainty in these matrix elements for every mixture occupancy examined are listed in Table S1 (Supporting Information). The Fickian diffusion matrix predicted by the SSK approximation for each occupancy that was considered is compared with the MC result in Figure 2. The overall statistics for this comparison are summarized in Table 2. The values of the matrix elements span almost 2 orders of magnitude when the entire range of occupancies is considered. It is clear from Figure 2 and Table 2 that the SSK approximation is very accurate in this case. Most of the predicted matrix elements lie within 5% of the MC values, and almost all lie with 10% of the MC values. Example 2: Heterogeneous Lattice With Moderate Site Trapping. This example is the first one in which the site energies of the lattice are inhomogeneous. As mentioned above (and defined in detail in Table 1), the trap sites bind both diffusing species with species 2 having a stronger preference for the 15% of the lattice covered with trap sites than species 1. The binary Fickian diffusivities computed using KMC for this model are listed in Table S2 (Supporting Information). The single-component diffusion coefficients for both species are shown in Figure 3. The magnitude of the diffusivities for


Langmuir, Vol. 22, No. 8, 2006 3711

Figure 2. A comparison of the mixture diffusion coefficients predicted by the SSK method, DSSK ij , with the actual mixture diffusion coefficients, Dij, for Example 1. The range of mixture compositions and occupancies used for this comparison is described in the text. Table 2. Statistics Summarizing the Comparison between Mixture Diffusion Coefficients Predicted by the SSK Method, DSSK ij , and the Actual Mixture Diffusion Coefficients, Dij, for Example 1a DSSK 11 /D11 DSSK 12 /D12 DSSK 21 /D21 DSSK 22 /D22 a

mean

max.

min.

0.961 ( 0.051 0.982 ( 0.035 1.023 ( 0.081 1.000 ( 0.024

1.017 at B θ ) (0.2,0.1) 1.041 at B θ ) (0.3,0.1) 1.289 at B θ ) (0.1,0.2) 1.048 at B θ ) (0.1,0.8)


Figure 4. The same as Figure 2, but showing mixture diffusion data for Example 2. Table 3. The Same as Table 2, but for Example 2 DSSK 11 /D11 DSSK 12 /D12 DSSK 21 /D21 DSSK 22 /D22

mean

max.

min.

0.884 ( 0.078 0.916 ( 0.057 1.198 ( 0.137 1.217 ( 0.148



small loadings as the trap sites that retard the diffusion of trapped particles become filled, forcing additional particles to diffuse among the nontrapping sites. The single-component diffusivities for species 2 in this example are well described by

Ds(θ) ) 0.2 × (0.2589 + 0.4470θ - 0.0971θ2 -

The uncertainty in the mean is one standard deviation.

4.6523θ3 + 6.9460θ4 - 2.9024θ5) (18) and

D0(θ) ) 0.2 × (0.2589 + 0.5113θ + 0.7184θ2 6.2212θ3 + 7.8604θ4 - 3.1276θ5) (19)

Figure 3. The same as Figure 1, but showing single-component diffusivities for Example 2.

species 1 is reduced relative to Example 1 by the presence of trap sites (cf. Figure 1), but the overall shape of the singlecomponent diffusivities is only weakly affected by the traps. The self- and corrected diffusion coefficients for species 1 are well described by

Ds(θ) ) 0.79508 - 1.1446θ + 0.23827θ2 + 0.11127θ3 (16) and

D0(θ) ) 0.79508 - 0.7370θ + 0.0581θ2

(17)

The effect of site heterogeneity is more apparent in the singlecomponent behavior of species 2. Both the self- and corrected diffusivities for species 2 increase with increasing loading for

The comparison between the predictions of the SSK approximation and the Fickian diffusion matrix elements computed using MC for this example is shown in Figure 4. Comparing Figures 2 and 4 indicates that the SSK approximation is less accurate for this example than for the homogeneous lattice of Example 1. Table 3 summarizes the accuracy of the SSK approximation for Example 2. The approximate method systematically underestimates D11 and D12 (by ∼10% on average and by ∼30% in the worst examples) while systematically overestimating D21 and D22 (by ∼20% on average and by more than 50% in the worst cases). Example 3: A Second Heterogeneous Lattice With Moderate Site Trapping. The third example has trap sites with the same energies as those in Example 2, but has 40% of the lattice covered by trap sites. Table S3 (Supporting Information) lists all of the binary Fickian diffusion coefficients for this example. Figure 5 shows the single-component diffusivities for this case. The single-component diffusivities for species 1 were fitted using

Ds(θ) ) 0.5927 - 0.8023θ + 0.0753θ2 + 0.11353θ3 (20) and

Ds(θ) ) 0.5927 - 0.4291θ - 0.3265θ2 + 0.1638θ3

(21)

3712 Langmuir, Vol. 22, No. 8, 2006

Sholl

Figure 5. The same as Figure 1, but showing single-component diffusivities for Example 3.

Figure 6. The same as Figure 2, but showing mixture diffusion data for Example 3.

Table 4. The Same as Table 2, but for Example 3 DSSK 11 /D11 DSSK 12 /D12 DSSK 21 /D21 DSSK 22 /D22

mean

max.

min.

0.839 ( 0.112 0.865 ( 0.091 1.21 ( 0.177 1.199 ( 0.136



For species 2, the analogous fitting functions were

Ds(θ) ) 0.2 × (0.1158 - 0.0037θ - 0.2455θ2 + 0.1334θ3) (22) and

D0(θ) ) 0.2 × (0.1158 + 0.0146θ - 0.0911θ2 - 0.0393θ3) (23) Table 4 and Figure 6 compile information comparing the SSK predictions for this example with the binary diffusivities computed using MC. The accuracy of the SSK predictions is similar to that seen above for Example 2. D11 and D12 are underestimated by the SSK approximation, with DSSK 11 /D11 < 1 for every coverage examined. D21 and D22 are overestimated by ∼20% on average, with errors >50% occurring for some coverages. At every coverage, DSSK 22 /D22 > 1. Example 4: Heterogeneous Lattice with Strong Site Trapping. In this final example, the effects of site heterogeneity are stronger than those in the other three examples above. The site energies were chosen so that species 1 is very strongly trapped on the 15% of trap sites present on the lattice, but species 2 has no energetic preference to reside at the trap sites. The parameters defining this model are listed in Table 1, and the binary Fickian diffusivities are listed in Table S4 (Supporting Information). The single-component diffusivities of each species are shown in Figure 7. The single-component diffusivities of species 2 are identical to those in Example 1 above. The single-component diffusivities of species 1 are strongly nonmonotonic as functions of coverage because of the strength of the trap sites for this species. At low coverages, almost all particles from species 1 reside in trap sites, so diffusion is dominated by the slow hopping out of these sites. At moderate coverages, the particles that are not trapped can diffuse relatively rapidly on the remaining nontrap sites, so the overall diffusion coefficients increase. At high coverages, the diffusion coefficients decrease as usual. Finding empirical functions to fit the single-component diffusivities of species 1 was more difficult than finding ones

Figure 7. The same as Figure 1 but showing single-component diffusivities for Example 4.

for the examples already discussed. The curves shown in Figure 7 have the form

{

D(θ) ) f1(θ) for θ e 0.25 2 2 (1 - x )f1(θ) + x f2(θ) for θ e 0.25 for 0.25 < θ < 0.35 f2(θ) for θ g 0.35 (24) with x ) (θ - 0.25)/0.1. For the self-diffusivity of species 1,

f1(θ) ) 0.0433 - 0.3792θ + 9.2261θ2 - 18.7784θ3

(25)

and

f2(θ) ) 0.2039 + 0.5574θ - 1.5671θ2 + 0.8059θ3

(26)

For the corrected diffusivity, the fitting functions were similar:

f1(θ) ) 0.0433 - 0.3131θ + 8.6917θ2 - 16.4668θ3 (27) and

f2(θ) ) 0.06455 + 1.3493θ - 2.2950θ2 + 0.8811θ3 (28) It can be seen in Figure 7 that these empirical functions are not entirely satisfactory for 0 < θ < 0.10; in particular, the slope of the diffusivities is described incorrectly in part of this region. In the comparisons below, however, only values of the single-


Langmuir, Vol. 22, No. 8, 2006 3713

Figure 8. The same as Figure 2, but showing mixture diffusion data for Example 4. Table 5. The Same as Table 2, but for Example 4 DSSK 11 /D11 DSSK 12 /D12 DSSK 21 /D21 DSSK 22 /D22

mean

max.

min.

2.897 ( 2.159 2.842 ( 2.144 1.423 ( 0.303 1.196 ( 0.054



Table 6. The Same as Table 5, except that Only Coverages from Example 4 with θ1 > 0.3 Are Included DSSK 11 /D11 DSSK 12 /D12 DSSK 21 /D21 DSSK 22 /D22

mean

max.

min.

1.353 ( 0.156 1.331 ( 0.148 1.189 ( 0.053 1.199 ( 0.028



component diffusivities with θ g 0.10 were used, so this shortcoming of the fitting functions plays no role in these comparisons. Figure 8 compares the predictions of the SSK approximation with the actual binary diffusivities for this example. It is clear from Figure 8 that the SSK approximation is considerably less accurate for this example than for the other examples presented above. All the diffusion matrix elements at every coverage examined are overestimated by the SSK approximation, with D11 and D12 being represented much less accurately than D21 and D22. Some statistics representing this comparison are listed in Table 5. D22 is predicted by the SSK approximation with moderate accuracy; on average, this matrix element is overestimated by ∼20%, and none of the values are overestimated by more than 30%. The predictions for D21 are somewhat less accurate, with an average error of ∼40% and maximum errors exceeding 100%. Both D11 and D12 are overestimated, on average, by almost a factor of 3, although the large standard deviations in Table 5 indicate that the results are not distributed symmetrically about the mean values. In the worst cases, the matrix elements D11 and D12 are overestimated by more than a factor of 8. Because this is the first example to date where the SSK approach is quite inaccurate, it is worthwhile examining the causes of this inaccuracy. The performance of the SSK approximation for this example is considerably worse for low coverages of species 1 than for other situations. To make this point quantitatively, Table 6 shows the statistics for DSSK ij /Dij for the data in Figure 8 after excluding the data with θ1 e 0.3. The SSK results for this restricted set of occupancies are considerably more accurate than when the full range of occupancies is used.

Figure 9. A comparison between D0,1(θ1 + θ2) (filled symbols) and }1(θ1,θ2) (open symbols) using MC data for species 1 from Example 1 (squares) and Example 4 (circles).

The SSK approach to predicting binary diffusivities relies on three independent assumptions. First, it is assumed that the MS diffusivities }i(θ1,θ2) in the binary mixture can be predicted by the MS diffusivities of species i as a single component at the same total loading. This approximation is expressed in eq 14, where it should be noted that the single-component corrected and MS diffusivities are the same quantity.13 The validity of this approximation is explored in Figure 9 for species 1 in Examples 1 and 4. Although the approximation is not exact for Example 1, it is quite accurate. This observation is consistent with the fact that the SSK approach gives accurate predictions for Example 1 (cf. Figure 2). The situation for Example 4 is quite different: the approximation can only be considered to be accurate for a few of the mixture occupancies shown in Figure 9. The qualitative reasons for this inaccuracy are straightforward to understand when the different characteristics of the two diffusing species during single-component diffusion are considered. The other two approximations made in the SSK approach are encapsulated in eq 15. The first assumption is that the mixture self-exchange coefficients, }corr ii , can be estimated by using the mixture’s total occupancy in eq 13. The second assumption is that, by using these mixture self-exchange coefficients, the mixture exchange coefficients, }corr ij , can be estimated using eq 15. The first of these assumptions is analogous to the treatment of }i(θ1,θ2) via eq 14. Since it was just demonstrated that this approach is not at all accurate for Example 4 (see Figure 9 and the preceding paragraph), it is not surprising that this approach is also inaccurate for describing }corr ii in the mixture. A more interesting question is whether the correlation defined in eq 15 is accurate if the exact mixture self-exchange coefficients are used. This procedure cannot be used in practice if the goal is to make predictions from single-component data, but it is a useful way to probe the validity of the correlation in eq 15. If the mixture self-diffusivities, Ds,i(θ1,θ2) and mixture MS diffusivities defined in eq 5 are known, then these quantities are related to }corr ii without approximation by13,51

θ2 θ1 1 1 + corr ) + corr Ds,i(θ1,θ2) }i(θ1,θ2) }ii (θ1,θ2) }ij (θ1,θ2)

(29)

with j * i. To apply this expression, I used the binary Fickian diffusivities from KMC simulations to define }i(θ1,θ2) and }corr 12 (θ1,θ2) via eq 7. By combining these MS diffusivities with (51) Krishna, R.; Paschek, D. Phys. Chem. Chem. Phys. 2002, 4, 1891.

3714 Langmuir, Vol. 22, No. 8, 2006

Figure 10. A comparison of the values of }corr 12 (θ1,θ2) determined from KMC simulations, with the results predicted using the correlation defined by eq 15 using the values of }corr ii (θ1,θ2) computed using KMC. Results are shown for all mixture occupancies simulated for Examples 1 and 4.

values of Ds,i(θ1,θ2) computed from KMC simulations, eq 29 was used to define }corr ii (θ1,θ2). Figure 10 compares the actual (θ values of }corr ,θ ) for all mixture occupancies of Examples 1 2 12 1 and 4 with the results that are obtained from eq 15 using the exact values of }corr ii (θ1,θ2). The interpolation formula, while not exact, performs quite accurately for both lattice models. It could even be argued that the accuracy of this interpolation method is slightly better for Example 4 than it is for Example 1, a surprising result given that the mixture diffusion in Example 1 is considerably simpler in character than that in Example 4. I reiterate that the results in Figure 10 use the self-exchange coefficients determined from mixture KMC simulations. This type of calculation is not possible if only single-component data is available. Figures 9 and 10 strongly suggest that if improvements in the SSK method are to be found, these improvements should focus on eq 14 and its analogue for the MS self-exchange diffusion coefficients.

Sholl

mixture diffusion in nanoporous materials that are relatively energetically homogeneous had all shown that this method gave accurate predictions.13-17,20 By using lattice models of mixture diffusion for a broad range of occupancies and compositions, the results presented here provide a stringent test on the accuracy of the SSK method (or any other similar approach). The SSK method was found to be very accurate for energetically homogeneous examples, but the predicted mixture diffusivities were rather inaccurate in examples with strong energetic heterogeneity. It is dangerous to make statements about any approximate method of the sort tested here being “very accurate” or “rather inaccurate” because the level of accuracy required in any specific application is of course dependent on the context in which the diffusivities are to be used. Nevertheless, the strong deviations seen between the predicted diffusivities and the actual binary diffusivities for the example with the strongest energetic heterogeneity (see Example 4 above) indicate that the SSK method as described here is not suitable for all conceivable materials. This study suggests several worthwhile directions for future work. It would of course be valuable to develop mixing theories for diffusion that perform better than the SSK method. Alternatives to the type of correlations used in the SSK method have been explored for the prediction of binary self-diffusion coefficients;39 it would be interesting to pursue these or other methods for the prediction of binary Fickian diffusion coefficients. This direction is based on the inability of the SSK method to treat strong energetic heterogeneity in a quantitative way. An alternative point of view is that there are clearly many examples in which the SSK method does work quantitatively. It would be valuable to understand the characteristics of single-component diffusion that would suggest when this method would (or would not) give acceptable results. Finally, all of the data presented here has been for relatively artificial lattice models. Performing fully atomistic simulations of molecular adsorption and diffusion in materials of practical interest with spatially heterogeneous energy landscapes could give useful insight into whether the strongly heterogeneous examples studied here are typical of real materials.

Conclusion

Acknowledgment. This work was supported by the NSF (CTS-0413027).

The aim of this paper has been to examine the accuracy of the SSK method for predicting macroscopic mixture diffusivities from single-component information in examples where diffusion occurs on a material defined by an energetically heterogeneous range of binding sites. Previous applications of this method to

Supporting Information Available: Binary Fickian diffusivities computed from KMC simulations of Examples 1-4. This material is available free of charge via the Internet at http://pubs.acs.org. LA053405B

Testing Predictions of Macroscopic Binary Diffusion

Recommend Documents