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 “ $\alpha$ -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.

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 “ $\alpha$ -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{\alpha }^{0.53}({\mu }_0\overline{u}{)}^{0.67}{E}^{\mathrm{\prime }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/R1+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 “ $\alpha$ -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 $(\mathrm{\partial }\ln h_c/\mathrm{\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:

$$\mu ={\mu }_0{\left(1+{\displaystyle \frac{{\alpha }_0}{q}p}\right)}^q={\mu }_0{\left(1+{\displaystyle \frac{\alpha *}{q-1}p}\right)}^q,q>1$$

(2)This equation is accurate for the PAO blend at 70 °C to a pressure of 500 MPa using ${\mu }_0=$ 24 mPa·s, ${\alpha }_0=13.9\mathrm{GP}{\mathrm{a}}^{\text{-}1}$ , $\alpha *=12.0\mathrm{GP}{\mathrm{a}}^{\text{-}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

$$\eta =\mu {\left[1+{\left({\displaystyle \frac{\tau }{G}}\right)}^2\right]}^{\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{\alpha }^{0.47}({\mu }_0\overline{u}{)}^{0.692}{E}^{{}^{\mathrm{\prime }}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=\mathrm{\partial }\ln h_c/\mathrm{\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 $(\alpha */{\alpha }_0)$ is constant and the exponent, $Q_{\alpha }$ , is independent of the definition of $\alpha$ . On average, the non-Newtonian film thickness is 72% of the Newtonian case. Adjusting the value of ${\alpha }_0$ from 13.9 to 7 GPa⁻¹ would cause the average effect of shear dependence to be reconciled by the change in “ $\alpha$ -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={\displaystyle \frac{\ln \left[{\displaystyle \frac{h_c(R_2)}{h_c(R_1)}}\right]}{ln\left({\displaystyle \frac{R_2}{R_1}}\right)}}$$

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^{\mathrm{\prime }}}$ ; however the range of operating parameters is different from Pan and Hamrock and these exponents are very small, meaning that $E^{\mathrm{\prime }}$ 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_{{\mu }_0}$ and $Q_{\overline{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{\alpha }^{0.469}({\mu }_0\overline{u}{)}^{0.644}{E}^{{}^{\mathrm{\prime }}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 “ $\alpha$ -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^{\mathrm{\prime }}$ =223 GPa and R = 5 mm, entrainment velocity, $\overline{u}=$ 7 m/s, and Hertz pressure, $p_H=$ 2 GPa. Consider also an EHD rig with $E^{\mathrm{\prime }}$ =121 GPa and R = 10 mm, entrainment velocity, $\overline{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“ $\alpha$ -value” be adjusted from ${\alpha }_0=13.9\mathrm{to}9.4\mathrm{GP}{\mathrm{a}}^{\text{-}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 “ $\alpha$ -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 “ $\alpha$ -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 “ $\alpha$ -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.