Impact of Molecular Architecture on the High ... - ACS Publications

7774

J. Phys. Chem. B 2000, 104, 7774-7783

Impact of Molecular Architecture on the High-Pressure Rheology of Hydrocarbon Fluids Loukas I. Kioupis and Edward J. Maginn* Department of Chemical Engineering, UniVersity of Notre Dame, Notre Dame, Indiana 46556 ReceiVed: March 14, 2000

Molecular dynamics simulations are conducted on three C18 poly-R-olefin isomers under extreme conditions typical of traction fluids or lubricants under elastohydrodynamic lubrication conditions. The viscosity, selfdiffusivity, and rotational relaxation times of the molecules are computed at pressures ranging from atmospheric to as high as 1 GPa. The dynamics of all three isomers are slowed as pressure increases, but a highly branched isomer shows a more dramatic reduction in mobility with pressure than does a linear or singly branched isomer. In particular, the viscosity of the highly branched molecule exhibits a much larger increase with pressure than does the viscosity of the other isomers, indicating that the highly branched molecule should exhibit more favorable traction properties than the other isomers. An explanation for the differences in dynamic properties between the isomers is given in terms of a reduction in liquid void volume coupled with the greater backbone stiffness of the highly branched molecule. A free volume analysis is conducted and shown to provide a better means of correlating the pressure dependence of diffusivity and viscosity than commonly used engineering models.

1. Introduction In a previous paper1 we examined the temperature dependence of the dynamics of poly-R-olefin (PAO) hydrocarbon molecules commonly found in synthetic motor oil base stock.2 The focus was to understand how molecular architecture impacts the viscosity number, an important lubricant performance measure that quantifies the susceptibility of a lubricant to viscosity reduction with increasing temperature. It was shown that the number and location of branch points in PAO isomers have a large effect on the temperature dependence of diffusivity and viscosity. This behavior was explained in terms of intramolecular mobility differences between the various isomers. In the present work, we examine two equally important issues. First, we wish to understand how dynamic quantities such as the viscosity change with pressure for different PAO isomers. Second, we want to compare the magnitude of the change in dynamics with pressure across the different isomers. In so doing, we hope to gain insight into what structural features can be “tuned” to achieve a desirable level of pressure dependence for the viscosity and other dynamic properties. An understanding of how lubricants perform under a wide range of pressures is important from a practical standpoint. Real lubricants operate frequently under conditions of high pressure, temperature, and shear rate. The pressures encountered in the elastohydrodynamic lubrication regime (EHL) often exceed 1 GPa,3 and the shear rates can be orders of magnitude higher (. 106 s-1) than those that can be achieved with modern rheometers.4 The dynamic behavior of lubricants under these extreme conditions is generally unknown and cannot be predicted easily. The viscosity of the fluid at these conditions is the key element in determining whether the lubricant can sustain a sufficient fluid film between moving parts. Most of the viscosity data available in the literature are at ambient conditions. The viscosity of fluids can change dramatically with pressure, however; hence, their performance as lubricants under * Corresponding author. Tel.: (219) 631-5687. Fax: (219) 631-8366. E-mail: [email protected]. http://www.nd.edu/∼ed.

operating conditions can be very different from that observed at atmospheric pressure. Fluids may also solidify at high pressures, which can result in very poor lubrication of the machine elements. To predict the film thickness with EHL theory, accurate rheological and equation of state models are needed for the pressure, temperature, and shear rate dependence of the viscosity. Ideally, these models should also be sensitive to differences in the structure, shape, size, and specific interactions of the lubricant molecules. The development of such models will require a more fundamental understanding of the link between molecular structure and the pressure dependence of lubricant properties. Another class of fluids for which pressure-viscosity characteristics are important are those used in automotive automatic transmissions. An example is the traction fluids used in traction drive devices, such as the continuously variable transmission (CVT).5 CVTs are being developed as fuel-efficient alternatives to conventional automatic transmissions. They transmit mechanical power via the shear resistance of a fluid film. The key property of these fluids is the traction coefficient, defined as µ ) F/N, where F is the tangential or shear force and N is the normal load applied on the film. The traction force F can be controlled in CVTs by changing the normal load N. Because high normal loads (corresponding to gigapascal pressures) are necessary for high-torque vehicles such as trucks, small regions of steel in the transmission can experience plastic deformation. This is due to the fact that conventional CVT fluids do not have large enough traction coefficients to achieve the required shear forces without the imposition of unreasonably large normal loads. A more complete fundamental understanding of what gives a fluid a large traction coefficient, while still maintaining other desirable features (such as lubrication performance over wide ranges of pressure, temperature, and sliding velocity, a sufficient viscosity to prevent roller-roller contact in the transmission, and fluidity at low temperatures) is clearly needed. The objective of the current work is to investigate the effect that molecular architecture has on the performance of fluids used in the above applications. We will examine the extent to which

10.1021/jp000966x CCC: $19.00 © 2000 American Chemical Society Published on Web 07/20/2000

High-Pressure Rheology of Hydrocarbon Fluids the degree of branching and backbone stiffness impact the dynamic properties of these fluids as a function of pressure using molecular dynamics (MD) simulations. There are two reasons for choosing MD to examine this problem. First, it is extraordinarily difficult to perform rheological experiments under the extreme conditions of pressure, temperature, and shear rate encountered in many applications. MD provides a viable and relatively simple alternative for studying these systems. Second, MD simulations can probe the behavior of liquids at the atomistic level, thus providing a fundamental understanding of how fluids “work” at these conditions. Such understanding cannot be obtained directly from experiments or from continuum mechanics models. This information can be useful from a formulation standpoint. The suitability of new molecules that have not been either tested or synthesized can easily be evaluated using MD. Such efforts can help direct and target the formulation of new materials having molecular structures with specific features. Numerous molecular simulation studies of hydrocarbon molecules at equilibrium or under shear conditions have been conducted. We refer the reader to two previous papers from our group1,6 for an extensive list of these works. Some previous MD studies that deal specifically with the issues addressed in the current work should be mentioned. Chynoweth et al.7 performed MD simulations of hexadecane at conditions typical of those found in EHL applications. The impact of pressure and temperature on the rheology of hexadecane at low and high shear rates was examined. Yamano et al.8 conducted an MD study of benzene and cyclohexane confined between two sliding surfaces. They computed the traction performance of these two molecules and found that cyclohexane exhibits a higher traction coefficient. The cyclohexane is able to transmit force to a greater degree than benzene through a mechanism in which the “chair” conformations of cyclohexane interlock to a greater degree than do the planar benzene molecules. The effect of molecular architecture on the frictional properties of hydrocarbon fluids has been studied by Tamura et al.9 The authors performed MD simulations on cyclohexane, n-hexane, and isohexane thin films confined between two solid Fe(001) surfaces. Finally, simulations that studied the conformational properties of realistic PAO molecules have been performed by Lahtela et al.10,11 Traction fluids have been studied experimentally as well. Of particular relevance to the current work are the studies of Hata and Tsubouchi,12,13 in which they measured the traction properties of several model compounds and developed structureproperty relations. Molecular features such as stiffness, polarity, size or shape, and hydrogen bonding were analyzed and discussed with respect to the traction properties of fluids. 2. Simulation Details Three C18 PAO isomers were simulated (see Figure 1): a highly-branched PAO molecule (4,5,6,7-tetraethyldecane, designated “HB”), a star PAO molecule (7-butyltetradecane, designated “S”), and a linear molecule (n-octadecane, designated “L”). As discussed in ref 1, these molecules represent three types of PAO fluids that can be produced from the trimerization of 1-hexene with different degrees of isomerization. Nonequilibrium molecular dynamics (NEMD) simulations were carried out to compute the viscosity and study the fluids under shear, while equilibrium molecular dynamics (EMD) simulations were used to calculate several other dynamic and static properties of the fluids at equilibrium conditions. The methods, equations of motion, and algorithmic details are

J. Phys. Chem. B, Vol. 104, No. 32, 2000 7775

Figure 1. Three different molecular structures simulated. Their threedimensional structures are also shown.

described in previous papers1,6 and have also been treated extensively in several excellent references.14-17 The simulations were performed in the isothermal-isobaric ensemble at 200 °C and at pressures ranging from 0.1 MPa up to 1 GPa. These conditions are typical of lubricants under EHL conditions, which are compressed within a very short period of time and subjected to a significant increase in temperature. Even at the highest pressures, all three isomers are liquids under these conditions.18 The united atom transferable potentials for phase equilibria (UA-TraPPE) force field20,21 was used to model the interactions of the hydrocarbon molecules. A fully-flexible representation of the TraPPE model was used by accounting for bond-stretching and bond-bending potentials using harmonic functions.1 We have used the UA-TraPPE model in past work and found it to be relatively good at describing the dynamics of linear and branched hydrocarbons. Although it tends to yield dynamics that are slightly faster than what is observed experimentally (i.e., viscosities are slightly too low, diffusivities are slightly too high), it does a remarkably good job of capturing trends as well as relatiVe differences between molecules. To test the accuracy of the UA-TraPPE model at high-pressure conditions, NEMD simulations of n-octane at 75 °C were performed and the results compared with experimental data.22 The computed Newtonian viscosities are shown in Figure 2. As expected, the computed viscosity is slightly lower than the experimental values, but the trend in viscosity with pressure is captured remarkably well. We therefore anticipate that the viscosities calculated for the PAO isomers may be slightly lower than what would be obtained experimentally, but that the trend with pressure will be captured faithfully. A Nose´-Hoover chain thermostat of length M ) 5, combined with a Nose´-Hoover barostat, was employed to maintain conditions of constant temperature and pressure.1,17 The thermostat controls the atomic kinetic energy of the system measured with respect to the streaming velocity. We note that Travis et al.23 observed differences in the non-Newtonian viscosity when an “atomic” or a “molecular” (i.e., translational) thermostat is used with the SLLOD equations of motion. A profile unbiased thermostat (PUT)24,25 is more appropriate to use in order to ensure the correct behavior at high shear rates. However, the limitations of conventional atomic thermostats are significant only at very high shear rates (i.e., for reduced shear rates γ˘ τ* g 2, where τ* ) (mσ/2)1/2). It has been shown that the linear (i.e., Newtonian) regime is insensitive to the thermostating method. For the relatively small shear rates used in

7776 J. Phys. Chem. B, Vol. 104, No. 32, 2000

Kioupis and Maginn 3. Results and Discussion 3.1. Shear Viscosity. NEMD simulations were performed at different shear rates, and pressures of 0.1, 200, and 400 MPa. The computed viscosities are shown in Figure 3 as a function of shear rate. The results are fit with the Carreau model28

η(γ˘ ) )

Figure 2. Comparison of simulation with experimental data22 for the viscosity of n-octane at 75 °C and at different pressures.

this study (γ˘ τ* < 1), it is generally acceptable to use atomic thermostats,26,27 and we have adopted this methodology for simplicity. As in ref 1, the reversible reference system propagator algorithm (rRESPA)17 was used to integrate the equations of motion. For the present study, a small time step of δt ) 1 fs and a large time step of ∆t ) 4 fs were used to integrate the fast and slow modes, respectively. The length of the EMD simulation runs varied from 1.5 ns to about 3 ns depending on the pressure. The NEMD lengths ranged from 0.4 ns at high shear rate and low pressure to 5 ns at low shear rate and high pressure. The EMD simulations for the L and S liquids were performed using 60 molecules in the periodic box. For the highly branched EMD simulations, 50 molecules were used, but the results were averaged over two different runs at each state point. All NEMD runs were performed using 50 molecules in each box, for all different molecules and state points.

η0 [1 + (λγ˘ )2]p

(1)

where the parameters η0, p, and λ have the following physical meaning. In the limit of zero shear rates (γ˘ f 0) the viscosity plateaus off to its Newtonian value η0. At high shear rates (γ˘ f ∞) the viscosity decreases (shear-thinning), and the viscosity versus shear rate curve approaches a power law η ∝ γ˘ -2p, where 2p is the slope of the log η vs log γ˘ curve. The parameter λ is a characteristic time constant, which approximately equals the inverse of the shear rate corresponding to the transition from the Newtonian to non-Newtonian regime. The parameters from the fit are listed in Table 1. A number of interesting qualitative observations can be made with regard to these results. The shear viscosity of the HB molecule is higher than the other two PAOs at all shear rates. This is consistent with calculations performed on these isomers at lower temperatures and a pressure of 2 MPa.1 All three isomers show significant shear-thinning, but the S molecule exhibits a somewhat greater resistance to shear-thinning than do either of the other two molecules. At low shear rates, the S molecule has a lower Newtonian viscosity than the L molecule. The situation is reversed at high shear rates, owing to excessive shear-thinning of the L molecule caused by its greater tendency to align with the flow field. This behavior was also observed in previous simulations of L and S molecules,1 as well as in the simulations of Evans and co-workers for linear and star tridecane molecules,29 and in experimental findings for star and linear unentangled polymers.30 The Newtonian viscosities η0 are shown as a function of pressure in Figure 4. The viscosity of all three fluids increases significantly with pressure, owing to fluid densification. Interestingly, the viscosity of the L and S molecules exhibits a similar dependence on pressure, while the HB fluid shows a much

Figure 3. Shear viscosity (η) versus shear rate (γ˘ ) from NEMD simulations at 200 °C and pressures of 0.1, 200, and 400 MPa. The lines are fits using the Carreau model (eq 1).

High-Pressure Rheology of Hydrocarbon Fluids


TABLE 1: Parameters of the Carreau Model (Eq 1) from Fitting the NEMD Results PAO molecule

P (MPa)

η0 (cP)

λ (ps)

p

highly branched (HB)

0.1 200 400 0.1 200 400 0.1 200 400

0.396 ( 0.032 2.347 ( 0.211 9.184 ( 1.010 0.289 ( 0.027 1.271 ( 0.083 2.819 ( 0.240 0.365 ( 0.022 1.468 ( 0.103 2.961 ( 0.296

21 ( 2 113 ( 10 540 ( 60 18 ( 2 44 ( 3 94 ( 8 35 ( 2 82 ( 6 163 ( 16

0.301 ( 0.024 0.244 ( 0.022 0.262 ( 0.029 0.233 ( 0.022 0.268 ( 0.017 0.258 ( 0.022 0.224 ( 0.013 0.273 ( 0.019 0.252 ( 0.025

star (S) linear (L)

Figure 4. Pressure dependence of the Newtonian shear viscosity η0.

greater increase in viscosity with pressure. The dependence of viscosity on pressure is commonly correlated using the empirical Barus equation3,5

η0 ) η1 exp(RPP)

(2)

where η1 is the viscosity of the fluid at atmospheric pressure and RP is an empirical parameter known as the pressureviscosity coefficient. This equation is known to be applicable only over a limited range of pressures and is unreliable when used to predict the viscosity at high pressure using values of η1 and RP determined at moderate pressures. The results in Figures 2 and 4 confirm this. Attempts at developing more sophisticated models that account for the pressure and temperature dependence of viscosity have been made.3,5,31-33 Nevertheless, the Barus equation remains the most commonly used model for estimating the pressure dependence of viscosity, mainly because it utilizes only a single parameter RP, which can be compared readily for different lubricants. Large pressure-viscosity coefficients are desirable for traction fluids, since this correlates well with the traction coefficient. For example, the viscosity of the HB molecules at 200 MPa is about the same as the viscosities of the L and S molecules at 400 MPa. Therefore, the shear force F ) Aη0γ˘ transmitted with the HB fluid at a normal load of 200 MPa should be about the same as the force transmitted via the L and S fluids at normal loads of 400 MPa, assuming that γ˘ , the contact area A, and the interactions of the fluids with the roller's surface are the same. This means that the HB fluid should have a higher traction coefficient µ ) F/N than either of the other two isomers. Using

a highly branched fluid, we can achieve better traction at smaller normal loads and avoid the problems associated with the application of high loads in CVTs (e.g., plastic deformation of steel surfaces). In a previous publication,1 we saw that one of the distinguishing features of the HB molecule is that it is considerably stiffer than the other two PAO isomers. The correlation between backbone stiffness and traction properties of fluids has been discussed in a review of molecular structures of traction fluids by Hata and Tsubouchi.12 In addition, Muraki34 has reported that the traction coefficient of hydrocarbon fluids increases with increased degree of branching. Using the current MD simulation results, a discussion on why stiffness and intramolecular mobility are so important in determining traction performance is given in section 3.4. 3.2. Self-Diffusivity. The Newtonian viscosity of fluids is strongly related to other dynamical properties, such as the selfdiffusivity and rotational relaxation rates. It is therefore useful to examine the way in which these quantities also change with pressure. Self-diffusion coefficients (Ds) computed as a function of pressure using EMD simulations are presented in Table 2. The self-diffusivity is a measure of the translational dynamics of the molecules and is obtained by calculating the long-time slope of the mean-square displacement of the center of mass of each chain and applying the Einstein equation for diffusion.14

d 1 Ds ) lim 〈[r(t) - r(0)]2〉 6 tf∞dt

(3)

The self-diffusivity is much easier to compute than is the viscosity. As a result, simulations over a broad range of pressures were conducted to more accurately establish the dependence of Ds on pressure. Figure 5 is a semilog plot of these results. The linear molecule (L) has the largest selfdiffusion coefficient at all pressures, whereas the highly branched (HB) molecule has the smallest diffusivity. The star (S) PAO falls somewhere in between. This trend was previously observed in simulations at conditions close to ambient pressure and explained1 in terms of intramolecular mobility differences arising from the degree of branching of the isomers. The same factor accounts for the differences observed at high pressure. From Figure 5, it can also be seen that the self-diffusivity decreases dramatically with pressure for all the molecules, as one expects, given that density increases with pressure. Interestingly, the dependence of Ds on pressure is especially strong for the HB molecule, where Ds drops by almost 2 orders of magnitude in moving from 0.1 MPa to 1 GPa. The linear and star molecules exhibit a less drastic dependence on pressure.35 This behavior is consistent with the results from the previous section, where it was shown that the viscosity of the HB molecule was much more dependent on pressure than the other isomers. As was also the case with viscosity, the pressure dependence of Ds cannot be described adequately with a simple exponential model of the form Ds ) A exp(-RPP).


Kioupis and Maginn

TABLE 2: Density (G), Square Radius of Gyration (R2g), Overall Torsion Angle Rotational Diffusivity (DO), Rotational Relaxation Time (τ1), and Self-Diffusion Coefficient (Ds) Values Computed from EMD Runsa PAO HB

S

L

a

P (MPa)

F (g/cm3)

R2g (Å2)

Dφ (10-1 rad2/ps)

τ1 (ps)

Ds (10-9 m2/s)

τrun (ns)

0.1 100 200 400 600 1000 0.1 100 200 400 600 1000 0.1 100 200 400 600 1000

0.6994(9) 0.7951(5) 0.8420(3) 0.9018(3) 0.9427(3) 1.0005(4) 0.6723(10) 0.7700(5) 0.8190(3) 0.8797(2) 0.9215(2) 0.9812(2) 0.6723(7) 0.7662(3) 0.8147(4) 0.8760(3) 0.9177(2) 0.9771(3)

10.21(1) 10.15(1) 10.15(1) 10.11(1) 10.11(2) 10.07(3) 18.53(1) 18.51(2) 18.49(2) 18.41(2) 18.36(3) 18.35(3) 27.98(5) 28.03(9) 27.92(6) 27.70(6) 27.72(11) 27.56(8)

0.162(7) 0.140(6) 0.132(5) 0.131(8) 0.113(11) 0.086(17) 2.08(6) 1.92(8) 1.79(8) 1.79(9) 1.54(9) 1.25(11) 2.85(8) 2.69(11) 2.67(12) 2.44(12) 2.35(14) 1.91(18)

18(1) 43(1) 77(2) 215(9) 676(47) 5538(1108) 21(1) 40(1) 52(2) 140(6) 169(8) 383(25) 39(1) 74(3) 123(6) 216(11) 350(21) 775(70)

2.22(9) 0.80(3) 0.39(2) 0.13(1) 0.044(5) 0.006(1) 2.94(9) 1.35(5) 0.70(3) 0.40(2) 0.20(1) 0.07(1) 3.80(11) 1.73(7) 1.11(5) 0.58(3) 0.33(2) 0.17(2)

1.5 2 2.5 3 3.15 3.15 1.5 1.5 2 3 3 3 1.5 1.5 2 3 3 3

Numbers in parentheses represent statistical uncertainties in the last reported digits.

Figure 5. Pressure dependence of the self-diffusion coefficient Ds.

Figure 6. Pressure dependence of the rotational relaxation time τ1.

3.3. Rotational Relaxation. The rotational relaxation rate of the longest axis of a molecule gives an indication of how fast the molecules “tumble” and “reorient”. A relaxation time τ1 associated with this process can be computed from the decay of the rotational autocorrelation function 〈e1(0)‚e1(t)〉, where e1 is the unit vector of the longest axis of the molecule. The longest axis of the molecule has the same direction as the principal axis of an ellipsoid inertially equivalent to the molecule. The unit vector e1 is found from the eigenvector corresponding to the smallest eigenvalue of the averaged tensor of inertia of each molecule.29 The relaxation time τ1 was computed via the following expression

molecules. That is, τ1 is expected to be largest for those molecules having the longest axis of rotation, as measured by the square radius of gyration, R2g. The basic picture is that “shorter” molecules are able to rotate faster whereas long molecules tumble sluggishly in a liquid. This picture is consistent with the low-pressure results shown in Figure 6. At 0.1 MPa, the L molecule has the largest R2g (see Table 2) and the largest value of τ1. The HB molecule, which has the smallest R2g, also has the smallest τ1 at low pressure. Exactly the same trend was observed in ref 1 at lower temperatures and pressures. Interestingly, this simple picture of τ1 scaling with R2g does not hold up at higher pressures. The rotational relaxation time of the HB molecule increases dramatically with pressure, exceeding τ1 of the L molecule at about 400 MPa. The magnitude of the increase in the rotational dynamics is similar to what was observed for the self-diffusivity. The τ1 observed at 1 GPa is 2 orders of magnitude larger for the HB, and about 1 order of magnitude larger for the L and S molecules, than their τ1 values at 0.1 MPa. The rotational relaxation dynamics often are related to the critical shear rate at which shear-thinning first occurs. Above

τ1 )

∫0∞ 〈e1(0)‚e1(t)〉 dt ) ∫0∞ e-(t/τ)

β

dt

(4)

where the autocorrelation function has been fitted to the Kohlraush-Williams-Watt (KWW) expression.36 The values of the computed relaxation times τ1 are listed in Table 2. The dependence of τ1 on pressure is shown in Figure 6. Differences in rotational relaxation times between molecules are usually explained in terms of the relative size of the



this critical shear rate, a rearrangement of the structural configuration of the fluid must take place to relieve stress. This frequently results in an increased intermolecular alignment of the molecules with the flow field and a resulting decrease in the apparent shear viscosity. Thus “slower” molecules are expected to shear thin at lower shear rates because they cannot viscous-dissipate fast enough the excess energy imposed due to shearing. The results in Figure 3 and values of λ in Table 1 indicate that shear-thinning occurs at progressively lower shear rates as pressure increases. This is particularly true for the HB molecule and is consistent with the observation that τ1 increases with pressure for all the isomers. The results in the three previous sections raise a number of interesting questions regarding the dependence of the dynamics of PAO molecules on pressure. Why are the dynamics of the HB molecule more dependent on pressure than the L and S molecules? How does the difference in the molecular architecture of the molecules affect their dynamics at high pressures? Is intramolecular mobility still as important in governing the dependence of dynamic properties on pressure as it was for temperature?1 To answer these questions we turn to an examination of intramolecular dynamics and free volume. 3.4. Intramolecular Flexibility. In a previous MD study,1 it was found that species having low intramolecular mobility had viscosities that exhibited a much greater dependence on temperature than did molecules that had a high degree of intramolecular mobility. To see if this can also be used to explain trends in viscosity with pressure, the intramolecular mobility of each isomer is quantified by computing a rotational diffusion coefficient Dφ, defined as1,37

1 d Dφ ) lim 〈[φ(t) - φ(0)]2〉 2 tf∞ dt

(5)

The rotational diffusion coefficient is a measure of the rate with which torsion angles change. Low values of Dφ indicate slower rotation around carbon-carbon bonds and increased stiffness in the structure. The dependence of Dφ on pressure was determined via EMD simulations. Results are listed in Table 2 and plotted in Figure 7. It can be seen that, at a given pressure, Dφ is much smaller for the HB molecule than for the other two isomers. This indicates that the HB molecule is much “stiffer” than the other two molecules. This is consistent with a previous study,1 where it was also found that Dφ for torsion angles close to a branch point are considerably smaller than those corresponding to torsion angles far from a branch point. Thus, the HB molecule exhibits low values of intramolecular flexibility owing to the high concentration of branch points and high steric hindrance in its structure. The star PAO exhibits smaller Dφ than the linear PAO, owing to the one branch point in its structure. Although there is a slight dependence of Dφ on pressure, these results show that intramolecular mobility is only a weak function of pressure. Thus relative changes in the intramolecular mobility with pressure alone cannot explain why the viscosity of the HB molecule exhibits a much greater dependence on pressure than the other two isomers. Another possible explanation for the fact that the viscosity of the HB molecule is more dependent on pressure than the viscosity of the L or S isomers is due to fluid density differences. Very small changes in density have a dramatic impact on dynamical properties of fluids. As shown in Table 2, the densities of these fluids change very much with pressure. The fluids are compressed by about 45% as pressure increases from 0.1 MPa to 1 GPa. If the density of the HB fluid increased more with pressure than the other fluids, then this

Figure 7. Pressure dependence of the overall torsion angle diffusivity Dφ.

would explain the trends in viscosity with pressure. The simulation results show, however, that the relatiVe change in the density of each fluid with pressure is roughly the same. Density differences by themselves do not explain the viscosity trends. While changes in intramolecular mobility and density with pressure cannot by themselVes explain why the viscosity of the HB molecule is more sensitive to pressure than the other isomers, we postulate that a combination of these two factors is responsible. This postulate is based upon the following physical picture. As the density increases with pressure, the voids in the liquid become quite small and few in number. The motion of the molecules in the liquid relies upon molecules either making “jumps” between these voids or displacing neighboring molecules into these void regions. This behavior is much like the motion of a penetrant in a polymer matrix. In either case, species that have an inherently low intramolecular flexibility (such as the HB molecule) will find it increasingly difficult to make these jumps or accommodate the movement of adjacent molecules as the fluid density increases. These molecules become “trapped” as the pressure (i.e., density) increases. On the other hand, species that have a high degree of backbone flexibility, such as the L and S molecules, will be able to more readily adopt conformations that favor molecular rearrangement. Even at high compressions, flexible molecules can still rearrange their shape more effectively and diffuse or “wiggle” through the dense medium easier than stiff molecules. This qualitative picture can explain why the dynamic properties of the “stiff” HB molecule have a greater pressure dependence than those of the L and S molecule. It can also serve to explain why we expect the HB molecule to act as a better traction fluid. In addition, the rigid arms of the HB molecule can interlock with other species under shear, enabling the molecule to act somewhat like a gear, transmitting force to adjacent molecules. These qualitative ideas are consistent with experimental observations that branched or stiff molecules are superior traction fluids. To provide a more quantitative basis for these observations, we turn to a free volume analysis. 3.5. Free Volume Analysis. In 1959 Cohen and Turnbull38 developed a model for viscous flow based on the concept of


Kioupis and Maginn

TABLE 3: Average Cell Volume W, Average Occupied van der Waals Volume in the Cell W0, Average Free Volume of the Cell Wf, and average van der Waals Volume of a Single Chain Wm, Computed from EMD Runs at Different Pressures P PAO molecule

P (MPa)

V (Å3/molecule)

V0 (Å3/molecule)

Vf ) V - V0 (Å3/molecule)

Vm (Å3/molecule)

highly branched (HB)

0.1 100 200 400 600 1000 0.1 100 200 400 600 1000 0.1 100 200 400 600 1000

603.1 530.5 500.6 467.7 447.4 421.6 627.3 547.7 514.9 479.6 457.8 429.9 627.3 550.4 517.6 481.6 459.7 431.6

317.4 317.9 317.7 316.4 315.0 314.4 324.7 323.6 323.3 322.2 321.4 320.2 325.0 323.8 323.6 323.0 321.4 320.3

285.7 212.6 182.9 151.3 132.4 107.2 302.7 224.2 191.6 157.3 136.4 109.8 302.4 226.7 194.0 158.6 138.3 111.3

318.8 318.5 318.3 318.0 317.8 317.2 325.2 325.1 324.9 324.6 324.4 324.1 325.4 325.3 325.1 324.9 324.7 324.3

star (S)

linear (L)

free volume of a liquid of hard spheres. The theory assumes that diffusion occurs when a void larger than some minimum value V* forms in the liquid owing to density fluctuations, and the molecule jumps into the void. On the basis of simple considerations of statistical redistribution of free volume, the following relation for the self-diffusivity can be derived

Ds ) gR*u exp(-γV*/Vf)

(6)

where g is a geometric factor, R* is approximately equal to the molecular diameter, u is the gas kinetic velocity u ) (3kBT/ m)1/2, γ is a constant introduced to correct for possible overlap of the free volume (0.5 < γ < 1), V* is the minimum void volume needed for a jump, and Vf is the free volume. Using the Stokes-Einstein relation

Ds )

kBT 6πRη0

(7)

(R is the radius of the spherical molecule, and η0 the viscosity), an expression for the viscosity can be obtained

η0 ) A′T1/2 exp(γV*/Vf)

(8)

where A′ is a constant. For conditions of constant temperature, the above relation is equivalent to the empirical equation derived earlier by Doolittle39,40

η0 ) A exp(BV0/Vf)

(9)

where V0 represents the van der Waals or occupied volume of the molecule, and Vf is the free volume in the fluid. In this case, we will assume that V* ) BV0 ) γV* approximately represents the void volume barrier that a molecule needs to overcome in order to take an appreciable diffusive step. For simple liquids, the constant B is of the order unity. This means that the void volume that needs to open up in a fluid of simple molecules is approximately equal to the molecular volume. The key point in applying the above models is the definition of free volume Vf. Unfortunately, Vf has never been defined in an exact way. Initially, Vf was considered to be the average fluid volume V per molecule minus the van der Waals volume V0 (Vf ) V - V0). In practice, Vf most often represents the difference between the liquid specific volume and either the volume of a glasslike condensed phase30 or the volume of the liquid

extrapolated to zero temperature. Alternatively, V0 can be treated as a third adjustable parameter.41 The advantage of a molecular simulation study is that these volume terms may be computed directly and the predictions of the theory rigorously tested. For this study, we consider Vf ) V - V0, where V is the volume of the simulation box per molecule (the liquid molar volume), and V0 is the occupied volume of the molecules in the box. The occupied volume is calculated analytically using the method of Dodd and Theodorou42 as the sum of the van der Waals volumes of all atoms minus the overlapping volume. The Lennard-Jones parameter σii is used as the atomic diameter of each atom i. This is by no means the only legitimate way to define the occupied volume. Unfortunately, there are many ways to define a molecular volume in addition to the van der Waals volume, including the Voronoi volume, the excluded volume, and the volume enclosed by the surface defined by rolling a spherical probe on the van der Waals surface of the molecules. However, we believe that the differences between the various representations of volume should have little impact on the relative changes of the volume with pressure, which is what we seek here. For an interesting discussion on the volume of macromolecules, we refer the reader to ref 43. The average values of V, V0, and Vf computed directly from the simulations are shown in Table 3. Also shown is the average molecular volume of each chain Vm. The molecular volume, Vm, was also computed using the analytic procedure of ref 42. The difference is that in the computation of V0 both intramolecular and intermolecular overlapping of the spheres were accounted for, whereas in the calculation of Vm each chain was treated separately; overlapping between united atom spheres that belong to different chains are not accounted for. From a comparison of the values in Table 3, we see that the free volume decreases significantly with increasing pressure. This is driven not by changes in V0 or Vm, which are essentially constant with pressure. Rather, the free volume is reduced by the elimination of void volume between molecules (i.e., a reduction in the liquid volume). The computed viscosity and self-diffusivity results are plotted against the ratio V0/Vf in Figures 8 and 9 to test the free volume model. The data were fit using equations of the form of eq 9, with A and B as adjustable parameters. For the self-diffusivity (Figure 9), the fit is good for pressures up to 600 MPa for the S and L molecules, and up to 400 MPa for the HB fluid. Therefore, the straight lines shown in Figure 9 were only fit to the results in these pressure ranges. The straight lines shown,


Figure 8. Viscosity η0 vs occupied to free volume ratio V0/Vf. The straight lines are fits using the equation η0 ) A exp(BV0/Vf).

Figure 9. Self-diffusion coefficient Ds vs occupied to free volume ratio V0/Vf. The straight lines are fits using the equation Ds ) A exp(-BV0/Vf).

however, were extrapolated to 1 GPa pressures for comparison purposes. Despite the slight deviation at high pressure and the relative simplicity of the model, the viscosity and diffusivity appear to be captured quite well by the model. This suggests that free volume models such as that in eq 9 are better able to capture the pressure dependence of viscosity and diffusivity than simple Barus-type models. The free volume model also gives physical insight into the transport mechanisms in these fluids. The parameters B from Figure 8 for the L, S, and HB molecules are 2.19, 2.34, and 3.17, respectively. From Figure 9 the slope B for L, S, and HB is -1.93, -2.02, and -2.91, respectively. The fact that B > 1 indicates that voids having a volume larger than the volume of molecules are needed for diffusion and that the behavior of the fluids deviates from that of simple liquids. The parameter B

J. Phys. Chem. B, Vol. 104, No. 32, 2000 7781 also gives a measure of the relative dependence of the dynamics on changes in the fraction of occupied to free volume V0/Vf of the fluids. The fact that BHB > BS > BL indicates that the highly branched molecule is more sensitive to changes in free volume and therefore to the density (since V0 = const) than S, which in turn is only slightly more sensitive to free volume changes than L. This explains the stronger dependence of the HB dynamics on pressure. Despite the fact that the HB molecule has the smallest R2g, it requires the largest void volume V* ) BV0 for a jump. This is directly related to its backbone stiffness and shape. The HB molecule can only move or rotate as a whole owing to its low flexibility. Thus, in order to take an appreciable diffusive step, it has to sweep out a larger amount of neighboring molecules. In other words, a larger hole must open up in order for the HB molecule to jump. Since large void volumes occur infrequently, particularly at high pressures, the HB molecule has the smallest diffusivity and slowest rotation rate. Consequently, its dynamic properties are impacted the most by reductions in free volume. The other two PAO fluids can more easily diffuse and penetrate smaller voids because their greater backbone flexibility enables them to change configurations and readjust their shape faster. The self-diffusivity results in Figure 9 show a deviation from the free volume model at the highest pressures (close to 1 GPa). There could be several reasons for this. First, the simulation results are most uncertain at the highest pressures because of the slow dynamics. Thus the deviations could just be the result of inaccuracies in the computed diffusivity. More likely is that the model itself breaks down at these extreme pressures. The free volume equations described here were developed for hard sphere liquids. The inadequacy of this model for large and complex molecules is likely to become more evident at high pressures, where a significant decrease in free volume and the increased intermolecular overlapping makes other mechanisms of diffusion more significant. In support of this, we simulated the self-diffusivity for smaller alkanes (i.e., octane isomers) at pressures up to 1 GPa. It was found that eq 9 fit the results over the entire pressure range, without the high-pressure deviations seen with the C18 isomers. This suggests that more sophisticated models that account for the shape of the molecules may be more useful for describing the way diffusivity and viscosity vary with pressure for complex and large molecules. One such model that appears to capture some of the qualitative aspects of the results presented here has been proposed by Mauritz, Storey, and co-workers.44-51 This model is based on free volume theory but assumes that large molecules diffuse with small jumps along the three principal axes of inertia. Depending on the shape and size of the molecule, diffusion along each axis of inertia of the molecule can occur at different rates. This is because different cross sections of the same molecule have different shapes and areas and require different sized holes into which they can jump. For example, n-paraffin molecules, such as the L molecule, will diffuse by motion along their chain contour, owing to their smaller cross section in this direction. Similarly, molecules that can be approximated as prolate ellipsoids will diffuse faster along their longest axis. On the other hand, the more globular HB molecule should require about the same size holes in all different directions. This characteristic allows linear molecules to diffuse more effectively under high compressions (i.e., high pressures) than branched molecules, and therefore exhibit less dependence on pressure. The theory in refs 44-51, however, assumes that molecules are rigid in their minimum configuration and does not account for differences in flexibility between different molecules and


Kioupis and Maginn an equivalent ellipsoid directly from the simulations. Thus, if eq 10 is valid, we can compute the viscosity from self-diffusivity values that are usually obtained with less computational expense. 4. Conclusions

Figure 10. Self-diffusion coefficient Ds vs 1/Rη0 (see text). The straight lines are fits using the Stokes-Einstein relation.

their ability to change shapes. Molecular simulations should provide effective means to test these and related theories. Finally, the validity of the Stokes-Einstein relation (eq 7) is tested in Figure 10 for all the systems studied and for pressures up to 400 MPa. The values of the self-diffusivity Ds are plotted against the quantity 1/Rη0, where η0 is the Newtonian viscosity and R is the radius of a sphere of equal volume to the molecule. The radius R is computed from the expression Vm ) 4/3πR3 using the values of the molecular volume Vm from Table 3. It is seen that the self-diffusion coefficient scales linearly with 1/Rη0 as predicted by the Stokes-Einstein relation. A plot of Ds vs 1/η0 (not shown) also reveals that the self-diffusion coefficient is inversely proportional to viscosity. This is because the equivalent radius R is almost independent of pressure (see values of Vm from Table 3). The Stokes-Einstein relation was used in Cohen-Turnbull free volume theory to correlate self-diffusivity with viscosity. It is shown here that this simple relation can be used as a qualitative model for the longer C18 hydrocarbons as well. However, the slopes in Figure 10 are different than the value of kBT/6π ) 0.346 cP Å3/ps (for T ) 473 K), and the model cannot be used for quantitative predictions from firstprinciples. This indicates that the diffusion of the C18 hydrocarbons cannot be simply approximated to that of a sphere immersed in the fluid. This is particularly evident for the linear C18 molecule. For assymmetric molecules a simple correction can be used on the basis of Perrin’s extension of Stokes equation to ellipsoids of revolution.52

Ds )

kBT 6πFRη0

(10)

In this modified expression the molecules can be approximated as ellipsoids. The coefficient F is a shape factor that characterizes the deviation of molecular shape from sphericity. From the value of the slope in Figure 10 we find that F ) 1.73 for the L molecule. Perrin has derived analytical expressions for F as a function of the semiaxes of prolate and oblate ellipsoids of random orientation.52 For a prolate ellipsoid the value of F ) 1.73 corresponds to an axial ratio a/b = 13. It would be interesting to explore this further and compute the semiaxes of

The viscosity and dynamic behaviors of three C18 poly-Rolefin isomers were examined at conditions of high pressure and temperature, typical of those found in many elastohydrodynamic lubrication applications. The viscosity of the fluids was computed as a function of shear rate using the NPT-Sllod algorithm. In addition, EMD simulations in the NPT ensemble were performed to compute several dynamic and static properties of the fluids under equilibrium conditions. It was shown that the viscosity of the fluids increases dramatically with increasing pressure. Similarly, the selfdiffusivity decreases and the rotational relaxation time increases by an order of magnitude or more in moving from 0.1 MPa to 1 GPa pressure. Pressure slows the dynamics of the molecules significantly, owing to a large compression of the fluid. The HB molecule had the highest viscosity and the lowest diffusivity of all three isomers. Furthermore, these quantities exhibited a much stronger dependence on pressure for the HB molecule than for the L or S molecules. This behavior was explained in terms of a free volume model in which the stiffer branched molecule is less efficient at moving in a fluid with reduced free volume (as occurs under extreme pressure) than the more mobile S and L molecules. The results suggest that stiff, branched molecules such as the HB isomer would exhibit superior traction properties. The free volume models of Cohen-Turnbull and Doolittle were found to correlate the viscosity and diffusivity results much better than the commonly used Barus equation. This suggests a viable strategy for estimating the viscosity and diffusivity of complex molecules under extreme conditions. The free volume of a fluid can be computed from a short simulation at the conditions of interest. Easy to obtain ambient results (either experimental or simulation) can be extrapolated to the extreme conditions using the model and computed free volumes. More sophisticated models can be developed that are expected to more accurately capture the pressure dependence of viscosity and diffusivity. Acknowledgment. E.J.M. thanks the National Science Foundation for a CAREER award (CTS-9701470). We thank the U.S. Army Research Office (Grant DAAG55-98-1-0091) for the computers on which many of these calculations were performed, and the Mobil Foundation for partial financial support. Dr. Michael Greenfield of Ford Motor Company is acknowledged for several helpful discussions. References and Notes (1) Kioupis, L. I.; Maginn, E. J. J. Phys. Chem. B 1999, 103, 10781. (2) Wu, M. M. High Viscosity Index Synthetic Lubricant Compositions. U.S. Patent 4,827,064, May 2, 1989. (3) Hamrock, B. J. Fundamentals of Fluid Film Lubrication; McGraw Hill: New York, 1994. (4) Spearot, J. A., Ed. High-Temperature, High-Shear Oil Viscosity. Measurement and Relationship to Engine Operation; ASTM (STP 1068): Philadelphia, PA, 1989. (5) Stachowiak, G. W.; Batchelor, A. W. Engineering Tribology; Elsevier: Amsterdam, 1993. (6) Kioupis, L. I.; Maginn, E. J. Chem. Eng. J. 1999, 74, 129. (7) Chynoweth, S.; Coy, R. C.; Michopoulos, Y. J. Eng. Tribol. 1995, 209, 243. (8) Yamano, H.; Shiota, K.; Miura, R.; Katagiri, M.; Kubo, M.; Stirling, A.; Broclawik, E.; Miyamoto, A.; Tsubouchi, T. Thin Solid Films 1996, 281-282, 598.

High-Pressure Rheology of Hydrocarbon Fluids (9) Tamura, H.; Yoshida, M.; Kusakabe, K.; Young-Mo, C.; Miura, R.; Kubo, M.; Teraishi, K.; Chatterjee, A.; Miyamoto, A. Langmuir 1999, 15, 7816. (10) Lahtela, M.; Pakkanen, T. A.; Nissfolk, F. J. Phys. Chem. 1995, 99, 10267. (11) Lahtela, M.; Pakkanen, T. A.; Nissfolk, F. J. Phys. Chem. A 1997, 101, 5949. (12) Hata, H.; Tsubouchi, T. Tribol. Letters 1998, 5, 69. (13) Tsubouchi, T.; Hata, H. Jpn. J. Tribol. 1996, 41, 463. (14) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (15) Evans, D. J.; Morriss, G. P. Statistical Mechanics of Nonequilibrium Liquids; Academic Press: London, 1990. (16) Cummings, P. T.; Evans, D. J. Ind. Eng. Chem. Res. 1992, 31, 1237. (17) Martyna G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (18) Although we were unable to locate experimental solidification pressures for the PAO isomers studied here, we are confident that they are in the liquid state for all the conditions examined. We base this on the fact that 9-n-octylheptadecane goes through a solidification at about 0.5 GPa at 37 °C, but this transition moves to greater than 1 GPa at 100 °C.19 Likewise, n-hexadecane has solid transitions at about 88 MPa at 37 °C, and 0.46 GPa at 100 °C.19 Given that the simulation temperature is 200 °C, it is reasonable to expect that these C18 isomers will be liquids even at 1 GPa pressure. (19) Dr. Peter Gordon, Mobil Technology Company, private communication. (20) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569. (21) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508. (22) Harris, K. R.; Malhotra, R.; Woolf, L. A. J. Chem. Eng. Data 1997, 42, 1254. (23) Travis, K. P.; Daivis, P. J.; Evans, D. J. J. Chem. Phys. 1995, 103, 1109. (24) Travis, K. P.; Daivis, P. J.; Evans, D. J. J. Chem. Phys. 1995, 103, 10638. (25) Bagchi, K; Balasubramanian, S.; Mundy, C. S.; Klein, M. L. J. Chem. Phys. 1996, 105, 11183. (26) Cui, S. T.; Cummings, P. T.; Cochran, H. D. J. Chem. Phys. 1996, 104, 255. (27) Balasubramanian, S.; Mundy, C. S.; Klein, M. L. J. Chem. Phys. 1996, 105, 11190.

J. Phys. Chem. B, Vol. 104, No. 32, 2000 7783 (28) Hieber, C. A.; Chiang, H. H. Polym. Eng. Sci. 1992, 32, 931. (29) Daivis, P. J.; Evans, D. J.; Morriss, G. P. J. Chem. Phys. 1992, 97, 616. (30) Ferry, J. D. Viscoelastic Properties of Polymers; John Wiley & Sons Inc.: New York, 1970. (31) Sorab, J.; Vanarsdale, W. E. Tribol. Trans. 1991, 34, 604. (32) Yasutomi, S.; Bair, S.; Winer, W. O. J. Tribol. ASME Trans. 1984, 106, 291. (33) Yasutomi, S.; Bair, S.; Winer, W. O. J. Tribol. ASME Trans. 1984, 106, 304. (34) Muraki, M. Tribol. Intern. 1987, 20, 347. (35) The dynamics of these molecules at 1 GPa are quite slow, such that true diffusive behavior was not observed over the time scale of the simulations, despite the fact that the mean-square displacement was linear in time. We therefore regard the diffusivities at 1 GPa with some caution. (36) Smith, G. D.; Yoon, D. Y. J. Chem. Phys. 1994, 100, 649. (37) Runnebaum, R.; Maginn, E. J. J. Phys. Chem. B 1997, 101, 6394. (38) Cohen, M. H.; Turnbull, D. J. Chem. Phys. 1959, 31, 1164. (39) Doolittle, A. K. J. Appl. Phys. 1951, 22, 1471. (40) Doolittle, A. K.; Doolittle, D. B. J. Appl. Phys. 1957, 28, 901. (41) Hogenboom, D.; Webb, W.; Dixon, J. A. J. Chem. Phys. 1967, 46, 2586. (42) Dodd, L. R.; Theodorou, D. N. Mol. Phys. 1991, 72, 1313. (43) Paci, E.; Velikson, B. Biopolymers 1997, 41, 785. (44) Mauritz, K. A.; Storey, R. F.; George, S. E. Macromolecules 1990, 23, 441. (45) Mauritz, K. A.; Storey, R. F. Macromolecules 1990, 23, 2033. (46) Coughlin, C. S.; Mauritz, K. A.; Storey, R. F. Macromolecules 1990, 23, 3187. (47) Coughlin, C. S.; Mauritz, K. A.; Storey, R. F. Macromolecules 1991, 24, 1526. (48) Coughlin, C. S.; Mauritz, K. A.; Storey, R. F. Macromolecules 1991, 24, 2113. (49) Storey, R. F.; Mauritz, K. A.; Carter, M. L.; Wang, D. Polym. Prepr. (Am. Chem. Soc., DiV. Polym. Chem.) 1991, 32, 130. (50) Wang, D.; Mauritz, K. A.; Storey, R. F. Polym. Prepr. (Am. Chem. Soc., DiV. Polym. Chem.) 1991, 32, 269. (51) Wang, D.; Mauritz, K. A. J. Am. Chem. Soc. 1992, 114, 6785. (52) Tanford, C. Physical Chemistry of Macromolecules; John Wiley & Sons: New York, 1961; p 327.

Impact of Molecular Architecture on the High ... - ACS Publications

Recommend Documents