ARTICLE pubs.acs.org/IECR
Multicomponent MaxwellStefan Diffusivities at Infinite Dilution Xin Liu,†,‡ Andre Bardow,†,‡ and Thijs J. H. Vlugt*,‡ † ‡
Lehrstuhl f€ur Technische Thermodynamik, Schinkelstrasse 8, 52062 Aachen, Germany Process & Energy Laboratory, Delft University of Technology, Leeghwaterstraat 44, 2628 CA Delft, The Netherlands
bS Supporting Information ABSTRACT: Diffusion plays an important role in (bio)chemical processes. It is usually difficult to obtain MaxwellStefan diffusivities from experiments as well as molecular simulation. Therefore, predictive models based on easily measurable quantities are highly desired. The Vignes equation is commonly used to describe the concentration dependence of MaxwellStefan diffusivities. In mixtures containing at least three components, the generalized Vignes equation requires the value of the quantity ^xijkf1, which describes the friction between components i and j when both are infinitely diluted in component k. Over the past decades, several empirical models were proposed for estimating ^xijkf1, and all of these are lacking a sound theoretical basis. In this study, we show that ^xijkf1 actually exists (i.e., its value does not depend on the molar ratio xi/xj), and we derive an analytical expression for ^xijkf1 that is based on the linear response theory and the Onsager relations. We find that ^xijkf1 can be expressed in terms of binary and pure-component self-diffusivities and integrals over velocity cross-correlation functions. By neglecting the latter terms, we obtain a convenient predictive model for ^xijkf1. Molecular dynamics simulations are used to validate the assumptions made in this model. The following test systems are considered: a ternary system consisting of particles interacting using WeeksChandlerAndersen (WCA) interactions and the ternary systems n-hexanecyclohexanetoluene and ethanolmethanolwater. Our results show that (1) for the WCA system, as well as the n-hexanecyclohexanetoluene system, neglecting the integrals over velocity crosscorrelation functions results in accurate predictions for ^xijkf1; (2) for the WCA system, our model prediction is superior, compared to the existing models for ^xijkf1; (3) in the ethanolmethanolwater system, the integrals over velocity cross-correlation functions cannot be neglected, because of the presence of hydrogen bonds (models for predicting ^xijkf1 in this system will require detailed information on the collective motion of molecules); and (4) our model may provide an explanation why the MaxwellStefan diffusivity describing the friction between adsorbed components in a porous material is usually very large.
1. INTRODUCTION Diffusion plays an important role in (bio)chemical processes. The MaxwellStefan (MS) theory provides a sound theoretical basis and is therefore preferred for the modeling of diffusion.13 Unfortunately, it is difficult to obtain MS diffusivities from experiments or simulations. In experiments, Fick diffusivities are measured,4,5 which must be converted to MS diffusivities via a thermodynamic factor. This introduces large uncertainties.6,7 Extracting MS diffusivities from molecular simulations may lead to extensive CPU-time requirements.8,9 Therefore, predictive models for the concentration dependence of MS diffusivities based on easily measurable quantities are highly desired. The Vignes equation is often recommended to predict the concentration dependence of MS diffusivities.1,3 In ternary mixtures, the generalized Vignes equation is given as10 x f 1 xj
^ij ¼ ð^xiji f 1 Þxi ð^ijj
Þ ð^xijk f 1 Þxk
^ijxk f 1 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x f 1 xi f 1 ^ijj ^ij
ð2Þ
o KT (Kooijman and Taylor11): ^ijxk f 1 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^xikk f 1 ^jkxk f 1
ð3Þ
o VKB (Krishna and van Baten12):
ð1Þ
in which ^ij is the MS diffusivity that describes the friction between components i and j, and xi is the mole fraction of component i. The terms ^ijxif1 and ^ijxjf1 describe the friction between a diluted component and the solvent; therefore, these are equal to their value in a binary mixture of i and j only. The generalized Vignes equation requires the value of ^xijkf1, which describes the friction between components i and j when r 2011 American Chemical Society
both are infinitely diluted in a third component k. This quantity is not easily accessible in experiments, since no direct measurement is possible.6 In the past decades, several empirical models were proposed for estimating ^xijkf1 from binary diffusion data: o WK (Wesselingh and Krishna10):
^ijxk f 1 ¼ ð^ikxk f 1 Þxi =ðxi þ xj Þ ð^jkxk f 1 Þxj =ðxi þ xj Þ ð4Þ
Received: December 16, 2010 Accepted: February 24, 2011 Revised: February 18, 2011 Published: March 15, 2011 4776
dx.doi.org/10.1021/ie102515w | Ind. Eng. Chem. Res. 2011, 50, 4776–4782
Industrial & Engineering Chemistry Research o DKB (Krishna and van Baten12): xj xi ^ijxk f 1 ¼ ^ xk f 1 þ ^ xk f 1 xi þ xj ik xi þ xj jk
ARTICLE
ð5Þ
where Γij is the thermodynamic factor, which is defined as ! D ln γi Γij ¼ δij þ xi ð8Þ Dxj p, T , xk, k6¼i, j
ð6Þ
Here, δij is the Kronecker delta and γi denotes the activity coefficient of component i. MD simulations can be used to compute the MS diffusivities ^ij directly from local fluctuations. First, the so-called Onsager coefficients (Λij) can be obtained directly from the motion of the molecules in MD simulations,12
o RS (Rehfeldt and Stichlmair13,14): x f1
^xijk f 1 ¼ ð^ikxk f 1 ^jkxk f 1 ^ijj
^ijxi f 1 Þ1=4
These empirical models can be classified into three categories: (1) the friction between the two diluted components is taken into account and the friction between the diluted components and the solvent is neglected (i.e., the WK model); (2) the friction between the diluted components and the solvent is taken into account and the friction between the diluted components themselves is neglected (i.e., the KT, VKB, and DKB models); and (3) a combination of the first two categories (i.e., the RS model). Since there is no physical basis for any of these models, it is unclear which one to use for a specific system. It is important to note that errors introduced by modeling the concentration dependence of MS diffusivities using the generalized Vignes equation (eq 1) may be either be reduced or enhanced by a particular choice for a model for ^xijkf1.15 In the remainder of the manuscript, we will not consider this issue further and, instead, will focus on how ^xijkf1 can be predicted. In this study, we show that ^xijkf1 is a well-defined quantity that is independent of the ratio xi/xj. Based on linear response theory and the Onsager relations, we propose a new model for estimating ^xijkf1 from easily measurable self-diffusivities. Unlike the WK, KT, VKB, DKB, and RS models, our model has a sound theoretical basis. The key assumption is that velocity cross-correlations between particles are negligible. We perform molecular dynamics (MD) simulations to test the quality of our new model. It turns out that the model works very well for systems consisting of nonassociating molecules. For the watermethanolethanol system, the earlier mentioned cross-correlations between particle velocities cannot be neglected, because of the collective motion of molecules, originating from networks of hydrogen bonds. This paper is organized as follows. In section 2, we explain how MS diffusivities can be obtained from MD simulations. We show that ^xijkf1 does not depend on the molar ratio xi/xj and we derive a model for estimating ^xijkf1 using easily measurable pure-component and binary self-diffusivities. Details about the MD simulations are addressed in section 3. In section 4, we test our model for the following ternary systems: (1) a system without attractions (i.e., particles interact using a WCA potential);16 (2) the n-hexanecyclohexanetoluene system, in which molecules interact using Lennard-Jones (LJ) interactions; and (3) the ethanolmethanolwater system, in which electrostatic interactions play an important role. We compare the predictions of ^xijkf1 with the results from MD simulations. Our conclusions are summarized in section 5.
2. OBTAINING MAXWELLSTEFAN DIFFUSIVITIES FROM MD SIMULATIONS The MS formulation uses chemical potential gradients as driving forces for mass transport. For liquid mixtures at constant temperature and pressure, the MS equation can be expressed as17 n1
Λij ¼
10 *0 N i @ ðrl, i ðt þ mΔtÞ rl, i ðtÞÞA@
∑
l¼1
k¼1 k6¼ i
t
Nj
∑
k¼1
1+ ðrk, j ðt þ mΔtÞ rk, j ðtÞÞA
ð9Þ An alternative but equivalent expression for obtaining Λij is *N + Z Nj i 1 ¥ Λij ¼ dt vl, i ð0Þ vk, j ðtÞ ð10Þ 3N 0 l¼1 k¼1
∑
∑
In these equations, N is the total number of molecules in the simulation, m the number of time steps, Δt the integration time step used in the simulations, rl,i(t) the position of the center of mass of the lth molecule of component i at time t, and vl,i(t) its velocity. Note that the matrix [Λ] is symmetric; that is, Λij = Λji. Second, a matrix [Δ] is defined with elements that follow from the Onsager coefficients Λij.12 For a ternary system, they are given by Λ11 Λ13 Λ21 Λ23 Λ31 Λ33 x1 Δ11 ¼ ð1 x1 Þ þ x1 x3 x1 x3 x1 x3 ð11Þ
Δ12
Λ12 Λ13 Λ22 Λ23 Λ32 Λ33 x1 ¼ ð1 x1 Þ þ x2 x3 x2 x3 x2 x3
ð12Þ
Δ21
Λ21 Λ23 Λ11 Λ13 Λ31 Λ33 x2 ¼ ð1 x2 Þ þ x1 x3 x1 x3 x1 x3
ð13Þ Δ22 ¼ ð1 x2 Þ
Λ22 Λ23 Λ12 Λ13 Λ32 Λ33 x2 þ x2 x3 x2 x3 x2 x3
ð14Þ Third, the matrix [Δ] can be inverted to obtained the matrix [B]: ½B ¼ ½Δ1
ð15Þ
Finally, the MS diffusivities follow directly from the elements of the matrix [B]:12
n
xi Jk xk Ji ∑ Γij rxj ¼ ∑ c^ j¼1
1 1 lim m f ¥ 6N mΔt
^13 ¼
ð7Þ
ik
4777
1 x2 B12 B11 þ x1
ð16Þ
dx.doi.org/10.1021/ie102515w |Ind. Eng. Chem. Res. 2011, 50, 4776–4782
Industrial & Engineering Chemistry Research ^12 ¼
ARTICLE
1 x1 þ x3 B11 B12 x1
ð17Þ
1 x1 B22 þ B21 x2
ð18Þ
^23 ¼
^ij ¼ f ðxi , xj , xk , ½ΛÞ
ð19Þ
Because the resulting expression is rather complicated, it is provided as eq S1 in the Supporting Information. Using this expression, we consider the situation in which components i and j are infinitely diluted in component k. It is possible to find a convenient expression for the term Λii using the GreenKubo formulation of the Onsager coefficients (eq 10). Since component i is infinitely diluted in component k, we can neglect correlations in the velocities between two distinct molecules of component i. Therefore, *N + Z Ni i 1 ¥ dt vl, i ð0Þ vg , i ðt Þ Λii ¼ 3N 0 g ¼1 l¼1
∑
1 ¼ 3N
Z
dt
¼
Z
Z
Ni 3N
∑ vl, i ð0Þ vl, i ðtÞ
dt * dt ¥
+
Ni
Ni
∑ ∑ vl, i ð0Þ vg, i ðtÞ
+
l¼1 g ¼1 g 6¼ i
0
0
Z
Ni
*
¥
¥
∑
l¼1
0
1 þ 3N
1 3N
*
¥
∑
Ni Nj 3N
^xijkf1.
Using eqs 1118, it is possible to write ^ij as a function of xi, xj, xk, and all elements of the matrix [Λ]: 2.1. Deriving an Expression for
in which Ckk and C*kk account for self- and cross-correlations of the velocities of molecules of type k, respectively. The value of Ckk is independent of the total number of molecules. For Λij with i 6¼ j (i.e., correlations between unlike molecules), we can write + Nj Z ¥ * Ni 1 dt vl, i ð0Þ vk, j ðtÞ Λij ¼ 3N 0 l¼1 k¼1
Ni
∑ vl, i ð0Þ vl, i ðtÞ
+
l¼1
dt vi, 1 ð0Þ vi, 1 ðt Þ ¼ xi Cii
ð22Þ
Z
In the limit of infinite dilution (which requires the thermodynamic limit), eqs 22 and 21 will be exact. For the solvent k, the expression for Λkk is + Z ¥ * Nk Nk 1 dt vl, k ð0Þ vg , k ðtÞ Λkk ¼ 3N 0 g¼1 l¼1 * + Z Nk 1 ¥ ¼ dt vl, k ð0Þvl, k ðtÞ 3N 0 l¼1 *N + Z Nk k 1 ¥ þ dt vl, k ð0Þvg , k ðtÞ 3N 0 l ¼ 1 g ¼ 1, g6¼ l
∑
∑
∑
∑ ∑
xk Ckk þ xk 2 NCkk
ð22Þ
dt Æv1, i ð0Þv1, j ðtÞæ ¼ Nxi xj Cij
ð23Þ
0
By inserting eqs 2223 into eq S1 (see the Supporting Information) and setting xj = axi (with a being a positive constant), we take the limit that xi f 0. The resulting expression is ^xijk f 1 ¼
Cii Cjj Ckk þ Cx
ð24Þ
with Cx ¼ NðCij Cik Cjk þ Ckk Þ
ð25Þ
Equation 24 clearly shows that ^xijkf1 is independent of the ratio a = xj/xi. This was already suggested in the study by Kooijman and Taylor 20 years ago;11 however, in that study, no formal proof was provided. The fact that ^xijkf1 should be independent of xi/xj suggests that eqs 4 and 5 violate Onsager’s reciprocal relations at infinite dilution. Note that all terms in eq 25 originate from velocity cross-correlations (i.e., correlations of the velocities of different molecules). Because the Onsager coefficients Λ are intensive properties, the terms NC*, kk NCij, NCik, and NCjk will approach finite values when N f ¥. Therefore, the term Cx converges to a finite value when N f ¥. 2.2. Obtaining a Physical Model for ^xijkf1. In general, the self-diffusivity of component i in a medium can be computed from a GreenKubo relation18 Z 1 ¥ Di, self ¼ dt Ævi ð0Þvi ðtÞæ ð26Þ 3 0 in which Di,self is the self-diffusivity of component i, and vi(0) and vi(t) are the velocities of a molecule of component i at time 0 and time t, respectively. From eqs 22 and 26, it follows that
0
in which Cii is a constant that is not dependent on the total number of molecules. Similarly, we can write, for the other component j that is diluted in component k, ð21Þ Λjj xj Cjj
¥
∑
Cii ¼ Di, self
ð27Þ
Therefore, eq 24 can be written as ^xijk f 1 ¼
Di, self Dj, self Dk, self þ Cx
ð28Þ
in which Di,self, Dj,self, and Dk,self are taken from a system in which xk f 1. It is natural to assume that integrals of velocity autocorrelations are much larger than integrals of velocity cross-correlations, suggesting that Cx is small, compared to Dk,self. Using this assumption, eq 28 then becomes ^ijxk f 1 ¼
^xikk f 1 ^xjkk f 1 Di, self Dj, self ¼ Dk, self Dk, self
ð29Þ
This equation allows us to predict the ternary diffusivity ^xijkf1 based on the binary diffusion coefficients ^xijkf1, ^xijkf1 and the pure component self-diffusivity Dk,self. In the remainder of the paper, we will verify the predictions of eqs 28 and 29 using MD simulations. It is worthwhile to note that eq 29 reduces to the KT model (eq 3) in the case that Dk,self = (^xijkf1 ^xijkf1)1/2. 4778
dx.doi.org/10.1021/ie102515w |Ind. Eng. Chem. Res. 2011, 50, 4776–4782
Industrial & Engineering Chemistry Research
3. DETAILS OF MOLECULAR DYNAMICS SIMULATIONS The predictions of eqs 28 and 29 are tested for the following ternary systems: (1) A simple system in which only repulsive interactions are considered; here, particles interact using a WCA potential.16 The three components only differ in their molar mass. (2) The n-hexanecyclohexanetoluene system, in which particles interact using LJ potentials. (3) The ethanolmethanolwater system in which both LJ and electrostatic interactions are considered. For all systems, diffusivities were obtained from MD simulations in the NVE ensemble (constant number of particles, constant total energy, and constant volume). Self- and MS diffusivities are computed using eqs 1018, as well as eq 26. The details of extracting correlation functions and mean-squared displacements from the MD trajectories are described in refs 15,18, and 19. Errors in computed diffusivities are calculated by performing at least five independent MD simulations for each system (i.e., simulations using different initial positions of the molecules) and analyzing their differences. 3.1. WCA System. Simulations were carried out in a system containing 400 particles that interact using a WCA potential.16 The WCA potential is constructed by truncating and shifting the LJ potential at 21/6σ. A linked-cell algorithm is applied to improve the efficiency.18 For convenience, we express all quantities in reduced units18 by setting the LJ parameters σ and ε as units for length and energy. The mass of the lightest component is set as a mass unit (i.e., M1 = 1). We investigate the influence of the mass of the solvent molecules (here, M3) and the total number density (F). 3.2. n-Hexane-Cyclohexane-Toluene. The united-atom model is used in which CHx groups are considered as single interaction sites. For n-hexane and cyclohexane, the SHAKE algorithm is used to constrain the distance between two neighboring atoms.20 Bond bending and torsion potentials describe the interaction between three and four consecutive atoms, respectively. Toluene is treated as a rigid molecule. The nonbonded interactions are described by LJ potentials. The LorentzBerthelot mixing rule is applied to obtain the LJ parameters for the interaction of unlike atoms.21 LJ interactions are truncated and shifted at 12 Å. The equations of motion were integrated using the velocity Verlet algorithm with a time step of 1 fs. The simulations were performed at a total energy and density corresponding to a pressure of 1 atm and a temperature of 298 K. The force field parameters for this system are listed in Table S1 of the Supporting Information, as well as in ref 15. 3.3. EthanolMethanolWater. In the ethanolmethanol water system, intermolecular interactions are described by LJ and electrostatic interactions with parameters listed in Table S2 of the Supporting Information. LJ interactions are truncated and shifted at 12 Å. Electrostatic interactions are handled by the Ewald summation, using a relative precision of 104. Ethanol and methanol are described using TraPPE-UA model,22 in which the CHx groups are considered as united atoms. A flexible simple point-charge (SPCFW) model is used to describe water.23 The simulations were performed at a total energy and density corresponding to a pressure of 1 atm and a temperature of 298 K. 4. RESULTS AND DISCUSSION 4.1. WCA System. In Table 1, a comparison is made between the MS diffusivities obtained from MD simulations and the
ARTICLE
Table 1. Computed MS Diffusivities Compared to the Predictions of eqs 28 and 29 for Ternary Mixtures in Which Particles Interact Using a WCA Potentiala MS Diffusivity, ^12 Prediction MD simulation eq 28 ADb eq 29 ADb Number Densityc 0.1
3.294
3.212 2%
3.341 1%
0.2
1.411
1.401 1%
1.441 2%
0.5
0.318
0.310 2%
0.315 1%
0.7
0.144
0.145 1%
0.149 3%
0.8
0.094
0.092 2%
0.097 4%
Mass of Species 3, F = 0.2 d
5
1.411
1.401 1%
1.441 2%
50
2.692
2.686 0%
2.560 5%
100
3.348
3.344 0%
3.288 2%
500
6.238
6.237 0%
6.619 6%
2500
17.609
17.609 0% 18.923 7%
41.794
41.792 0% 44.051 5%
5000 Mass of Species 3, F = 0.5 e
5 10
0.318 0.270
0.310 2% 0.267 1%
0.315 1% 0.275 2%
50
0.218
0.217 0%
0.227 5%
100
0.172
0.172 0%
0.161 7%
f
Ratio of x1/x2, with x3 = 0.95 1/3
0.318
0.312 2%
0.312 2%
1
0.318
0.310 2%
0.315 1%
3
0.318
0.312 2%
0.326 2%
1/3 1
0.318 0.318
0.316 1% 0.320 1%
0.320 1% 0.325 2%
3
0.319
0.322 1%
0.321 1%
g
Ratio of x1/x2, with x3 = 0.97
a
Equations 28 and 29 are parametrized using self-diffusivities and velocity cross-correlations. Details of the self-diffusivities and velocity cross-correlations are listed in Tables S3 and S4 of the Supporting Information. Here, we consider the case that components 1 and 2 are diluted in component 3. All reported quantities are in reduced units. The statistical errors in the computed diffusivities are
