## Abstract

Realistic prediction of the characteristics of the elastohydrodynamic lubrication (EHL) contact requires consideration of the appropriate constitutive equation for the lubricant. In many applications, the lubricant exhibits a shear-thinning behaviour which significantly affects the film thickness. In this paper, we present a generalized formulation that can efficiently treat shear-thinning fluids with provision for compressibility in the EHL line contact. Specifically, the Carreau model and the sinh-law model are investigated. An extensive set of numerical solutions and comparison with experiments reveal that the Carreau equation properly captures the film thickness behaviour under both rolling and sliding conditions.

## 1. Introduction

The lubrication regime of many machine elements with non-conformal contacts such as rolling element bearings, gears, etc. is generally classified as elastohydrodynamic lubrication (EHL). In EHL, the non-conformance of the bodies results in high contact pressures so that deformation of the surfaces plays an important role. The origin of the theory dates back to 1949 when a paper by Grubin & Vinogradova (1949) set the stage for the EHL theory by incorporating both the elastic deformation of the surfaces and the pressure–viscosity characteristics of the lubricant. Subsequently, Dowson & Higginson (1959) developed a well-known expression for the minimum film thickness in the EHL line contact, which was verified by Crook's experiments (1961). Progress in analytical treatment of EHL continued in many different directions from the optimization of numerical simulations (Okamura 1982) to the development of algorithms that could handle higher loads (Houpert & Hamrock 1986).

In the 1970s, it was also realized that the non-Newtonian character of the fluid must be taken into account in order to close the gap between theory and experimentally obtained film thickness and traction data. Trachman & Cheng (1972) used the hyperbolic relationship in their traction analysis. Their results were in good agreement with the data by Johnson & Cameron (1967). Johnson & Tevaarwerk (1977) proposed a nonlinear Maxwell model which included the sinh law. Salehizadeh & Saka (1991) developed a thermal non-Newtonian EHL model under the pure rolling condition with the sinh-law constitutive equation. Wang *et al.* (1991) extended the thermal EHL problem with the sinh-law model. Sui & Sadeghi (1991) also solved the thermal EHL problem by using the sinh-law model. Bair & Winer (1979) used primary measurement to show that, in EHL applications, lubricants exhibit a limiting shear stress and proposed a new constitutive equation. Khonsari & Hua (1993) developed a generalized non-Newtonian EHL model by using the Bair–Winer constitutive equation, which yielded close agreement in the prediction of traction coefficient compared with experimental measurements.

The generalized Newtonian model considers only shear-dependent viscosity (shear-thinning) and ignores other aspects of non-Newtonian response such as time dependence, normal stress differences and enhanced resistance to extensional deformation, which are all a part of non-Newtonian flow. The Carreau viscosity model is one of the viscosity functions available in the generalized Newtonian flow application. Mongkolwongrojn *et al.* (2005) calculated the performance characteristics of the EHL line contact by varying the Carreau viscosity model parameters.

In this paper, we focus our attention on the EHL line contact of shear-thinning fluids. An appropriate formulation and numerical solution algorithm are provided, which can effectively treat the problem at high loads and operating speeds. Two types of constitutive relationships are examined: the Carreau model and the sinh-law model. Results are compared with the film thickness experiment of Dyson & Wilson (1965–66). It is shown that the Carreau model is capable of accurately predicting the behaviour of EHL films under both rolling and sliding contacts.

## 2. Theory

Elastohydrodynamic theory requires simultaneous solution of the modified Reynolds equation with provision for the generalized Newtonian constitutive equation and the elasticity equation for the surface deformation to obtain the pressure distribution and the lubricant film thickness. Before discussing the details of the formulation, we first present a simplified method for predicting film thickness based on an inlet zone analysis.

### (a) Inlet zone analysis

Bair & Khonsari (1996) developed an EHL inlet zone analysis to predict the film thickness for a pure rolling case for any constitutive equation. Bair (1998) extended this simple numerical method to include the sliding velocity. It can be shown that the equilibrium equation for thin-film lubrication simplifies to(2.1)Integrating equation (2.1), the shear stress across the film is given by(2.2)where *τ*_{1} is the shear stress on the lower surface. The constitutive law for a generalized Newtonian fluid is(2.3)Integrating equation (2.3) yields the local fluid velocity(2.4)where *u*_{1} is the velocity of the lower surface. Conservation of mass requires that(2.5)The central film thickness *h*_{0} must be found such that the isoviscous pressure is close to the asymptotic value at the Hertzian boundary. This requirement effectively ends the Poiseuille component of flow at the edge of the Hertz region and the pressure boundary used in this theory is (Blok 1963)(2.6)where *α* is the pressure–viscosity coefficient. Equations (2.1)–(2.6) enable one to set up an iteration scheme for predicting the central film thickness. To begin the numerical solution, the central film thickness, *h*_{0}, and the pressure gradient, , are guessed. In the inner loop, the program iterates on the shear stress on the lower surface *τ*_{1} until equations (2.1) and (2.2) are satisfied. In the second loop, the program iterates on the pressure gradient until equation (2.5) is satisfied. Then, in the outer loop, the program iterates on the central film thickness until the pressure boundary (2.6) is satisfied.

### (b) Generalized Reynolds equation

The appropriate Reynolds equation for an EHL line contact is (Dowson 1962)(2.7)The Reynolds equation (2.7) is applicable to any generalized Newtonian constitutive equation. The equivalent viscosity *η* appears in the integral coefficients *F*_{0}, *F*_{1} and *F*_{2}, and they are defined as(2.8)Integrating equation (2.7), we obtain(2.9)where the rolling velocity *u* and the sliding velocity *s* are defined as and , respectively. The slide–roll ratio is defined as . The dimensionless parameters are defined as(2.10)Parameter is the semi-width of the Hertzian contact and *w* is the load capacity per unit length. The equivalent radius *R* and elastic modulus *E* are and , respectively.

The dimensionless form of the Reynolds equation now reads(2.11)where *A* is a constant to be determined and . The integral coefficients in dimensionless form are(2.12)

The boundary conditions for the pressure distribution are(2.13)where *X*_{min} and *X*_{end} are the minimum values of *X* in the mesh and the outlet meniscus distance, respectively.

### (c) Velocity distribution

The velocity and the velocity gradient across the film are as follows:(2.14)Defining , the dimensionless velocity gradient is(2.15)where and .

### (d) Elastic deformation

The elastic deformation in a semi-infinite solid in a state of plane strain can be expressed as(2.16)Letting , the expression for the dimensionless elastic deformation reads(2.17)The elastic deformation is determined by assuming that the pressure is described by a second-degree polynomial in the interval [*X*_{j−1},*X*_{j+1}]. Therefore, the elementary deformation can be computed easily. The elementary deformation represents the deformation at node *i* due to pressure defined in the interval [*X*_{j−1},*X*_{j+1}]. The total elastic deformation at node *i* is the sum of all the elementary deformations (Houpert & Hamrock 1986), i.e.(2.18)Using equation (2.17), the elementary deformation is given by(2.19)Inserting the pressure gradient of a second-degree polynomial into equation (2.19) and integrating the equation, the elementary deformation can be expressed in the following form:(2.20)where(2.21)Using *Z*_{min}=*X*_{i}−*X*_{j−1} and *Z*_{max}=*X*_{i}−*X*_{j+1}, the parameters *M*_{1}, *M*_{2} and *K* used in equation (2.21) are defined as(2.22)The parameters *a*_{1}–*a*_{6} are(2.23)Inserting equation (2.20) into equation (2.18), the elastic deformation is given by(2.24)where if *j* is an even number. If *j* is an odd number, is the sum of in the interval [*X*_{j−2},*X*_{j}] and d*D*_{i,j−1} in the interval [*X*_{j},*X*_{j+2}]. For the first and last nodes, *D*_{i,1}=d*D*_{i,1} and *D*_{i,N}=d*D*_{i,N}. Finally, the film thickness is given by(2.25)*H*_{00} is an unknown parameter to be determined, which is the sum of the central film thickness, *H*_{0}, and the constant in equation (2.24).

### (e) Load condition

The dimensionless constant load condition is given by(2.26)Inserting the pressure distribution of a second-degree polynomial into equation (2.26), we can define the coefficients d*C*_{j,j−1}, d*C*_{j,j} and d*C*_{j,j+1} in the interval [*X*_{j−1},*X*_{j+1}], such that(2.27)where(2.28)Finally, the load condition can be written as(2.29)The weight factors *C*_{j} are defined as *C*_{j}=d*C*_{j,j} if *j* is an even number. If *j* is an odd number, *C*_{j} is the sum of d*C*_{j,}_{j+1} in the interval [*X*_{j−2},*X*_{j}] and d*C*_{j,}_{j−1} in the interval [*X*_{j},*X*_{j+2}]. In the interval [*X*_{N},*X*_{end}], the weight factor is corrected appropriately as suggested by Houpert & Hamrock (1986).

### (f) Location of *X*_{end}

*X*_{N} is the nearest node to *X*_{end}, such that *X*_{N}<*X*_{end}<*X*_{N+1}. Note that the pressure and its gradient is zero at *X*_{end}. *X*_{end} is calculated from the Reynolds equation(2.30)where the subscript c represents the location where the pressure gradient is zero in the contact zone.

### (g) Density

To obtain more accurate results, experimental measurements for the density are used for the simulations. The variation of density can be expressed by the Tait equation (Bair 2007)(2.31)For polydimethylsiloxane, the constants are =10.19 and *K*_{0}=0.85 GPa (Sachdev *et al.* 1998).

### (h) Numerical formulation

The system of the governing equations is solved by the Newton–Raphson method. The problem has (*N*+1) unknowns: (*N*−1) pressures at nodes 2 to *N*, and two constants *A* and *H*_{00} (Houpert & Hamrock 1986). The mesh size across the film is uniform, while the mesh size in the moving direction increases exponentially, leading to fine computational grids at the contact region and an improved computing time. The Newton–Raphson algorithm requires a linear system of (*N*+1) equations. *N* equations from the Reynolds equation and one equation from the load condition,(2.32)where(2.33)Therefore, a linear system of (*N*+1) equations is(2.34)The Gauss elimination method is used to solve this linear system. The partial differential coefficients of the Jacobian factors are(2.35)In the equations of (2.35), four partial differential coefficients , , and should be determined(2.36)In equation (2.36), adopting the Carreau model, the two partial differential coefficients and are(2.37)where(2.38)In equation (2.37), the two partial differential coefficients and are(2.39)andwhere(2.40)The partial differential coefficient can be determined from the first equation in (2.37) and (2.39). Letting , we have(2.41)where *M* is the number of computational nodes placed across the film. Similarly, from the second equation in (2.37) and (2.39), the partial differential coefficient can be determined. Letting , we have(2.42)The pressure gradient is calculated using the two-node formula(2.43)

The pressure is described by a second-degree polynomial between and . Therefore, the pressure gradient at is(2.44)

The solution algorithm begins by using the Hertzian pressure as an initial guess. Values for the two unknowns *A* and *H*_{00} are also assumed. Next, the film thickness is computed with the pressure distribution, and the viscosity *μ*_{1} is computed from equation (2.53). Then, the effective viscosity is determined iteratively using the constitutive equation. Finally, the linear system of (2.34) is solved for the pressure distribution using the Gaussian elimination method. The entire process is repeated until the pressure distribution converges. The above derivations and numerical solution algorithm can effectively treat any generalized Newtonian constitutive equation. In what follows, we focus our attention on the hyperbolic sine (sinh) law and the Carreau model for the characterization of a shear-thinning fluid.

### (i) The sinh law in rheology

The experimental origin of the sinh law,(2.45)is in the thermal softening that results from viscous heating in capillary viscometers, where the apparent shear rate is and the wall stress is . Before analytical solutions were obtained for viscous heating, these severely nonlinear flow curves were mistaken for non-Newtonian constitutive behaviour (Gruntfest 1965). This interpretation was strengthened by the appearance, in 1936, of Eyring's theory of stress-assisted thermally activated flow (Eyring 1936), which predicted a relationship between rate and stress in the form of a sinh law. It was thought that the nonlinearity in equation (2.45) was a manifestation of shear-thinning or thixotropy with 1<*τ*_{0}<30 kPa (Norton *et al.* 1941; Hahn *et al.* 1958) obtained from the experiment. The Eyring theory did not provide an estimate of *τ*_{0}. By 1948, Alfrey had revealed that the mechanics of the stress-assisted portion of the Eyring theory were unsound. Investigations of viscous heating in capillary flow (Gerrard *et al.* 1965; Sukanek & Lawrence 1974) clearly showed that flow curves of a sinh-law form are produced by thermal softening when a critical value of the Nahme–Griffith (Brinkman) number is exceeded. Even Eyring was active in refuting the constitutive interpretation of the simple sinh law (Ree *et al.* 1958). Ree & Eyring (1955) responded with a new model for shear-thinning and thixotropy employing a series of *N* inverse sinh terms (flow units), which with a large number *N* of flow units could reproduce complete experimental isothermal flow curves(2.46)Following these events, the sinh law all but disappeared as a description of shear-thinning (while it continues to be useful for the structural viscosity of thixotropic materials). The exception, of course, has been the field of EHL. Traction curves obtained from full-film lubricated concentrated contacts at intermediate pressures have the form of a sinh law,(2.47)where *C*_{E} is a constant for a given load; *τ* is the product of the traction coefficient and the average pressure; and the ‘Eyring stress’ is 0.8<*τ*_{0}<13 MPa (Hirst & Moore 1980). This contact behaviour is not disputed. Invoking Eyring's original theory, the EHL community has maintained that this traction response is evidence of the constitutive behaviour of the sinh-law form. In order to support this claim, however, it must be assumed that the piezoviscous response becomes very weak at pressures of the Hertzian zone, a proposition that has been refuted by experimental measurements for 80 years (Bridgman 1926). It is noteworthy that equation (2.47) can be useful for dry contacts as well, with 30<*τ*_{0}<90 MPa. The effective viscosity for the sinh-law model used in the theory is(2.48)

### (j) The Carreau model

In contrast with the usual EHL description of shear-dependent viscosity, experimental measurement (Tanner 2000), kinetic theory (Bird *et al.* 1987) and non-equilibrium molecular dynamics simulation (Kioupis & Maginn 1999) have clearly shown that the low shear behaviour of a generalized Newtonian liquid can be described with a shear-independent viscosity, *μ*, followed by a power-law regime at high shear, with , where *n* is the power-law index or exponent. The critical shear rate at the transition from limiting low shear response to the power-law regime occurs at , where *λ* is the rotational relaxation time of a molecule (Cui *et al.* 1998) for monodisperse, low-molecular-weight liquids. The Einstein–Debye relation gives(2.49)for the rotational relaxation time (Paluch *et al.* 2003). In equation (2.49), *M* is the molecular weight and *R*_{g} is the gas constant. While it may be appreciated that the Ree–Eyring model in equation (2.46) is more general than others and capable of describing the shear dependence of the viscosity of many-component mixtures, the requirement of as many as five flow units with as many as nine parameters (Kim *et al.* 1960) makes this model unwieldy for EHL calculations. Experimental measurements and molecular simulations favour the Carreau law (Carreau 1972) for simplicity and accuracy,(2.50)Since laboratories have been generating flow curves for more than 50 years, it is not surprising that there are many models with properties similar to the Carreau law (Bair & Qureshi 2003). In this simulation, the Carreau model is expressed by the following equation (Bair & Khonsari 2006):(2.51)where *μ*_{1} is the first Newtonian viscosity at low shear rates and *μ*_{2} is the second Newtonian viscosity at high shear rates. Parameter is a critical stress. Letting *Γ*=*G*_{cr}/*E*, the dimensionless Carreau model reads(2.52)For polydimethylsiloxane, which is the lubricant used in simulations that resulted in figures 2–7, the equation for the viscosity variation with the pressure (Bair & Qureshi 2003) is given by(2.53)

## 3. Results and discussion

### (a) Verification

We begin by providing a comparison between the EHL theory and the inlet zone analysis described in §2*a*. Figure 1 shows the central film thickness using the Carreau model plotted as a function of the critical stress. The input data used in this simulation are *E*=211 GPa, *R*=20 mm, *b*=0.39 mm, *μ*_{0}=0.08 Pa s and *α*=15 GPa^{−1}. In figure 1, the variation of the first Newtonian viscosity with the pressure is used instead of equation (2.53). *μ*_{0} is the absolute viscosity of the lubricant at atmospheric pressure. The dimensionless viscosity is , where *G*^{*}=*bG*/4*R* and *G*=*Eα*. Also shown are the film thickness predictions presented by Bair & Khonsari (1996) based on the inlet zone analysis, which is very simple and predicts the central film thickness effectively. However, it predicts a single film thickness rather than a profile since the computation domain is confined to the inlet region. The results show that the predicted central film thickness is very close to the film thickness computed by the present theory. The film thickness predicted using the inlet analysis consistently underpredicts the film thickness compared with those of the present theory. At low to moderate values of shear stress, the central film thickness increases with increasing critical stress. At high values of critical stress (approx. 1 MPa under current operating conditions), the increment in the film thickness becomes small since the shear-thinning effect decreases. The film thickness is larger with higher *n* due to the higher effective viscosity. As the critical stress increases, the film thickness with *n*=0.75 approaches that with *n*=0.5, since the shear-thinning effect becomes small, and therefore, the film thickness with a greater critical stress than approximately 3 MPa approaches the Newtonian predictions, as expected.

Figure 2 shows a comparison of the dimensionless film thickness profiles with Dyson & Wilson's experimental data (1965–66) under pure rolling and sliding conditions. Also shown are theoretical isothermal models (Dowson & Higginson formula and Newtonian, Carreau and sinh-law models). The lubricant used in Dyson & Wilson's experiments (1965–66) is polydimethylsiloxane. The input data used in this simulation are *E*=224 GPa, *R*=19.1 mm, *n*=0.5 and *μ*_{0}=0.86 Pa s, leading to the dimensionless parameter *W*=2.58×10^{−5}. The density and the molecular weight for polydimethylsiloxane are *ρ*=965 kg m^{−3} and *M*=28 kg gmol^{−1}, respectively. Using the relationship *G*_{cr}=*ρR*_{0}*T*/*M*, the critical stress is determined as *G*_{cr}=80 kPa at room temperature.

Let us turn our attention to the well-established Dowson & Higginson formula for the central film thickness(3.1)The results for the Newtonian model are in good agreement with the Dowson & Higginson formula at relatively low rolling speeds. At higher rolling speeds, the present analysis predicts a slightly lower film thickness than the generalized expression of equation (3.1).

The results in figure 2 reveal that the Carreau model also predicts that the film thickness increases nearly in a power-law fashion with increasing rolling speed, and that the film thickness is lower at a higher slide–roll ratio, as expected. The results for the Carreau model at *G*_{cr}=80 kPa are in very good agreement with the experimental data by Dyson & Wilson for both rolling (*Σ*=0) and sliding (*Σ*=0.5) conditions.

The same algorithm can be used to simulate any constitutive equation. The results for the sinh-law model at *τ*_{0}=80 kPa are presented in figure 2 to compare with the Carreau model. The results show that at *u*=5 cm s^{−1}, the film thickness for the pure rolling case matches well with the measurement. However, the rate of the reduction in the film thickness with increasing rolling speed is relatively large in comparison with the experimental data. It is also shown that the effect of the sliding for the sinh-law model is improperly large in comparison with the measurements or that for the Carreau model. Therefore, it seems that the sinh-law model does not lead to the proper treatment of film thickness. Figure 3 clearly shows that the sinh-law model could not capture the behaviour of the shear-thinning fluids at any *τ*_{0} in comparison with the experimental measurements.

The limitation of the Reynolds equation in the EHL regime, where the viscosity varies with pressure, is associated with the situations where the cross-film pressure gradient may be significant and secondary flows result, even for the flow between parallel plates (Bair *et al.* 1998; Rajagopal & Szeri 2003). It has been shown that the Reynolds equation adequately captures the mechanics of fluids only when *τα*<1 (Bair *et al.* 1998). In figure 2, the maximum shear stress, *τ*_{max}=3 MPa, occurs for the sliding case at *u*^{*}=2000 cm s^{−1}. A conservative value of *α* is 10^{−8} Pa^{−1} and, therefore, *τα*=0.03<1. In this paper, most of the simulations pertain to *τα*<1 and the Reynolds equation is valid.

### (b) Pressure profile

In what follows, we present the results of a series of parametric performances of shear-thinning fluids using the Carreau model. The property values used in these simulations are the same as those used in the prediction of Dyson & Wilson's experiments. Figure 4 shows the comparison of the pressure profiles between the Newtonian and Carreau models. In this simulation, the sliding velocity is assumed to be zero. The magnitude of the pressure spike for the Carreau model is significantly smaller than that for the Newtonian model due to the shear-thinning effect leading to a lower effective viscosity. The film thickness for the Carreau model is also predicted to be smaller than that for the Newtonian model. In the contact region, the pressure is very high and, therefore, the viscosity is an order of magnitude greater than its value at atmospheric pressure. As a result, the film thickness in the contact region is relatively flat such that the pressure gradients have a realistic value.

### (c) Effect of load

Figure 5 shows the variation of the pressure with applied load using the Carreau model. The film thickness decreases with increasing load, as expected. The magnitude of the pressure spike near the end of the contact drops significantly with increasing load, and in the limit when the load is very large, the pressure approaches a parabolic-shape Hertzian contact pressure and the spike disappears. Under this condition, the inlet meniscus moves to *X*=−1. These trends are easily captured by the computational algorithm derived in §2*h*.

The film thickness shows a relatively large reduction near the outlet region, which corresponds to the peak pressure. Under a high loading condition, the Hertzian load is a close approximation to the actual load. A large pressure gradient must exist near the outlet region for the pressure to drop to the ambient value. This large pressure gradient at a comparatively low pressure can be achieved only by a reduction in the film thickness (Cameron 1966).

### (d) Effect of speed

Figure 6 shows the variation of the minimum film thickness with applied load and the slide–roll ratio. The minimum film thickness decreases rapidly at a low load and shows a relatively small reduction at a high load. The results show that the minimum film thickness decreases with increasing slide–roll ratio due to the lower effective viscosity associated with the shear-thinning effect. At a relatively low load, the effect of sliding velocity is small.

Figure 7 shows the variation of pressure with rolling speed. The slide–roll ratio is fixed at *Σ*=0.1. The film thickness decreases with decreasing rolling speed, as expected. The diminution in film thickness near the end of the contact increases when the film thickness increases. The location of the peak pressure is shifted to the centre as the rolling speed increases. As the rolling speed decreases, the magnitude of the pressure spike near the end of the contact drops and the inlet meniscus moves to *X*=−1.

## 4. Concluding remarks

An EHL line contact model for simulating the behaviour of shear-thinning fluids is developed using a generalized Reynolds equation. The algorithm is efficient and the numerical solution converges at a high applied load as well as at high values of rolling–sliding ratios. Owing to the shear-thinning effects, the peak pressure and the film thickness are smaller than those for the Newtonian models. The evaluation of the critical stress and power-law exponent for the non-Newtonian lubricant is also important for the accurate prediction of the pressure and the film thickness. Two constitutive equations are considered: the sinh-law model and the Carreau model. The sinh-law model does not capture the behaviour of the EHL line contact when compared with experimental measurements. In contrast, the Carreau model effectively predicts the film thickness that closely matches experiments and properly characterizes the behaviour of shear-thinning lubricants.

## Acknowledgments

This research was in part supported by Caterpillar, Inc. Two of the authors (J.Y.J. and M.M.K.) gratefully acknowledge the financial support. The authors are solely responsible for the results and conclusions presented.

## Footnotes

- Received June 3, 2007.
- Accepted September 7, 2007.

- © 2007 The Royal Society