Ind. Eng. Chem. Res. 2003, 42, 7025-7036
7025
Characterizing Isoparaffin Transport Properties with Stokes-Einstein Relationships Peter A. Gordon* Corporate Strategic Research, ExxonMobil Research & Engineering, 1545 Route 22 East, Annandale, New Jersey 08801
We present extensive molecular dynamics simulation results for self-diffusion and viscosity for isoparaffins between C6 and C16 and examine the ability of several forms of Stokes-Einstein relationships to link the transport properties in a coherent framework. Based on a hydrodynamic model for drag past prolate ellipsoids in the slip limit, much of the variation between diffusion and viscosity can be explained as a function of molecular structure. For pure species, we demonstrate how one can exploit the linearity of the Stokes-Einstein relationship to make quantitative, extrapolative predictions of viscosity from computer simulations based on the selfdiffusion coefficient and systematic characterization of the molecule’s gyration ellipsoid. Introduction In a previous study1 we demonstrated that viscosity computed from equilibrium molecular dynamics (EMD) simulations using the TRAPPE-UA intermolecular potential for isoparaffins2 appeared capable of quantitatively describing the impact of chain branching on viscosity for a variety of low molecular weight isomers. In addition to several examples illustrating this point, we have computed transport property data for nearly 80 different isomers between carbon number 6 and 16, at a variety of temperatures and pressures. The motivation behind generating such a large set of transport data was to develop a deeper understanding of molecular structure-rheology relationships over a wide range of molecular chain architectures. The utility of a molecular model that is capable of accurately describing the relationship between branching and viscosity should be readily apparent. With such a model, one can explore both fundamental questions regarding the structure-transport properties of complex fluids and questions of direct relevance to lubricant basestock composition. For instance, for a given molecular weight, what molecular structure(s) represent upper and lower bounds on viscosity? Do the structures that represent these bounds change as a function of temperature and/or pressure? Can properties for which little data exists, such as pressure-viscosity coefficients and low-temperature viscosity, be better understood in terms of molecular structure? The answers to these types of questions should allow one to identify ‘ideal’ molecular structures for a given performance attribute, if such ideals exist, or tailor the composition of a lubricant basestock for improved properties at the extremes of low temperature and high pressure required in many lubrication applications. A complicating feature of this approach is the computational burden associated with computing viscosity. With equilibrium molecular dynamics, the time required to obtain reasonably converged estimates of viscosity rises quickly with increasing molecular weight, making EMD viscosity calculations on long chain molecules very time-consuming with today’s computing power. This * To whom correspondence should be addressed. Tel.: (908) 730-2546. Fax: (908) 730-3031. E-mail: Peter.A.Gordon@ exxonmobil.com.
Figure 1. Self-diffusivity plotted against the product of viscosity and radius of gyration for selected isoparaffins between C6 and C16 at T ) 293 K. For molecules of similar shape, defined by φeff, the characteristic Stokes-Einstein slope is nearly equivalent. Similar trends are observed for data at 373 K as well.
poor scaling stems from the slowly converging stress autocorrelation function, which requires longer and longer simulations to obtain reasonable statistics. By contrast, properties of individual molecules, such as translational and rotational diffusion, converge much more rapidly in a simulation. It would then be quite desirable to be able to understand the relationship between viscosity and other, more simply computed properties. If such an understanding could be developed, viscosity could be inferred through calculation of these simple quantities. In the course of analyzing the isoparaffin transport data from previous work, a strong correlation between self-diffusion, viscosity, and molecular shape was observed that appeared to be valid for a wide variety of chain architectures and is depicted in Figure 1. When self-diffusivity is plotted against the product of viscosity and the radius of gyration, similarly shaped molecules, as characterized by the molecule’s gyration ellipsoid and effective aspect ratio, φeff, are found to obey the same characteristic slope, regardless of chain branching.
10.1021/ie030512x CCC: $25.00 © 2003 American Chemical Society Published on Web 11/18/2003
7026 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
Before exploring in detail the full meaning of this plot, we present these data, describing its meaning qualitatively, to underscore our motivation. We seek to develop a unified model that can both capture molecular structure-viscosity relationships and make extrapolative predictions of viscosity from more easily computed molecular properties such as self-diffusion and molecular shape. Figure 1 is highly suggestive of the idea that such a unifying model can be developed. In this paper, we employ Stokes-Einstein relationships to link translational diffusion to viscosity. We explore several different formulations of such relationships, examining the influence of boundary conditions employed in the hydrodynamic theory and the effect of incorporating molecular shape in the model. It is found that a hydrodynamic description for creeping flow past prolate ellipsoids undergoing isotropic translation solved in the slip limit provides a very good description of the relationship between molecular shape, diffusion, and viscosity. This is demonstrated through the rationalization of transport data from extensive simulations of several C8 and C12 isoparaffins. The power of the approach is demonstrated by making quantitative extrapolations of viscosity to low temperature and high pressure using solely self-diffusion and geometric considerations. This paper is organized as follows. We begin with a brief introduction to the development of Stokes-Einstein equations, emphasizing differences in the resulting equations from different hydrodynamic boundary conditions and incorporating asphericity into the model. We then outline several computational procedures by which we can characterize molecular geometry and shape. We compare the ability of these Stokes-Einstein formulations to rationalize extensive transport data for n-octane and then examine the models’ ability to rationalize transport data for a wide variety of isoparaffin molecular structures at the same statepoint. Finally, we focus on sets of simulations for single species, conducted over a wide range of temperature and pressure, to test the range of validity of the model and to demonstrate the ability to make extrapolative predictions of viscosity.
The frictional resistance in eq 1 can be related to the diffusivity of a solute particle moving in a solvent through Einstein’s theory of Brownian motion. In the case of a bulk liquid, the solute molecule can be thought of as a tagged particle that is otherwise identical to the solvent molecules. The motion of this particle is denoted as the self-diffusivity, D, and is identical to the macroscopic (Fickian) transport diffusivity in the limit of zero concentration. Combining eq 1 with Einstein’s equation, we arrive at an equation for the self-diffusivity of the form
Background
with CBC ) 6π in the stick boundary condition limit. The geometric correction factor K| is solely a function of the aspect ratio of the prolate, φ ) a/b, and in the stick limit is given by the expression
Stokes-Einstein Relationships. A central problem in low Reynolds number hydrodynamic analysis is the estimation of the drag on an object subjected to a creeping (low velocity) flow field in some particular orientation.3,4 Stokes’ law expressing the drag force F past a sphere moving with velocity U in a fluid of viscosity η is
f)
F ) CBCRη U
(1)
The friction coefficient, f (also called the frictional resistance), is often used in place of drag force, where it is understood that the direction of the drag force opposes that of the direction of motion. The constant CBC depends on the boundary conditions applied at the surface of the sphere. In the no-slip or ‘stick’ limit, the fluid velocity on the surface is zero. In the slip limit, the boundary conditions applied at the surface require that (i) the fluid velocity component normal to the surface be zero, and (ii) no tangential stresses exist. In these two limits, CBC is given as 6π and 4π, respectively.
D)
kT kT ) f CBCRη
(2)
which Einstein deduced from the conditions of dynamic equilibrium experienced by particles suspended in a liquid.5 This equilibrium is a balance between forces caused by osmotic pressure differences from local concentration gradients with the frictional drag induced by the particle motion interacting with the solvent. For spherical particles, the friction coefficient of eq 1 has a simple form. The friction coefficient can also be analyzed for translation of aspherical objects. In this case, however, the resistance to translation will be a function of the object’s shape and its orientation with respect to the velocity field. For axisymmetric objects, for instance, one would expect differences in the drag force resisting translation parallel or perpendicular to the symmetry axis. In the stick limit, the friction coefficient has been derived for a variety of axisymmetric particle shapes, including prolate and oblate ellipsoids, disks, and cylinders.3,6 Development of corresponding expressions in the slip limit is generally more complicated, but, where available for comparison, these expressions exhibit qualitatively different behavior. For instance, for a prolate ellipsoid with major semiaxis a and equivalent minor semiaxes b, translation parallel to the long axis results in a drag coefficient given by
f)
K|(φ) )
[
F ) CBCbηK//(φ) U
(
(3)
) (x )]
2φ2 - 1 4 φ + - 2 coth-1 3 φ -1 (φ2 - 1)3/2
-1
φ
φ2 - 1
(4)
For translation of a prolate ellipsoid perpendicular to its long semiaxis, K| in eq 3 is replaced by K⊥
K⊥(φ) )
[
(
)
]
8 φ 2φ2 - 3 ln(φ + xφ2 - 1) + 2 (3/2) 3 φ2 - 1 (φ - 1)
-1
(5)
For a molecule in a bulk liquid undergoing Brownian motion, the motion of an individual particle is governed by the random (isotropic) impact of neighboring solvent molecules. The net resistance to this motion, KAV, can be found by averaging the resistance over all possible flow directions.3 This averaging can be decomposed along any set of orthogonal body-fixed coordinates; in
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 7027
Figure 2. Average friction coefficient for prolate ellipsoids undergoing isotropic translation: effect of boundary conditions. ([) Prolate ellipsoids in the stick limit, (9) prolate ellipsoids in the slip limit, (]) sphere with equivalent surface area, stick limit, (0) sphere with equivalent surface area, slip limit.
the case of prolate ellipsoids, it is natural to express this is as
[
] [
1 1 1 1 1 1 2 1 ) + + ) + KAV 3 K⊥ K⊥ K| 3 K⊥ K|
]
limit appears to be a more accurate description of the molecular motion of these objects. Note that this applies to systems where the solvent and solute are of roughly the same molecular dimensions. Linking Hydrodynamics and Molecular Geometry. Previous studies by several authors9-12 have employed the unit vectors of the principal moments of inertia as a means to characterize the orientational motion of molecules. This interpretation is particularly useful for branched molecules, where identification of a single vector quantity to track molecular conformation, such as the end-to-end vector, is not straightforward. During the course of analyzing isoparaffin simulations between C4 and C26, however, it was also noted that the average shape of the gyration ellipsoids for these structures closely approximates that of a prolate ellipsoid. Thus, we can attempt to link the gyration ellipsoid representation with the hydrodynamic models given by eqs 3-8, in both the stick and slip limits. The radius of gyration lends itself to a natural geometric interpretation of the molecular motion. For a molecule composed of a set of bonded interaction sites with masses mI, the squared radius of gyration is related to the on-diagonal components of the moments of inertia
(6a) R2gyr )
with the average frictional resistance given as
fAV ) CBCbηKAV(φ)
2
K|(φ) )
1 2 + exp(-0.3φ) 2 3
(7) (8)
The average (reduced by ηb) resistance a prolate ellipsoid experiences in the stick and slip limits is plotted as a function of φ in Figure 2. In the limit of spherical particles (φf1), the stick and slip limits approach their respective spherical counterparts of 6π and 4π, respectively. The qualitative differences in the friction coefficient obtained in the stick and slip limit are apparent. With increasing aspect ratio, the frictional resistance in the stick limit increases monotonically. In the slip limit, however, we see the frictional resistance levels off and actually modestly decreases. The frictional resistance encountered by a translating sphere of equivalent surface area to an ellipsoid of aspect ratio φ (and unit minor semiaxis, b) is also shown in the stick and slip limits. The slip boundary condition is still a monotonically increasing function of the surface area in the spherical case. The apparent decoupling of the frictional resistance of prolate ellipsoids from the longaxis of the object suggests a short time preferential diffusion direction along the path of least resistances the long axis of the prolate. Computer simulations of hard prolates in the bulk7 and as solutes in spherical solvents8 have demonstrated convincingly that the slip
∑
(9) mi
sites
(6b)
In the limiting case of no-stick or ‘slip’ boundary conditions, no simple analytic expression exists analogous to eqs 4 and 5. However, Tang and Evans7 have developed an approximate solution for hard prolate ellipsoids of aspect ratio 1 e φ e 10. Over this range, eq 6b applies with CBC ) 4π, and the geometric factors are approximately given as
K⊥(φ) ) 1 + 0.54(φ - 1) - 0.011(φ - 1)2
Tr(I)
and the inertia tensor element IRβ (R,β vary over x, y, and z directions) is defined as
IRβ )
mi[δRβ(|r bi - b r com|2) - (ri,R - rcom,R)(ri,β ∑ sites rcom,β)] (10)
The sum in eq 10 is over interaction sites in the molecule, and b ri is the position vector of the ith interaction site. The arbitrary reference point is taken to be the molecule center of mass, b rcom. The inertia tensor may be represented geometrically by its inertia ellipsoid equation
∑R ∑β IRβ(rbcom)xRxβ ) γ
(11)
where xR (R ) 1,..3) refers to the x, y, and z directions, respectively. The quantity γ is an arbitrary positive constant with dimensions [Mass×Length4]. Again, this equation is referenced to the center of mass of the molecule. The moment of inertia tensor may be diagonalized through solution of the associated eigenvalue problem. The three eigenvalues found, λk (k ) 1,3), are the ondiagonal elements of the moment of inertia tensor. The eigenvalues are directed along b ek, the corresponding normalized eigenvectors. Equations 9 and 11 can be expressed in terms of the eigenvalues and eigenvectors as 3
R2gyr )
∑ λk k)1 2
∑ mi
sites
(12)
7028 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
and
λ1|e b1|2 + λ2|e b2|2 + λ3|e b3|2 ) γ
(13)
Note that eq 13 has the form of an ellipsoid with semiaxes a, b, and c, and volume Vellipse given as
a ) xγ/λ1, b ) xγ/λ2, etc. Vellipse )
4 πγ3/2 (14) 3 λλλ x123
It is written in terms of the reference frame of the gyration ellipsoid; expansion of the eigenvectors, themselves expressed in a space-fixed coordinate system, yields the equation for the gyration ellipsoid in the space-fixed frame. To complete the characterization of the gyration ellipsoid, the quantity γ must be specified in a systematic way. There are several ways this can be accomplished. One method is to consider the gyration ellipsoid of a body of total mass M and uniform density, bounded by the surface
y2 z2 x2 + 2 + 2 )1 2 aunif bunif cunif
(15)
where the subscript of the semiaxes refers to the uniform mass density of the body. The moments of inertia for this body are related to the semiaxes as13
Ixx )
M 2 M M (b + c2), Iyy ) (a2 + c2), Izz ) (a2 + b2) 5 5 5 (16)
If we impose the requirement that this ellipsoid is ‘inertially equivalent’ to that of the discrete mass distribution of sites in a molecule, then the set of λk can be equated with the principal moments of inertia of eq 16. With this assumption, eq 16 can be rearranged for the semiaxes of the inertially equivalent ellipsoid (INRTEQV) of uniform mass distribution
5 (λ + λ3 - λ1) a2unif ) 2M 2
(17)
with similar formulas for the other semiaxes obtained by permutation of indices. Another means of characterizing γ is to equate the volume of the inertially equivalent ellipsoid above to the ellipsoid volume defined by eq 14. This volume-matched ellipsoid (VMTCH-INRTEQV) is given by eq 13 with 2/3
γ ) (aunifbunifcunif) [λ1λ2λ3]
1/3
elegant method to analytically determine the true volume of a collection of overlapping spheres, based on a geometric analysis of the planes formed by intersecting spheres. If we define this measure of molecular volume as VSPHPL, we can match this to the ellipsoid volume of eq 13 for the discrete mass distribution of the molecule and define the semiaxes of the gyration ellipsoid in an analogous manner to eq 19. In this case,
γSPHPL )
(4π3 V
SPHPL
)
2/3
[λ1λ2λ3]1/3
(20)
(18)
and semiaxes
adisc ) xγ/λ1, etc
Figure 3. Representations of molecular volume for TRAPPE-UA model of 3,3,5-trimethylheptane. Top: Molecular volume is given by the space filled by overlapping Lennard-Jones united atom spheres. Bottom: Volume approximated by radius of gyration ellipsoid. The atomic radii have been reduced for visual clarity. The lighter translucent ellipsoid is the inertially equivalent gyration ellipsoid corresponding to an object of uniform mass density (eqs 15 and 17). The ellipsoid shaded purple is the gyration ellipsoid with discrete mass distribution volume matched to the inertially equivalent ellipsoid (eqs 13 and 18).
(19)
A third method of estimating molecular volume can be employed by directly characterizing the space occupied by the molecule. A typical molecular model is composed of a set of bonded spheres, each with a characteristic radius. An estimate of the molecular volume can be obtained from computation of the total volume of the assemblage of spheres. Bonded spheres in the models employed in this work contain a high degree of overlap, and this overlap must be accounted for in order to obtain an accurate estimate of the true molecular volume. Dodd et al.14 have developed an
We employ the label for the molecular volume, ‘SPHPL’, to refer to the sphere-planes algorithm14 used in characterizing the volume of overlapping spheres and use the TRAPPE-UA model values of σ, the interaction site diameters, for this calculation. We denote this ellipsoid characterization as VMTCH-VSPHPL. An example of the different representations of molecular volume employed here is shown in Figure 3. The molecular volume of 3,3,5-trimethylheptane, with atomic radii given by the Lennard-Jones parameters of the TRAPPE-UA model employed in this work, is depicted in the top portion of the figure. Characterization of the same molecule through the two prolate ellipsoid representations using both the inertially equivalent ellipsoid of uniform density (eq 17) and the volume-matched ellipsoid (eqs 18 and 19) are shown in the bottom portion
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 7029 Table 1. Summary of Hydrodynamic Models and Molecular Geometry Analysis Conventions Used in Comparing Stokes-Einstein Characterizations of Transport Data in This Work title INRTEQV VMTCH-INRTEQV VMTCH-VSPHPL EQVSPH
description
shape
inertially equivalent ellipsoid of uniform density ellipsoid volume matched to inertially-equivalent ellipsoid ellipsoid volume matched to volume of overlapping spheres equivalent sphere of volume ) volume of overlapping spheres
prolate ellipsoid prolate ellipsoid prolate ellipsoid sphere
of the figure. We note that the volume of the gyration ellipsoid volume does not necessarily reflect the actual volume occupied by the molecular model but should provide a self-consistent means of applying the ellipsoid model to molecules of different shapes. In fact, StokesEinstein relationships are often used in the inverse form; assuming a particular hydrodynamic model is valid, such as that for spherical particles represented by eq 2, the effective molecular size of the particle can be deduced from transport data. If we wish to test the validity of the Stokes-Einstein equation, we need an independent method of assessing molecular volume and/ or shape. The above approaches offer methods of doing so, but these are by no means unique. While each of the above descriptions fully characterize the gyration ellipsoid, we must approximately map the ellipsoid to a prolate shape in order to make use of the hydrodynamic connections described by eqs 3-8. The shape of the prolate ellipsoid is characterized by the parameter φeff, the prolate aspect ratio, given as
φeff )
x
a1 a12 ) a2a3 bav
(21)
where a1 is the longest semiaxis of the gyration ellipsoid, and bav is the geometric mean of the minor semiaxes, a2 and a3. The aspect ratio characterizes the deviation from a spherical shape and varies from 1 (spherical) to ∞ (needle). We also note that the effective aspect ratio, φeff, is completely determined by the principal moments of inertia and is unaffected by the methods by which we scale the ellipsoid volume. Thus, given the prescription for the geometric characterization, we can obtain a time-averaged shape and volume of the gyration ellipsoids by averaging over each molecule in the system during the course of a simulation. The values of bav and φeff obtained from simulations can then be used to characterize the geometry and frictional resistance described in eq 6b and its variants. Finally, we consider a simpler representation of the molecules, which assumes a spherical particle shape with the volume given by the spheres-planes algorithm. This representation, which we denote EQVSPH, is simply characterized by the sphere radius
RSPH )
(
)
3VSPHPL 4π
1/3
(22)
This characterization is analogous to estimating the molecular volume via the fluidity analysis of Hildebrand.15,16 By extrapolating the linear, low-density region of a plot of fluidity () 1/viscosity) vs molar volume along an isotherm, the limiting molar volume, I0 corresponding to infinite viscosity is obtained. Dymond17 suggested equating this state to a hexagonally close-
hydrodynamic boundary conditions
molecular geometry characterization
stick (eqs 4-6) slip (eqs 6-8) stick (eqs 4-6) slip (eqs 6-8) stick (eqs 4-6) slip (eqs 6-8) stick (eq 2, CBC ) 6π) slip (eq 2, CBC ) 4π)
eqs 15 and 17 eq 13, γ from eq 18 eq 13, γ from eq 20 RSPH from eq 22
packed array of hard spheres of diameter σ. The sphere diameter can be obtained from the relationship
x2I0 1 ≈ 1.384 ) φhcp N σ3
(23)
A
where NA is Avogadro’s number. While not strictly appropriate for aspherical molecules, this method does offer a systematic way to extract molecular volume from experimental data. In addition, for small molecules with little conformational flexibility, the effective diameter is expected to be only weakly dependent on temperature. The fluidity analysis for characterizing molecular has been employed in a number of studies and has demonstrated that eqs 1 and 2 hold over a very wide range of temperatures and pressures for several simple molecular fluids, including benzene, tetramethylsilane,18,19 cyclohexane,20 methylcyclohexane,21 pyridine,22 and perfluorocyclobutane.23 A complicating feature of this method is the requirement for molar volume data over a range of pressures. To summarize, the mass distribution corresponding to the collection of atoms in a molecule can be visualized as a gyration ellipsoid, a spheroid with direct physical correspondence to shapes whose drag forces (undergoing Brownian motion) can be described by hydrodynamics. Different approximations for systematically characterizing molecular shape and volume can be imagined, and we test each model’s ability to rationalize transport property data of a variety of molecular species. The approaches employed in subsequent sections, along with acronyms describing their features, are summarized in Table 1. Simulation Methods To test various representations of the Stokes-Einstein formulas, EMD simulations were conducted to obtain transport properties (self-diffusivity and viscosity) for a number of isomers of octane and dodecane in the bulk over a wide range of temperatures and pressures. In addition, simulations of nearly 80 isoparaffin systems were carried out at 293 and 373 K and 0.1 MPa. Linear and branched paraffins were modeled using the TRAPPE-United Atom potential parametrized by Martin et al.2 In this representation, hydrogen atoms are lumped with their bonded carbon neighbors to form an effective interaction site. While this model generally underpredicts viscosity relative to experimental values, several studies have validated that it correctly captures many of the essential relationships between structural variation and viscosity.1,24,25 The general protocols employed in carrying out these calculations are the same as reported in earlier work.1,26
7030 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
Green-Kubo relationships27,28 were employed to compute the transport coefficients of interest. Self-diffusivity is related to the integral of the velocity autocorrelation function
D)
1 3Nmol
Nmol
∫0 ∑ 〈vi(t0)‚vi(t)〉dt i ∞
(24)
and viscosity is computed through integration of the stress tensor autocorrelation function
η)
V
∞ 〈 Pst (t )Pst (t)〉dt ∫ 0 ∑ Rβ 0 Rβ 10kT
(25)
Rβ
where the summation Rβ is over the six independent elements of the second-order symmetrized traceless pressure tensor. The molecular pressure formalism29,30 is employed in all calculations in this work. The eigenvectors associated with the principal moments of inertia (eq 13) are used to track molecular orientation and rotational motion and to characterize the size and shape of the gyration ellipsoid. The eigenvector associated with the smallest eigenvalue, denoted as b ei,min, corresponds to the longest axis of the molecule; an orientational autocorrelation function can be defined from this vector quantity as
oacf(t) )
1 Nmol
Figure 4. Viscosity (hollow symbols) and self-diffusivity (filled symbols) of bulk n-octane over a range of temperatures and pressures computed via equilibrium molecular dynamics. Data are plotted for 3 isotherms: (]) T ) 293 K, (0) T ) 353 K, (4) T ) 410 K. Self-diffusivity isotherms are plotted on the secondary axis. Uncertainties in the viscosity calculations are shown for the 293 K isotherm.
Nmol
∑i 〈ebi,min(t0)‚ebi,min(t)〉
(26)
The orientation autocorrelation function generally exhibits stretched exponential decay with an exponent slightly less than unity. We can extract a characteristic rotational relaxation time through the KohlrauschWilliams-Watts equation31
τrot(t) )
∫0∞oacf(t)dt ) ∫0∞exp(-(t/τ)β)dt
(27)
This characteristic time is a useful measure of the total simulation time required in order to estimate viscosity via eq 25. The quantities associated with characterizing the molecular shape and volume were obtained through thermodynamic averaging over each molecule’s configuration as sampled during the simulation at the temperature and pressure of interest. For the gyration ellipsoids, the semiaxes and the aspect ratio were computed for each molecule at fixed intervals during the course of a simulation run. For the effective sphere representation, the equivalent spherical radius is the quantity sampled. Results Single Species - n-Octane: Tests of Different Formulations of Stokes-Einstein Equations. In Figure 4 the viscosity and self-diffusivity of bulk noctane obtained via EMD are plotted as a function of pressure for 3 different isotherms. Representative uncertainties in the viscosity calculations are shown for the 293 K isotherm for reference. The expected exponential-like increase in viscosity with pressure is observed, particularly at the lowest temperature. Similarly, self-diffusion decreases with increasing pressure and decreasing temperature.
Figure 5. Stokes-Einstein plot of n-octane transport data using the model of inertially equivalent prolate ellipsoids of uniform density (INRTEQV). The symbols correspond to the hydrodynamic model applied with different stick (O) and slip (0) boundary conditions. The constant CBC is 6π for the stick model, and 4π for the slip model. The dashed line of slope ) 1 represents predicted scaling of each Stokes-Einstein model.
Figure 5 depicts the Stokes-Einstein relationship with the inertially equivalent (INRTEQV) prolate ellipsoid model applied to the n-octane transport data of Figure 4. Both the stick (CBC ) 6π) and slip (CBC ) 4π) boundary condition models are shown. In these representations, the transport data fall onto a single curve, and the linearity of the relationship appears to hold over the entire range of pressures and temperatures considered. The primary difference between the two models is the value of the average resistance, Kav. The aspect ratio of the n-octane molecule varies only modestly over the range of T and P considered; φeff varies from 4.0 at 233 K to 3.33 at 410 K. Furthermore, the average aspect ratio shows almost no variation with pressure along an isotherm. The corresponding values of Kav range from 2.3 to 2.8 in the stick model and 1.33-1.35 for the slip model. Thus, the lower slope obtained in the stick model in Figure 5 is a consequence of the larger estimated drag
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 7031
Figure 6. Stokes-Einstein plot of n-octane transport data using the model of prolate gyration ellipsoids with discrete mass distribution, volume matched to the inertially equivalent ellipsoid (VMTCH-INERTEQV). Symbols and lines are the same as Figure 5.
Figure 7. Stokes-Einstein plot of n-octane transport data using the model of prolate gyration ellipsoids with discrete mass distribution, volume matched to the volume determined from method of overlapping spheres (VMTCH-VSPHPL). Symbols and lines are the same as Figure 5.
force predicted for the translating molecule. We note that if the Stokes-Einstein model were quantitative, the data should fall on the dotted line in the figure (denoted slope ) 1). The slope of the data represented by the slip model is higher than expected and lower than expected for the stick model. These deviations appear to be nearly equal in magnitude. Figure 6 shows the application of the prolate ellipsoid model volume matched to the inertially-equivalent ellipsoid (VMTCH-INRTEQV) to the same n-octane transport data. Qualitatively, the same behavior is similar to that shown in Figure 5. We note that the approximation made in mapping the gyration ellipsoid to a prolate (eq 21) is somewhat different in the two cases. As an example, for n-octane at 293 K and 0.1 MPa, the INRTEQV model yields average gyration ellipsoid semiaxes of 5.96, 1.49, and 0.68 Å, for , , and , respectively. The VMTCH-INRTEQV model yields 4.28, 1.18, and 1.14 Å for 〈a〉, 〈b〉, and 〈c〉. The ellipsoid shape of the discrete massdistribution model is more closely matched to the prolate shape (where 〈b〉 ≈ 〈c〉) than the model which smears the mass distribution of the molecule uniformly over the ellipsoid volume. Figure 7 depicts the VMTCH-VSPHPL prolate ellipsoid model with both stick and slip boundary condi-
Figure 8. Stokes-Einstein plot of n-octane transport data using the equivalent sphere model, with volume determined from method of overlapping spheres (EQVSPH). Symbols and lines are the same as Figure 5.
tions. This characterization is identical to that shown in Figure 6, except that the volume of the ellipsoid is determined through the overlapping spheres method. The slope of the data represented by the slip model very nearly matches the dotted line in the figure, while the stick model still deviates significantly, underpredicting the slope of the line. We note that, in general, the SPHPL volume is much larger than the gyration ellipsoid volume determined in the previous figures. For example, for n-octane at 293 K and 0.1 MPa, we obtain VSPHPL ) 149 Å3/molecule, while the volume obtained for the ellipsoids in the INRTEQV and VMTCHINERTEQV are ≈ 25 Å3/molecule. In the representations of Figure 7 it is apparent that the slip description yields a better representation of the data, particularly at lower values of viscosity (high T, low P). However, we observe a departure from the slip model as viscosity increases. This may be due simply in part to the larger uncertainties associated with our viscosity data at extreme conditions (as indicated by the error bars associated with the points on the curves). However, this departure may be linked to an increase in translational-rotational coupling with increasing density; as free volume decreases, more energy dissipation occurs through rotational motion rather than through translation. This effectively lowers the expected diffusion relative to the slip model prediction. While the transport data in this figure can still be well represented by a straight line, this sort of deviation from the slip model prediction has been regarded by others as a failure of the Stokes-Einstein equation to correlate transport data for pure fluids.21,22,32 In Figure 8 we depict the EQVSPH representation of the data, where the molecule is characterized by an effective particle radius. In this case, the stick and slip representations differ only by a constant. We find that the slip model is closer to the expected slope ) 1 line than the stick model. The effective radius of the n-octane ‘particle’ at 293 K and 0.1 MPa is 3.3 Å and is found to be fairly insensitive to both temperature and pressure. This value is close to that obtained from the Hildebrand fluidity analysis15 (eq 23) applied to the data for the 293 K isotherm, which yields a value of approximately 3.1 Å. Stick vs Slip - Tests across Many Isomers. To compare the Stokes-Einstein models across different molecular architectures, we have generated the trans-
7032 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
Figure 9. INRTEQV model applied to transport data of 78 isoparaffins between C6 and C16 at T ) 293 K and 373 K, and P ) 0.1 MPa. Symbols and lines are the same as Figure 5.
Figure 10. VMTCH-INRTEQV model applied to transport data of 78 isoparaffins between C6 and C16 at T ) 293 K and 373 K, and P ) 0.1 MPa. Symbols and lines are the same as Figure 5.
port properties of 78 different isoparaffin isomers between C6 and C16 in the bulk at 0.1 MPa and 2 temperatures, 293 K and 373 K. In the following comparisons, we are primarily interested in the ability of each of these representations to collapse the transport data from a variety of different molecular architectures onto a single curve. Application of the INRTEQV model applied to our data is shown in Figure 9 in both the stick and slip representations. With the stick boundary condition model, the same general trend is observed as with the n-octane transport data of Figure 5, but there is a fair degree of scatter. By comparison, most of the data in the slip prolate model collapses more uniformly onto a line. We consider this to be strong evidence of the superior representation of the slip prolate model. Deviations in the Stokes-Einstein slopes are similar to that of Figure 5, again a consequence of equating the molecular volume to the prolate ellipsoid volume. While the average slope of the slip model data is larger than that predicted by the model, it would appear that the scaling of the prolate ellipsoid volume and shape capture a good portion of the transport property variation with molecular architecture. Figure 10 depicts the transport data with the VMTCHINRTEQV model in the stick and slip limits. While the main features of the plot are essentially similar to that of Figure 9, the semiaxes obtained more closely approximate a prolate ellipsoid shape than the INRTEQV
Figure 11. VMTCH-VSPHPL model applied to transport data of 78 isoparaffins between C6 and C16 at T ) 293 K and 373 K, and P ) 0.1 MPa. Symbols and lines are the same as Figure 5.
model. This may be a fine distinction, but we judge this model to be preferable because this formalism of the more natural mapping to the prolate ellipsoid geometry required in application of eq 21. It is also worth noting that for both the stick and slip models in Figures 9 and 10, the long, linear or near-linear paraffins such as n-hexadecane, 8-methylpentadecane, and 2-methylpentadecane, lie somewhat below the rest of the data. Among the species considered in our study, these molecules possess the greatest potential for intramolecular conformational rearrangements. While not accounted for in the model, it appears that the increased molecular flexibility of these species results in enhanced diffusion and/or increased viscosity relative to shorter, less flexible molecules. Figure 11 shows the data interpreted through the VMTCH-VSPHPL model representation. While the slip model again shows improvement over the stick model, there is comparatively more scatter in the slip model compared with the representations of Figures 9 and 10. A consequence of the spheres-planes volume approach of this model is that the molecular volume does not vary significantly among isomers at a given carbon number. Since the aspect ratio of the prolate ellipsoids in this representation is equivalent to that of Figure 10, it appears that the hydrodynamic volume employed in this approach modestly is the source of the increase in scatter. For completeness, we show in Figure 12 the EQVSPH model for the isoparaffin transport data for stick boundary conditions. The slip model differs only by a constant CBC. The scatter in the curve again indicates that this simplified molecular characterization is too coarse to capture the relationship between transport data with molecular architecture. Application of Stokes-Einstein Models: Prediction of Pure Species Viscosity. Figures 9 and 10 suggest that the drag model for the prolate ellipsoid with slip boundary conditions, combined with a reasonable estimate of the molecular volume, offer a realistic representation of the relationship between self-diffusivity and viscosity. From a practical standpoint, however, the strong degree of linearity observed between D and η in any of the representations depicted in Figures 5-8 allows for extrapolation of transport properties to extremes in temperature and pressure. Assuming we can quantify the observed linearity between D and η, we can make predictions of viscosity
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 7033
Figure 12. EQVSPH model applied to transport data of 78 isoparaffins between C6 and C16 at T ) 293 K and 373 K, and P ) 0.1 MPa. As the slip model differs only by a constant, only the stick limit model is plotted. Figure 14. Viscosity vs temperature predicted for (0) 2,2,3trimethylpentane and (4) 2,2,4-trimethylpentane from the StokesEinstein relationships of Figure 13. The viscosity is estimated from the computed self-diffusion coefficient and the VMTCH-INRTEQV model characterization of the molecule’s size and shape. Corresponding filled symbols represent Green-Kubo viscosity calculations, along with estimated uncertainties.
Figure 13. VMTCH-INERTEQV model in the slip limit for prolate ellipsoids for (- - -) 2,2,3-trimethylpentane and (‚‚‚‚‚) 2,2,4-trimethylpentane. Lines represent best fits to transport data over a range of temperatures and pressures. Points correspond to additional simulations performed at lower temperatures and atmospheric pressure: (2) 2,2,3-trimethylpentane and (0) 2,2,4trimethylpentane. See text for additional details.
at a given temperature and pressure with knowledge of the self-diffusion coefficient and the molecular volume. Figure 13 demonstrates this idea for 2 C8 isomers, 2,2,3-trimethylpentane and 2,2,4-trimethylpentane, using the VMTCH-INRTEQV model in the slip boundary condition limit. The dashed curves represent the fitted lines to a set of transport data collected at three temperatures (293, 353, 410 K) and a range of pressures (between atmospheric and 600 MPa). The slopes of these curves are both approximately 1.5. Simulations were then performed at lower temperatures, down to 183 K. Assuming the linear relationship holds, we can use the self-diffusion coefficient, and the prolate ellipsoid size and shape, to predict the viscosity. Knowledge of the self-diffusion coefficient enables us to read off the predicted value of the abscissa, from which the only unknown is the viscosity. This read-across is depicted in the figure for 2,2,3-trimethylpentane at 233, 203, and 183 K. The Green-Kubo viscosity calculations of these states, along with several other low temperature states
Figure 15. Stokes-Einstein plot using VMTCH-INERTEQV model in the slip limit for (4) n-dodecane and (0) 5-n-propylnonane. The open symbols represent the data from which the Stokes-Einstein lines were parametrized. Corresponding filled symbols represent additional simulations at state points beyond the original set of simulations. Error bars for points on the abscissa are derived from the viscosity calculations.
for 2,2,4-trimethylpentane, are labeled as points on the curve. Note the excellent agreement with the viscosity computed via molecular dynamics and the StokesEinstein prediction made solely from diffusivity and molecular shape. Figure 14 plots the predicted viscosity as a function of temperature for these isomers as well as the Stokes-Einstein predictions. It is important to emphasize that we have made an extrapolative quantitative prediction of the model fluid viscosity, up to 110 K below the temperature range with which we determined the Stokes-Einstein slope. A similar example is shown in Figure 15 for ndodecane and 5-n-propylnonane. In this case, the open symbols represent the data from which the best fit
7034 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003
Figure 16. Predicted viscosity vs pressure for five isomers of C12 paraffins at 293 K. Viscosity is estimated from self-diffusion coefficients in the same manner as described for Figure 14. The different symbols correspond to (]) n-dodecane, (0) 5-n-propylnonane, (×) 2,4-dimethyldecane, (4) 2-methyl-5-ethylnonane, and (*) 2,6-dimethyl-4-isopropylheptane.
Stokes-Einstein slopes were obtained, which were 1.14 and 1.51, respectively. For n-dodecane, the range of state points covered by these data were isotherms at 293 K (0.1-150 MPa), 358 K (0.1-300 MPa), and 423 K (0.1-450 MPa). For 5-n-propylnonane, the data covers 293 K (0.1-250 MPa), 358 K (0.1-300 MPa), and 423 K (0.1-600 MPa). Additional simulations at higher pressures were conducted, and several points are shown as filled symbols in the figure. These points agree quite well with the Stokes-Einstein lines, even though they represent an extrapolation of the curves. Following the same methodology, we depict in Figure 16 the predicted viscosity of several dodecane isomers at high pressure via the corresponding Stokes-Einstein relationships for prolate ellipsoids. The rapid rise in viscosity of the more highly branched molecules is notable, particularly for the extremely branched example of 2,6-dimethyl-4-isopropylheptane. These predictions agree well with viscosity computed via molecular dynamics for the state points for which we have data. Most importantly, the method enables predictions of viscosity at state points that would be prohibitively expensive to compute directly via equilibrium molecular dynamics. The preceding predictions illustrate the utility of this approach as a scoping and predictive tool for the molecular design of lubricants. Relatively minor molecular rearrangements can give rise to unexpectedly large differences in viscosity at low temperature and high pressure. For extremes in pressure, it has been demonstrated experimentally that fluid viscosity at high pressure is quite sensitive to molecular structure, which has a direct consequence on lubrication properties.33 The ability to explore the effect of such structural changes on rheological response at extreme conditions with molecular simulation tools is an approach we expect to see applied with increasing frequency in the future. Discussion Computational Requirements. Although the above examples may appear at first blush to be a roundabout
means of estimating viscosity, the real potential of the methodology becomes clear upon considering the scaling of the computational expense of the viscosity calculation. Reasonable convergence of the stress autocorrelation function integral necessary to compute viscosity via eq 25 generally requires fairly long simulation times. This time can be approximated as Nητrot, where τrot is the rotational relaxation time of the individual molecules in the system (eq 27), and Nη is an integer multiple of this time scale needed for convergence of eq 25, typically between 200 and 300 for the isoparaffin systems considered here. With increasing carbon number, τrot increases, necessitating longer total simulation times to achieve the same degree of convergence of the viscosity integral. In addition, the time required to integrate the equations of motion through a single time step is dominated by the calculation of forces among the N atomic sites in the system; this calculation time scales roughly as NlogN. Combining these effects, we can estimate the time required to obtain a reasonable estimate of the viscosity for united atom model paraffins as
Tη ) K0Nητrot(T,P)[NmolClog(NmolC)]
(28)
where τrot is a function of the state conditions of the system and the particular molecule of interest, Nmol is the number of molecules in the simulation, C# is the carbon number of the paraffin, and K0 is a constant of proportionality related to the attributes of the hardware configuration and simulation software. The self-diffusion coefficient and rotational relaxation time, on the other hand, are related to the velocity autocorrelation function and orientational autocorrelation function, respectively. Both of these quantities are averaged over the Nmol molecules in the system; thus, every time step, Nmol new pieces of information are collected in the evaluation of eqs 24 and 26. These integrals converge after just a few multiples of τrot. Denoting ND such that NDτrot is the simulation time required for convergence of the ‘N’ particle autocorrelation functions, these properties are generally obtained reliably for ND ) 5-10. With this simple analysis, we now consider the estimated time required to compute viscosity and selfdiffusivity for a system of 120 n-alkane molecules interacting with the TRAPPE-UA potential at 373 K and 0.1 MPa. We have computed the rotational relaxation time of these systems as a function of carbon number up to C26, and the scaling is well-approximated as
1 5 τrot(ps) ) C2 - C 3 4
(29)
In Figure 17 we depict the estimated time required for the viscosity computation as a function of carbon number. In this figure we have assumed Nη ) 300 and ND ) 10. Given a benchmark of 1.3 h for the reference n-hexane system, the scaling of the expected viscosity and diffusion calculations is compared. We see that, for the same degree of convergence of the viscosity integral, n-C35 will require a factor of 103 more effort than the n-C6 system. While the n-C35 viscosity calculation becomes prohibitively expensive, exceeding 1 month per datum, the diffusion calculation is achievable at a small fraction of this cost.
Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 7035
remaining scatter in Figures 9 and 10 and make predictions for viscosity for any isoparaffin. Conclusions
Figure 17. Estimated time requirements vs carbon number for transport property calculations of n-alkane systems at 373 K and 0.1 MPa. The estimates assume 120 molecules in a system. Assuming a benchmark calculation for n-hexane takes 1.3 h, the symbols depict estimated time required for (9) viscosity and (2) diffusion calculations.
Possible Improvements. The linear relationships observed in the Stokes-Einstein plots of 13 and 15 demonstrate how we can make extrapolative predictions of viscosity from self-diffusion and geometric considerations, once we have determined the characteristic slope of the plot. However, for the pure species in which we have characterized transport properties over a wide range of conditions, we still observe some variation in the slope of the lines as a function of molecular structure. Across many isomers, we see strong evidence that the Stokes-Einstein formalism for prolate ellipsoids in the slip limit is capable of capturing a large portion of the structural variation observed among isomers. While clearly a superior representation to the models of spherical particles, there is still some variation in the data plotted in Figures 9 and 10. The key obstacle to overcome is to develop an understanding of the small differences in the Stokes-Einstein slope obtained as a function of molecular structure. One obvious deficiency that is neglected in the StokesEinstein approach is the incorporation of molecular flexibility. Hydrodynamic equations for drag in creeping flow are developed for rigid bodies; yet as carbon number increases, we expect that dihedral angle rotations in the molecule will contribute to increasing conformational flexibility. Increased molecular flexibility should manifest itself in a lower overall slope in the Stokes-Einstein plots for pure species. This is observed, for example, among the C12 isomers in Figure 15, where the slope of the more flexible n-dodecane molecule is approximately 35% lower than 5-n-propylnonane. To make predictions for viscosity of any given isoparaffin, we must quantify the effect of molecular flexibility on the resulting slopes obtained in the previous figures. We are currently undertaking simulations of additional pure species over an extended range of temperature and pressure in order to quantify the variation in the Stokes-Einstein slope for the prolate ellipsoid model as a function of structure. Given this additional data, it should be possible to incorporate an additional parameter based on molecular flexibility to reduce the
In this paper, we have examined the potential of Stokes-Einstein relationships to link translational diffusion to viscosity. We have investigated several formulations of such models, testing to see how well various hydrodynamic approximations can represent molecular motions in liquids in the bulk. In particular, the hydrodynamic drag for prolate ellipsoids in the slip limit appears capable of explaining much of the structural variation among isoparaffins in transport properties when placed in the Stokes-Einstein formalism. We have demonstrated that one can exploit the linearity of these types of models to make quantitative extrapolative predictions of viscosity from computer simulations from the more rapidly converging quantities of self-diffusion and gyration ellipsoids. Further work to quantify the role of molecular flexibility on the resulting StokesEinstein slopes is ongoing and should yield a universal equation capable of making viscosity predictions for chain molecules of arbitrary architecture. Acknowledgment The author would like to thank Larry Dodd (for permission to use) and Michael Greenfield and Ed Maginn (for copies of) for the molvol code used in the analytical determination of the volume of overlapping spherical assemblies of particles. Thanks also to ExxonMobil Research & Engineering for support of this work and permission to publish the results. Literature Cited (1) Gordon, P. A. Modeling Structure Property Relationships in Synthetic Basestocks. In ACS Symposium Series: Dynamics and Friction at Submicron Confining Systems, American Chemical Society: 2002; in press. (2) Martin, M. G.; Siepmann, J. I. Novel Configurational-Bias Monte Carlo Method for Branched Molecules. Transferable Potentials for Phase Equilibria. 2. United-Atom Description of Branched Alkanes. J. Phys. Chem. B 1999, 103, 4508. (3) Happel, J.; Brenner, H. Low Reynolds Number Hydrodynamics; Prentice-Hall: Englewood Cliffs, NJ, 1965. (4) Lamb, H. Hydrodynamics, 6th ed.; Dover Publications: New York, 1945. (5) Einstein, A. Investigations on the Theory of Brownian Movement; Dover Publications: New York, 1926. (6) Tirado, M. M.; Garcia de la Torre, J. Translational Friction Coefficients of Rigid, Symmetric Top Macromolecules. Application to Circular Cylinders. J. Chem. Phys. 1979, 71, 2581. (7) Tang, S.; Evans, G. T. A Critique of Slip and Stick Hydrodynamics for Ellipsoidal Bodies. Mol. Phys. 1993, 80, 1443. (8) Ravichandran, S.; Bagchi, B. Anisotropic Diffusion of Nonspherical Molecules in Dense Liquids: A Molecular Dynamics Simulation of Isolated Ellipsoids in the Sea of Spheres. J. Chem. Phys. 1999, 111, 7505. (9) Mondello, M.; Grest, G. S. Molecular Dynamics of Linear and Branched Alkanes. J. Chem. Phys. 1995, 103, 7156. (10) Mondello, M.; Grest, G. S. Viscosity Calculations of nalkanes by Equilibrium Molecular Dynamics. J. Chem. Phys. 1997, 106, 9327. (11) Daivis, P. J.; Evans, D. J. Computer Simulation Study of the Comparative Rheology of Branched and Linear Alkanes. J. Chem. Phys. 1992, 97, 616. (12) Zhang, Y.; Venable, R. M.; Pastor, R. W. Molecular Dynamics Simulations of Neat Alkanes: The Viscosity Dependence of Rotational Relaxation. J. Phys. Chem. 1996, 100, 2652. (13) Desloge, E. A. Classical Mechanics; John Wiley & Sons: New York, 1982.
7036 Ind. Eng. Chem. Res., Vol. 42, No. 26, 2003 (14) Dodd, L. R.; Theodorou, D. N. Analytical Treatment of the Volume and Surface Area of Molecules Formed by an Arbitrary Collection of Unequal Spheres Intersected by Planes. Mol. Phys. 1991, 72, 1313. (15) Hildebrand, J. H.; Lamoreaux, R. H. Fluidity: A General Theory. Proc. Natl. Acad. Sci. U.S.A. 1972, 69, 3428. (16) Hildebrand, J. H. Motions of Molecules in Liquids: Viscosity and Diffusivity. Science 1971, 174, 490. (17) Dymond, J. H. Corrected Enskog Theory and the Transport Coefficients of Liquids. J. Chem. Phys. 1974, 60, 969. (18) Parkhurst, H. J., Jr.; Jonas, J. Dense liquids. II. Effect of density and temperature on viscosity of tetramethylsilane and benzene. J. Chem. Phys. 1975, 63, 2705. (19) Parkhurst, H. J., Jr.; Jonas, J. Dense liquids. I. Effect of density and temperature on self-diffusion of tetramethylsilane and benzene-d6. J. Chem. Phys. 1975, 63, 2698. (20) Jonas, J.; Hasha, D.; Huang, S. G. Density effects of transport properties in liquid cyclohexane. J. Phys. Chem. 1980, 84, 109. (21) Jonas, J.; Hasha, D.; Huang, S. G. Self-diffusion and Viscosity of Methylcyclohexane in the Dense Liquid Region. J. Chem. Phys. 1979, 71, 3996. (22) Fury, M.; Munie, G.; Jonas, J. Transport Processes in Compressed Liquid Pyridine. J. Chem. Phys. 1979, 70, 1260. (23) Finney, R. J.; Fury, M.; Jonas, J. Density and Temperature Dependence of Self-diffusion and Shear Viscosity of Perfluorocyclobutane in the Dense Fluid Region. J. Chem. Phys. 1977, 66, 760. (24) Kioupis, L. I.; Maginn, E. J. Molecular Simulation of Polya-olefin Synthetic Lubricants: Impact of Molecular Architecture on Performance Properties. J. Phys. Chem. B 1999, 103, 10781.
(25) Kioupis, L. I.; Maginn, E. J. Impact of Molecular Architecture on the High-Pressure Rehology of Hydrocarbon Fluids. J. Phys. Chem. B 2000, 104, 7774. (26) Gordon, P. A. Influence of Simulation Details on Thermodynamic and Transport Properties in Molecular Dynamics of Fully Flexible Molecular Models. Mol. Sim. 2003, 29, 479. (27) Haile, J. M. Molecular Dynamics Simulation; John Wiley & Sons: New York, 1992. (28) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (29) Cui, S. T.; Cummings, P. T.; Cochran, H. D. The Calculation of the Viscosity from the Autocorrelation Function Using Molecular and Atomic Stress Tensors. Mol. Phys. 1996, 88, 1657. (30) Allen, M. P. Atomic and Molecular Representations of Molecular Hydrodynamic Variables. Mol. Phys. 1984, 52, 705. (31) Williams, G.; Watts, D. C. Non-Symmetrical Dielectric Relaxation Behavior Arising from a Simple Empirical Decay Function. Trans. Faraday Soc. 1970, 66, 80. (32) Walker, N. A.; Lamb, D. M.; Adamy, S. T.; Jonas, J.; DareEdwards, M. P. Self-diffusion in the Compressed, Highly Viscous Liquid 2-ethylhexylbenzoate. J. Phys. Chem. 1988, 92, 3675. (33) Bair, S. The High-Pressure Rheology of Some Simple Model Hydrocarbons. Proc. Instn. Mech. Engrs. Part J: J. Eng. Tribology 2002, 216, 139.
Received for review June 16, 2003 Revised manuscript received September 30, 2003 Accepted October 9, 2003 IE030512X
