Effects of a typical shear dependent viscosity on analytical elastohydrodynamic lubrication film thickness predictions: A critical issue for the classical approach

Abstract

Many engineering estimates of the film thickness in the concentrated contacts of real machines have come from extrapolations of measurements of elastohydrodynamic lubrication film thickness performed in glass on steel elastohydrodynamic rigs. Such estimates are likely to have large errors due to shear dependence of viscosity. The classical film thickness formulas employed have not been validated except for some Newtonian reference liquids at room temperature because the real pressure-viscosity response measured in viscometers has been ignored. A blend of polyalpha olefin base oils with mild shear-thinning has been employed in a line contact calculation to assess the effects of shear-dependent viscosity on the power-law exponents of the classical film thickness formula. The exponents on the pressure-viscosity coefficient and on the elastic modulus of the solids are not sensitive to the non-Newtonian effect. The exponents on ambient pressure viscosity and velocity are slightly reduced by shear-thinning. The exponents on pressure and scale are substantially increased by the shear dependence. The usual practice of measuring film thickness in an elastohydrodynamic rig to obtain an effective “ $α$ -value” to use in a classical Newtonian film thickness formula will overstated the film thickness of a liquid which is shear dependent in the elastohydrodynamic lubrication inlet.

Keywords

Classical elastohydrodynamic lubrication quantitative elastohydrodynamic lubrication film thickness shear-thinning rheology

Introduction

The pressure and shear dependence of lubricating oils have been measured at high pressure for > 50 years.¹ Despite these available data, the classical approach to elastohydrodynamic lubrication (EHL) has not used the real viscosity from high-pressure viscometers. Instead, fictional accounts of work by Barus, Roelands, and Eyring have been employed to justify analysis with viscosities which cannot be measured in viscometers.² Barus³ measured the effective viscosity of a solid and his equation was linear, not exponential. Roelands⁴ ignored data that did not fit his correlation and it is, therefore, not used by laboratories that make accurate measurements. Eyring and co-authors⁵ offered a theory to explain shear dependence observed in early viscometers, a shear dependence later found to be caused by viscous heating. The validity of the assumption of a Newtonian inlet zone has been obscured by the practice of adjusting the “ $α$ -value” to make the Newtonian formula agree with a film measurement.

Many engineering estimates of the film thickness in the concentrated contacts of real machines have come from extrapolations of measurements of EHL film thickness performed in glass on steel elastohydrodynamic (EHD) rigs. There are standards which recommend this practice.^6,7 The formula⁸ often chosen is, for circular contacts,

h_{c} = 1.55 α^{0.53} (μ_{0} \bar{u})^{0.67} E^{' 0.061} R^{0.33} p_{H}^{- 0.201}

(1)

The traditional dimensionless groups are useful for idealized liquids which have a weak compressibility independent of temperature and are Newtonian to the highest shear stress within the EHL inlet zone. The Hamrock & Dowson⁸ formula (H&D) above is written here in dimensional form in terms of pressure rather than load. Contacts of differing scale,

R = 1 / (1 / R 1 + 1 / R_{2})

where R₁ and R₂ are the local radii of curvature of the contacting solids at the contact point, should be compared at the same pressure, not load. For example, an increase in scale at constant load will greatly reduce the pressure and may take the contact out of the full EHL regime and into the isoviscous rigid regime. This formula (1) assumes a Newtonian inlet zone. Extrapolations from EHD rigs to real machines, employing effective “

α

-values,” are being used without investigation of the implications of shear dependence in either situation, the EHD rig or the machine. Here, the effects of shear dependence on the exponents in a typical power-law film thickness formula are investigated through isothermal non-Newtonian EHL simulations.

Early inlet zone calculations⁹ indicated that the exponent on R in H&D (0.33) is only a lower bound for Newtonian response; the exponent for real lubricants may be much larger. This has been experimentally validated.¹⁰ It was found that $(\partial \ln h_{c} / \partial \ln R)$ may be as large as 0.64.¹⁰ It was also found experimentally¹¹ that the exponent on $p_{H}$ can be as large in magnitude as −0.49.

These facts indicate that many engineering estimates of the film thickness in the concentrated contacts of real machines have significant errors. It is proposed that rolling line contact simulations be used to assess the effect of mild shear dependence on each of the H&D parameters for one example of the shear dependence of a blend of different viscosity base oils. The results of this work have more significance because of the recent publication¹² of an experimental investigation of the effect of blending base oils on micropitting in sliding contacts where the film thickness is critical for interpretation but cannot be directly measured. In that work, the film thickness was estimated with a classical formula from measurements in a different contact, although it was understood that the Newtonian assumption is not valid.¹²

Rheology of the example blend

The shear dependence that affects EHL film thickness occurs at relatively low shear stress which can be applied in some types of high-pressure viscometers. Ambient pressure viscometers are not useful because of shear cavitation.¹³ Shear thinning at inlet stress can be observed in polymer-thickened oils, in high molecular weight base oils,¹⁴ and in base oil blends of very different molecular weights.¹³ Here, a blend of two PAO (polyalpha olefin) base oils is used as a model lubricant. It is a PAO-4 with 25% by weight of PAO-600.¹⁵ An interesting property of base oil blends is that for a small concentration of the higher viscosity component, the Newtonian limit is lower than either of the components. See the blend of PAO-4 and PAO-40 by Bair and Qureshi.¹³

The pressure dependence of a blend of PAO-4 with 25% by weight of PAO-600 at 70 °C is plotted in Figure 1. McEwen¹⁶ proposed the following model which is widely employed by laboratories with high-pressure viscometers:

μ = μ_{0} {(1 + \frac{α_{0}}{q} p)}^{q} = μ_{0} {(1 + \frac{α *}{q - 1} p)}^{q}, q > 1

(2)

This equation is accurate for the PAO blend at 70 °C to a pressure of 500 MPa using

μ_{0} =

24 mPa·s,

α_{0} = 13.9 GP a^{- 1}

α * = 12.0 GP a^{- 1}

, and

q = 7.3

Figure 1.

The pressure dependence of a blend of polyalpha olefin (PAO)-4 with 25% by weight of PAO-600 at 70 °C.

The shear dependence of a blend of PAO-4 with 25% by weight of PAO-600 is plotted in Figure 2. The shear dependence can be modeled with a modified Carreau⁹ form

η = μ {[1 + {(\frac{τ}{G})}^{2}]}^{\frac{1 - \frac{1}{n}}{2}}

(3)

Figure 2.

The shear-dependent viscosity of a blend of polyalpha olefin (PAO)-4 with 25% by weight of PAO-600.

G = 30 kPa and n = 0.82. This is not a severe case of shear thinning. For a polymer-thickened oil, G will be smaller. For a high molecular weight base oil, G will be larger and n will be much smaller. There is likely to be a second transition in the shear dependence due to the shear-thinning of the PAO-4 at shear stress greater than 10 MPa; however, this would not affect the film thickness since it is greater than the stress in the inlet.

The equation of state that is frequently used in EHL simulations, the Dowson and Higginson version of the Hayward equation with fixed compressibility,⁸ will be employed. While this is only applicable for mineral oil at room temperature and low pressure, it is retained here for consistency.

Pan and Hamrock¹⁷ found for the compressible central film thickness of a line contact,

h_{c} = 2.154 α^{0.47} (μ_{0} \bar{u})^{0.692} E^{^{'} 0.110} R^{0.308} p_{H}^{- 0.332}

(4)

Considering the mild shear dependence above, for each of the six parameters on the right-hand side of equation (4),

X_{i}

, the effective exponent,

Q_{i} = \partial \ln h_{c} / \partial l n_{i}

and how it depends on the other five parameters can be calculated.

Numerical simulations

The numerical model employed in this work is based on the full-system finite element approach.¹⁸ Steady-state smooth isothermal line contacts are considered. The generalized Reynolds equation of Yang and Wen¹⁹ is used in resolving the hydrodynamic pressure field that forms between the solid surfaces. Compared to the standard Reynolds equation, it allows for incorporating viscosity variations with shear across the film thickness through cross-film integral terms. The shear dependence is described by equation (3). The hydrodynamic pressure is applied as a normal boundary load to the surface of an equivalent solid domain with reduced material properties, such that it would accommodate the total elastic deformation of the two solids. This deformation is obtained by applying the linear elasticity (i.e. Navier) equations to a square solid domain of sufficiently large dimensions, such that a half-space approximation would be attained. Given the infinite nature of the contact along the roller principal axis, a plane strain configuration is adopted. A load balance equation is used to ensure that the correct external load is applied to the contact, by monitoring the value of the rigid body separation. This equation balances the hydrodynamic pressure generated within the lubricant film with the external load applied to the contact. Finally, a shear stress equation is used to determine shear stress variations within the lubricating film. These are used for evaluating the shear response of the lubricant. All equations are written and solved in dimensionless form, with the pressure, elastic deformation, and shear stress fields being discretized using second-order non-regular, and non-structured Lagrange finite elements. Triangular elements are used for discretizing the two-dimensional (2D) solid domain, while line elements are used for discretizing the one-dimensional (1D) hydrodynamic domain. Projections of the 2D triangular elements of the solid domain over the hydrodynamic domain are used as the mesh of the latter to avoid any unnecessary (and computationally expensive) mapping operations between the two meshes. A sufficiently fine mesh is employed to guarantee grid-independent solutions (i.e. discretization errors are reduced to machine precision), with mesh size being the finest in the vicinity of the contact domain where the highest resolution is required. Mesh size is increased gradually while moving towards the peripheral areas. As for the rigid body separation, no discretization is required since it is a scalar. All equations are assembled into a single monolithic algebraic system of equations, and they are solved simultaneously, guaranteeing fast convergence rates. Given the highly nonlinear character of the generalized Reynolds equation (mainly due to the viscosity-pressure dependence), a damped-Newton procedure²⁰ is employed for the resolution. For highly loaded contacts, special stabilized finite element formulations²¹ are required for the treatment of the generalized Reynolds equation, to cure any non-physical numerical oscillations that would arise in the pressure field. As for the treatment of the free cavitation boundary that arises at the contact outlet, a penalty method is employed, as proposed by Wu.²² The numerical model description provided here is only a brief summary. For more details about governing equations and their associated boundary conditions, geometrical and computational domains, finite element formulations, meshing, convergence criteria, etc., interested readers are referred to Habchi.¹⁸

Results and discussion

There have been many accurate predictions of EHL film thickness from the real pressure and shear dependence of viscosity measured in viscometers. Some examples may be found in references.^{10,11,14,23,24} This experience gives confidence in the accuracy of the simulations here.

A total of 128 cases were simulated as shown in Table A1, 64 each for the Newtonian cases and another 64 for non-Newtonian. The combinations of the six parameters from equation (4) are given in the table along with the resulting central film thicknesses. Because q has been held constant, the ratio $(α * / α_{0})$ is constant and the exponent, $Q_{α}$ , is independent of the definition of $α$ . On average, the non-Newtonian film thickness is 72% of the Newtonian case. Adjusting the value of $α_{0}$ from 13.9 to 7 GPa⁻¹ would cause the average effect of shear dependence to be reconciled by the change in “ $α$ -value,” as can be seen from the results of Table A1. This is the procedure of classical EHL. However, the exponents of the classical formulas are not applicable to the shear-dependent cases.

The exponents have been calculated by finite differences for the Newtonian cases in Table A2 from the film thicknesses in Table A1. For example,

Q_{R} = \frac{\ln [\frac{h_{c} (R_{2})}{h_{c} (R_{1})}]}{l n (\frac{R_{2}}{R_{1}})}

The average and the maximum and minimum values are given at the bottom of the table. The average values of the exponents are in good agreement with those in equation (4). There seems to be poor agreement of the

Q_{E^{'}}

; however the range of operating parameters is different from Pan and Hamrock and these exponents are very small, meaning that

E^{'}

has little effect on the film thickness.

The exponents have been calculated for the non-Newtonian cases in Table A3 from the film thicknesses in Table A1. The average values of the exponents are not in good agreement with those in equation (4). There is little difference in the average exponents for the pressure-viscosity coefficient between the non-Newtonian results and equation (4). The average exponents for entrainment velocity and ambient viscosity, $Q_{μ_{0}}$ and $Q_{\bar{u}}$ , are slightly smaller than the Newtonian values, 0.64 versus 0.69. The average exponents for the composite radius and pressure, $Q_{R}$ and $Q_{p_{H}}$ , being 0.356 and −0.395, are significantly greater in magnitude than the Newtonian values, 0.308 and −0.332, respectively.

If the average exponents are inserted into equation (4) the result is

h_{c} = k α^{0.469} (μ_{0} \bar{u})^{0.644} E^{^{'} 0.069} R^{0.356} p_{H}^{- 0.395}

(5)

The leading coefficient, which is no longer dimensionless, was found, by minimizing the standard deviation, to be k = 17.7 N^0.151m^−0.302. The standard relative deviation is 2%. The ordinary dimensionless groups do not have parameters for shear-dependent viscosity or, incidentally, for compressibility. Equation (5) is only applicable for the shear dependence described by equation (3) with G = 30 kPa and n = 0.82.

Instead of using a high-pressure viscometer, an effective “ $α$ -value” is often extracted from an EHD rig. The consequences of not understanding the effect of shear dependence, even the mild shear-thinning of a blend of base oils, can be appreciated by examining a practical example. Consider a pair of machine elements forming a concentrated contact lubricated by the PAO blend considered here. The elements are steel with $E^{'}$ =223 GPa and R = 5 mm, entrainment velocity, $\bar{u} =$ 7 m/s, and Hertz pressure, $p_{H} =$ 2 GPa. Consider also an EHD rig with $E^{'}$ =121 GPa and R = 10 mm, entrainment velocity, $\bar{u} =$ 0.5 m/s, and Hertz pressure, $p_{H} =$ 0.5 GPa. The central film thickness measured in the EHD rig will be 87.6 nm, given by equation (5). However, the Newtonian assumption used in the Pan and Hamrock equation (4) requires that the“ $α$ -value” be adjusted from $α_{0} = 13.9 to 9.4 GP a^{- 1}$ to recover the measured film thickness. This is a typical adjustment of the pressure-viscosity coefficient which is then employed in a classical film thickness formula. Using this effective “ $α$ -value” in the Pan and Hamrock equation (4) predicts the film thickness to be 299 nm in the real machine. However, equation (5) predicts the film thickness to be 226 nm in the real machine. Thus, the Newtonian formula (4) overstates the thickness by 32% for a blend of base oils. The adjustment of the “ $α$ -value” is insufficient to accurately estimate film thickness.

Conclusion

Engineering estimates of film thickness in real machines are likely to have substantial errors because the classical film thickness formulas assume that the response in the inlet zone is Newtonian. This assumption has persisted because the pressure and shear dependence of viscosity measured in viscometers has been ignored. These errors are exacerbated if the pressure-viscosity coefficient is derived from a film thickness measurement using the same film thickness formula.

A blend of PAO base oils with mild shear-thinning has been employed in a rolling line contact to assess the effects of shear-dependent viscosity on the power-law exponents of the classical film thickness formula. The exponents on the pressure-viscosity coefficient and the elastic modulus of the solids are not sensitive to the non-Newtonian effect. The exponents on ambient pressure viscosity and velocity are slightly reduced by shear-thinning. The exponents on pressure and scale are substantially increased by the shear dependence. The usual practice of measuring film thickness in an EHD rig to obtain an effective “ $α$ -value” to use in a classical Newtonian film thickness formula will understate the film thickness of a liquid which is shear-dependent in the EHL inlet because of the sensitivity of the exponents to shear-thinning. It seems that the only accurate estimates come from a full simulation with viscometer data. The effect of sliding has not been considered here.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Scott Bair

Wassim Habchi

Appendix

References

Novak

Winer

. Some measurements of high pressure lubricant rheology. J Lubr Technol 1968; 90: 580–590.

Bair

. The rheological assumptions of classical EHL: what went wrong? Tribol Int 2019; 131: 45–50.

Barus

. ART. X.–Isothermals, isopiestics and isometrics relative to viscosity. Am J Sci, Third Series 1893; XLV: 87–96.

Roelands

CJA

. Correlational aspects of the viscosity-temperature-pressure relationship of lubricating Oils , PhD Thesis, University of Technology, Delft. 1966.

Ree

Eyring

. Relaxation theory of transport problems in condensed systems. Ind Eng Chem 1958; 50: 1036–1040.

American Gear Manufacturers Association . Effect of lubrication on gear surface distress. AGMA 925-A03, 2003.

Aerospace Information Report AIR5433 Revision C . Lubricating characteristics and typical properties of lubricants used in aviation propulsion and drive systems. Warrendale, PA: SAE Inc., 2019.

Hamrock

Dowson

. Ball bearing lubrication: the elastohydrodynamics of elliptical contacts. NewYork: Wiley, 1981.

Bair

Winer

. A simple formula for EHD film thickness of non-newtonian liquids. In: Dowson

Taylor

Childs

THC

Dalmaz

Berthier

Flamand

Georges

Lubrecht

(eds) Proceedings of the 23rd Leeds-Lyon symposium on tribology, elastohydrodynamics 1996. Tribology series. 32. Amsterdam: Elsevier, 1997, pp.235–241.

10.

Krupka

Bair

Kumar

, et al. An experimental validation of the recently discovered scale effect in generalized Newtonian EHL. Tribol Lett 2009; 33: 127–135.

11.

Krupka

Kumar

Bair

, et al. The effect of load (pressure) for quantitative EHL film thickness. Tribol Lett 2010; 37: 613–622.

12.

Ueda

Wong

Spikes

. Influence of dumbbell base oil blends on micropitting. Tribol Int 2023; 185: 108578.

13.

Bair

Qureshi

. Time-temperature-pressure superposition in polymer thickened liquid lubricants. J Tribol 2014; 136: 2.

14.

Bair

Vergne

Marchetti

. The effect of shear-thinning on film thickness for space lubricants. Tribol Trans 2002; 45: 330–333.

15.

Bair

Flores-Torres

. The viscosity of polyalphaolefins mixtures at high pressure and stress. J Tribol 2019; 141: 021802.

16.

McEwen

. The effect of variation of viscosity with pressure on the load-carrying capacity of the oil film between gear-teeth. J Inst Pet 1952; 38: 646–672.

17.

Pan

Hamrock

. Simple formulas for performance parameters used in elastohydrodynamically lubricated line contacts. J Tribol 1989; 111: 246–251.

18.

Habchi

. Finite element modeling of elastohydrodynamic lubrication problems. Chichester, UK: Wiley, 2018, ISBN: 978-1-119-22512-6.

19.

Yang

Wen

. A generalized Reynolds equation for non-newtonian thermal elastohydrodynamic lubrication. ASME J of Tribol 1990; 112: 631–636.

20.

Deuflhard

. Newton methods for nonlinear problems, affine invariance and adaptive algorithms. Germany: Springer, 2004.

21.

Habchi

Eyheramendy

Vergne

, et al. Stabilized fully-coupled finite elements for elastohydrodynamic lubrication problems. Adv Eng Softw 2012; 46: 4–18.

22.

. A penalty formulation and numerical approximation of the Reynolds-Hertz problem of elastohydrodynamic lubrication. Int J Eng Sci 1986; 24: 1001–1013.

23.

Habchi

Bair

Qureshi

, et al. A film thickness correction formula for double-Newtonian shear-thinning in rolling EHL circular contacts. Tribol Lett 2013; 50: 59–66.

24.

Liu

Wang

Bair

, et al. A quantitative solution for the full shear-thinning EHL point contact problem including traction. Tribol Lett 2007; 28: 171–181.