Article pubs.acs.org/IECR
Nonequilibrium Linear Response Theory: Application to Onsager− Stefan−Maxwell Diffusion Charles W. Monroe,*,† Dean R. Wheeler,‡ and John Newman§ †
Department of Chemical Engineering, University of Michigan, Ann Arbor, Michigan 48109-2135, United States Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602, United States § Department of Chemical and Biomolecular Engineering, University of California, Berkeley, California 94720-1462, United States ‡
ABSTRACT: The original regression hypothesis of Onsager contains more physics than the fluctuation−dissipation theorem. A discrepancy in a prior derivation of the reciprocal relation for Stefan−Maxwell coefficients through fluctuation theory is resolved. Along with a demonstration of the correct reciprocal relation, it is shown how to use the regression hypothesis to analyze data from an equilibrium simulation and determine values of the transport coefficients.
■
INTRODUCTION Some years ago, Monroe and Newman published in this journal1 how to obtain a reciprocal relation for transport properties by applying Onsager’s regression hypothesis2,3 to a complete macroscopic description of the transport system and enforcing microscopic reversibility in the result. A generally asymmetric reciprocal relation for multicomponent diffusion was obtained, which still provided the expected number of relationships among transport coefficients.1 In contrast, Wheeler and Newman used methods of statistical mechanics to show that +ij = +ji for diffusion coefficients in the same Onsager− Stefan−Maxwell transport equation,4 reinforcing the conventional understanding of the reciprocal relation as a symmetry of transport coefficients. The resolution of this discrepancy, which is the primary purpose of this paper, is a correction of Monroe and Newman’s calculation of the initial values of fluctuation correlations. In addition, it is shown below how the theory using the regression hypothesis1 yields more general transport relations that express the consequences of the macroscopic transport theory and which could permit analysis of the results of molecular-dynamics simulations, particularly to predict the decay of fluctuations and to measure the macroscopic transport coefficients. Examination of differences between these simulations and the predictions of the macroscopic theory can elucidate processes that the macroscopic theory ignores, which may affect molecular simulations in the early stages. The simplest such process to conceive is the relaxation time for diffusion due to inertial effects neglected when the transport laws are taken to be instantaneous. After it is shown that the corrected results of the regression hypothesis are reliable, the two methods of analysis are discussed in more detail, along with philosophical comments about nonequilibrated systems that are close to equilibrium.
manifestations of the behavior of the same system, but the hypothesis requires some definition of how it is to be implemented. There may be existing at the same time an experimental system, a macroscopic description of it, and a molecular simulation of it, say by the Monte Carlo method or by molecular dynamics. The last ones are, of course, imperfect approximations of the first. One way to implement the regression hypothesis on the macroscopic level is to express the physics in a more-or-less complete mathematical description, including in the present case material balances and transport laws for an isotropic, isothermal, isobaric, multicomponent solution. These laws will be taken to have the form of the Onsager−Stefan−Maxwell equation, i.e., the extended Stefan−Maxwell equation that uses gradients of the component chemical potentials μi as diffusion driving forces.5 In addition to the balances and flux constitutive laws, there are various equations of state that express how properties relate to thermodynamic variables. These include how μi, +ij , and the total molar solution concentration cT depend on the temperature T, pressure p, and component mole fractions yi. State equations also express thermodynamic constraints, for instance, that the mole fractions sum to unity. Paper 11 develops this framework, with many restrictions, such as the absence of chemical reactions and no explicit account for molecular vibrations. All of the governing equations are applied to a one-dimensional slab of material with thickness L and are linearized about a particular homogeneous equilibrium state. (Equilibrium values are indicated with a superscript ∞.) The fundamental dependent variables, which in this case include the mole fractions and component velocities vi, are expanded spatially in Fourier series, or equivalently in Fourier transforms. These Fourier relations define αim(t) and
CONSEQUENCES OF THE REGRESSION HYPOTHESIS The regression hypothesis of Onsager is that microscopic fluctuations decay in the same manner as macroscopic variations. This is intuitively obvious because it describes two
Special Issue: Scott Fogler Festschrift
■
© XXXX American Chemical Society
Received: September 30, 2014 Revised: January 14, 2015 Accepted: January 14, 2015
A
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
fraction basis, γi, is defined by the constitutive law μi = μ0i + RT ln(γiyi), in which μ0i represents the chemical potential of pure component i at T and p. The Stefan−Maxwell transport matrix M has entries
βim(t) as time-dependent Fourier coefficients for the principal variables of mole fractions by ∞
yi(1) =
⎛ mπ z ⎞ ⎟ L∞ ⎠
∑ αim(t ) cos⎜⎝
m=1
(1)
Mij =
and component velocities by ∞
vi(1)
=
∑ m=1
⎛ mπ z ⎞ βim(t ) sin⎜ ∞ ⎟ ⎝ L ⎠
R ij =
(3)
A regular perturbation about an equilibrium state is performed by replacing all of the dependent variables in the governing system with forms like this. It is understood that when dependent variables are multiplied, the order of the product is (1) (1) (2) additive: for example, y∞ = O(1), y(1) = O(2), y(1) = i vj i vj i ∇μj (3) (1) (1) (1) (3) O , yi vj cT = O , etc. The equations of state are also replaced with their Taylor expansions about equilibrium to identify how changes in the fundamental state variables (T, p, and a set of independently variable yi) contribute at each order of approximation. Terms like y(1) solve the linearized system i obtained by discarding quantities of O(2) and higher from the governing system. Terms like yi(2) arise from the most significant nonlinearities neglected in the first-order analysis; they solve the system at O(2) obtained by discarding quantities of O(3) and higher while incorporating the solutions for all quantities at O(1). Like y(1) i , the higher-order functions can be expanded in Fourier series such as ∞
yi(2) =
⎛ mπ z ⎞ ⎟ L∞ ⎠
∑ αim(2)(t ) cos⎜⎝
m=1
and
j≠n
j≠n
(4)
dαjm dτ
(5)
To ensure that the bases for composition and excess flux are minimal and, consequently, that property matrices are invertible, paper 1 introduces the (n − 1) × (n − 1) matrices Q and R, which arise when applying thermodynamic and kinematic constraints. Thermodynamic quantities that convert the forces driving diffusion from −cTyi∇μi for chemical potentials to −∇yi for mole fractions are included in the nonsingular thermodynamic matrix Q, whose entries are Q ij = δij +
⎛ ∂ ln γ ⎞∞ i⎟ ⎟ ∂ y ⎝ j ⎠ T ,p,y
Mii = −∑ Mik
and
k≠i
yi∞ +∞ ij
−
yi∞ +∞ in
(7)
R ii = −
yi∞ +∞ in
if i ≠ j −
∑ k≠i
yk∞ +∞ ik
dC = −PC dτ
(8)
(9)
where entries of the correlation matrix C(τ) are defined following Onsager,1,2 as the ensemble averages Cil(τ ) = ⟨αim(τ0 + τ ) αlm(τ0)⟩
(10)
This is the basic equation, embodying the Onsager regression hypothesis, and applies both to equilibrium ensembles and to nonequilibrium ensembles relaxing toward equilibrium in the absence of external forces. Correlation-evolution equation (9) can be used in a general manner. If one can calculate in a simulation, equilibrium or nonequilibrium, both C and dC/dτ at the same time, one can calculate P from
yi∞⎜⎜
k≠j,n
if i ≠ j
These relations arise naturally when reducing the numbers of fluxes and driving forces from n to n − 1 by the method in paper 1. Equation 5 is the translation of the macroscopic description of the system into a form that shows explicitly how the n − 1 independent mole fractions vary. The time derivative comes from the component material balance. The wavelength of the disturbance is incorporated into the time variable τ = (mπ/ L∞)2t. Because the different wavelengths do not interfere in the linearized system, the form of τ can be ignored until one wants to deal with general disturbances. Macroscopic equation (5) makes concrete the linear transport system that is needed to initiate Casimir’s regression process6 for Onsager−Stefan−Maxwell diffusion. Define the transport matrix P = −R −1 Q. (The second law of thermodynamics assures that for a stable single-phase system all properties approach their equilibrium values. Consequently, P is diagonalizable, has positive eigenvalues, and is invertible.) For an isothermal, isobaric multicomponent diffusion system in which the independent mole fractions are adopted as the basis of intensive state variables, Casimir’s general transport matrix equals −P. To the best of our knowledge, that general matrix has not previously been expressed in terms of +ij and γi. Composition Correlations. The regression hypothesis says that microscopic systems decay in the same way as macroscopic ones. Here this is taken to mean that spatially averaged correlations of fluctuations evolve with respect to time by following an equation similar to eq 5:
(The superscript is appended to αim only for higher-order terms.) Paper 1 details the above manipulations for a diffusion system containing n components, with elimination of various (Fourier) variables until only the mole fractions remain. Thus, the dynamical laws of the linearized system can be encapsulated in the expression
∑ Q ijαjm = ∑ R ij
+ij
where R is the gas constant. Note that M is an n × n matrix and necessarily singular. Entries of the nonsingular transport matrix R derive from those of M:
(2)
In these expansions, z indicates the position within the slab. To express the formal linearization procedure in more detail, adopt the notation yi = yi∞ + yi(1) + yi(2) + O(3)
RTc Tyy i j
(6)
in which δij is the Kronecker delta, equal to 1 if i = j and 0 otherwise. The activity coefficient of component i on a moleB
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
P=−
dC −1 C dτ
C(∞) − C(0) = −P
(11)
dC −1 −1 C Q dτ
V(∞) − V(0) = −YPY −1
(12)
δij yi∞
+
mπ L∞
∑ Yij j≠n
dαjm dτ
∫0
∞
V (τ ) d τ
(13)
R−1QC0 = (R−1QC0)T
1 yn∞
P = C0 (
C dτ )−1
(23)
∫0
∞
V dτ )−1Y
(24)
The reciprocal relation can also be derived from either of these formulas. Algebraic manipulation of eq 23 gives P−1C0 =
∫0
∞
C(τ ) dτ
(25)
Because C is symmetric for an equilibrium ensemble, the integral is also symmetric, and therefore the quantity on the left is symmetric, according to
(16)
P−1C0 = (P−1C0)T = (C0)T P−T (17)
(26) T
T T
where the transpose identity (AB) = B A has been used. Multiplication through the equalities from the left by P and from the right by PT, followed by additional applications of the transpose identity, shows that
Substitution into eq 9 and algebraic manipulations show that (18)
PC0 = (PC0)T
−1
Because P and YPY are similar matrices, the decay constants for V should be the same as those for C. Rearrangement yields dV −1 V Y dτ
∞
P = Y −1V 0(
where
dV = −YPY −1V dτ
∫0
and eq 21 rearranges to
Note that the use of relative component velocities in eqs 13 and 15 leads to physical relationships that are invariant in moving reference frames. Formation of the velocity correlation V(τ), according to the regression hypothesis, gives
Vij(τ ) = ⟨(βim − βnm)|τ0+ τ (βjm − βnm)|τ0 ⟩
(22)
as shown in paper 1. The reciprocal relation is independent of the wavelength chosen. The integral relations are also simplified for equilibrium ensembles, providing simple equations that can be used to calculate P. Equation 20 becomes
(14)
⎛ mπ ⎞ 2 V = ⎜ ∞ ⎟ YPCPTYT ⎝L ⎠
(21)
EQUILIBRIUM ENSEMBLES For an ensemble under equilibrium conditions, the principle of microscopic reversibility requires that C and V be symmetric. The correlations also decay to zero as τ approaches ∞ according to eqs 9 and 18 because there are no external forces. Finally, the initial values are specified by a thermodynamic potential appropriate for describing the ambient conditions, such as the entropy or Gibbs energy of the system.1 These initial values for equilibrium ensembles are denoted C0 and V0 to distinguish them from the initial values for nonequilibrium systems and to emphasize the symmetry property. Necessarily under equilibrium conditions, eq 9 leads to the reciprocal relation
With this definition in hand, one can show that RTc∞ T R = M′Y, where M′ is the nonsingular matrix formed by deleting the nth row and the nth column of M. This identity substantiates the indicial expressions of Rij given in eq 7 and paper 1. From eq 5 and the definition of P, the Fourier-amplitude representation of the velocity difference is found to be mπ βim − βnm = ∞ ∑ ∑ YijPjkαkm L j≠n k≠n (15)
P = − Y −1
(20)
■
Here the sum of mole fractions has been used to eliminate αnm, and Y is a nonsingular symmetric matrix defined as
Yij =
C (τ ) d τ
Under nonequilibrium conditions, matrices C(0) and V(0) may be asymmetric, but C(∞) and V(∞) will still vanish in the absence of external forces.
Because the governing physics is based on macroscopic laws, these formulas, and similar results below, are best applied in the phenomenological period, after inertial and other early effects have ceased to be important. Velocity Correlations. Statistical mechanics often uses velocity correlations rather than the mole-fraction correlations described above. Intermediate equations in paper 1, coming from the material balances, relate the velocity and mole-fraction Fourier coefficients. Subtraction of the nth transformed material balance from the ith, so that the velocity of species i is referred to that of species n, gives βim − βnm = −
∞
Integration of eq 18 shows that for the velocity correlations
and consequently R−1 =
∫0
(27)
This is the reciprocal relation from eq 22, obtained in a more cumbersome manner by the integral relationship. Similar analysis of the velocity-correlation integral (24) establishes that YP−1Y−1V0 is symmetric. Substitution of eq 16 shows this symmetry to be equivalent to the reciprocal relation in eqs 22 and 27. Expansion of the Gibbs Function. The probabilities of the initial fluctuations are determined by the total Gibbs free energy in an isothermal, isobaric system. The Gibbs function G
(19)
This equation can be used to extract P or R from a simulation in a manner similar to the one involving C. Correlation Integrals. Molecular simulations also typically use a Green−Kubo integral of correlations over time. A similar equation results by integration of eq 9 from 0 to ∞: C
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research of a fluctuation state can be expressed in terms of local species chemical potentials, G=
∫V cT ∑ yi μi dV
G − G∞ =
i
= (28)
i
∫V cT ∑ yi (μi − μi∞) dV
Its value in a homogeneous equilibrium state, G∞, is
=
∫V
∞
c T ∑ yi (μi − μi∞) dV + Δ i
∫V ∑ cT∞yi∞μi(1) dV ∞
i
∞
G =
∑ niμi
∞
=
i
∫V
c∞ ∞ T
∑
yi∞μi∞
i
+
dV
+
∫V
∞
c T∞yi∞ dV =
∫V cTyi dV
∫V cT ∑ yi μi∞ dV i
∫V cT ∑ yi (μi − μi∞) dV i
∞
(34)
In an isothermal, isobaric system, the chemical potentials depend on the composition only. Taylor expansion about the equilibrium state and insertion of the perturbed mole-fraction distributions defined in eq 3 yield μi − μi
∞
⎛ ∂μ ⎞∞ = ∑ ⎜⎜ i ⎟⎟ yj(1) + j ≠ n ⎝ ∂yj ⎠
⎛ ∂μ ⎞∞ ∑ ⎜⎜ i ⎟⎟ yj(2) j ≠ n ⎝ ∂yj ⎠
⎛ ∂ 2μ ⎞∞ ∑ ∑ ⎜⎜ i ⎟⎟ yj(1)yk(1) + O(3) j ≠ n k ≠ n ⎝ ∂yj ∂yk ⎠
1 + 2
(35)
(30)
Also, the isothermal, isobaric Gibbs−Duhem equation requires that
∑ i
⎛ ∂μ ⎞∞ i⎟ ⎟ y ∂ ⎝ j ⎠T , p , y
yi∞⎜⎜
=0 (36)
l≠j,n
and ⎛ ∂ 2μ ⎞∞ ∑ yi∞⎜⎜ i ⎟⎟ ⎝ ∂yj ∂yk ⎠T , p , y i
(31)
Subtraction of this relationship from eq 28 yields an integral for the Gibbs-energy fluctuation associated with composition fluctuations, G − G∞ =
∫V ∑ cT∞yi∞μi(2) dV + O(3) i
Recognizing that the chemical potential of species i in the equilibrium state, μ∞ i , is also a constant, one can use eq 30 to replace the integral over V∞ in eq 29 with an integral over the fluctuating volume V, writing the equilibrium free energy instead as G∞ =
∞
i
(29)
where ni indicates the total number of moles of species i. Computation of the initial correlations requires expansion of the difference G − G∞ in terms of the fluctuation amplitudes, taking care to discard all terms of order (3) and higher. The discrepancy between the (asymmetric) reciprocal relation derived in paper 1 and the (symmetric) reciprocal relation derived by Wheeler and Newman4 is resolved by revisiting the expansion of G − G∞. Even in a fluctuation state away from equilibrium, ni is the integral of cTyi over the volume V and is constant, so that ni =
∫V ∑ (cT(1)yi∞ + cT∞yi(1))μi(1) dV
l≠j,k ,n
⎛ ∂(μ − μ ) ⎞∞ k n ⎟ = −⎜⎜ ⎟ y ∂ ⎝ ⎠T , p , y j
l≠j,n
(37)
The first sum on the right of eq 35 establishes μ(1) i ; the second and third sums comprise μ(2) i . Substitution of these quantities into eq 34 and application of the identities from eqs 36 and 37 lead to the result
(32)
c∞ G−G = T 2
Note that this subtraction eliminates the reference states for chemical potential, μ0i , from the integrand. Expressing the freeenergy fluctuation in this way has the advantage that μi − μ∞ i is a small quantity, and consequently only two terms in the expansions of cT and yi are needed. Equation 32 is rigorous, but an integral over the equilibrium volume V∞ suffices to account for all terms of order (2). Let the error in the Gibbs energy due to neglect of the difference between V and V∞ be defined as a quantity Δ,
∞
⎛ ∂(μ − μ ) ⎞∞ ∑ ⎜⎜ i n ⎟⎟ yi(1)yj(1) dV + O(3) ∂yj V∞ ⎠ i,j≠n ⎝
∫
(38)
Insertion of the Fourier series from eq 1 and integration give a relationship in terms of the Fourier amplitudes, n G − G∞ = T 4
⎛ ∂(μ − μ ) ⎞∞ ∑ ∑ ∑ ⎜⎜ i n ⎟⎟ αimαjm + O(3) ∂yj ⎠ m=1 i≠n j≠n ⎝ ∞
(39)
Δ=
∫V −V
c ∑ yi (μi − μi∞) dV ∞ T i
where nT = ∑i ni is the total number of moles in the system. This is the desired expression of the free-energy fluctuation as a quadratic form of the independent fluctuation amplitudes, neglecting all quantities of order (3) and higher. Definition of K and Relationship to Q and Y. Let Kij represent the quadratic coefficients in the expansion of G with respect to the independent composition variables, i.e.,
(33)
∞ Here μi − μ∞ is of order (2); i is of order (1), and V − V therefore, Δ should be of order (3). Expansion of the right side of eq 32 up to terms of order (3) thus gives
D
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research G − G∞ = kT ∞
∞
This result is independent of the wavelength index m. The contributions to G − G∞ are thus nominally the same for each value of m. To avoid obtaining an infinite result, it is necessary to reason physically that the Fourier series in eqs 1 and 2 must be terminated when L∞/m reaches atomic dimensions. Insertion of eqs 40 and 47 into eq 46 yields an integral that can be solved explicitly by diagonalization of K. In matrix form, the result is
∑ ∑ ∑ K ijαimαjm + O(3) (40)
m=1 i≠n j≠n
where the divisions by temperature and Boltzmann’s constant k make Kij dimensionless. In agreement with eq 39, let ⎛ ⎞∞ n T ⎜ ∂(μi − μn ) ⎟ K ij = ⎟ ∂yj 4kT ∞ ⎜⎝ ⎠T , p , y
C0 =
k≠j,n
⎛ ⎞∞ n T ⎜ ∂ 2(G /n T) ⎟ = 4kT ∞ ⎜⎝ ∂yi ∂yj ⎟⎠ T ,p,y
(41)
R−1QK−1 = K−TQTR−T
In paper 1, the higher order of Δ (cf. eq 33) and the identity from eq 37 were not accounted for properly, leading to the inclusion of some second-order derivatives of total concentration in the quadratic expansion of G − G∞. These issues are resolved by using eq 41 to define the entries of matrix K rather than eq 47 of paper 1. K is symmetric, K = KT (or Kij = Kji), by Maxwell relations. In terms of activity coefficients on a mole-fraction basis, its entries become ∞ ⎡ ∂ ln(γi /γn) ⎤⎥ n TNA ⎢ δij 1 K ij = + + ⎥⎦ 4 ⎢⎣ yi yn ∂yj
RTQ−TKT = KQ−1R = (KQ−1R)T
⎛ ∂ ln γ ⎞∞ ⎛ ∂ ln γ ⎞∞ n⎟ k⎟ ∞⎜ ⎟ = − ∑ yk ⎜ ∂y ⎟ = 1 − y ∂ ⎝ ⎝ k≠n j ⎠ j ⎠
YR = (YR)T
∑ Q kj k≠n
∑ YikQ kj
NT N YQ = KT = T (YQ)T 4 4
(44)
(45)
Thus, K is extensive, proportional to the total number of molecules in the system. Determination of C0 for the Gibbs Ensemble. The initial correlations are computed according to Cij0 = A m
∞
∞
∫−∞ ∫−∞ αimαjme−G/kT
∞
dαim dαjm
(46)
where Am is a normalization factor. The coefficients Kij are extensive, as mentioned earlier. For large enough systems, terms of third and higher order contribute negligibly to the integral in eq 46, and ⎛ G∞ ⎞ det(K) A m = exp⎜ ∞ ⎟ ⎝ kT ⎠ π n − 1
(51)
DISCUSSION Harmony with the conventional statement of the reciprocal relation is achieved because the relation is reduced to eq 51, which becomes Mij = Mji or +ij = +ji for Onsager−Stefan− Maxwell diffusion. For many years there were three proofs of the reciprocal relation: (1) regression hypothesis,2,3,6 valid but based on an unusual driving force, one equivalent to the derivative of entropy with respect to fluctuation amplitudes; (2) minimum dissipation principle, stated by Onsager5 in 1945 based on the 1931 papers but refuted by Coleman and Truesdell;7 (3) pedagogic proof, based on Newton’s third law of motion but stated to be unreliable by Truesdell.8 For Onsager−Stefan−Maxwell diffusion, we now have Monroe and Newman’s detailed implementation of the first proof (as corrected here). The method could be applied to any set of fluxes and driving forces. Also, there is Wheeler and Newman’s extension of the fluctuation−dissipation theorem of Johnson9 and Nyquist.10 It is difficult to assess the validity of one theory versus another when various approaches lead to the same result. The method of Monroe and Newman, although it included an error in the initial correlations that is corrected in the present paper, is unassailable as an implementation of Onsager’s original procedure. It follows directly from the regression hypothesis, including the principle of microscopic reversibility and the determination of the magnitude of initial correlations on the basis of the Gibbs function. The only other known direct application of this procedure appears to be that of the Onsager papers involving heat conduction in an anisotropic crystal2 and the more general multicomponent derivation involving a driving force of (essentially) the derivative of
where the definition of Yij from eq 14 has been used and NT = nTNA represents the total number of molecules. In isothermal, isobaric, multicomponent transport systems, it may be concluded that K=
RY −1 = (RY −1)T
■
allows the elimination of γn, with the result
k≠n
or
(Recall that Y is symmetric.) Because both Q and K have disappeared, the same reciprocal relation applies independently of the thermodynamic ideality of the solution. The conventional reciprocal relation +ij = +ji results from eq 51 because RY−1 ∝ M′.
(43)
⎞ N Q ∑ kj⎟⎟ = T 4 ⎠ k≠n
(50)
Substitution of eq 45 into eq 50 shows that the reciprocal relation reduces to
(42)
yn∞⎜⎜
(49)
By matrix manipulations similar to those used earlier (after eq 25), one obtains
because R/k equals Avogadro’s number N A . Another application of the Gibbs−Duhem equation, in the form
NT ⎛⎜ Q ij 1 + ∞ 4 ⎜⎝ yi∞ yn
(48)
Resolution of the Discrepancy. Insertion of eq 48 into the reciprocal relation in eq 22 yields
k≠i ,j,n
K ij =
1 −1 K 2
(47) E
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research entropy with respect to αim.3 In other words, the method was not shown to apply to Onsager−Stefan−Maxwell diffusion. An interesting question here is whether Onsager knew a more convincing argument but neglected to add any additional condition that would have satisfied the refutation of Coleman and Truesdell. On the other hand, Onsager may simply have jumped to the conclusion that because his forces and fluxes, which lead to a symmetric transport matrix, add up to the entropy dissipation when multiplied pairwise, then any set of forces and fluxes that thus add to the entropy dissipation would lead to a symmetric transport matrix. In another paper,11 Monroe and Newman presumed the validity of eqs 45 and 48 (which are substantiated for the first time here) and demonstrated what conditions are required of “proper” flux/force pairs describing multicomponent diffusion. Any set of driving forces related to the set of independent −∇yi by an invertible linear transformation, which preserves the dissipation and does not combine any fluxes with the forces, yields a symmetric transport matrix. Such restrictions would likely have appeared natural to Onsager. It should be mentioned that there is a simple pedagogic argument (reproduced in ref 12) based on Newton’s third law of motion that action equals reaction. Applied to multicomponent diffusion, it suggests that symmetric retarding forces between pairs of species in transit would directly relate to pairwise diffusivities. Even though it leads to the correct result, namely, +ij = +ji , this is also regarded as an invalid argument.8 Wheeler and Newman4 did not set out to prove a symmetric reciprocal relation (i.e., it was taken for granted) but instead to derive a Green−Kubo-type relation for generating Onsager− Stefan−Maxwell diffusion coefficients from molecular-dynamics simulations. Linear-response theory suggests that an observable can be expressed as an integral over its past history.13 In the derivation of Wheeler and Newman,4 the initial state is at equilibrium, constant external forces are applied beginning at t = 0, the external forces create a nonequilibrium ensemble, and a theoretical thermostat removes the dissipated energy, so that a steady state can be reached. They express the linear response of the velocity difference Ji = vi − vn as ⟨Ji (t )⟩1 − ⟨Ji (0)⟩ =
V kT
∫0
t
reversibility, so that Lij is also symmetric, and this leads to the conclusion that +ij = +ji . (The matrix L is the negative inverse of the truncated, invertible Stefan−Maxwell matrix M′; if L is symmetric, then so is M′.) Thus, the derivation of eq 54 became a proof of the symmetry of the reciprocal relation, one that relies on linear-response theory13 and the principle of microscopic reversibility. The derivation of Wheeler and Newman can be regarded as complementary to that of Monroe and Newman. They both lead to the same result but by different routes. It is interesting to note that the initial value of the correlation does not appear in the proof of Wheeler and Newman. It is a clear participant in the symmetry proof of Onsager. Because the system is linear, initial values are necessary. For example, the solution of eq 9 requires initial values; for the equilibrium ensemble, C0 was used, but for a nonequilbrium system, the initial value could be asymmetric, in which case the values evolving over time would generally continue to be asymmetric. The value for the initial correlation in the equilibrium system is crucial for getting the proper reciprocal relation by the method of Monroe and Newman. The two methods also differ in the manner of imposing driving forces. The derivation of eq 54 involves no concentration gradients and relies on a uniform external body force for each species to create the fluxes. A companion nonequilibrium method14 likewise produces diffusion coefficients using a thermodynamic driving force. In contrast, the regression-hypothesis method uses a gradient in concentration or mole fraction as the driving force and would produce diffusion coefficients like those measured experimentally. It is analogous to the restricted-diffusion method, a technique whereby one-dimensional diffusion is generated in an enclosed volume of finite length containing a chemical mixture of two or more species, with measurement of the longest time constant.15,16 Thermodynamic matrix Q can be used to convert diffusion coefficients from one driving-force basis to the other; Q can be measured or directly calculated in a Monte Carlo or equilibrium molecular-dynamics simulation.17 Computing independent Q and transport matrices from the differing methods allows a consistency check to be made. Real molecular systems exhibit a Knudsen or inertial or prephenomenological regime for short observation times, in which an instantaneous transport relation such as the Onsager− Stefan−Maxwell law is not followed.4 During the course of collisions between molecules, the system transitions to linearrelaxation behavior. (Typically the duration of this transition is on the order of picoseconds for molecular liquids at 300 K.) When using molecular simulations to study linear transport coefficients, one must therefore ensure that correlation functions (e.g., eq 10, 17, or 54) are computed to times sufficiently long that prephenomenological behavior is not dominating the observed response. Lastly, one might ask whether an objection similar to that of Coleman and Truesdell7 could be raised with respect to transformed driving forces. That is, first write the transport relation in the form
⟨Ji (t0) ∑ Jj (0) Xj⟩1 dt0 j≠n
(52)
where Xj is a constant external force associated with −cTyj∇μj. The subscript 1 on angle brackets indicates that the ensemble average is taken under the perturbation of the external force. Because the system is initially equilibrated, ⟨Ji(0)⟩ = 0, although this does not mean that Jj(0) is zero in the correlation in the integral. Because Xj does not change with the position and periodic boundaries are employed, there is no overt variation of the concentration with the position. If a macroscopic transport relation
⟨Ji (∞)⟩1 =
∑ LijXj (53)
j≠n
is written, then Wheeler and Newman show that Lij =
V kT
∫0
∞
⟨Ji (t0) Jj (0)⟩ dt0
Xj = −∑ Mji′Ji
(54)
i≠n
which, in the limit of vanishing external perturbations Xj, can be computed in an equilibrium ensemble in the manner of a Green−Kubo formula. They noted that the correlation in the integrand is symmetric by the principle of microscopic
(55)
which is equivalent to eq 53, given the inverse relationship between M′ and L. Next take Mij = Mji. Construct a new set of forces Xj′ = Xj − ∑iWjiJi, where W is a nonzero skew-symmetric F
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research matrix, that is, Wij = −Wji or W = −WT. The properties now satisfied are as follows: (1) The entropy-production rate is given by either set: Tg =
boundary conditions can also be treated. Almost any quantity can be expanded in small quantities to yield a linear system. Monroe and Newman do this in detail to substantiate the reciprocal relations for Onsager−Stefan−Maxwell diffusion. However, more is contained in the regression hypothesis than just the reciprocal relation. Equation 9 can be used to follow the regression of fluctuations in an equilibrium ensemble or a nonequilibrium ensemble with no external forces, and eqs 11, 19, 23, and 24 show how one could calculate the transport coefficients in a method complementary to the Green−Kubotype integral in eq 54. A forthcoming paper by Nichols and Wheeler will illustrate such an equilibrium calculation for a binary mixture, showing that both complementary methods yield consistent results from the same molecular-dynamics simulation.
∑ XiJi = ∑ (Xi′ + ∑ WikJk )Ji i≠n
i≠n
=
k≠n
∑ Xi′Ji (56)
i≠n
(2) Forces and fluxes are still linearly related through Xj′ = −∑ (Mji′ + Wji)Ji = −∑ Mji″Ji i≠n
i≠n
(57)
(3) The new transport matrix is not symmetric, since Mji″ = Mji′ + Wji ≠ Mij′ + Wij
■
(58)
Coleman and Truesdell thus showed that forces and fluxes chosen according to inspection of the dissipation sum can be either symmetric or asymmetric and that there is no physical content in such an exercise. Would the method of Wheeler and Newman be subject to this problem? It has been asserted that the methods of modern statistical mechanics have been justified by the fluctuation− dissipation theorem,13 although not necessarily in the context of multicomponent diffusion. It can also be argued that the alternative forces defined above, in the manner of Coleman and Truesdell, are contrived because they mix forces and fluxes and that the forces used in the derivation of Wheeler and Newman are clearly defined and related to the driving force of the gradient of the chemical potential.
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS C.W.M. was supported by the U.S. National Science Foundation (Grant CBET 1253544). D.R.W. was partially supported by the Advanced Battery Materials Research (BMR) program of the U.S. Department of Energy.
■
LIST OF SYMBOLS
Matrices
■
C = (n − 1) × (n − 1) composition−fluctuation correlation matrix defined in eq 10 K = (n − 1) × (n − 1) thermodynamic matrix derived from quadratic expansion of G around equilibrium, as in eq 40; K = 1/4NTYQ; for an equilibrium ensemble of sufficient size, C0 = 1/2K−1 L = inverse of the transport matrix −M′; see eq 53 M = n × n singular Stefan−Maxwell matrix; see eq 7 M′ = (n − 1) × (n − 1) nonsingular Stefan−Maxwell matrix formed by deleting the nth row and the nth column from M; −1 at the equilibrium composition M′ = RTc∞ T RY M″ = (n − 1) × (n − 1) constructed asymmetric transport matrix; see eq 57 P = (n − 1) × (n − 1) transport matrix P = −R−1Q, similar to an inverse matrix like L but based on the −∇yi driving force Q = (n − 1) × (n − 1) thermodynamic matrix that converts independent diffusion driving forces from −cTyi∇μi to −∇yi; see eqs 6 and 45 R = (n − 1) × (n − 1) modified resistance matrix M′ after Fourier coefficients for component n (yn and vn) are eliminated; see eq 8; at the equilibrium composition R = M′Y/RTc∞ T V = (n − 1) × (n − 1) velocity-correlation matrix defined in eq 17 W = (n − 1) × (n − 1) skew-symmetric transport matrix; see eq 58 Y = (n − 1) × (n − 1) nonsingular compositional matrix ∞ ∞ defined in eq 14; note that (Y−1)ij = δijy∞ i − yi yj
CONCLUDING REMARKS In treating a nonequilibrium linear system, departures from equilibrium are proportional to an external driving force or arise from a composition gradient. The Onsager−Stefan−Maxwell equations assume an instantaneous transport relation between forces and fluxes and therefore neglect the prephenomenological behavior of real molecular systems. Molecular simulations can be used to observe and ensure the onset of linear relaxation behavior. One must verify that sufficient sampling is performed outside the prephenomenological time domain when linear transport coefficients are being computed. The external-driving-force method leads to separate but equivalent algorithms for obtaining transport coefficients in equilibrium and nonequilibrium simulations, under the assumption of linear transport.3,14 The regression-hypothesis method (as corrected here) yields the proper reciprocal relation among transport coefficients1 and likewise can be used to calculate transport coefficients from a molecular-dynamics simulation of an equilibrium ensemble or a nonequilibrium ensemble without external forces. This method yields P rather than M or R. Because it is based on the regression hypothesis, this method can predict the decay of fluctuations in an equilibrium or a nonequilibrium ensemble. Onsager’s regression hypothesis is not contained in the fluctuation−dissipation theorem and provides an alternative approach to transport analysis. After one has the basic governing laws, macroscopic theory provides a guide to calculating (on a linear system) for an arbitrary disturbance. Take its Fourier transform and superpose the results. Initial concentration profiles or situations imposed through the
Roman Letters
cT = total molar concentration, mol·m−3 Cij = entry in the correlation matrix C, dimensionless
G
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
+ij = Stefan−Maxwell coefficient of i through j, m2·s−1 g = entropy-generation rate, J·K−1·m−3·s−1 G = Gibbs function, J Ji = excess velocity of component i relative to n, m·s−1 k = R/NA, Boltzmann’s constant, 1.3806 J·K−1 Kij = entry in the thermodynamic matrix K, dimensionless; see eqs 41 and 44 L = fluctuating length of a one-dimensional slab of material, m Lij = entry in the inverted truncated Stefan−Maxwell matrix L, m5·s−1·J−1 Mij = entry in the Stefan−Maxwell matrix M, J·s·m−5 M′ij = entry in the nonsingular truncated Stefan−Maxwell matrix M′, J·s·m−5 M″ij = entry in the matrix M′ + W, J·s·m−5 n = number of components in the system, dimensionless ni = number of moles of component i, mol nT = total number of moles, nT = ∑ini, mol NA = R/k, Avogadro’s number, 6.0221 × 1023 mol−1 NT = total number of molecules, dimensionless O = order symbol p = ambient pressure (assumed to be constant), J·m−3 Pij = entry in the transport matrix P, m2·s−1 Qij = entry in the thermodynamic matrix Q, dimensionless R = kNA, universal gas constant, 8.3145 J·mol−1·K−1 Rij = entry in the transport matrix R, s·m−2 t = time, s T = absolute temperature (assumed to be constant), K vi = velocity of component i, m·s−1 V = system volume, m3 Vij = entry in the velocity-correlation matrix V, m2·s−2 Wij = entry in the antisymmetric transport matrix W, J·s·m−5 Xi = force driving diffusion of component i, J·m−4 Xi′ = transformed driving force, J·m−4; see eq 56 yi = mole fraction of component i, dimensionless Yij = entry in the compositional matrix Y, dimensionless z = position in a one-dimensional slab of material, m
(4) Wheeler, D. R.; Newman, J. Molecular Dynamics Simulations of Multicomponent Diffusion. I. Equilibrium Method. J. Phys. Chem. B 2004, 108, 18353−18361. (5) Onsager, L. Theories and Problems of Liquid Diffusion. Ann. N.Y. Acad. Sci. 1945, 46, 241−265. (6) Casimir, H. B. G. On Onsager’s Principle of Microscopic Reversibility. Rev. Mod. Phys. 1945, 17, 343−350. (7) Coleman, B. D.; Truesdell, C. On the Reciprocal Relations of Onsager. J. Chem. Phys. 1960, 33, 28−31. (8) Truesdell, C. Mechanical Basis of Diffusion. J. Chem. Phys. 1962, 37, 2336−2344. (9) Johnson, J. B. Thermal Agitation of Electricity in Conductors. Phys. Rev. 1928, 32, 97−109. (10) Nyquist, H. Thermal Agitation of Electric Charge in Conductors. Phys. Rev. 1928, 32, 110−113. (11) Monroe, C. W.; Newman, J. Onsager’s shortcut to proper forces and fluxes. Chem. Eng. Sci. 2009, 64, 4804−4809. (12) Newman, J.; Thomas-Alyea, K. E. Electrochemical Systems, 3rd ed.; John Wiley and Sons, Inc.: Hoboken, NJ, 2004; pp 297−298. (13) Evans, D. J.; Morriss, G. P. Statistical Mechanics of Nonequilibrium Liquids; Academic Press: New York, 1990. (14) Wheeler, D. R.; Newman, J. Molecular Dynamics Simulations of Multicomponent Diffusion. II. Nonequilibrium Method. J. Phys. Chem. B 2004, 108, 18362−18367. (15) Harned, H. S.; French, D. M. A Conductance Method for the Determination of the Diffusion Coefficients of Electrolytes. Ann. N.Y. Acad. Sci. 1945, 46, 267−281. (16) Newman, J.; Chapman, T. W. Restricted Diffusion in Binary Solutions. AIChE J. 1973, 19, 343−348. (17) 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.
Greek Letters
αim = Fourier amplitude of yi(1) at wavenumber m, dimensionless βim = Fourier amplitude of v(1) at wavenumber m, m·s−1 i δij = Kronecker delta, 1 if i = j and 0 otherwise Δ = error in the Gibbs function due to neglect of V − V∞, J γi = activity coefficient of component i on a mole-fraction basis, dimensionless τ = scaled time variable, τ = (mπ/L∞)2t, s·m−2 μi = chemical potential of component i, J·mol−1 Superscripts and Subscripts
= initial value (as in C0) or reference state (as in μ0i ) = value in a homogeneous equilibrium state (1) = value obtained with a linearized governing system (2) = value obtained with a second-order governing system (3) = value obtained with a third-order governing system = initial value (as in τ0) 0 0
∞
■
REFERENCES
(1) Monroe, C. W.; Newman, J. Onsager Reciprocal Relations for Stefan−Maxwell Diffusion. Ind. Eng. Chem. Res. 2006, 45, 5361−5367. (2) Onsager, L. Reciprocal Relations in Irreversible Processes. I. Phys. Rev. 1931, 37, 405−426. (3) Onsager, L. Reciprocal Relations in Irreversible Processes. II. Phys. Rev. 1931, 38, 2265−2279. H
DOI: 10.1021/ie503875c Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
