# «David J. Keffer and Parag Adhangale Department of Chemical Engineering University of Tennessee 1512 Middle Drive Knoxville, TN 37996-2200 ...»

The Composition Dependence of Self and Transport Diffusivities

from Molecular Dynamics Simulations

David J. Keffer and Parag Adhangale

Department of Chemical Engineering

University of Tennessee

1512 Middle Drive

Knoxville, TN 37996-2200

dkeffer@utk.edu

Abstract

Using molecular dynamics simulations we determine the composition dependence of the

self-diffusivity and transport diffusivity of a methane/ethane mixture at high pressure. The

transport diffusivity is generated from the self-diffusivities using the Darken equation. While obtaining the transport diffusivity via the Darken equation is an approximation, it has great computational advantage over obtaining it from phenomenological coefficients, in terms of both speed and statistical accuracy. We perform a careful analysis of the molecular dynamics simulations and show that it is possible to reproduce the results in the microcanonical, canonical, and isobaric-isothermal ensembles. We demonstrate that in order to capture the sensitive dependence of the diffusivities on composition, it is necessary to run simulations with larger systems and for longer durations than is typical. We report the trends in the diffusivities as a function of composition, temperature, pressure, and density. We modify an existing empirical correlation, which when combined with a corresponding states chart, is capable of quantitatively reproducing the simulated diffusivity dependence on composition, temperature, pressure, and density. Finally, we quantify the effect that the choice of equation of state used to evaluate the thermodynamic factor in the Darken equation has on the transport diffusivity.

## I. INTRODUCTION

The estimation of diffusivities using theoretical, experimental, and simulation approaches has been of practical interest for many decades. Kinetic theory can predict the self-diffusivity of pure (single-component) dilute gases quantitatively.1 It can also predict the self-diffusivity of pure dense fluids quantitatively, given that values of the collision integral are available. Many empirical correlations to determine the self-diffusivity as a function of temperature, pressure, and molecular weight have also been developed.2,3 Additionally, plots based on corresponding states, which are fitted to experimental and simulation data, can predict self-diffusivities of dense fluids semi-quantitatively.3 However, for many systems there are no theories or correlations that can accurately predict the self-diffusivity. Frequently, the best and sometimes the only recourse is to turn to experiment. Typically in the determination of self-diffusivities, the experiments fall into two categories: laboratory experiments and computer simulations. We classify computer simulations as a type of experiment, since the only theory that necessarily feeds into the simulations is Newton’s second law.Before we can continue our introduction, we must distinguish between self-diffusivities and transport diffusivities. We provide a qualitative description here and a mathematical definition later. A self-diffusivity is a measure of the mobility of fluid molecules in the absence of a driving force for diffusion, e.g. a composition gradient. Self-diffusivities arise due to the Brownian (random-walk) motion of the molecules. Self-diffusivities are defined for pure fluids as well as for each component in a mixture.

Transport diffusivities, also known as Fick or Fickian diffusivities, appear in Fick’s law.

They relate a mass or mole flux to a driving force, such as a gradient in the concentration of a particular species. Transport diffusivities are typically not defined for pure fluids. For a multicomponent fluid containing nc components, the transport diffusivity is typically represented as an ncxnc matrix, although not all of the elements of this matrix are independent.

Transformations of this matrix are possible, generally to an (nc-1)x(nc-1) matrix, where each element of the upper triangular component of the matrix is independent. Transport diffusivities arise due to gradients in the system.

When one chooses to evaluate Fick’s law, one requires transport diffusivities rather than self-diffusivities. In addition to self-diffusivities, kinetic theory can also predict transport diffusivities. It is easy to use kinetic theory to predict the transport diffusivities of low-pressure gases and substantially more difficult to predict them for dense fluids. Empirical correlations and corresponding-states plots, while intended for self-diffusivities of pure fluids, have also been modified to estimate transport diffusivities, at least for binary fluids. However, these modifications do not predict the transport diffusivity (or the self-diffusivity) as a function of the composition of the mixture. In other words, there is no dependence of the self-diffusivities and transport diffusivity on mole fractions.

In this work, we employ molecular dynamics simulations to estimate the self-diffusivities and transport-diffusivities of a high pressure binary mixture of methane and ethane. We specifically determine the composition dependence of the diffusivities. We describe how existing correlations and charts can be easily modified to give quantitative agreement with simulation.

In Section II. Background, we provide a review of the pertinent literature. In Section III.

Simulation, we provide a detailed discussion of our simulation techniques, including our efforts to minimize statistical error. In Section IV. Results and Discussion, we present the results of our simulations. We also test the validity of the modified correlations we have developed. Finally, in the last section, we present our conclusions.

## II. BACKGROUND

II.A. Simulation II.A.1. Self-Diffusivities from Simulation It is a standard procedure to estimate self-diffusivities in pure components or mixtures using equilibrium molecular dynamics (MD) simulations.4-6 One can use either of two equivalent methods to obtain a self-diffusivity for each component in the mixture. The first method is to use Einstein’s relation, which relates the mean square displacement to the observation time viawhere Dself,α is the self-diffusivity of component α, d is the dimensionality of the system, τ is the observation time, and rα,i is the position of the ith molecule of component α. The brackets indicate an ensemble average over time, t, and over all molecules of component α. There is no ambiguity in defining the self-diffusivity. Moreover, there are no practical obstacles in the determination of a self-diffusivity from a molecular dynamics simulation, aside from the ubiquitous concerns that the simulation include a sufficient number of molecules and be of adequate duration so as to ensure acceptable statistical accuracy of the results.

A second method to determine the self-diffusivity is to use the velocity auto-correlation function in a Green-Kubo relation.4-6

where vα,i is the velocity of the ith molecule of component α. Again, there is no practical impediment to implementation, except to be aware that this integral contains a long-time tail.

II.A.2. Transport Diffusivities from Simulation There are several ways that one can measure the transport diffusivity of a mixture in an MD simulation. First, one can perform a non-equilibrium MD simulation. For example, Heffelfinger et al. developed and used the Dual Control-Volume Grand-Canonical Molecular Dynamics (DCV-GCMD) simulation to measure transport diffusivities.7 Using this technique one established a gradient in chemical potential. One measures the flux. One can then use Fick’s law to back out the transport diffusivity. Maginn et al. also used non-equilibrium MD simulations to generate transport diffusivities.8 The drawback to this method is its computational expense.

A second method to determine the transport diffusivity is to run equilibrium MD simulations, as described above, and to invoke the Darken equation.9 The Darken equation can be written as

where kB is Boltzmann’s constant, T is temperature, aα is the activity of component α, and xα is the mole fraction. The derivation of the Darken equation is available elsewhere.9 Although the Darken equation was published in 1948 to describe diffusion in binary alloys, it has been widely used for a variety of systems in the subsequent half century. There are numerous approximations involved in the Darken equation and there exist decades of literature attempting to quantify the magnitude of error resulting from those approximations.10,11 Carman summarizes his findings, “…the Darken type of equation and probably any equation based upon a simple premise cannot be expected to give a highly accurate correlation with experimental data in most systems.”10 Nevertheless, the Darken equation is still widely used in both experimental and simulated systems to compute transport diffusivities from self-diffusivities.12 The appeal of the Darken equation is that it requires only self-diffusivities, which can be obtained in an unambiguous manner and with some statistical significance from a single equilibrium MD simulation.

A third method of obtaining the transport diffusivities is to adopt the formalism of linear irreversible thermodynamics (LIT), in which one can relate transport diffusivities to phenomenological coefficients. These phenomenological coefficients in turn can be obtained from equilibrium MD simulations.8,13,14 The paper of Sanborn and Snurr is particularly cited as being especially revealing and reader-friendly.14 Again the phenomenological coefficients can be obtained from equivalent Green-Kubo or Einstein-like relations. The disadvantage with this method lies in the fact that the phenomenological coefficients display a far greater susceptibility to statistical noise than do the self-diffusivities. Therefore, in order to obtain a statistically reliable diffusivity from the phenomenological coefficients, one must average the results of many simulations. Previous researchers have used an average over 10 to 14 14.

Sometimes the transport diffusivity is called the “mutual diffusion coefficient”.15,16 This mutual diffusion coefficient is a transport diffusivity for a binary system relative to a local center of volume frame of reference. It effectively is using linear irreversible thermodynamics to obtain the transport diffusivity. As such, numerous simulations are required to obtain a reasonable standard deviation. Jolly and Bearman averaged 27 estimates.15 Aside from these methods listed above for obtaining the transport diffusivity, there are other techniques. For example, Ali et al. employed Mode Coupling Theory to estimate transport diffusivities of binary Lennard-Jones systems.17 The full advantage of this approach has not yet been explored.

Of these methods, the Darken equation is the least rigorous, but has the advantage that it requires the least computation time because self-diffusivities take less simulation effort to compute than transport diffusivities. The reason lies in the fact that in a simulation of N molecules, each molecule contributes independent information to the self-diffusivity at every time step, effectively providing N data points to which the self-diffusivity will be fit. However, the transport diffusivities examine system properties, like center-of-mass position, which provide only 1 piece of information per time step. To obtain the transport diffusivity from phenomenological coefficients to the same degree of accuracy as one can obtain self-diffusivities from a single simulation would require N simulations. In this work, we use N=10,000 particles per simulation. Clearly, one cannot perform 10,000 simulations per data point. In this work we determined transport diffusivities from (i) the self-diffusivity via the Darken relation and (ii) the phenomenological coefficients via LIT.

Regardless of whether one uses the Darken equation or phenomenological coefficients to obtain the transport diffusivity, one must still calculate the same thermodynamic factor. There are a variety of ways to do this. Here we discuss two methods. First, if one has an equation of state that describes the simulated fluid relatively well, then one can simply use the equation of state to evaluate the activity and the necessary partial derivative of the activity, which is the thermodynamic factor in the Darken equation. If one chooses this approach, then one does not require additional simulation time and maintains the advantage of the Darken equation.

Alternatively one can use molecular simulations to determine the thermodynamic factor.

For example, there are methods to determine the chemical potential of a molecular-level simulation. One common method is called Widom’s particle insertion method.4,18 In this method, one attempts to regularly insert particles into the simulation. Depending upon the ease with which the particles are accepted, one can determine the chemical potential. This method only works well for gases. As the fluid becomes more dense, the acceptance probability is so low that the method typically becomes computationally impractical.

To obtain the derivative of the chemical potential with respect to mole fraction, as needed for the thermodynamic factor in Darken’s equation, one could then use a centered-finite difference formula. If the system of interest was at T, p, and mole fraction, xα, then one would simulate the system at T, p, and mole fraction, xα+∆x and at T, p, and mole fraction, xα-∆x.

Widom’s particle insertion method would be performed in all three simulations. Clearly these simulations would be most easily accomplished in the isobaric-isothermal ensemble. Other researchers used Grand Canonical Monte Carlo (GCMC) simulations to determine the thermodynamic factor.14 In this work, we are simulating bulk fluids. We will demonstrate that the Lennard-Jones equation of state predicts the simulated thermodynamic properties well.19 We will use the Lennard-Jones equation of state to determine the thermodynamic factor.

II.B. Theory, Empirical Correlations and Corresponding States Plots In this section, we make no attempt to completely review the rich and lengthy history of the theoretical estimation of diffusivities. Instead, we cite several specific theories and empirical correlations that we will test, modify, and compare with our simulation data.