Non-Newtonian Thermo-Elastohydrodynamics and Sub-Surface Stress Field of High-Performance Racing Spur Gears

Abstract

Meshing teeth pairs of involute spur gears often form the final drive of high-performance motorsport transmissions. They are subject to high normal and shear loading. Under transient conditions pertaining to a meshing cycle, the contact conditions alter from the onset of teeth pair engagement through to maximum normal loading, followed by contact separation. Sliding motion only ceases instantaneously at the pitch point. The regime of lubrication remains mostly in non-Newtonian thermo-elastohydrodynamic conditions. The results show that a starved inlet boundary is attained throughout most of the meshing cycle which leads to the diminution of the pressure spike at the exit from the contact conjunction. The reversing sub-surface shear stresses are the main source of the onset of any inelastic deformation, which is dominated by the primary pressure peak in compliance with the Hertzian maximum pressure. The shear stress field is supplemented by an induced field due to the presence of the pressure spike. Under starved conditions this secondary stress field is diminished. The combined solution of elastohydrodynamics with a thermal network model, non-Newtonian lubricant traction, and sub-surface stress evaluation provides for a comprehensive solution not hitherto reported in the literature.

1. Introduction

Spur gears are used in many engineering applications, including some automotive transmission systems. These are highly stressed components subjected to high normal contact pressures, friction and generated contact temperatures. Therefore, gears are often the limiting factor in the durability, operational integrity and performance of the systems of which they are a part. Their harsh operating environment often contributes to inefficiency, power loss, emissions, and noise, vibration and harshness (NVH). Unfortunately, some of these required key attributes are contradictory. For example, to reduce NVH, good tractive action is required, which is contrary to the desired operational efficiency that requires reduced friction. Nevertheless, some degree of traction is essential for the transmission of power [1,2]. The aim is to optimize the extent of traction in order to attain a high degree of efficiency. Therefore, the refinement of all gearing systems is a key attribute that requires detailed predictive analysis [3].

An early tribo-dynamic approach for a wavy-surfaced disc was presented by Mehdigoli et al. [4]. They showed agreement with the experimental findings of Dareing and Johnson [5] who noted that under elastohydrodynamic conditions lubricant viscous damping is rather negligible. Any attenuation in NVH is actually the result of generated contact friction. Later, Li and Kahraman [2] developed a tribo-dynamic model to predict friction throughout a meshing cycle of a spur gear pair. Their study coupled gear dynamics with the tribology of meshing teeth for a representative involute tooth profile. Of course, deviations from an involute shape, such as tip relief and crowning [6,7,8,9,10,11], as well as the surface topography of the meshing teeth, affect the contact conditions [12,13].

Meshing gears are subjected to highly loaded concentrated contacts, resulting in two main forms of failure: wear and fatigue spalling [14]. The former is often caused by lubricant starvation, resulting in the formation of a thin lubricant film of insufficient thickness [15,16,17,18]. The latter occurs because of high sub-surface stresses which can result in inelastic deformation in the bulk of the solid when sufficiently high stresses coincide with any existing pores, inclusions or flaws. Near-surface and sub-surface stresses are the intrinsic causes of failure in pitting and spalling [19]. Clearly, such failures affect the durability and structural integrity of mechanical components such as gears and bearings [14,20].

The onset of fatigue spalling determines the useful life of gearing pairs [21,22,23], so an accurate prediction of sub-surface shear stresses is essential. For bearings and gears, the reversing orthogonal shear stresses are the main determining factor for useful life [14,24]. Huber and Fuchs [25] were the first to present methods of evaluation of the sub-surface stress field. Lyman [26] emphasized the importance of cyclic reversing orthogonal shear stresses in rolling contacts, whilst Poritsky [27] investigated their effect in gearing contacts and locomotive wheels. Johnson [28,29] reviewed earlier investigations of contact mechanics and the fatigue of solids. Johns-Rahnejat and Gohar [30] provided analyses of 3D sub-surface stress fields for point contact conditions, presenting sub-surface stress fields for elastohydrodynamic (EHD) conditions. Under EHD, the secondary pressure peak (i.e., pressure spike) at the contact exit induces localized sub-surface fields of its own quite near the surface, which can lead to premature failure. Houpert et al. [31] studied the effect of pressure spike on rolling bearing fatigue. To capture these effects and those due to asperity interactions of rough counter faces [32] a suitably detailed numerical analysis is required as highlighted by Johns-Rahnejat [30,33]. Recently, a generic semi-analytical model was presented by Johns-Rahnejat et al. [34] which can be adopted for the evaluation of sub-surface stress fields for all gear contacts with varying degrees of conformity. This approach is used here for the case of lubricated contact of spur gear pairs.

The paper presents thermo-elastohydrodynamics of a highly loaded spur gear pair of a racing vehicle, where the lubricant is subjected to high shear, thus resulting in non-Newtonian lubricant behaviour. Realistic applied loads, meshing geometry and contact kinematics are necessary for tribological studies. These are obtained by lubricated loaded tooth contact analysis (LLTCA), highlighted in detail by Oglieve et al. [35]. Unlike the usual TCA [36,37,38] LLTCA takes into account the effect of friction in sliding lubricated contact, as well as the measured profile of contacting surfaces including their real rough surface topography and any applied micro-geometrical form modifications [35]. The current paper also incorporates the assessment of contact integrity through the determination of the sub-surface stress field, including the effect of generated elastohydrodynamic lubrication (EHL) pressure spikes and friction. Therefore, the paper presents a novel comprehensive integrated analysis, not hitherto reported in the literature.

2. Methodology

2.1. Lubricated Loaded Tooth Contact Analysis (LLTCA)

The LLTCA approach is used to obtain the instantaneous contact geometry and sliding velocity of a pair of meshing involute teeth. A comprehensive description of LLTCA is provided by Oglieve et al. [35] and Sivayogan et al. [39]. Here, only a brief description is provided. Essentially, the surface geometry of the meshing teeth pairs is measured using a coordinate measuring machine (CMM) with a high measurement sensitivity of ±1.5 μm. This action enables the inclusion of realistic surface anomalies. Under steady-state conditions, cyclic meshing assumes that all the mating teeth pairs undergo the same normal loading and surface shear. Thus, to reduce the burden of computation, only 3 consecutive teeth pairs in a meshing cycle are measured. Finite element analysis is used to determine the contact mechanics (deformation) and kinematics (surface velocities) of contacting teeth flanks within a meshing cycle. This is the usual approach for the TCA of meshing teeth, as described in [36,37,38]. With LLTCA, the information obtained through TCA is supplemented by lubricated viscous friction. For the high-performance transmission, originally studied by Elisaus et al. [9] and the subject of the current analysis, high contact loads and traction lead to the non-Newtonian shear of the lubricant. Under these conditions, friction is obtained as:

$$F_f=\mu w$$

(1)

where $w$ is the contact load obtained as the integrated instantaneous elastohydrodynamic pressure distribution (Equation (41)).

The average (Pascal) pressure is obtained as:

$$\overline{p}=\frac{W}{A}$$

(2)

where the area of contact of an involute gear teeth pair A is a rectangular strip (Figure 1) obtained through Hertzian theory as:

$$A=2aL$$

(3)

where L is the length of the rectangular contact strip (assumed to be the width of the tooth flank along the line of action) and a is the contact semi-half-width (Figure 1):

$$a={\left(\frac{4wR}{\pi L{E}^{*}}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$

(4)

where R is the reduced radius of a non-conforming pair of contacting ellipsoidal solids, such as a pair of involute teeth, and $E^{*}$ is the effective Young’s modulus of elasticity of the contacting pair (plane strain modulus for the same material pair):

$$\frac{1}{R}=\frac{1}{R_p}+\frac{1}{R_w}$$

(5)

where the radii of curvature of the meshing pinion and wheel teeth pairs at any instant of time are obtained following the standard methodology highlighted by Merritt [40].

$$E^{*}=\frac{E}{1-{\nu }^2}$$

(6)

where E is the Young’s modulus of elasticity and $\nu$ is the Poisson’s ratio.

For the non-Newtonian behaviour of the lubricant at high load and shear, the coefficient of friction is given as [41]:

$$\mu =0.87\alpha {\tau }_0+1.74\frac{{\tau }_0}{\overline{p}}\ln \left[\frac{1.2}{{\tau }_0{h}_c}{\left(\frac{2k_l{\eta }_0}{1+9.6\zeta }\right)}^{1/2}\right]$$

(7)

where:

$$\zeta = \frac{4k_l}{\pi h_c/R}{\left(\frac{2\overline{p}}{E^{*}{k}_{s}R{\rho }^{\prime }{c}^{\prime }U}\right)}^{\frac{1}{2}}$$

(8)

where $\overline{p}$ and ${\tau }_0$ are the average contact pressure (Equation (2)) and characteristic shear stress of the lubricant, respectively. $U$ is the speed of lubricant entrainment into the contact (Equation (10)). The mean pressure piezo-viscosity coefficient of the lubricant is obtained as:

$$\alpha =\ln \left(\eta /{\eta }_0\right)/\overline{p}$$

(9)

Density and the specific heat capacity of the gear material are denoted by ${\rho }^{\prime }$ and $c^{\prime }$, respectively. $K_s$ and $K_l$ are the thermal conductivity of the solid and lubricant, respectively, and $h_c$ is the central contact lubricant film thickness obtained through elastohydrodynamic lubrication analysis.

The speed of lubricant entrainment into the contact at any instant within a meshing cycle is required. This is given as:

$$U= \left(u_p+u_w\right)/2$$

(10)

The tooth contact analysis provides the surface speed of the meshing teeth pair with the progressive pinion angle $\phi$ as (see ) [40]:

$$u_p=R_p{\omega }_p\left(\sin \phi +\frac{2l}{D_p^{\prime }}\right)$$

(11)

$$u_w=R_w{\omega }_w\left(\sin \phi +\frac{2l}{D_w^{\prime }}\right)$$

(12)

where ${\omega }_p$ and ${\omega }_w$ are the angular velocities of the pinion and wheel gear, and the pitch diameter for the pinion and wheel are shown by $D_p^{\prime }$ and $D_w^{\prime }$ in Figure 2.

2.2. Elastohydrodynamic Lubrication (EHL)

The general form of Reynolds equation is:

$$\frac{\partial }{\partial x}\left[\frac{\rho h^3}{6\eta }\left(\frac{\partial p}{\partial x}\right)\right]+\frac{\partial }{\partial y}\left[\frac{\rho h^3}{6\eta }\left(\frac{\partial p}{\partial y}\right)\right]=\frac{\partial \left(\rho hU\right)}{\partial x}+\frac{\partial \left(\rho hV\right)}{\partial y} + \frac{\partial \left(2\rho h\right)}{\partial t}$$

(13)

where p, h, η and ρ are pressure, film thickness, dynamic viscosity and density. U and V are surface velocities along the direction of entraining motion, x, and in the side-leakage direction, y. Olver [15] showed that spur gears generally run starved. Therefore, it can be assumed that the effect of lubricant side-leakage is negligible (i.e., V = 0), so a one-dimensional (1D) EHL solution would suffice. Furthermore, a quasi-static analysis of the meshing cycle is carried out, where the effect of squeeze film motion $\partial \left(2 \rho h\right)/\partial t$ is neglected. This is a common approach [12,16,42] when gear dynamic response or NVH problems originating from gear contact, such as the gear rattle reported in [43], are not of primary concern. In such cases convergence is sought in line with instantaneous quasi-static equilibrium. The current study’s main focus and contribution is in the area of the analysis of the sub-surface stress field. A more comprehensive approach would require the effect of transience through the inclusion of a squeeze film term as an extension of the current analysis [44]. With the stated assumptions, Reynolds equation for this case simplifies to:

$$\frac{\partial }{\partial x}\left[\frac{\rho h^3}{6\eta }\left(\frac{\partial p}{\partial x}\right)\right]=\frac{\partial \left(\rho hU\right)}{\partial x}$$

(14)

The elastic film shape is given as:

$$h=h_0+s\left(x\right)+\delta \left(x\right)$$

(15)

where $h_0$ is the minimum film thickness, s is the geometric profile of the teeth pair conjunction, and $\delta$ is the localised elastic deflection. It is apparent that the aspect ratio of the contact footprint would be high for spur gear teeth pair conjunctions. Therefore, for the case of an undeformed line contact, the profile of an equivalent solid near a semi-infinite elastic half-space becomes:

$$s\left(x\right)=\frac{x^2}{2R}$$

(16)

The radius of curvature $R$ of the equivalent solid is given by Equation (6). The general form of the localised elastic deflection is given by the elasticity potential equation as [14]:

$$\delta \left(x,y\right)=\frac{4}{\pi E^{*}}\int p\left(x^{\prime }\right)\ln \left(x-x^{\prime }\right)dx\prime$$

(17)

The variation of lubricant viscosity with pressure is given as [45]:

$$\eta = {\eta }_0\exp \left\{\ln \left(\frac{{\eta }_0}{{\eta }_r}\right)\left[{\left(1+\frac{p}{p_r}\right)}^Z-1\right]\right\}$$

(18)

where the ambient viscosity of the lubricant, ${\eta }_0$, is temperature-dependent and is found using the Vogel’s equation [46]:

$${\eta }_0=a_v\exp \left(\frac{b_v}{T-c_v}\right)$$

(19)

where $a_v$, $b_v$ and $c_v$ are lubricant-specific constants found through the measurements of viscosity at three different temperatures and through subsequent curve-fitting [14].

The constants ${\eta }_r$ and $p_r$ in Equation (18) are 6.31 × 10⁻⁵ and 1.9609 × 10⁸, respectively, and the piezo-viscosity index, Z, is obtained as:

$$Z=\frac{{\alpha }_0{p}_r}{\ln \left(\frac{{\eta }_0}{{\eta }_r}\right)}$$

(20)

where the pressure–viscosity coefficient, ${\alpha }_0$, is usually obtained through experimental measurements.

With the transmission gears of high-performance racing vehicles, as in the current study, high shear highly loaded contacts are subject to non-Newtonian traction. Therefore, although the Newtonian Reynolds equation is used here, it is essential to take into account the shear thinning of the lubricant as also shown in [39,47]. This is in order to provide a more accurate prediction of traction in the contact. Approaches to solve the modified Reynolds equation, including non-Newtonian effects, are described in other works such as [48,49,50]. Therefore, although the non-Newtonian Reynolds equation is not considered in this study, with the use a modified version of Havriliak and Negami [51] it is possible to accommodate the non-Newtonian lubricant behaviour by adjusting the viscosity of the lubricant as discussed in Paouris et al. [47]. This yields an effective viscosity as:

$${\eta }_{eff}=\frac{\eta }{F\left(\lambda \right)}$$

(21)

where $F\left(\lambda \right)$, the non-Newtonian function, and the shear rate $\dot{\gamma }$ are, respectively, given as:

$$F\left(\lambda \right)= {\left[1+{\left(\lambda \dot{\gamma }\right)}^{{\alpha }_{HN}}\right]}^{{\beta }_{HN}}$$

(22)

and

$$\dot{\gamma }=\frac{\mathsf{\Delta }U}{h_c}$$

(23)

where $\mathsf{\Delta }U=\left|u_p-u_w\right|$. Furthermore, $\lambda$,${\alpha }_{HN}$, and ${\beta }_{HN}$ are lubricant-specific parameters obtained experimentally [51]. The shear rate is calculated using the lubricant film thickness at the centre of the contact domain. Note that Sivayogan et al. [39] showed that under the current operating conditions some regions of the meshing cycle can be shear-independent. In these regions $F\left(\lambda \right)=1$, resulting in lubricant viscosity being solely dependent on the generated pressures.

The effect of pressure and temperature upon lubricant density is also included in the analysis as [52]:

$$\rho = {\rho }_0\left(1+\frac{0.6\times {10}^{-9}p}{1+1.7\times {10}^{-9}p}\right)\left[1-0.65\times {10}^{-3}\left({\theta }_e-{\theta }_0\right)\right]$$

(24)

where the density at ambient temperature and pressure is denoted by ${\rho }_0$. Effective lubricant contact and bulk temperatures are denoted by ${\theta }_e$ and ${\theta }_0$, respectively.

It is important to determine the lubricant temperature in the EHL contact as well as the flash temperature rise of the contacting surfaces. A more detailed method would be to solve the energy equation which provides information about the temperature distribution within the contact [42,53,54]. However, this approach would significantly increase the computation times. For the current analysis, it is essential to use the actual viscosity of the lubricant in the contact to better predict the generated traction. The viscosity in the contact is normally lower than that of the bulk lubricant viscosity measured at the sump due to the increased contact temperature as the result of generated friction. For this reason, a thermal network model similar to that originally presented by Morris et al. [55] and used for meshing teeth of spur gears in [39] is used. The thermal network model lacks in providing detailed information about the distribution of temperature in the contact, but it is a computationally powerful tool in providing effective viscosity values as experienced in the contact. The thermal network approach is presented in detail in [39,55], but a brief description is provided here for the sake of the completeness of the TEHL model. Morris et al. [55] noted that some of the friction-generated heat flows to the bounding contacting surfaces, and some is taken away by the flow of lubricant through the contact. In all cases, the flow of heat needs to overcome several thermal resistive barriers, as shown in Figure 3. Therefore, an analytical approach based on a thermal network model is suitable.

The rate of heat generation through contact friction is:

$$\dot{Q}=F_f\mathsf{\Delta }U$$

(25)

For highly loaded EHL contacts a combination of thin lubricant film and low lubricant flow rate through the contact means that the rate of convective heat transfer, ${\dot{Q}}_{cv}$, is almost negligible [14,56]. Most of the heat is conducted away through the bounding contacting solids (i.e., the meshing teeth). However, convective heat transfer through lubricant flow is retained for the sake of completeness of the analysis:

$${\dot{Q}}_{cv}=\dot{m}{c}_p\left({\theta }_e-{\theta }_0\right)$$

(26)

where ${\theta }_0$ is the inlet bulk temperature of the lubricant and ${\theta }_e$ is the effective lubricant contact temperature.

The inlet lubricant temperature, ${\theta }_0$, is calculated using solid body temperatures, with inlet heating constituting heat convection from the solid surfaces into the entrant lubricant at the contact inlet, thus:

$${\theta }_0=\frac{{\theta }_{s\mathrm{p}}{u}_{\mathrm{p}}+{\theta }_{s\mathrm{w}}{u}_{\mathrm{w}}}{u_{\mathrm{p}}+u_{\mathrm{w}}}$$

(27)

The lubricant mass flow rate in Equation (26) is calculated as [14,57]:

$$\dot{m}= -\frac{\mathsf{\rho }{h}^3}{12{\eta }_{eff}}\left(\frac{\partial p}{\partial x}\right)+ U\left(\frac{\mathsf{\rho }\mathrm{h}}{2}\right)$$

(28)

${\dot{Q}}_p$ and ${\dot{Q}}_w$ are the proportions of generated heat conducted away through the contacting surfaces. They are obtained as:

$${\dot{Q}}_i= \frac{{\theta }_e-{\theta }_{si}}{R_{si}}$$

(29)

where $i =\mathrm{p}, \mathrm{and} \mathrm{w}$ refers to the contacting surfaces.

It is clear that effective contact and surface temperatures are required for the TEHL. The resistances to heat flow in the thermal network model of Figure 3 are used to determine ${\theta }_e$ and ${\theta }_{si}$:

$$R_{si}=R_{l_i}+ R_{f_i}, i =\mathrm{p},\mathrm{w}$$

(30)

where $R_{l_i}$ and $R_{f_i}$ are the thermal resistances due to the lubricant film and the flash temperature rise of the solid surfaces. For heat generated at the centre of the contact and conveyed equally to the mating surfaces, the thermal resistance of the formed lubricant film thickness becomes [47,58]:

$$R_{l_i}=\frac{h_T}{2k_l\mathrm{A}}$$

(31)

and for resistance due to the flash temperature of the contacting surfaces [56]:

$$R_{f_i}=\frac{1.06S_{f_i}}{k_s\mathrm{A}}$$

(32)

where [47]:

$$S_{f_i}=\sqrt{\frac{2k_{s}a}{{\rho }_s{c}^{\prime }\mathsf{\Delta }U}}$$

(33)

The equivalent thermal resistance for lubricant mass flow rate is [55]:

$$R_e=\frac{1}{\dot{m}{c}_p}$$

(34)

The heat flow balance (generated and conveyed through convection and conduction) can be stated as:

$$\dot{Q}={\dot{Q}}_1+ {\dot{Q}}_2+{\dot{Q}}_{cv}$$

(35)

Replacing the rates of heat flow from the above relations in Equation (35) leads to the following relationship for the average effective contact temperature:

$${\theta }_e=\left[\dot{Q}+ \left(\frac{{\theta }_{s1}}{R_{s1}}+\frac{{\theta }_{s2}}{R_{s2}}+\frac{{\theta }_0}{R_e}\right)\right]{\left(\frac{1}{R_{s1}}+\frac{1}{R_{s2}}+\frac{1}{R_e}\right)}^{-1}$$

(36)

where the surface temperatures are obtained as:

$${\theta }_{si}={\theta }_{ini}+ \mathsf{\Delta }{\theta }_{si}$$

(37)

where ${\theta }_{ini}$ is the assumed initial surface temperature of the contacting surfaces at the start of a meshing cycle. For the spur gear pairs of the dry sump high-performance transmission system studied in the current paper, this initial surface temperature is at the environmental sump temperature of $40–50 °\mathrm{C}$. The surface temperature rise $\mathsf{\Delta }{\theta }_{si}$ is obtained as:

$$\mathsf{\Delta }{\theta }_{si}= \frac{S_{f_i}}{k_{s}A}{\dot{Q}}_i$$

(38)

The simultaneous solution of Equations (36)–(38) yields the effective contact temperature, ${\theta }_e$ and the surface temperatures. The contact temperature is used to update the lubricant rheological state, viscosity and density.

2.3. Method of Solution

The unknowns in the TEHL (thermo-elastohydrodynamics) analysis are pressure, $p$, dynamic viscosity, $\eta$, density,$\rho$, deflection, $\delta$, film thickness, $h$, and temperature, $\theta$. The number of unknowns equates to the number of equations at each instant of meshing. However, due to the non-linear nature of the relationships, it is necessary to use an iterative numerical technique to solve the problem at any instant of time within a meshing cycle. The Effective Influence Newton–Raphson (EIN) method, which is well-established for the solution of EHL problems, especially at high loads [39,59,60,61], is used within the following procedure:

• Parameters from LLTCA are input at the start of a meshing cycle;
• An initial guess is made for the lubricant film thickness at the centre of the contact;
• The computational domain is set with an inlet length of $12.42a$ and contact exit position of $4.42a$ where $a$ is the contact semi-half-width (Equation (4)). The number of elements used in the direction of lubricant entrainment, x, is 2051;
• Iterative pressure residuals are found using EIN iterations and the pressures are updated using the recursive expression:

  $$p^n=p^{n-1}+\mathsf{\Omega }\mathsf{\Delta }{p}^n$$

  (39)

  where n denotes the iteration step and $\mathsf{\Omega }$ is the under-relaxation factor, typically: ${10}^{-2}$-${10}^{-1}$;
• The pressure convergence criterion used is:

  $${\displaystyle \sum }_i{\displaystyle \sum }_j\left|\frac{p_{i,j}^n-p_{i,j}^{n-1}}{p_{i,j}^n}\right|\le {10}^{-3}$$

  (40)
• The lubricant reaction is obtained as:

  $$w_p=\int pdx$$

  (41)
• The instantaneous equilibrium condition is sought using the load convergence criterion:

  $$\left|\frac{w-w_p}{w}\right|\le {10}^{-3}$$

  (42)
• If the stated equilibrium condition is not satisfied, the film thickness is updated through the modification of the undeformed gap as:

  $$h_0=h_0{\left(\frac{w_p}{w}\right)}^{\varsigma }$$

  (43)

  where $\varsigma$ is the damping (load relaxation) factor. In the current analysis, the value of $\varsigma ={10}^{-3}$ is used. Subsequently, density and viscosity are updated using the converged pressures. Then, steps 4 to 8 are repeated until step 7 is satisfied;
• Once the film thickness is determined, the thermal network model is used to obtain the temperature of the lubricant as well as the flash temperatures of the contacting surfaces;
• The lubricant temperature is used to adjust the lubricant density and dynamic viscosity.

2.4. Sub-Surface Stress Field

The combination of generated contact pressures and surface traction in meshing gear teeth pairs induces sub-surface stresses which are commonly responsible for inelastic deformation, as noted in Section 1. For the case of the concentrated counter-formal contact of ellipsoidal solids of revolution, as in the instantaneous contact of involute spur gear teeth, the classical Hertzian theory may be used to determine the sub-surface stresses when the contacting solids can be considered as smooth frictionless semi-infinite elastic half-spaces [14,29]. The Hertzian assumption regarding friction may be relaxed to include the effect of contact traction [29]. However, the semi-infinite and dry elastostatic nature of the contact remains crucial for the use of Hertzian theory. The result of the classical approach fails to comply with certain practical situations, such as lubricated EHL pressure distribution, and has contact conditions which abrogate the semi-infinite assumption, such as increased contact conformity or the presence of a thin layer of coating [62]. A lot of gear meshing conditions fall into these exceptions as they are normally lubricated, increasingly use hard, wear-resistant coatings, or have a higher degree of conformity, such as circular arc, bevel and Novikov gears. Therefore, a more generic solution is preferred. In the case of involute spur gears of high-performance transmissions, lubricated TEHL contacts occur.

For the general case of 2D sub-surface stresses in the bulk of a semi-infinite elastic solid, a parabolic pressure distribution may be assumed which closely approximates the 2D elliptical Hertzian pressure profile, as noted by many authors such as Love [63] and Muskhelishvili [64].

Johns-Rahnejat et al. [34] provided an analytical generic solution for sub-surface stresses generated in an elastic solid under a parabolic pressure distribution. However, EHL contact pressure distribution does not follow a simple analytical expression. Johns-Rahnejat et al. [34] spread the arbitrary pressure distribution over the contacting face of an elastic solid of revolution approximated by a parabola, as shown in Figure 4. Here, the solid of revolution is a long roller, representing an instantaneous geometry of a gear tooth during a meshing cycle. The generated EHL pressure obtained in Section 2.3 is approximated by a parabolic distribution. This is also augmented by tangential traction due to generated friction (Equation (1)). Provided that the contact undergoes shear loading, a random traction distribution will also act upon the elastic solid. The contact conjunction can be discretized into multiple small elements. The sub-surface stresses due to normal pressure and tangential traction acting on strip elements can be evaluated at any position within the bulk of the solid positions: $\left(x,z\right)$.

The instantaneous approximate parabolic pressure distribution can be stated as [34]:

$$p=p_m\left[1-{\left(\frac{\beta }{{\beta }_m}\right)}^2\right]$$

(44)

where $\beta$ is any location within the arc of contact, measured from the position of the maximum pressure, $p_m$, and $2{\beta }_m$ is the arc length of the parabolic pressure distribution.

A pressure element $p$ (acting in the z-direction) and tangential traction $q$, acting along the x-direction on the surface of the roller, induce sub-surface stresses at any arbitrary position $\left(X,Z\right)$. in the bulk of the elastic solid:

$${\sigma }_x=-\frac{2p}{\pi }\left[\frac{\left(\frac{D}{2}-z\right)x^2}{{\left[x^2+{\left(\frac{D}{2}-z\right)}^2\right]}^2}+\frac{\left(\frac{D}{2}+z\right)x^2}{{\left[x^2+{\left(\frac{D}{2}+z\right)}^2\right]}^2}-\frac{1}{D}\right]-\frac{2q}{\pi }\left[\frac{x^3}{{\left[x^2+{\left(\frac{D}{2}-z\right)}^2\right]}^2}+\frac{x^3}{{\left[x^2+{\left(\frac{D}{2}+z\right)}^2\right]}^2}-\frac{\sqrt{\frac{D}{z}-1}}{D}\right]$$

(45)

$${\sigma }_z=-\frac{2p}{\pi }\left[\frac{{\left(\frac{D}{2}-z\right)}^3}{{\left[x^2+{\left(\frac{D}{2}-z\right)}^2\right]}^2}+\frac{{\left(\frac{D}{2}+z\right)}^3}{{\left[x^2+{\left(\frac{D}{2}+z\right)}^2\right]}^2}-\frac{1}{D}\right]-\frac{2q}{\pi }\left[\frac{{\left(\frac{D}{2}-z\right)}^{2}x}{{\left[x^2+{\left(\frac{D}{2}-z\right)}^2\right]}^2}+\frac{{\left(\frac{D}{2}+z\right)}^{2}x}{{\left[x^2+{\left(\frac{D}{2}+z\right)}^2\right]}^2}-\frac{\sqrt{\frac{D}{z}-1}}{D}\right]$$

(46)

$${\tau }_{xz}=+\frac{2p}{\pi }\left[\frac{{\left(\frac{D}{2}-z\right)}^{2}x}{{\left[x^2+{\left(\frac{D}{2}-z\right)}^2\right]}^2}-\frac{{\left(\frac{D}{2}+z\right)}^{2}x}{{\left[x^2+{\left(\frac{D}{2}+z\right)}^2\right]}^2}\right]+\frac{2q}{\pi }\left[\frac{\left(\frac{D}{2}-z\right)x^2}{{\left[x^2+{\left(\frac{D}{2}-z\right)}^2\right]}^2}-\frac{\left(\frac{D}{2}+z\right)x^2}{{\left[x^2+{\left(\frac{D}{2}+z\right)}^2\right]}^2}\right]$$

(47)

where D is the instantaneous radius of the elastic solid. In the current analysis, it is the diameter of the equivalent long elastic cylinder (see Equation (6)): $D=2R$.

For any arbitrary pressure distribution and tangential traction, the contact region can be discretized into multiple $\mathrm{d}\beta$ strips. The sub-surface stresses can then be evaluated for each discretized element and the overall effect of all such strips upon the sub-surface stress field is determined. As there are slight variations in the angle $\beta$ for each elemental strip $\mathrm{d}\beta$, it is more convenient to express the sub-surface stresses in the general coordinates $X,Z$. The transformations between the pressure co-ordinates $x,z$ and the general sub-surface coordinates $X,Z$ are:

$$z=Z\cos \beta +X\sin \beta , x=-Z\sin \beta +X\cos \beta$$

(48)

In addition to the coordinate system transformation, the sub-surface stresses ${\sigma }_x$, ${\sigma }_z$ and ${\tau }_{xz}$ should be expressed in the general plane of stresses in terms of ${\sigma }_X$, ${\sigma }_Z$ and ${\tau }_{XZ}$. These transformations are:

$${\sigma }_X={\sigma }_z{\sin }^2\beta +{\sigma }_x{\cos }^2\beta +{\tau }_{xz}\sin 2\beta$$

(49)

$${\sigma }_Z={\sigma }_z{\cos }^2\beta +{\sigma }_x{\sin }^2\beta -{\tau }_{xz}\sin 2\beta$$

(50)

$${\tau }_{XZ}=\frac{1}{2}\left({\sigma }_z-{\sigma }_x\right)\sin 2\beta +{\tau }_{xz}\cos 2\beta$$

(51)

and:

$$p=p\left(\beta \right), \mathrm{d}p=\left(\frac{D}{2}\right)p\mathrm{d}\beta$$

(52)

$$q=q\left(\beta \right), \mathrm{d}q=\left(\frac{D}{2}\right)q\mathrm{d}\beta$$

(53)

Implementing the above transformation from the surface co-ordinates (pressure and traction) for each element $d\beta$, the corresponding sub-surface stresses are obtained as shown in Equations (A1) to (A3) in the Appendix A.

In the current analysis, the average friction obtained through the use of Equations (1), (7) and (8) is distributed equally to all elemental strips $d\beta$. However, Equations (A1)–(A3) are generic and any tractive distribution can be applied. The overall sub-surface stresses at any location of $X,Z$ within the bulk of the elastic contacting solid are obtained by the summation of the effect of the entire elemental pressure-tractive elements $d\beta$ as:

$${\sigma }_X={\displaystyle \sum }_{\beta }d{\sigma }_X, {\sigma }_Z={\displaystyle \sum }_{\beta }d{\sigma }_Z, {\tau }_{XZ}={\displaystyle \sum }_{\beta }d{\tau }_{XZ}$$

(54)

3. Results and Discussion

Table 1. Mechanical and thermal properties of the surfaces and lubricant.

Parameter	Value	Unit
Modulus of elasticity of gear material, $E$	206	GPa
Poison ratio of solid, $\nu$	0.3	-
Density of the solid, ${\rho }^{\prime }$	7800	kg/m3
Thermal conductivity of solid, $k_s$	46.7	W/m·K
Specific heat capacity of the solid, $c^{\prime }$	470	J/kg·K
Dynamic viscosity of lubricant at 40 °C, ${\eta }_0$	0.03034	Pa.s
Pressure–viscosity coefficient at 40 °C, ${\alpha }_0$	1.67 × 10−8	1/Pa
Thermal conductivity of lubricant, $k_l$	0.137	W/m·K
Specific heat capacity of lubricant, $c_p$	1670	J/kg·K
Characteristic shear stress, ${\tau }_0$	2	MPa
Havriliak–Negami non-Newtonian model parameters:
${\alpha }_{HN}$	0.7	-
${\beta }_{HN}$	1	-
${\lambda }_{HN}$	7.9 × 10−8	s

lists the parameters used in the current analysis. The involute spur gear pair considered comprise the final drive in high-performance racing transmission. The gears were made of steel with super-finished teeth and a composite surface average roughness of the order of 0.1 µm. Meshing teeth pair contacts were subjected to high loads and shear, resulting in the thermo-elastohydrodynamic regime of lubrication (TEHL), as discussed below. The thermo-mechanical data for the gears are provided in Table 1. Low viscosity lubricant was used to reduce frictional power loss. The lubricant rheological data are also provided in Table 1.

A number of teeth pairs are in simultaneous mesh at any instant of time within a meshing cycle. For any given teeth pair contact, the path of action is divided into 37 points. LLTCA was carried out with elastohydrodynamic analysis and tractive action at each of these points, leading to generated contact temperature. Figure 5 shows the variations in the governing operating parameters of contact load, sliding velocity, shear stress and generated contact temperature in a meshing cycle. The contact load is the share of the overall load per meshing teeth pair, which increases from the initial engagement to the region of high contact load and subsequently drops as separation commences. The sliding velocity decreases with progressive engagement until the pitch point, where the meshing teeth are in an instantaneous pure rolling condition. Here, the relative sliding velocity of meshing surfaces momentarily diminishes. Thereafter, with the gradual separation of a mating pair, the sliding velocity increases until complete disengagement. Shear stress is directly proportional to the sliding velocity under Newtonian shear. However, this is not the case here, as can be seen in Figure 5b. The lubricant remains in non-Newtonian traction, with the shear stress exceeding its characteristic shear stress, marking the limit of Newtonian shear (see Figure 5b). The limiting shear stress of lubricants at atmospheric conditions was extensively measured by Jacobson [65], who found this to be in the range of 1 MPa -5 MPa. In the current study, a value of 2 MPa was chosen, in line with the data used by Elisaus et al. [9]. It should be noted that, in gear teeth contacts, the shear stress is mainly dominated by the contact pressure effect. Therefore, the effect of characteristic limiting shear stress at atmospheric conditions is rather trivial [66]. The constant value indicated for the characteristic shear stress corresponds to that under ambient conditions. In fact, its value varies with pressure and temperature. Therefore, it is shown in the figure merely to demonstrate the dominant non-Newtonian tractive behaviour in the case studied here. As expected, the generated temperature followed the shear and the sliding velocity variations. Its diverging correspondence with the contact load variation (Figure 5a) indicates the dominance of viscous shear heating over any compressive heating, which is only significant during the gradual engagement and separation of meshing surfaces.

Since the gears considered in this study are for a high-performance application, the typical composite surface roughness values are often well below 0.1 μm. Thus, with the film thickness values encountered in this study, the probability of asperity contact is expected to be negligible. Figure 5c shows the variations in central film thickness for a meshing cycle in the form of the Stribeck film ratio ($\lambda =h_c/\sigma$). As can be seen from this figure, the central film thickness remains generally well above the threshold for the full film regime of lubrication. Therefore, it is expected that the effect of asperity contact on the predicted results would be insignificant. At relatively lower sliding speeds, however, the effect of asperity contact can become significant and should be considered [54].

The variation in temperature in Figure 5a follows the shear stress and sliding speed as shown in Figure 5b, depending on the lubricant response to these parameters. For instance, in the region between B and C in Figure 5a, both shear stress and sliding speed decrease, leading to a reduction in temperature. In the region C to D, the shear stress continues its downward shallow slope, whilst the sliding velocity rises. Therefore, a rise in temperature is observed. The variation in lubricant temperature in the regions A to B and from D onwards is also affected by shear stress and the sliding velocity. Thus, the fluctuations in the generated temperature follow the stated trends.

Four distinct locations for teeth pair contact are highlighted in

Table 2. Operating conditions at key points in a meshing cycle.

Location, Grid Point	Input Parameters			Output Parameters
	Sliding Velocity (m/s)	Contact Temperature (°C)	Equivalent Radius of Curvature (mm)	Max. Shear Stress (MPa)	Max. Pressure (Gpa)	EHL Pressure Spike (Gpa)
A: Initial approach	10.699	107.8	0.608	6.94	0.277	--
B: Max. Hertz pressure	6.329	147.4	4.713	34.77	2.872	1.556
C: Pitch point	1.896	56.4	7.599	27.98	2.271	1.834
D: Onset of separation	10.093	70.0	8.156	13.22	1.550	1.550

. These are: (i) the point of approach (point A: initial engagement of a teeth pair); (ii) the point of maximum generated contact pressure (point B); (iii) the pitch point (point C) where there is the pure rolling of mating pairs; and (iv) the point of contact separation (point D: disengagement of a meshing pair).

Figure 6a–d show the generated TEHL pressure distribution, the corresponding film thickness profile, viscous shear stress and sub-surface orthogonal reversing shear stress distribution for the parts of a typical meshing cycle described by points A–D in Table 2. In a, the generated pressures are below 0.5 GPa with the film thickness shape corresponding to undeformed hydrodynamic conditions. In fact, Sivayogan et al. [39], studying the same pair of involute gears operating under similar conditions and using the Greenwood chart [67], showed that at the onset of engagement the conditions pertained to iso-viscous-elastic conditions, sometimes referred to as soft EHL. The average surface viscous shear stress was far in excess of its characteristic value of 2 Mpa shown in Figure 5b. The alternating orthogonal sub-surface shear stress contours show that the material at the leading edge of the contact (inlet) was under compression (${\tau }_{zx}<0$) whilst the trailing edge (outlet) was in tension (${\tau }_{zx}>0$). In the rolling and sliding contact of a meshing teeth pair, represented by an instantaneous equivalent roller contacting an elastic half-space, the generated sub-surface orthogonal shear stresses alternate between the leading and trailing edges of the contact. Although their magnitude is less than the maximum shear stress, their cyclic nature and double amplitude variation of approximately 0.5 of the maximum EHL primary pressure peak (same as the maximum Hertzian pressure) is considered to be the determining factor in any sub-surface inelastic deformation. Their symmetric cyclicity is more clearly observable under hard EHL conditions, as can be seen from Figure 6b–d. This is not the case for Figure 6a because of rather low pressures. In other cases (b–d), any asymmetry in the sub-surface stress field is because of the secondary pressure peak (EHL pressure spike) which induces a sub-surface stress field of its own. These localized sub-surface stresses are shown in the zoomed insets of Figure 6b–d. The asymmetry is more pronounced with a high-pressure spike, as in the case of Figure 6d.

The presence of a significant pressure spike usually shifts the region of maximum shear stresses towards the contact surface, which in the presence of any surface anomaly such as inclusions or pores can induce micro-pitting. This is often observed in high-performance transmissions. The pressure spike occurs essentially as the result of rapid change in lubricant viscosity. In the high-pressure (Hertzian) region of the contact area the lubricant acts similarly to an amorphous solid with high viscosity. As it is drawn towards the contact exit the viscosity is dramatically reduced. This is because the prevailing conditions there tend to be atmospheric. This rapid change causes a pressure spike. With thermal conditions the lubricant viscosity within the contact is somewhat attenuated due to heat. Therefore, under TEHL, the viscosity effect is not as dramatic as that under isothermal conditions (reduced pressure spike). This is shown, for example, for the case of hypoid gear transmissions under non-Newtonian conditions in [42]. The pressure spike increases in magnitude with the increased speed of the entraining motion of the lubricant, $U$, into the contact and decreases with increased contact load $w$. This has been noted experimentally by Gohar and Safa [57] for the case of line contact EHL and by Johns-Rahnejat and Gohar [68] for the case of point contact. This is evident in the results of Figure 6, as the contact load is higher for cases (b) and (c) with lower sliding velocity. In these cases, the pressure spike has diminished. This is partly due to progressive starvation at the inlet meniscus. In the usual numerical method of solution for EHL contacts, any generated negative pressures are set to zero, i.e., $p<0\to p=0$. This is commonplace in all such solutions, resulting in the inlet meniscus moving inwards from the idealized initially assumed: $p\to 0 at x\to -\infty$ towards the edge of the high-pressure Hertzian region. This constitutes the progressive starvation of the contact. In fact, when $p=0 at x=-a$ (the edge of Hertzian rectangular strip), the contact would nominally be under a dry Hertzian condition, with no pressure spike at the trailing edge of the contact (exit). It should be noted that these conditions occur as the result of an instantaneous equilibrium in the iterative procedure, not due to any imposed inlet boundary condition. There are solutions noted for imposed inlet boundary starvation, such as that originally proposed by Hamrock and Dowson [69].

Figure 7 shows the ratio of pressure spike in the vicinity of contact outlet to the maximum primary pressure peak ($p_s/p_0$) for a typical meshing cycle. It is interesting to note that the ratio remains less than unity for the main part of the meshing cycle where the load is higher (Figure 5a). Under these conditions the pressure spike is reduced, particularly with any inlet starvation, which is the case here. At the onset of initial contact, the conditions are hydrodynamic as already noted; therefore, there is no pressure spike. As the separation of a teeth pair commences at the end of a meshing cycle, the increasing sliding speed and reducing load results in flooded inlet conditions and the ratio of $p_s/p_0$ exceeds unity. The pressure spike can be significantly larger than the primary pressure (maximum Hertzian pressure). In fact, the results under EHL impact conditions show the ratio well in excess of unity, as recorded experimentally by Safa and Gohar [70] and verified theoretically by Al-Samieh and Rahnejat [71]. Using a thin-film miniature pressure transducer, Johns-Rahnejat [33] resolved the rapidly transitory pressure spike under various contact operating loads and speed conditions and showed that the pressure spike increases with speed and diminishes with load for flooded inlets. Starved inlets lead to the significant diminution of the pressure spike. The starvation is signified by the inlet distance to the edge of the high-pressure (Hertzian) contact. This is given by the ratio $\overline{x_i}=x_i/a$. Tipei [63] equated the potential flow gradient (without a surface tension effect) with the Poiseuille flow at the inlet to determine the point at which the flow is inbound towards the Hertzian contact region. This equality of potential and Poiseuille flow gradients leads to the film ratio of $h_i/h_c =9$, which leads to the following expression [63]:

$$\frac{{\overline{{\mathrm{x}}_{\mathrm{i}}}}^2-1}{{\overline{{\mathrm{x}}_{\mathrm{c}}}}^2-1}=9$$

(55)

where $x_c$ is the position in the inlet region prior to the edge of the Hertzian zone where the film thickness becomes parallel. This depends on the load and speed of the entraining motion, which alter instantaneously during a meshing cycle. Tipei [72] showed that: $\overline{x_i}\approx -2$ for the onset of starvation. Figure 8 shows that this is the case for most of the meshing cycle, and thus the diminution of the pressure spike (i.e., $p_s/p_0 <1$ in Figure 7). In fact, a more detailed analysis (e.g., including the effect of surface tension, etc.) including observations of the inlet region was carried out by Birkhoff and Hays [73] who showed that the film ratio at the onset of starvation is actually $h_i/h_c =11.298$ for the cylindrical rolling condition (similar to instantaneous conditions in the current analysis). Therefore, the conditions are predominantly starved.

4. Concluding Remarks

This paper presents a detailed contact analysis of a high-performance spur gear pair, such as those used for the final drive in motorsport racing applications. Meshing teeth pair contacts are subject to high orthogonal and shear loading. Under these conditions, the meshing teeth pairs undergo a non-Newtonian thermo-elastohydrodynamic regime of lubrication with starved inlet boundaries. These conditions generate high contact pressures with a diminished EHL pressure spike at the contact outlet. Significant sub-surface reversing orthogonal shear stresses occur within the bulk of the contacting solids, which can cause inelastic deformation of the surfaces if they coincide with any pores, pits, inclusions or other flaws. This problem is exacerbated by the pressure spike at the outlet of the EHL conjunction, which induces a sub-surface stress field of its own as a near-surface effect. However, the results clearly show that the potency of the pressure spike is somewhat palliated with starved inlets, as shown by the representative analysis here. This merely implies that with the increasing trend towards use of dry sumps in transmissions [74] in order to conserve hydrocarbon resources and to palliate emissions, the diminution of the pressure spike would improve sub-surface structural integrity. However, it should be noted that starved inlets can also result in the diminution of lubricant film thickness, direct boundary interactions and potential wear. Therefore, the problem of gear lubrication is multi-variate, and an optimum solution may not be easily attained. It is also noteworthy that another role of lubrication is to cool the contacting solids, which can otherwise be subjected to thermo-structural instability. Nevertheless, multi-physics analyses, such as the one reported here, provide detailed findings to aid the development of high-performance modern transmissions. Further extensions of this work can include the investigation of coated gear surfaces [75,76], particularly with application in wind turbines, where residual stresses can affect the contact fatigue life [77]. An improvement to the current expounded methodology could be achieved by combining gear pair inertial dynamics with the tribology of meshing teeth pairs with the inclusion of transient effects, including the squeeze film effect.