Extrapolation of Rheological Properties for Lubricant Components with

Quantifying physical properties of Weyburn oil via molecular dynamics simulation. M. Uddin , D. Coombe , J. Ivory. Chemical Engineering Journal 2016 3...
1 downloads 0 Views 147KB Size
Download PDF
5828

Ind. Eng. Chem. Res. 2005, 44, 5828-5835

Extrapolation of Rheological Properties for Lubricant Components with Stokes-Einstein Relationships Peter A. Gordon† Corporate Strategic Research, ExxonMobil Research and Engineering, 1545 Route 22 East, Annandale, New Jersey 08801

In this paper we examine the scaling of transport properties obtained from molecular dynamics simulations on a number of model C16 isoparaffin systems and interpreted in a Stokes-Einstein framework. The properties span a wide range of temperatures and pressures. The scaling behavior can be nonlinear and shows variation among different isomers. Nonetheless, for a given molecule, the scaling behavior of transport property data appears to be well-behaved, and we propose two approximate methods by which viscosity can be estimated from the self-diffusion coefficient and geometric characterization of the molecule in terms of an inertially equivalent prolate ellipsoid. Both methods rely on high-temperature simulations, which provide important scaling information about transport properties with reasonable computational effort. We compare the accuracy of the predictions of these approaches with several examples of lube-ranged model compounds. Introduction In previous work1 on modeling structure-rheology dependencies, we have explored various formulations of Stokes-Einstein relationships in order to develop a formalism capable of rationalizing transport data obtained from equilibrium molecular dynamics (EMD) simulations. Different models describing frictional drag were considered, along with different interpretations of molecular shape and volume. The success of these formulations was judged on the basis of the ability to rationalize transport data of a wide range of low molecular weight (C6-C16) isoparaffins of differing chain architectures described with the TraPPE-UA2 intermolecular potential. It was found that a hydrodynamic drag model for creeping flow past prolate ellipsoids with slip boundary conditions at the surface was best able to scale the transport data obtained. The quality of this model is depicted in Figure 1, which plots self-diffusivity vs the friction coefficient for nearly 80 isoparaffin isomers at 293 and 373 K. In this comparison, an approximation to the frictional drag, detailed in this paper, must necessarily be employed. This model shows a modest degree of variation in the scaling of transport data that appears to have some residual molecularstructure dependence, but represents a significant improvement over other formalisms considered. By contrast, a description that exactly captures the frictional resistance would be capable of collapsing all the transport data in the figure onto a single unified curve. A Stokes-Einstein model that can collapse transport data of many different molecular structures onto a single curve could be exploited to estimate the zeroshear viscosity from the basis of the self-diffusion coefficient and geometric considerations of the molecule. This would enable predictions of the viscosity of large molecules, where direct calculations of the viscosity can be prohibitively expensive from a computational resource standpoint. Direct calculation of viscosity via † Phone: 908-730-2546. Fax: [email protected].

908-730-3031.

E-mail:

Figure 1. Transport data of selected TraPPE-UA model C6-C16 isoparaffins at P ) 0.1 MPa and 293 and 373 K interpreted through the Stokes drag coefficient past a prolate ellipsoid in the slip limit. The dashed line indicates a “best fit” of the data.

EMD becomes time-intensive for systems that exhibit long memory, or “slow” dynamics. Molecules typical of lubricant basestocks are in the C20-C50 (or higher) carbon number range and, when subjected to conditions typical of elastohydrodynamically lubricated contacts, reach pressures of 1 GPa or more. The combination of high molecular weight and extreme conditions can lead to very long memory effects, making direct viscosity calculations quite laborious. Given the computational burden associated with systems exhibiting slow dynamics, it is tempting to employ the best-fit line in Figure 1 as a “universal” curve describing the transport properties of molecular systems. However, when transport property data over a wide range of conditions are considered for a given molecule, we have found significant species-dependent variations in the scaling. This is demonstrated in Figure 2, which depicts the transport properties of n-hexadecane and 5-n-propyl-5-n-butylnonane, two C16 isoparaffins, over a range of conditions of T and P. The

10.1021/ie050156d CCC: $30.25 © 2005 American Chemical Society Published on Web 06/15/2005

Ind. Eng. Chem. Res., Vol. 44, No. 15, 2005 5829

Methods We begin with a brief review of the hydrodynamic model and the geometric characterization of molecules employed in the analysis of transport properties within the context of a Stokes-Einstein formalism. The StokesEinstein relationship, in a general sense, can be expressed as

kT ) fdrag ≈ CBCη × g(molecular size, shape, D boundary conditions) (1)

Figure 2. Diffusion vs Stokes drag coefficient for TraPPE-UA model (b) n-hexadecane, (0) 5-n-propyl-5-n-butylnonane, and the best fit line derived from Figure 1. The transport data cover three isotherms over a range of pressures: T ) 333 K, P ) {0.1 f 125 MPa}; T ) 393 K, P ) {0.1 f 550 MPa}; T ) 453 K, P ) {0.1 f 800 MPa}.

universal curve of Figure 1 appears to trace the transport properties of 5-n-propyl-5-n-butylnonane, but diverges from the n-hexadecane data. To develop approaches that can estimate the rheological properties of systems at extreme states, an alternative to that suggested by Figure 1 is needed. Despite the apparent structure-sensitive variations, the data in Figure 2 appear to exhibit a well-defined relationship for a given molecule. This led us to wonder about the range of validity of the transport data scaling. In particular, very high temperature state points can be examined for higher molecular weight species at conditions where the viscosity calculation converges with reasonable computational effort. These states may not be of interest, or even physically realistic (i.e., at temperatures above the expected decomposition point of such materials). If the scaling of the transport data holds over such an extended range, we should be able to use high-temperature results to make predictions of viscosity at conditions of temperature and pressure of interest. This paper is organized as follows. We begin by briefly reviewing the details of the Stokes-Einstein formalism used to characterize transport properties in this work, outlining the procedure by which we characterize molecular shape and volume. We devote considerable attention to the transport properties of five different systems of C16 isoparaffin isomers with varying degrees of branching, computed over a very wide range of temperatures and pressures. We examine a number of features of the scaling behavior of the transport data, identifying shape and density-dependent departures from hydrodynamic behavior. Based on these observations, we propose two distinct methods by which viscosity can be estimated from self-diffusion and geometric considerations alone. We compare these approaches for estimating viscosity of the C16 systems with direct Green-Kubo calculations. Finally, we present a number of tests of these viscosity estimation methods applied to representative synthetic lubricant basestock molecules, systems for which GreenKubo viscosity calculations are infeasible.

where D is the self-diffusion coefficient and fdrag is the Stokes drag coefficient. The form of drag depends on how the low Reynolds number hydrodynamic problem is formulated, but generally includes a constant, CBC, that depends on the boundary conditions applied on the surface of the object, the viscosity, η, and a description of the molecular shape and size. For creeping flow past a prolate ellipsoid with long axis a, short axis b, and aspect ratio φ ) a/b, g can be expressed as

g ) b × KAV

(2)

where KAV is the average resistance to translation (assuming isotropic translation) resolved along directions parallel (K|) and perpendicular (K⊥) to the long axis of the prolate.

[

1 2 1 1 ) + KAV 3 K⊥ K|

]

(3)

Using slip boundary conditions on the surface of the ellipsoid, Tang et al.3 have computed numerically the resistances parallel and perpendicular to the long axis of a prolate. For 1 e φ e 10, these resistances can reasonably be approximated as

K||(φ) )

1 2 + exp(-0.3φ) 2 3

K⊥(φ) ) 1 + 0.54(φ - 1) - 0.011(φ - 1)2

(4) (5)

and CBC in eq 1 is 4π. Precise definitions of the molecular shape and volume are open to interpretation for flexible molecular species. The approach that appeared to be the most promising in our previous work, and adopted here, is to employ the inertia ellipsoid as the basis of the characterization. We can express this through the diagonalized moment of inertia tensor of the molecule 3

λR|e bR | 2 ) γ ∑ R)1

(6)

where λR and b eR refer to the eigenvalue and eigenvector of the inertia tensor for the molecule in the R direction. Here it is assumed that the eigenvalues are presented in ascending order so that the longest axis of the molecule corresponds to λ1. The quantity γ is an arbitrary constant which scales the volume of the gyration ellipsoid; the length of the longest semiaxis of the ellipsoid is

x

γ λ1

5830

Ind. Eng. Chem. Res., Vol. 44, No. 15, 2005

We specify γ by comparing the inertia tensor of the molecule, with discrete mass points centered on atoms, to that of an inertially equivalent ellipsoid with uniform mass density. The latter object is bounded by the surface

(ax) + (by) + (zc) 2

2

2

)1

(7)

and would be inertially equivalent to the mass distribution of the molecule if its principal moments are

a2 )

5 (λ + λ3 - λ1) 2M 2

b2 )

5 (λ + λ3 - λ2) 2M 1

c2 )

5 (λ + λ2 - λ3) 2M 1

(8)

where M is the total mass of the molecule. With a, b, and c defined by eq 8, we can require that the volume of the uniform density ellipsoid given by eq 7 equals that of discrete mass distribution defined by the molecule (eq 6). This specifies γ in eq 6 as

γ ) (abc)2/3 (λ1λ2λ3)1/3

(9)

During a molecular dynamics simulation of a fluid at the conditions of interest, the gyration ellipsoid defined by eqs 6-9 are periodically evaluated for each molecule in the system. This produces complete geometric characterization of molecular size and shape. The mapping of the three-dimensional mass distributions to that of the axisymmetric prolate ellipsoid shape required to use eqs 5 and 4 is accomplished through the geometric mean of the small axis of the gyration ellipsoid,

xγ bAV ) (λ2λ3)0.25 φAV )

(λ2λ3)0.25

xλ1

(10)

(11)

The quantities bAV and φAV are time-averaged during a simulation and are used (in place of b and φ) to compute the average translational resistance in eqs 2-5. The molecular dynamics simulations have been conducted in LAMMPS4 and modifications therein. Simulation protocols have been detailed in previous work, and the interested reader is referred to these references for additional details.1,5,6 Except where noted, the intermolecular potential employed in our calculations is the TraPPE united atom model; details on the potential form and parameters can be found elsewhere.2 We note one minor modification to the potential for use in hightemperature simulations. In the united atom description for branched molecules employed in previous work, we have adopted the harmonic umbrella improper torsion potential suggested by Mondello et al.7 For a quartet of atoms i-j-k-l where atoms j, k, and l are all bonded to i,

Uχ ) Kχ(χ - χ0)2

(12)

where χ is the angle between the planes defined by atoms i-j-k and j-k-l, χ0 is the equilibrium angle, and

Figure 3. Stokes-Einstein plot for n-hexadecane. (0) Isotherms at state points shown in Figure 2, (b) state points at high temperature and a range of densities, and (- - -) best fit curve through high-temperature transport data. See text for additional details.

Kχ is the force constant. The purpose of this term is to overcome a physical defect inherent only to united atom models; it prevents single carbon branch points from an unphysical inversion. This potential is symmetric with respect to χ ) 0o. Test simulations at high temperatures showed that the potential barrier about χ ) 0 was insufficient to prevent inversions. For this reason, we have employed a modified form of the inversion potential, which seeks to maintain the harmonic character of the potential well and prevent inversion flips,

Uχ ) Kχ(χ - χ0)2 + Kχ2exp

(-|χ| ∆ )

(13)

The second term in eq 13 requires symmetry breaking of the improper torsion; this can be accomplished by monitoring the sign of the scalar triple product b rlk‚(r bkj ×b rij). Values of Kχ2 ) 20 kcal/(mol rad2) and ∆ ) 2o were found to be sufficient to prevent inversions across the χ ) 0 barrier at the highest temperature simulations undertaken here. Results and Discussion C16 Transport Properties: n-Hexadecane. Figure 3 depicts the Stokes-Einstein behavior of n-hexadecane over an extended range of states. In addition to the data shown in Figure 2, simulations at 980 K and a range of densities are plotted. Note that even at very high temperatures, most of these points line up with the data obtained at lower temperatures. In general, the hightemperature transport data are obtained with good precision at a fraction of the computational expense of the state points at ambient temperatures and elevated pressures. Of particular interest are the points on the far left of the plot that show noticeable deviations from the best fit of the data. These points correspond to states at a supercritical temperature of 980 K (note: Tc ≈ 725 K) and F ) 0.1 and 0.2 g/cm3, respectively, the lowest densities considered. One may explain the departures of the transport properties of these state points from the general StokesEinstein scaling of the rest of the data in terms of the conditions necessary for Brownian motion. A particle

Ind. Eng. Chem. Res., Vol. 44, No. 15, 2005 5831

Figure 4. Normalized force autocorrelation function n-hexadecane at selected states. (s) T ) 293 K, F ) 0.788 g/cm3, (-‚-‚-) T ) 980 K, F ) 0.5 g/cm3, (- - -) T ) 980 K, F ) 0.3 g/cm3, (‚‚‚) T ) 980 K, F ) 0.1 g/cm3.

undergoing Brownian motion may be considered to experience forces of two general types. The first is the frictional force due to the drag (with proportionality defined by drag coefficient f) exerted on the molecule by the rest of the fluid. The second is considered a rapidly fluctuating, random force whose origin is the constant molecular collisions exerted by the surrounding fluid. The Langevin equation summarizes the governing motion,8

A B (t) du b ) -fu b+ dt m

(14)

It is assumed that the random force, A B , varies much more rapidly than variations in b u, the molecular velocity. This implies that there are time intervals δt over which changes in the velocity are small while correlations in the random force have decayed. We can test this assumption by focusing on the appropriate autocorrelation functions of our (molecular) Brownian particle. We expect that the time scale of decay of the autocorrelation function of the net force acting on a molecule should be fast for a particle undergoing rapid collisions with its neighbors. This behavior is depicted in Figure 4 at several different state points for n-hexadecane. At room temperature and atmospheric pressure, the normalized force autocorrelation function exhibits a steep decrease at small times. The negative region centered around 0.125 ps is typical of the rapid “cage-like” collisions experienced by molecules with near-neighbors in a liquid structure. The correlation function is seen to decay rapidly to zero in a few picoseconds. By contrast, at 980 K and 0.1 g/cm3, the magnitude of the negative trough has decreased substantially, and the presence of a longer-time tail in the autocorrelation function emerges. At 980 K and intermediate density (0.3 g/cm3), the negative trough in the autocorrelation function is still present, but the decay to zero occurs on roughly the same time scale as the state point at 293 K. Referring back to Figure 3, the 980 K, 0.3 g/cm3 state point follows the general scaling of the transport data, but the low-density point falls significantly below the curve. It would appear that the presence of the long-time tail in the force autocorrelation function, indicative of more ballistic-like tra-

Figure 5. Stokes-Einstein plot of transport data of selected C16 isoparaffins. Symbols and axes labels the same as Figure 3.

jectories of molecules, is responsible for the observed deviations in the scaling of transport properties. C16 Isoparaffin Transport Properties: Variations among Isomers. Transport properties of four isoparaffin isomers of C16 are shown in Figure 5. The range of conditions covered by the data is similar to that of Figure 3, although we have omitted high-temperature low-density state points that show deviations from the general scaling of the data. For all cases, the hightemperature transport property data obey the same general scaling as the rest of the lower temperature results. This suggests that we may fit the high-temperature transport data to a simple function and use this relationship to extrapolate to extreme conditions. For each of the 5 C16 isomers considered, we have found that power law expressions of the form y ) A0xA1 are necessary to adequately fit the high-temperature transport data. Table 1 summarizes the power law fit parameters for each species. The power law fits are also plotted in Figures 3 and 5 and appear to possess good extrapolative power. Thus, at a given temperature and pressure, we can compute the self-diffusion coefficient and prolate ellipsoid geometric factors, and estimate viscosity through rearrangement of the power law equation,

ηSES )

( )

kT 1 4πbAVKAV A0D

1/A1

(15)

In what follows, we denote this approach the StokesEinstein scaling (SES) method of estimating viscosity. From the point of view of hydrodynamics, it is somewhat surprising that the data in Figures 3 and 5 scale with varying degrees of nonlinearity. We can learn something about the nature of these departures by considering previous studies on the transport properties of hard ellipsoids.9,10 In these works, molecular dynamics simulations were conducted for prolates of various aspect ratios over a range of densities, and the results were compared with predicted hydrodynamic behavior

5832

Ind. Eng. Chem. Res., Vol. 44, No. 15, 2005

Table 1. Summary of High-Temperature Transport Data Fits for Various C16 Isomersa SES parameters (eq 15)

SEP parameters (eq 17)

species

A0

A1

B0

B1

B2

B3

n-hexadecane 2,7,12-trimethyltridecane 6,7,8-trimethyltridecane 4,8-dimethyl-6-n-propylundecane 5-n-propyl-5-n-butylnonane

1.82 1.99 1.86 1.89 1.93

0.846 0.878 0.896 0.907 0.922

0.018 0.016 0.018 0.019 0.021

-0.032 -0.012 -0.022 -0.018 -0.027

0.21 0.15 0.19 0.17 0.12

-1.4 × 10-2 -4.3 × 10-3 -0.042 -0.045 -0.018

9-octyldocosaneb 9-octylheptadecane

1.32 1.66

0.928 0.890

0.098 0.074

0.087 0.086

0.007

SES parameters {Ai} (eq 15) correspond to power law fits to the Stokes-Einstein scaling plots of Figures 3 and 5; parameters {Bi} (eq 17) are the best fit parameters of the Stokes-Einstein product vs density 980 K as depicted in Figure 7. b For SKS intermolecular potential. a

Figure 6. Stokes-Einstein product vs reduced density of hard prolate ellipsoids taken from Tang et al.10 Points correspond to MD simulations of systems of differing elongations: (O) φ ) 2, (4) φ ) 3, and (0) φ ) 5. Dashed and solid lines of corresponding filled symbols indicate the hydrodynamics prediction in the slip and stick limits, respectively.

Figure 7. Stokes-Einstein product vs density of selected C16 paraffins interpreted through prolate ellipsoid mapping. Transport data shown along T ) 980 K isotherm. (O) n-Hexadecane, φAV ) 2.82, (4) 6,7,8-trimethyltridecane, φAV ) 2.12, and (0) 5-n-propyl5-n-butylnonane, φAV ) 1.22. Dashed and solid lines with filled symbols have similar meanings to those of Figure 6.

(in stick and slip limits). A direct comparison can be made by forming the Stokes-Einstein product (SEP),

product for several species in terms of the prolate ellipsoid mapping defined by eqs 6-11. Along an isotherm, the average asphericity, φAV, is essentially independent of density. We find the same qualitative trends in the Stokes-Einstein product with increasing density and with increasing average elongation; nhexadecane, with φAV ) 2.82, shows a more rapid increase in the Stokes-Einstein product with density than the nearly spherically shaped 5-n-propyl-5-nbutylnonane (φAV ) 1.22), whose product remains relatively flat up to the highest densities considered here. This helps explain the varying degrees of curvature observed in Figures 3 and 5. For liquidlike states of nearly spherical molecules such as 5-n-propyl-5-nbutylnonane, the Stokes-Einstein product remains fairly constant at liquidlike densities. Thus, kT/D is linearly proportional to ησ, which explains the linear scaling of the 5-n-propyl-5-n-butylnonane transport data in Figure 5. For the more elongated n-hexadecane, the Stokes-Einstein product rises with increasing density in the liquid region, implying kT/D ≈ ησ/CSEP(F), where CSEP(F) is an increasing function of density in the liquid region. This implies that the local slope of the transport data plotted in Figure 3 decreases with increasing density. Interestingly, the temperature dependence of the Stokes-Einstein product appears to be fairly weak. Figure 8 shows the SEP of 6,7,8-trimethyltridecane for several isotherms. The data at lower temperatures are subject to large uncertainties, but the general trends

SEP ≡

Dησ kT

(16)

where, following the convention of Happel and Brenner,11 σ ) 2(ab2)1/3 ≈ 2(〈a〉〈b〉2)1/3 is the equivalent diameter of the prolate ellipsoid. Some of the results of Tang et al. are depicted in Figure 6, which shows the Stokes-Einstein product for hard ellipsoids as a function of density for aspect ratios of 2, 3, and 5. At low density, this quantity decreases with increasing density. At intermediate densities, a minimum is reached, and the product then increases with increasing density. The onset of this increase occurs at lower density, and its apparent magnitude appears more pronounced for ellipsoids of increasing elongation. By contrast, hydrodynamics (in either the stick or slip limit) predicts that the Stokes-Einstein product should be independent of density. This systematic variation of the Stokes-Einstein product with density and elongation is generally attributed to the increased degree of coupling between molecular orientation and shear, which, at high density, results in an increase of the Stokes-Einstein product above that predicted by hydrodynamics or kinetic theory. These observations are consistent with our hightemperature i-C16 transport data interpreted as effective prolate ellipsoids. Figure 7 depicts the Stokes-Einstein

Ind. Eng. Chem. Res., Vol. 44, No. 15, 2005 5833 Table 3. Comparison of Viscosity Estimated by Stokes-Einstein Approaches to Those Computed by NEMD12 for SKS Model of 9-Octyldocosane

suggest that the SEP may be approximated as a function of density only. This would imply an alternative approach to estimate the viscosity; if temperature dependence is assumed to be negligible, we can fit CSEP(F) to the results of a series of simulations at high temperature and a range of densities. The viscosity can then be estimated at more extreme states as

kT CSEP(F)Dσ

F (g/cm3)

ηSES (cP)

ηSEP (cP)

ηNEMD (cP)

372 311

0.761 0.800

1.80 6.82

1.68 5.90

1.75 5.84

give reasonable estimates of the Green-Kubo derived viscosity; the SES slope estimation is generally within about 10%, and the SEP method is within 5%. At the statepoint 453 K and 800 MPa, the quality of the predictions are comparable. In general, the SES approach slightly overestimates the Green-Kubo result, while the SEP approach slightly underestimates this quantity. The worst result appears to be the SEP prediction of n-hexadecane at the high-pressure point, which underpredicts the Green-Kubo viscosity by approximately 24%. This may be due to violation of the assumption that the SEP is truly independent of density, particularly at very high densities. This point will be expanded upon in the following section. Lube-Ranged Materials. In this section, we test the SES and SEP viscosity estimation approaches on lubricant-ranged materials. For these cases, direct computation of viscosity from equilibrium molecular dynamics is considered infeasible. We first consider the temperature-viscosity behavior of 9-octyldocosane, a C30 tribranched isoparaffin. The viscometric properties of this molecule have been studied with nonequilibrium molecular dynamics by Moore et al.12 In their work, the intermolecular potential model developed for n-paraffins by Smit, Karaborni, and Siepmann13 (SKS), and subsequently extended to branched paraffins by Mondello and Grest7 was employed. Using the SKS model, we have performed hightemperature simulations at selected points ranging between 0.7 and 1.025 g/cm3. Fit parameters for the SES slope and SEP are given in Table 1. We compare the predicted viscosity at 372 and 311 K to the results obtained by Moore et al. in Table 3. At 372 K, both estimation methods are in very good agreement with the NEMD viscosity result of 1.75 cP. At 311 K, the SES prediction exceeds the NEMD result by 17%, while the SEP prediction is still within 1%. The poorer estimate of the SES prediction at 311 K is presumed to be due in part to the extrapolation required of the power fit curve beyond the range of the high-temperature data from which it was parametrized. In Figure 9 the viscosity of 9-octylheptadecane is plotted as a function of pressure at 373 K. The TraPPE-

Figure 8. Stokes-Einstein product vs density for 6,7,8-trimethyltridecane at several temperatures. (b) T ) 980 K, (0) T ) 453 K, (4) T ) 393 K, (×) T ) 333 K. Dashed line is best fit of 980 K isotherm.

ηSEP(T,F) )

T (K)

(17)

We have fit the 980 K isotherm data for each of the i-C16 isomers to the equation CSEP(F) ) B0/F + B1 + B2F + B3F2 over the density range 0.1 e 1.0 g/cm3. Parameters for each species are given in Table 1. To summarize this section, we have proposed two alternative methods to estimate the viscosity from selfdiffusion and molecular geometry based on observations regarding the scaling of high-temperature transport data. The so-called Stokes-Einstein scaling (SES) employs a simple power law model to high-temperature transport data. The second approach employs the Stokes-Einstein product (SEP), under the assumption that this quantity is solely a function of density. The scaling relationships are determined solely from hightemperature simulations of transport properties at 980 K and various densities. We compare the predictive capabilities of these approaches relative to the bruteforce Green-Kubo calculations of viscosity for two state points in Table 2. At 333 K and 0.1 MPa, both methods

Table 2. Comparison of Viscosity Estimated by Stokes-Einstein Approaches to Those Computed by Green-Kubo Formula for Isomers of Hexadecane at Selected State Points species

T (K)

P (MPa)

n-hexadecane 2,7,12-trimethyltridecane 6,7,8-trimethyltridecane 4,8-dimethyl-6-n-propylundecane 5-n-propyl-5-n-butylnonane

333 333 333 333 333

0.1 0.1 0.1 0.1 0.1

n-hexadecane 2,7,12-trimethyltridecane 6,7,8-trimethyltridecane 4,8-dimethyl-6-n-propylundecane 5-n-propyl-5-n-butylnonane

453 453 453 453 453

800 800 800 800 800

F (g/cm3)

tηa (ns)

tDa (ns)

ηGKb (cP)

0.759 0.757 0.776 0.766 0.785

59 49 56 26 48

9 6 6 6 12

0.91 ( 0.08 0.92 ( 0.10 1.05 ( 0.23 0.88 ( 0.23 1.39 ( 0.02

0.949 0.949 0.963 0.956 0.973

227 291 305 122 126

15 15 6 15 45

8.95 ( 0.24 12.0 ( 1.3 14.9 ( 1.0 16.0 ( 1.9 85 ( 19.4

ηSES (cP) 1.02 1.02 1.20 0.98 1.52 8.6 14.3 19.2 18.6 100.0

ηSEP (cP) 0.92 0.94 1.09 0.91 1.35 6.8 11.4 15.3 14.8 71.6

a t and t are the durations over which D and η were computed, respectively. b Uncertainties in Green-Kubo estimates of viscosity η D computed from three sub-block averages.

5834

Ind. Eng. Chem. Res., Vol. 44, No. 15, 2005

Figure 9. Viscosity vs pressure for 9-octylheptadecane at 373 K. (b) Experiment,14 (0) SES estimation method, (×) SEP estimation method, (2) NEMD simulation,15 (4) NEMD simulation results plotted at P-F state points computed from model. See text for additional details.

UA model was employed in our calculations. In addition to our predictions, the figure also depicts experimental data14 for this compound and NEMD predictions using the SKS model from McCabe et al.15 Compared to experiment, the viscosity estimated by the SES and SEP is underpredicted. This observation is consistent with direct viscosity calculations with the TraPPE-UA model on similar types of isoparaffins. As pressure increases, the two estimation methods produce more divergent predictions, with the SEP showing a larger degree of underprediction. We speculate that this may be due to some temperature dependence in the Stokes-Einstein product, which increases with increasing density. Further simulations will be done to confirm this suspicion. A useful measure of the overall rise in viscosity with pressure is R*, the generalized pressure-viscosity coefficient,16

R* )

[

η

0 dp ∫0∞ η(p)

]

-1

(18)

where η0 is the ambient pressure viscosity. Application of eq 18 to the data in Figure 9 yields R/expt ) 8.7 / GPa-1, RSES ) 9.0 GPa-1, and R/SEP ) 8.1 GPa-1. Thus, while the estimation methods consistently underpredict the viscosity relative to experiment, the overall pressure viscosity response is well-predicted. We also plot in Figure 9 nonequilibrium molecular dynamics simulation results for 9-octylheptadecane at high pressure from McCabe et al.15 We note two important differences between the calculations performed in that work compared with our calculations. First, the SKS potential model is slightly different from the TraPPE-UA model employed in our work. There is also a mismatch in the definition of NEMD state points compared with our results. This comparison is discussed in more detail in the Appendix. While it would be more instructive to directly compare our results with the same potential model, the behavior predicted is qualitatively similar to the NEMD predictions. As a final example, SES predictions of viscosity for squalane using the TraPPE-UA model are depicted in Figure 10. Again, at both temperatures considered, the SES approach underestimates the viscosity compared with available experimental data, but the values of R*

Figure 10. Viscosity vs pressure for squalane at several temperatures. Filled symbols are experimental data17 and open symbols are from EMD + Stokes-Einstein calculations. (9) T ) 313 K, (b) T ) 373 K, (2) NEMD simulation at 372 K.12

agree favorably. For squalane, R/expt ) 17.5 and 12.1 GPa-1 at 313 and 373 K, respectively. The simulations yield 16.7 and 11.3 GPa-1. In addition, the NEMD result of Moore et al. at 373 K is consistent with our predictions. Conclusions In this paper we have examined the transport properties of a number of model isoparaffins over a wide range of state conditions within a Stokes-Einstein framework. The scaling behavior can be nonlinear and shows variation among different C16 isomers. These observations indicate a weakness of the assumptions inherent in low Reynolds number hydrodynamics in describing the drag force acting on a particle. Nonetheless, for a given molecule, the scaling behavior of transport property data appears to be wellbehaved, and we have proposed two approximate methods by which viscosity can be estimated from the selfdiffusion coefficient and geometric characterization of the molecule in terms of an inertially equivalent prolate ellipsoid. The first of these methods employs fitting high-temperature transport data to determine the overall Stokes-Einstein scaling (SES). This curve can then be used to back out the viscosity at extreme states as given by eq 15. The second approach estimates the Stokes-Einstein product (SEP) at high temperature over a range of densities and assumes this quantity is independent of temperature. Viscosity can then be estimated from eq 17. The presumed advantage of these estimation methods is that they require quantities that can be computed more easily from a molecular dynamics simulation than viscosity, particularly for systems exhibiting slow dynamics such as lube-range molecules at low temperatures and/or high pressures. These approaches appear to have respective merits and drawbacks. The SEP method appears to be very accurate at moderate densities, but predictions at very high density (pressure) appear to yield systematic underpredictions. It is believed that this is due to the fact that the Stokes-Einstein product is in fact temperature-dependent and becomes more so with increasing density. Further work will investigate whether this defect can be remedied. The SES approach is not as quantitative overall, but does not appear to suffer the same density-dependent underpredictions of the SEP.

Ind. Eng. Chem. Res., Vol. 44, No. 15, 2005 5835

The SES method appears a more suitable choice for making predictions about pressure-viscosity behavior. Acknowledgment Thanks to ExxonMobil Research & Engineering for support of this work and permission to publish the results. Appendix State Point Definition for Comparison of 9-Octylheptadecane Simulations. We noted in the main text that the state points simulated in the NEMD and EMD calculations of 9-octylheptadecane are derived from different sources. The state points chosen by McCabe et al. correspond to PFT conditions from experimental measurements on the fluid. By contrast, we have determined state points for transport property calculations by first performing simulations in the isothermal isobaric ensemble to estimate the density of the fluid corresponding to the pressure of interest. Transport properties are then computed in the canonical ensemble at the equilibrium density obtained. This procedure was adopted because if we are to employ simulations in a predictive and exploratory mode, we cannot assume that we will have access to experimental equation of state data to use in determining the PFT state of the simulation. In addition, the former case assumes no error in the equation of state predicted by the molecular model, whereas the latter approach uses the state point as predicted by the model. Such differences may be quite significant at high pressure. To demonstrate this, we have computed the pressure of liquid 9-octylheptadecane using the SKS potential model at 372 K and at the densities simulated by McCabe et al. This allows us to directly compare the density vs pressures predicted by the SKS model and TraPPE model employed in our work. The values are Table 4. Comparison of High-Pressure State Points Obtained for 9-Octylheptadecane at 372 K model

F (g/cm3)

PMOLEC (MPa)

PLRC (MPa)

PTOTAL (MPa)

SKSa SKS SKS TraPPEb TraPPE TraPPE

0.943 0.964 0.980 0.9767 0.9946 1.0096

612 749 865 780 913 1030

-76.0 -79.4 -82.1 -72.6 -75.4 -77.7

536 670 783 707 838 952

a Cutoff radius of interactions r b cut ) 9.85 Å. Cutoff radius of interactions rcut ) 10.0 Å.

shown in Table 4. Although the two potential models are slightly different, the resulting F vs P data are fairly consistent. In addition, the viscosities obtained from that model are replotted in Figure 9 as open triangles at these pressures. The resulting shift in the data to lower pressures yields much improved agreement with the experimental data, and the pressure-viscosity coefficient remains essentially unchanged. This demonstrates the extreme sensitivity at high pressure to the state point chosen; at 372 K and at a given pressure, the corresponding density predicted by the TraPPE-UA model is approximately 2% higher than the experimental density. At the highest pressures considered here, deviations in density of this magnitude can lead to a factor of 2 difference in viscosity predictions. Literature Cited (1) Gordon, P. A. Ind. Eng. Chem. Res. 2003, 42, 7025-7036. (2) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4058-4517. (3) Tang, S.; Evans, G. T. Mol. Phys. 1993, 80, 1443-1457. (4) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1-19; www. cs.sandia.gov/ sjplimp/lammps.html. (5) Gordon, P. A. Modeling Structure-Property Relationships in Synthetic Basestocks. In ACS Symposium Series 882. Dynamics and Friction in Submicrometer Confining Systems; Braihman, Y., Drake, J., Family, F., Klafter, J., Eds.; American Chemical Society: Washington, DC, 2004. (6) Gordon, P. A. Mol. Simul. 2003, 29, 479-489. (7) Mondello, M.; Grest, G. S. J. Chem. Phys. 1995, 103, 71567165. (8) McQuarrie, D. A. Statistical Mechanics; Harper Collins Publishers: New York, 1976. (9) Bereolos, P.; Talbot, J.; Allen, M.; Evans, G. J. Chem. Phys. 1993, 99, 6087-6097. (10) Tang, S.; Evans, G. J. Chem. Phys. 1995, 102, 3794-3811. (11) Happel, J.; Brenner, H. Low Reynolds Number Hydrodynamics; Prentice-Hall: Englewood Cliffs, NJ, 1965. (12) Moore, J.; Cui, S.; Cochran, H.; Cummings, P. J. Chem. Phys. 2000, 113, 8833-8840. (13) Smit, B.; Karaborni, S.; Siepmann, J. I. J. Chem. Phys. 1995, 102, 2126-2140. (14) ASME. Pressure-Viscosity Report; Technical Report; American Society of Mechanical Engineers: New York, 1953; two volumes. (15) McCabe, C.; Cui, S.; Cummings, P. T.; Gordon, P. A.; Saeger, R. B. J. Chem. Phys. 2001, 114, 1887-1891. (16) Bair, S. Trib. Trans. 2000, 1, 91-99. (17) Bair, S. Proc. Inst. Mech. Eng. J: J. Eng. Trib. 2002, 216, 139-149.

Received for review February 8, 2005 Revised manuscript received April 29, 2005 Accepted May 10, 2005 IE050156D