Fick Diffusion Coefficients in Ternary Liquid Systems from Equilibrium

Article pubs.acs.org/IECR

Fick Diffusion Coefficients in Ternary Liquid Systems from Equilibrium Molecular Dynamics Simulations Xin Liu,†,‡ Ana Martín-Calvo,¶ Erin McGarrity,‡ Sondre K. Schnell,‡ Sofía Calero,¶ Jean-Marc Simon,§ Dick Bedeaux,∥ Signe Kjelstrup,∥,‡ André Bardow,†,‡ and Thijs J. H. Vlugt*,‡ †

Lehrstuhl für Technische Thermodynamik, RWTH Aachen University, Schinkelstraße 8, 52062 Aachen, Germany Process & Energy Laboratory, Delft University of Technology, Leeghwaterstraat 44, 2628 CA Delft, The Netherlands ¶ Physical, Chemical, and Natural Systems, University Pablo de Olavide Ctra. de Utrera, km. 1 Sevilla, Spain § Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR 5209 CNRS-Universite de Bourgogne, Dijon, France ∥ Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway ‡

S Supporting Information *

ABSTRACT: An approach for computing Fick diffusivities directly from equilibrium molecular dynamics (MD) simulations is presented and demonstrated for a ternary chloroform−acetone−methanol liquid mixture. In our approach, Fick diffusivities are calculated from the Maxwell−Stefan (MS) diffusivities and the so-called matrix of thermodynamic factors. MS diffusivities describe the friction between different molecular species and can be directly computed from MD simulations. The thermodynamic factor describes the deviation from ideal mixing behavior and is difficult to extract from both experiments and simulations. Here, we show that the thermodynamic factor in ternary systems can be obtained from density fluctuations in small subsystems embedded in a larger simulation box. Since the computation uses the Kirkwood−Buff coefficients, the present approach provides a general route toward the thermodynamics of the mixture. In experiments, Fick diffusion coefficients are measured, while previously equilibrium molecular dynamics simulation only provided MS transport diffusivities. Our approach provides an efficient and accurate route to predict multicomponent diffusion coefficients in liquids based on a consistent molecular picture and therefore bridges the gap between theory and experiment. and the elements often strongly depend on concentration.2 Also, in multicomponent systems, the elements of the matrix [D] are unrelated to their binary counterparts. The off-diagonal Fick diffusivities Dij with i ≠ j, can even be negative.1 The Maxwell−Stefan (MS) approach provides a physically based comprehensive framework for describing multicomponent mass transport.1,2,4,5 The approach relates the driving force for diffusion of component i (i.e. the chemical potential gradient ▽μi at constant temperature and pressure) to friction forces of molecules of component i with molecules of other components:

1. INTRODUCTION Chemical engineers often need to compute properties of systems undergoing multicomponent mass transfer. Diffusion is often rate-limiting for mass transport in liquids. Transport diffusion is a process by which mass transfer occurs due to a gradient in chemical potential or concentration as the driving force. The presence of three or more components in a liquid mixture introduces phenomena that cannot be described by the well-known first law of Fick for binary systems.1 A consistent methodology to accurately describe and predict mass transfer in multicomponent systems is therefore required. Two theories are commonly used to describe the mass transport: generalized Fick’s law and the Maxwell−Stefan approach. Generalized Fick’s law in an n-component system with respect to a molar reference frame is given by1−3

−

Ji = −ct ∑ Dij∇xj (1)

in which Ji is the diffusion flux, ct is the total molar concentration, Dij are Fick diffusivities, and ∇xj is the gradient in mole fraction of component j. Transformations between fluxes in other reference frames can be easily performed.1,2,4 From eq 1 it follows directly that the elements of the matrix of Fick diffusivities [D] depend on the labeling of the components. In an n-component system, transport diffusion is thus described by (n − 1)2 Fick diffusion coefficients. It is important to note that the matrix of Fick diffusivities [D] is not symmetric © 2012 American Chemical Society

n

∑ j = 1, j ≠ i

xj(ui − uj) Đij

(2)

in which R and T are the gas constant and absolute temperature, respectively. xj is the mole fraction of component j. The friction force between components i and j is proportional to the difference in average velocities of the components, (ui − uj). The MS diffusivity Đij can be considered as an inverse friction coefficient describing the magnitude of the friction between components i and j. MS diffusivities do not depend on a reference frame. The MS diffusivities are symmetric, Đij = Đji, and

n−1 j=1

1 ∇T , p μi = RT

Received: Revised: Accepted: Published: 10247

April 18, 2012 June 28, 2012 July 6, 2012 July 9, 2012 dx.doi.org/10.1021/ie301009v | Ind. Eng. Chem. Res. 2012, 51, 10247−10258

Industrial & Engineering Chemistry Research

Article

Table 1. Densities of Pure Chloroform, Acetone, and Methanol at 298 K, 1 atm ρ/(g·mL−1) component

this work (MD)

experiment

acetone methanol CHCl3

0.76 0.75 1.48

0.78a,b, 0.79c 0.79a,b 1.47d

a c

Experiments by Campbell et al.62 bExperiments by Noda et al.63 Experiments reported in ref 64. dExperiments reported in ref 58.

Figure 1. Densities of chloroform (1)−acetone (2) mixtures at 298 K, 1 atm. Circles are the computed densities by this work using MD simulations. Squares are the experimental densities measured by Karr et al.58

therefore, n(n − 1)/2 MS diffusivities Đij are sufficient to describe diffusion in an n-component system. However, it is difficult to obtain MS diffusivities from experiments as chemical potentials cannot be measured directly. In contrast, obtaining MS diffusivities from molecular dynamics (MD) simulations is a natural choice but often requires large amounts of CPU time.6−17 It is important to note that the MS diffusivities are also concentration dependent. In our earlier work, we developed the multicomponent Darken model for describing this concentration dependence of MS diffusivities.11 The multicomponent Darken model accurately predicts the MS diffusivities for ideal diffusing mixtures.11 As generalized Fick’s law and the MS theory describe the same physical process, their coefficients are related. In an n-component system in a molar reference frame, this relation is1,2,4,7 −1

[D] = [B] [Γ]

Figure 2. Thermodynamic factor Γ11 in the binary system chloroform (1)−acetone (2) at 298 K, 1 atm. Open symbols are the computed Γ11 in this work using MD simulations. Filled symbols are the computed Γ11 using COSMO-SAC. The solid line represents Γ11 calculated from the NRTL model fitted to experimental VLE data of ref 59.

∂ ln γi ∂xj

xi + Đin

n

∑

(3)

j = 1, j ≠ i

xj Đij

⎛1 1 ⎞⎟ Bij = − xi⎜⎜ − Đin ⎟⎠ ⎝ Đij

∂xj

− T ,p,j′

∂ ln γi ∂xn

T ,p,n′

(7)

D11 = Γ11Đ12

with i = 1, ..., (n − 1)

in which the thermodynamic factor Γ11 is given by

(4)

⎛ ∂ ln γ1 ⎞ Γ11 = 1 + x1⎜ ⎟ ⎝ ∂x1 ⎠T , p , Σ

with i , j = 1, ..., (n − 1) and i ≠ j (5)

(8) 1,2

(9)

Equation 3 clearly shows the gap between molecular simulations and experiments. In experiments, Fick diffusion coefficients are measured while equilibrium MD simulations usually provide MS diffusivities.9,19−23 The two worlds are related via the matrix of thermodynamic factors [Γ], but this is usually known only with relatively large uncertainties.18,24 To obtain the matrix of thermodynamic factors, the experimental vapor−liquid equilibrium data can be fitted using excess Gibbs energy models, e.g. Margules, van Laar, NRTL, etc.25 Several models may provide estimates of ln γi that give equally good fits to the vapor−liquid equilibrium data. However, the thermodynamic factor Γij involves the first derivative of the activity coefficient ln γi with respect to the composition. Errors of the size of 20% and larger are expected for this derivative.26

The elements of the so-called matrix of thermodynamic factors [Γ] are defined by1−3,7 ⎛ ∂ ln γ ⎞ i⎟ Γij = δij + xi⎜⎜ ⎟ x ∂ ⎝ j ⎠ T ,p,Σ

T ,p,Σ

∂ ln γi

Here, the prime symbols in the derivative evaluations indicate that the mole fractions of all other components are held constant. For any system in the limit of infinite dilution, the i→0 values of [Γ] are known, i.e., Γxiii→0 = 1, Γxij,i≠j = 0. For a binary 1,2 system, eqs 3−6 reduce to

in which [D] is the (n − 1) × (n − 1) matrix of Fick diffusivities. The elements of the matrix [B] are given by1,3,7,18 Bii =

=

(6)

in which δij is the Kronecker delta and γi is the activity coefficient of component i in the mixture. The symbol Σ indicates that the partial differentiation of ln γi with respect to mole fraction xj is carried out at constant mole fraction of all other components except the nth one, so that Σni=1 xi = 1 during the differentiation. The constrained derivative of eq 6 can be written as a function of unconstrained derivatives as follows: 10248

dx.doi.org/10.1021/ie301009v | Ind. Eng. Chem. Res. 2012, 51, 10247−10258


Article

Figure 3. Computed MS diffusivities Đ12 and Fick diffusivities D11 in the binary system chloroform (1)−acetone (2) at 298 K, 1 atm. Open symbols represent Fick diffusivities calculated using the computed thermodynamic factor Γ11 and MS diffusivities Đ12 from MD simulations. The error bars of computed diffusivities are smaller than the symbol size.

Figure 5. Thermodynamic factor Γ11 in the binary system chloroform (1)−methanol (2) at 298 K, 1 atm. Open symbols are the computed Γ11 by this work using MD simulations. Filled symbols represent Γ11 computed using COSMO-SAC. The solid line represents Γ11 calculated from the NRTL model fitted to experimental VLE data of ref 59.

Figure 4. Densities of chloroform (1)−methanol (2) mixtures at 298 K, 1 atm. Circles are the computed densities from this work using MD simulations. Squares are the experimental densities taken from ref 60.

Figure 6. Computed MS diffusivity Đ12 and Fick diffusivities D11 in the binary system chloroform (1)−methanol (2) at 298 K, 1 atm. Fick diffusivities calculated using the computed thermodynamic factor Γ11 and MS diffusivities Đ12 from MD simulations. The error bars of computed diffusivities are smaller than the symbol size.

Currently used approaches for computing Fick diffusivities from molecular simulation suffer from inconsistencies or other problems: • Direct calculation of Fick diffusivities using nonequilibrium MD (NEMD) requires significant efforts and very high concentration gradients.27−29 This approach is usually not accurate and quite impractical as the concentration dependence of Fick diffusivities is not easily captured. • Combining MS diffusivities obtained from equilibrium MD simulations with experimentally obtained equations of state or models for the excess Gibbs energy is inconsistent, as experiments and molecular models in principle provide different values for the thermodynamic factor.22,30,31 • Currently used molecular simulation techniques to determine the thermodynamic factor of liquid mixtures

are quite inefficient.21,32 To the best of our knowledge, these techniques have only been applied to binary mixtures. Recently, we developed a consistent methodology to calculate binary Fick diffusivities using MS diffusivities Đij and thermodynamic factors [Γ] from equilibrium MD simulations.33,34 Thermodynamic factors can be computed by studying density fluctuations in small subvolumes inside a larger simulation box and correcting for finite-size effects.32,35 In our earlier work, this approach was validated for the binary systems acetone− methanol and acetone−tetrachloromethane.33,34 In these simulations, all molecules were treated as rigid bodies interacting through Lennard-Jones (LJ) and electrostatic interactions. Our MD results for these systems quantitatively agree with the experimental data for binary Fick diffusion coefficients D11 as well for the thermodynamic factor Γ11.33,34 Our approach thus bridges the gap between diffusion experiments and molecular simulations. In this work, this approach is extended and validated 10249



Article

This paper is organized as follows. In section 2, we explain how to obtain the MS diffusivities and thermodynamic factors from equilibrium MD simulations. The details of the simulations are addressed in section 3. In section 4, we validate the methodology for computing Fick diffusivities for the ternary system chloroform−acetone−methanol. Our findings are summarized in section 5.

2. COMPUTATION OF DIFFUSION COEFFICIENTS AND THERMODYNAMIC FACTORS 2.1. Obtaining Diffusion Coefficients from MD Simulations. In equilibrium MD simulations, representative trajectories of a system consisting of interacting molecules are obtained.36−38 From these trajectories, transport properties can be computed.36 The self-diffusivity is related to the motion of individual molecules and can be calculated using the average displacements of molecules of type i as described by the Einstein equation:36

Figure 7. Densities of chloroform (1)−acetone (2)−methanol (3) mixtures at 298 K, 1 atm, computed from MD. The densities are plotted as a function of the mole fraction of one of the components, while keeping the mole fractions of the other components equal to each other.

for the ternary system chloroform−acetone-methanol. Our MD results show good agreement with experiments. Even with simple molecular models taken from standard classical force fields, i.e. excluding polarization effect, we shall see that it is possible to predict Fick diffusivities with reasonable accuracy. This may suggest that the use of more realistic/complex force fields is not needed for predicting transport diffusivities of typical small molecules.

N

Di ,self =

i 1 1 lim ⟨(∑ (rl , i(t + m ·Δt ) − rl , i(t ))2 )⟩ 6Ni m →∞ m ·Δt l = 1

(10)

Here, Ni is the total number of molecules of component i, m is the number of time steps and Δt is the time step used in MD simulations. rl,i(t) is the position of the l-th molecule of type i at time t.

Figure 8. Thermodynamic factor Γij in the ternary system chloroform (1)−acetone (2)−methanol (3) at 1 atm. Open circles are the computed values of Γij using MD simulations at 298 K. Filled circles are the computed values of Γij using COSMO-SAC at 298K. Solid lines represent Γij calculated from the NRTL model fitted to experimental VLE data at 303 K of ref 61. x1 is varied while keeping x2 = x3. 10250



Article

Figure 9. Thermodynamic factor Γij in the ternary system chloroform (1)−acetone (2)−methanol (3) at 1 atm. Open circles are the computed values of Γij using MD simulations at 298 K. Filled circles are the computed values of Γij using COSMO-SAC at 298 K. Solid lines represent Γij calculated from the NRTL model fitted to experimental VLE data at 303 K of ref 61. x2 is varied while keeping x1 = x3.

MS diffusivities Đij are calculated as follows: first, the Onsager coefficients Λij are calculated:7

improvement in this area but these simulations remain quite inefficient. 42,43 Kirkwood and Buff showed that density fluctuations in the grand-canonical ensemble are related to the integrals of radial distribution functions over volume, resulting in the following expression for the so-called Kirkwood−Buff (KB) coefficients for components i and j:40,41

N

Λij =

i 1 1 1 lim ⟨(∑ (rl , i(t + m ·Δt ) − rl , i(t ))) 6 m →∞ N m ·Δt l = 1

Nj

× (∑ (rk , j(t + m ·Δt ) − rk , j(t )))⟩ k=1

(11)

Gij = V

in which N is the total number of molecules and n is the number of components. The Onsager coefficients are constrained by ∑ni=1MiΛij = 0 in which Mi is the molar mass of component i.7 Furthermore, by definition Λij = Λji. Using the computed values for Λij and the mixture composition, MS diffusivities can be computed. The mean-squared displacement and Onsager coefficients were updated with different frequencies according to the order-n algorithm described in refs 36 and 39. The expressions relating the Onsager coefficients Λij to the MS diffusivities in binary, ternary, and quaternary systems can be found in refs 7 and 11. 2.2. Obtaining Thermodynamic Factors from MD Simulations. The elements of the matrix of thermodynamic factors [Γ] can be expressed as density fluctuations in the grandcanonical ensemble, as derived by Kirkwood and Buff40,41 in 1951. The natural method for obtaining these averages are grandcanonical Monte Carlo (GCMC) simulations.36 However, GCMC simulations of liquid mixtures at room temperature are quite challenging as the insertion and deletion of molecules is very inefficient for dense liquids.36 There is some recent

⟨NN i j⟩ − ⟨Ni⟩⟨Nj⟩

= 4π

⟨Ni⟩⟨Nj⟩

∫0

∞

−

δij ci

[gij(r ) − 1]r 2 dr

(12) (13)

In these equations, V is the volume and ci is the number density of component i, defined by ⟨Ni⟩/V. The brackets ⟨···⟩ denote an ensemble average in the grand-canonical ensemble. gij(r) is the radial distribution function for molecules of type i and j in the grandcanonical ensemble, and it is natural to define the distance r between two molecules as either the distance between their centers of mass or the distance between their central atoms. As the choice of ensemble is irrelevant for large systems, in practice gij(r) is computed in the NpT or NVT ensemble. In binary mixtures, the thermodynamic factor Γ11 is related to the KB coefficients Gij by40 Γ11 = 1 − x1

c 2(G11 + G22 − 2G12) x1c 2(G11 + G22 − 2G12)

(14)

In ternary systems, the elements of thermodynamic factor [Γ] are related to the KB coefficients Gij by41,44 10251



Article

Figure 10. Thermodynamic factor Γij in the ternary system chloroform (1)−acetone (2)−methanol (3) at 1 atm. Open circles are the computed values of Γij using MD simulations at 298 K. Filled circles are the computed values of Γij using COSMO-SAC at 298 K. Solid lines represent Γij calculated from the NRTL model fitted to experimental VLE data at 303 K of ref 61. x3 is varied while keeping x1 = x2.

1 Γ11 = − [−c 2c3G22 − c 2 + 2c 2c3G23 − c 2c3G33 − c3 η + c1(c 2G12 − c 2G22 − 1 + c 2G23 − c 2G13)] Γ12 = −

(16)

c2 (c1G11 − c1G12 − c3G12 − c1G13 + c3G13 + c1G23 η + c3G23 − c3G33)

Γ22 =

(15)

c1 (c 2G12 + c3G12 − c 2G13 − c3G13 − c 2G22 η

+ c 2G23 − c3G23 + c3G33) Γ21 =

The KB coefficients allow the computation of the thermodynamic factor from MD simulations and thus avoid the insertion and deletion of particles as needed in simulations in the grandcanonical ensemble. It is important to consider the convergence of the integrals in eq 13. For infinitely large systems, gij(r) → 1 for r → ∞. For finite systems with periodic boundary conditions, however, gij(r) does not converge to 1 for r → ∞ resulting in a divergence of Gij.41 Recently, several studies showed that the obtained value of Gij using eq 13 depends on the system size.32,45 In the Supporting Information, we show that the integrals of the radial distribution function over volume strongly vary with system size for a binary Weeks−Chandler−Andersen (WCA) fluid and for a binary system of acetone−methanol. To improve the convergence of Gij, a natural choice is to make the system larger and to use longer simulations. However, this is computationally inefficient and thus seriously hinders the direct estimation of KB coefficients from simulations. Recently, some effort has been made in computing the KB coefficients using different methods. Wedberg et al. accounted for finite size effects in simulations by introducing tail corrections to the correlation functions based on integral equation theory.46 These corrections are based on the hypernetted chain closure which is accurate at high densities. This method requires the optimization of two parameters, the joining radius and the tail approximation. Nichols et al. discussed the effects of truncation of the KB integrals.47 In addition, a method for estimating the KB integrals using the Fourier transformed atomic positions is described by these authors. This method is shown to account for the corners

(17)

1 (c1c3G11 + c1 − 2c1c3G13 + c1c3G33 + c3 η + c 2(c1G11 − c1G12 − c1G13 + 1 + c1G23))

(18)

in which η = c1 + c 2 + c3 + c1c 2Δ12 + c 2c3Δ23 + c1c3Δ13 1 − c1c 2c3(Δ12 2 + Δ232 + Δ132 − 2Δ13Δ23 − 2Δ12 Δ13 4 − 2Δ12 Δ23)

(19)

Δij = Gii + Gjj − 2Gij

(20)

and

10252



Article

Figure 11. Comparison of the computed MS diffusivities Đij and the predicted MS diffusivities Đij in the ternary system chloroform (1)−acetone (2)− methanol (3) at 298 K, 1 atm. x1 varies while keeping x2 = x3. Circles are the computed Đij using MD simulations. Triangles are the predicted Đij using multicomponent Darken equation (eq 22). The computed self-diffusivities from MD simulations are used to parametrize eq 22 and listed in Table S2 of the Supporting Information.

obtained by simple extrapolation to L → ∞. This approach allows for the use of much smaller systems. More details concerning this method can be found in the Supporting Information of this work. 2.3. Obtaining Thermodynamic Factors from COSMOSAC. An alternative method to predict the thermodynamic factors for a mixture is through a model for the activity coefficients γi.50 The COSMO-SAC theory is a quantum-mechanically based predictive model for the activity coefficients.51,52 In brief, the molecules are represented by a set of screening charge densities known as their σ-profiles. These profiles are generated by placing each molecule in a conducting or dielectric cavity and compute a charge distribution using density functional theory (DFT) which completely screens the molecular charge distribution from the outside. The resulting σ-profiles represent “fingerprints” of each of the molecules which can be used to model a solution by mixing them in their appropriate mole fractions. Since COSMO-SAC is an activity coefficient model, the dispersion is assumed to cancel in the reference fluid. The outcome of this modeling is the activity coefficients of each species as a function of their concentration. The thermodynamic factors Γij then follow from numerical differentiation using eqs 6 and 7. This predictive excess Gibbs energy model is used as a benchmark and reference in this work.

of the simulation box and is less sensitive to cutoff effects than the traditional method of eq 13. An empirical method for extrapolating the partial volume, compressibility, and activity coefficients for binary mixtures through the use of seven parameters was also presented. Mukherji et al. obtained the KB coefficients by coupling a small all-atom region to a very large coarse-grained reservoir. Exchanging molecules between the small system and the reservoir is allowed via a hybrid region. The difficulties of this approach are that the system should be very large and the thermodynamic force which ensures the thermodynamic equilibrium over the whole system should be calculated for each concentration.45 The simplest and in the past frequently used approach is to simply use a switching function to force the radial distribution function g(r) to converge to 1 for large distance r. However, it turns out that the final result depends on the choice of the switching function.45,48 Recently, we found that the average particle fluctuations in eq 12 can also be calculated by considering a large system in which smaller subvolumes are embedded. These subvolumes can exchange particles and energy with the simulation box and therefore a subvolume can be considered as a system in the grand-canonical ensemble. However, corrections should be taken due to the finite size effect originating from the boundaries of the subvolumes:32,35,49 Gij(L) = Gij∞ +

(constant) L

3. SIMULATION DETAILS In the ternary mixture chloroform−acetone−methanol, the OPLS force field was used for acetone and methanol.53,54 A nonpolarizable five-site model for chloroform was used.55 All components are treated as rigid molecules. The LJ potentials

(21)

In this equation, L is the linear length of subvolume in one dimension. The KB coefficient in the thermodynamic limit (G∞ ij ) can thus be 10253



Article

Figure 12. Comparison of the computed MS diffusivities Đij and the predicted MS diffusivities Đij in the ternary system chloroform (1)−acetone (2)− methanol (3) at 298 K, 1 atm. x2 varies while keeping x1 = x3. Circles are the computed Đij using MD simulations. Triangles are the predicted Đij using multicomponent Darken equation (eq 22). The computed self-diffusivities from MD simulations are used to parametrize eq 22 and listed in Table S2 of the Supporting Information.

molecules used in these simulations was typically 7000. Simulation runs of at least 5 ns were needed to obtain accurate values of the thermodynamic factors. The error bars of the computed thermodynamic factors are less than 3%. The COSMO-SAC computations were performed with the 2011.03 release ADF Theoretical Chemistry software.57 The COSMO surfaces were calculated using BP functionals and the TZP bases.57 The standard values of the parameters for the model were used in all calculations.52 To obtain accurate numerical derivatives, calculations for the logarithmic activities of each mixture were performed at mole fraction intervals of at most 0.05 in each component. The unconstrained derivatives from eq 7 were then calculated in MATLAB directly from these activities using a central differencing scheme. In the case of the ternary system, the derivatives were computed using a barycentric coordinate system. The constrained derivatives were then computed using eq 6.

describe the intermolecular nonbonded interactions which are truncated and shifted at 10 Å. The Lorentz−Berthelot mixing rules are applied to obtain the LJ interaction between unlike atoms.37 Electrostatic interactions are handled by Ewald summation using a relative precision56 of 10−5. The force field parameters as well as the parameters defining the geometries of the molecules can be found in the Supporting Information of this work. Three dimensional periodic boundary conditions consistent with a cubic box were applied to obtain properties corresponding to bulk systems. MD simulations were carried out at a temperature of 298 K and a pressure of 1 atm. Systems were first equilibrated in an NpT ensemble with a Nosé-Hoover thermostat and barostat using the time constants of 0.1 and 1 ps, respectively. The equations of motion were integrated using the leapfrog Verlet algorithm with a time step of 1 fs. The MS diffusion coefficients are obtained from equilibrium MD simulations in the NVT ensemble using the Nosé-Hoover thermostat at a density corresponding to a pressure of 1 atm.36 The box sizes were typically around (28 Å)3 resulting in a total number of molecules of the order of 200. The simulations for extracting diffusion coefficients were run for at least 100 ns to obtain MS diffusivities with an accuracy of around 5%. KB coefficients required for the calculation of thermodynamic factors were obtained from MD simulations in the NpT ensemble. Temperature and pressure were controlled using the Nosé-Hoover thermostat and barostat, respectively. In these simulations, the box size needs to be larger, see also section 2.2. The box size was typically fluctuating around (85 Å)3. The maximum number of

4. RESULT AND DISCUSSION To study the diffusion in the ternary mixture chloroform−acetone− methanol, three binary mixtures are relevant: chloroform−acetone, chloroform−methanol, and acetone−methanol. As ternary experimental data of diffusion coefficients are lacking for the ternary system chloroform−methanol−acetone, we validated our method in the relevant binary mixtures and predicted the ternary diffusion coefficients assuming that the quality of our method in ternary mixtures is similar to that of binary systems. The binary system acetone−methanol has been studied in our previous 10254



Article

Figure 13. Comparison of the computed MS diffusivities Đij and the predicted MS diffusivities Đij in the ternary system chloroform (1)−acetone (2)− methanol (3) at 298 K, 1 atm. x3 varies while keeping x1 = x2. Circles are the computed Đij using MD simulations. Triangles are the predicted Đij using multicomponent Darken equation eq 22. The computed self-diffusivities from MD simulations are used to parametrize (eq 22) and listed in Table S2 of the Supporting Information.

fusivity Đ12. This is reasonable since the computed mixture densities are lower than in the experiments and larger values of diffusion coefficients are thus expected. 4.2. Chloroform−Methanol. In the binary mixture chloroform−methanol, the computed mixture densities agree very well with the experiments;60 see Figure 4. Figure 5 shows thermodynamic factor Γ11 as a function of composition for the same binary mixture. COSMO-SAC accurately predicts Γ11 when the concentration of chloroform is low. As the concentration of chloroform increases, Γ11 computed using MD simulations is closer to the experimental data.59 Figure 6 shows the computed MS diffusivity Đ12 and Fick diffusivity D11 from MD simulations. The Fick diffusivity D11 is calculated using the computed MS diffusivity Đ12 and the thermodynamic factor Γ11. We observe that the MS diffusivity Đ12 increases as the concentration of chloroform increases. However, the increasing chloroform concentration drives the mixture away from ideal mixing behavior resulting in a reduced concentration dependence of Fick diffusivity D11. The Fick diffusivity D11 has a minimum value around a concentration of x1 = 0.7, while the MS diffusivity Đ12 displays a maximum in the same concentration range. The MS diffusivity Đ12 ranges from 3 × 10−9 to 6 × 10−9 m2·s−1 and the Fick diffusivity D11 ranges from 1.5 × 10−9 to 2.5 × 10−9 m2·s−1. 4.3. Chloroform−Acetone−Methanol. Figure 7 shows the computed densities in chloroform−acetone−methanol mixtures at 298 K, 1 atm. We observe a linear relation between mixture densities and composition. Figure 8−10 show the thermodynamic factors Γij in the ternary mixtures as a function of the

work.33,34 The computed Fick diffusivity D11 and thermodynamic factor Γ11 quantitatively agree with the experiments.33,34 4.1. Chloroform−Acetone. Karr et al. measured the densities of chloroform−acetone mixtures at 298 K, 1 atm.58 In Figure 1, the computed densities using MD simulations are compared to the experiments as a function of the composition. The computed mixture densities show a slightly stronger dependence on mixture composition suggesting a lower density of pure acetone and a higher density of pure chloroform. This is also shown in Table 1 in which computed pure-component densities are compared to the experimental ones. Figure 2 shows the thermodynamic factor Γ11 obtained from experiments and MD simulation in the binary mixture chloroform− acetone for various compositions. In ref 59, the vapor−liquid equilibrium (VLE) data are provided from which the activity coefficients γi can be calculated. By fitting these activity coefficients γi using an excess Gibbs energy model, i.e. the NRTL model, the thermodynamic factor Γ11 was obtained. Compared to the experiment, Γ11 obtained from COSMO-SAC shows good agreement. The computed Γ11 using MD simulations is somewhat lower. Using the computed MS diffusivity Đ12 and the thermodynamic factor Γ11 from MD simulations, we calculate the Fick diffusivity D11 as a function of composition, see Figure 3. A good agreement of Fick diffusivity D11 between experiments and MD simulations is seen. The very good agreement is partially due to error cancellation: the computed value of Γ11 is lower than in experiments implying larger values of the computed MS dif10255



Article

MS diffusivities Đij can not be measured directly in experiments. However, it is possible to predict MS diffusivities Đij using easily obtained self-diffusivities Di,self. For this purpose, we have proposed a multicomponent Darken model.11 In this model, we neglect the velocity cross-correlations between unlike molecules.11 The resulting multicomponent Darken equation is n

Đij = Di ,self Dj ,self

∑ i=1

xi Di ,self

(22)

in which Di,self is the self-diffusivity of component i in the mixtures. The computed self-diffusivities Di,self in the ternary system chloroform−acetone−methanol can be found in the Supporting Information of this work. The multicomponent Darken equation has been shown to allow for the most robust and accurate prediction of diffusion coefficients in ideal diffusing mixtures where the neglect of velocity cross-correlations is justified.10,11 Figures 11−13 show the computed MS diffusivities Đij and the predicted MS diffusivities Đij. From these figures, it can be seen that the concentration dependence of MS diffusivities Đij is not well captured by the multicomponent Darken equation. The MS diffusivities are either underestimated or overestimated by the multicomponent Darken equation suggesting that the correlations of different molecules are important in the chloroform−acetone−methanol system. The deviation of the predicted Đij using the multicomponent Darken equation from the computed Đij was already observed for the binary mixtures acetone−methanol and acetone−tetrachloromethane.33,34 The comparison to the multicomponent Darken equation demonstrates the potential of molecular simulations for the prediction of multicomponent Fick diffusion coefficients as the molecular models also allow the prediction of nonideal systems. Figure 14 shows the computed Fick diffusivities Dij using the computed MS diffusivities Đij and thermodynamic factors Γij as a function of composition. It can be observed that: (1) the diagonal Fick diffusivities are always positive and the off-diagonal Fick diffusivities may be negative; (2) the diagonal Fick diffusivities are about 1 order of magnitude larger than the off-diagonal Fick diffusivities suggesting the diffusion flux of component i mainly depends on its own concentration gradient while the concentration gradients of other components play a minor role. This behavior is in accordance with the bound placed on off-diagonal coefficients by the entropy production for ternary diffusion.4 It is important to note that the values of [Γ] and [D] depend on the labeling of the components. This problem does not apply for MS diffusivities. The computed ternary Fick diffusivities can be considered as true predictions as they directly originate from the force field.

Figure 14. Fick diffusivities Dij in the ternary system chloroform (1)− acetone (2)−methanol (3) at 298 K, 1 atm. Fick diffusivities Dij are calculated using the computed MS diffusivities Đij and thermodynamic factors Γij: (a) x1 varies while keeping x2 = x3; (b) x2 varies while keeping x1 = x3; (c) x3 varies while keeping x1 = x2.

5. CONCLUSION In this work, we present a consistent method for computing ternary Fick diffusivities from equilibrium MD simulations. For this purpose, MS diffusivities and thermodynamic factors are computed to calculate the matrix of Fick diffusivities. Our approach is applied to a ternary mixture chloroform−acetone− methanol. Even though a simple molecular model is used, the computed thermodynamic factors [Γ] are in close agreement with experiments. This findings suggests that the fluctuation method is well-suited for computing thermodynamic properties of mixtures in general since the framework relies on the Kirkwood−Buff coefficients which become now accessible in practical computations. Validation data for diffusion is only

composition. Data are reported for Γij as a function of xi while keeping xj = xk (with i ≠ j ≠ k). Oracz et al. measured the vapor− liquid equilibrium for this ternary mixtures at 303 K.61 This VLE data were converted to the matrix of thermodynamic factors Γij using the NRTL model.1,18 The ternary VLE data are very well described by the NRTL model. The computed Γij using MD simulation show quantitative agreement with the experimental data and the results from COSMO-SAC deviate from the experiments. We feel that the agreement is very good considering the fact that simple molecular interaction models were used in the MD simulations, i.e. polarization effects are not accounted for. 10256



Article

(9) Liu, X.; Vlugt, T. J. H.; Bardow, A. Maxwell-Stefan diffusivities in liquid mixtures: using Molecular Dynamics for testing model predictions. Fluid Phase Equilib. 2011, 301, 110−117. (10) Liu, X.; Bardow, A.; Vlugt, T. J. H. Multicomponent MaxwellStefan diffusivities at infinite dilution. Ind. Eng. Chem. Res. 2011, 50, 4776−4782. (11) Liu, X.; Vlugt, T. J. H.; Bardow, A. Predictive Darken equation for Maxwell-Stefan diffusivities in multicomponent systems. Ind. Eng. Chem. Res. 2011, 50, 10350−10358. (12) Liu, X.; Vlugt, T. J. H.; Bardow, A. Maxwell-Stefan diffusivities in binary mixtures of ionic liquids with DMSO and H2O. J. Phys. Chem. B 2011, 115, 8506−8517. (13) Schoen, M.; Hoheisel, C. The shear viscosity of a Lennard-Jones fluid calculated by equilibrium Molecular Dynamics. Mol. Phys. 1985, 56, 653−672. (14) Fernandez, G. A.; Vrabec, J.; Hasse, H. Self-diffusion and binary Maxwell-Stefan diffusion in simple fluids with the Green-Kubo method. Int. J. Thermophys. 2004, 25, 175−186. (15) Guevara-Carrion, G.; Vrabec, J.; Hasse, H. Prediction of selfdiffusion coefficient and shear viscosity of water and its binary mixtures with methanol and ethanol by molecular simulations. J. Chem. Phys. 2011, 134, 074508. (16) Guevara-Carrion, G.; Vrabec, J.; Hasse, H. On the prediction of transport properties of monomethylamine, dimethylether and hydrogen chloride by molecular simulation. Fluid Phase Equilib. 2012, 316, 46−54. (17) Guevara-Carrion, G.; Vrabec, J.; Hasse, H. Prediction of transport properties of liquid ammonia and its mixture with methanol by molecular simulation. Int. J. Thermophys 2012, 33, 449−468. (18) Taylor, R.; Kooijman, H. A. Composition derivatives of activity coefficient models for the estimation of termodynamic factors in diffusion. Chem. Eng. Commun. 1991, 102, 87−106. (19) Wheeler, D. R.; Newman, J. Molecular Dynamics simulations of multicomponent diffusion. 2. nonequilibrium method. J. Phys. Chem. B 2004, 108, 18362−18367. (20) Wheeler, D. R.; Newman, J. Molecular Dynamics simulations of multicomponent diffusion. 1. equilibrium method. J. Phys. Chem. B 2004, 108, 18353−18361. (21) Keffer, D. J.; Adhangale, A. The composition dependence of self and transport diffusivities from Molecular Dynamics simulations. Chem. Eng. J 2004, 100, 51−69. (22) Keffer, D. J.; Edwards, B. J.; Adhangale, P. Determination of statistically reliable transport diffusivities from Molecular Dynamics simulation. J. Non-Newtonian Fluid Mech. 2004, 120, 41−53. (23) Bardow, A.; Göke, V.; Koß, H. J.; Marquardt, W. Ternary diffusivities by model-based analysis of Raman spectroscopy measurements. AIChE J. 2006, 52, 4004−4015. (24) Medvedev, O. O.; Shapiro, A. A. Modeling diffusion coefficients in binary mixtures. Fluid Phase Equilib. 2004, 225, 13−22. (25) Poling, B.; Prausnitz, J.; O’Connell, J. The properties of gases and liquids, 5th ed.; McGraw-Hill: New York, 2004. (26) Medvedev, O. O.; Shapiro, A. A. Verifying reciprocal relations for experimental diffusion coefficients in multicomponent mixtures. Fluid Phase Equilib. 2003, 208, 291−301. (27) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. Transport diffusivity of methane in silicalite from equilibrium and nonequilibrium simulations. J. Phys. Chem. 1993, 97, 4173−4181. (28) Tsige, M.; Grest, G. S. Molecular Dynamics simulation of solvent polymer interdiffusion: Fickian diffusion. J. Chem. Phys. 2004, 120, 2989−2995. (29) Tsige, M.; Grest, G. S. Interdiffusion of solvent into glassy polymer films: a Molecular Dynamics study. J. Chem. Phys. 2004, 121, 7513−7519. (30) Keffer, D. J.; Gap, C. Y.; Edwards, B. J. On the relationship between Fickian diffusivities at the continuum and molecular levels. J. Phys. Chem. B 2005, 109, 5279−5288. (31) Zabala, D.; Nieto-Draghi, C.; de Hemptinne, J. C.; López de Ramos, A. L. Diffusion coefficients in CO2/n-alkane binary liquid mixtures by molecular simulation. J. Phys. Chem. B 2008, 112, 16610− 16618.

available for two binary subsystems. Here, MD results and experiments do agree well. Therefore, we expect that the computed Fick diffusivities should also be comparable with experiments. The presented approach allows for an efficient and consistent prediction of multicomponent Fick diffusion coefficients from molecular models.

■

ASSOCIATED CONTENT

S Supporting Information *

We provide (1) the force field parameters as well as the parameters defining the geometries of the molecules simulated in this work; (2) the integrals of radial distribution functions over volume as a function of the system size; (3) details concerning the fluctuation method to compute KB coefficients; (4) the selfdiffusivities in the ternary mixture chloroform−acetone− methanol at 298 K, 1 atm; (5) computed Kirkwood−Buff coefficients as a function of the cutoff radius used in the simulation. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was performed as part of the Cluster of Excellence “Tailor-Made Fuels from Biomass”, which is funded by the Excellence Initiative by the German federal and state governments to promote science and research at German universities. T.J.H.V., S.K.S., and S.K. acknowledge financial support from NWO-CW through an ECHO grant. This work was also sponsored by the Stichting Nationale Computerfaciliteiten (National Computing Facilities Foundation, NCF) for the use of supercomputing facilities, with financial support from the Nederlandse Organisatie voor Wetenschappelijk onderzoek (Netherlands Organization for Scientific Research, NWO).

■

REFERENCES

(1) Taylor, R.; Krishna, R. Multicomponent mass transfer, 1st ed.; Wiley: New York, NY, 1993. (2) Krishna, R.; Wesselingh, J. A. The Maxwell-Stefan approach to mass transfer. Chem. Eng. Sci. 1997, 52, 861−911. (3) Bardow, A.; Kriesten, E.; Voda, M. A.; Casanova, F.; Blümich, B.; Marquardt, W. Prediction of multicomponent mutual diffusion in liquids: model discrimination using NMR data. Fluid Phase Equilib. 2009, 278, 27−35. (4) Kjelstrup, S.; Bedeaux, D. Non-equilibrium thermodynamics of heterogeneous systems, 1st ed.; World Science Publishing Co. Pte. Ltd.: Singapore, 2008. (5) Kuiken, G. D. C. Thermodynamics of irreversible processes: applications to diffusion and rheology, 1st ed.; Wiley: England, 1994. (6) van de Ven-Lucassen, I. M. J. J.; Vlugt, T. J. H.; van der Zanden, A. J. J.; Kerkhof, P. J. A. M. Using Molecular Dynamics to obtain MaxwellStefan diffusion coefficients in liquid systems. Mol. Phys. 1998, 94, 495− 503. (7) Krishna, R.; van Baten, J. M. The Darken relation for multicomponent diffusion in liquid mixtures of linear alkanes: an investigation using Molecular Dynamics (MD) simulations. Ind. Eng. Chem. Res. 2005, 44, 6939−6947. (8) Fernandez, G. A.; Vrabec, J.; Hasse, H. Self-Diffusion and binary Maxwell-Stefan diffusion coefficients of quadrupolar real fluids from molecular simulation. Int. J. Thermophys. 2005, 26, 1389−1407. 10257



Article

(56) Smith, W.; Forester, T. R.; Todorov, I. T.; Leslie, M. The DLPOLY2 user manual; 2006, ftp://ftp.dl.ac.uk. (57) ADF 2011.03, SCM, Theoretical Chemistry; Vrije Universiteit, Amsterdam, The Netherlands, 2011. (58) Karr, A. E.; Bowes, W. M.; Scheibel, E. G. Analysis of binary and ternary mixtures. Anal. Chem. 1951, 23, 459−463. (59) Gmehling, J.; Onken, U. Vapor-liquid equilibrium data collection, 1st ed.; DECHEMA: Germany, 2005. (60) Wei, I. C.; Rowley, R. L. Binary liquid mixture viscosities and densities. J. Chem. Eng. Data 1984, 29, 332−335. (61) Oracz, P.; Warycha, S. Vapour-liquid equilibria. The ternary system methanol-chloroform−acetone at 303.15K. Fluid Phase Equilib. 1997, 137, 149−162. (62) Campbell, A. N.; Kartzmark, E. M. Thermodynamic and other properties of methanol + acetone, carbon disulphide + acetone, carbon disulphide + methanol, and carbon disulphide + methanol + acetone. J. Chem. Thermodyn. 1973, 5, 163−172. (63) Noda, K.; Ohashi, M.; Ishlda, K. Viscosities and densities at 298.15K for mixtures of methanol, acetone and water. J. Chem. Eng. Data 1982, 27, 326−328. (64) Timmermans, J. The physico-chemical constants of binary systems in concentrated solutions, 1st ed.; Interscience Publishers: New York, 1959.

(32) Schnell, S. K.; Liu, X.; Simon, J. M.; Bardow, A.; Bedeaux, D.; Vlugt, T. J. H.; Kjelstrup, S. Calculating thermodynamic properties from fluctuations at small scales. J. Phys. Chem. B 2011, 115, 10911−10918. (33) Liu, X.; Schnell, S. K.; Simon, J. M.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Fick diffusion coefficients of liquid mixtures directly obtained from equilibrium Molecular Dynamics. J. Phys. Chem. B 2011, 115, 12921−12929. (34) Liu, X.; Schnell, S. K.; Simon, J. M.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Correction to: Fick diffusion coefficients of liquid mixtures directly obtained from equilibrium Molecular Dynamics. J. Phys. Chem. B 2012, 116, 6070. (35) Schnell, S. K.; Vlugt, T. J. H.; Simon, J. M.; Bedeaux, D.; Kjelstrup, S. Thermodynamics of a small system in a μT resorvior. Chem. Phys. Lett. 2011, 504, 199−201. (36) Frenkel, D.; Smit, B. Understanding molecular simulation: from algorithms to applications, 2nd ed.; Academic Press: San Diego, 2002. (37) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids, 1st ed.; Oxford University Press: New York, 1987. (38) Rapaport, D. The art of Molecular Dynamics simulation, 2nd ed.; Cambridge University Press: Cambridge, 2004. (39) Dubbeldam, D.; Ford, D. C.; Ellis, D. E.; Snurr, R. Q. A new perspective on the order-n algorithm for computing correlation functions. Mol. Sim 2009, 35, 1084−1097. (40) Kirkwood, J. G.; Buff, F. P. The statistical mechanical theory of solution. J. Chem. Phys. 1951, 19, 774−777. (41) Ben-Naim, A. Molecular theory of solutions, 2nd ed.; Oxford university press: Oxford, 2006. (42) Shi, W.; Maginn, E. J. Continuous fractional component Monte Carlo: an adaptive biasing method for open system atomistic Simulations. J. Chem. Theory Comput. 2007, 3, 1451−1463. (43) Shi, W.; Maginn, E. J. Improvement in molecule exchange efficiency in Gibbs ensemble Monte Carlo: development and implementation of the continuous fractional component move. J. Comput. Chem. 2008, 29, 2520−2530. (44) Ruckenstein, E.; Shulgin, I. Entrainer effect in supercritical mixtures. Fluid Phase Equilib. 2001, 180, 345−359. (45) Mukherji, D.; van der Vegt, N. F. A.; Kremer, K.; Site, L. D. Kirkwood-Buff analysis of liquid mixtures in an open boundary simulation. J. Chem. Theory Comput. 2012, 8, 375−379. (46) Wedberg, R.; Connell, J. P. O.; Peters, G. H.; Abildskov, J. Pair correction function integrals: computation and use. J. Chem. Phys. 2011, 135, 084113. (47) Nichols, J. W.; Moore, S. G.; Wheeler, D. R. Improved implementation of Kirkwood-Buff solution theory in periodic molecular simulations. Phys. Rev. E 2009, 80, 051203. (48) Perera, A.; Zoranic, L.; Sokolic, F.; Mazighi, R. A comparative Molecular Dynamics study of water-methanol and acetone-methanol mixtures. J. Mol. Liq. 2011, 159, 52−59. (49) Schnell, S. K.; Vlugt, T. J. H.; Simon, J. M.; Bedeaux, D.; Kjelstrup, S. Thermodynamics of small systems embedded in a reservoir: a detailed analysis of finite size effects. Mol. Phys. 2012, 110, 1069−1079. (50) Taylor, R.; Kooijman, H. A. Composition derivative of activity coefficient models for the estimation of thermodynamic factors in diffusion. Chem. Eng. Commun. 1991, 102, 87−106. (51) Lin, S. T.; Sandler, S. I. A priori phase equilibrium prediction from a segment contribution solvation model. Ind. Eng. Chem. Res. 2002, 41, 899−913. (52) Hsieh, C. M.; Sandler, S. I.; Lin, S. T. Improvements of COSMOSAC for vapor-liquid and liquid-liquid equilibrium predictions. Fluid Phase Equilib. 2010, 297, 90−97. (53) Jorgensen, W. L. Optimized intermolecular potential functions for liquid alcohols. J. Phys. Chem. 1986, 90, 1276−1284. (54) Tummala, N. R.; Striolo, A. Hydrogen-bond dynamics for water confined in carbon tetrachloride-acetone mixtures. J. Phys. Chem. B 2008, 112, 10675−10683. (55) Gupta, R.; Chandra, A. Structural, single-particle and pair dynamical properties of acetone-chloroform mixtures with dissolved solutes. Chem. Phys. 2011, 383, 41−49. 10258


Fick Diffusion Coefficients in Ternary Liquid Systems from Equilibrium

Recommend Documents