Issue 
Mechanics & Industry
Volume 20, Number 6, 2019



Article Number  626  
Number of page(s)  16  
DOI  https://doi.org/10.1051/meca/2019064  
Published online  02 December 2019 
Regular Article
From Hertzian contact to spur gears: analyses of stresses and rolling contact fatigue
^{1}
Univ. Lyon, INSALyon, CNRS UMR5259, LaMCoS, 20 Avenue Albert Einstein, 69621 Villeurbanne, France
^{2}
Univ. Lyon, ECAM Lyon, INSA Lyon, LabECAM, 40 Montée St Barthélémy, Lyon, France
^{3}
Univ. Lyon, INSA Lyon, MATEIS UMR CNRS 5510, Bât. Blaise Pascal, 7 Avenue Jean Capelle, Villeurbanne, France
^{4}
SAFRAN Transmission Systems, 18 Boulevard Louis Seguin, Colombes, France
^{*} email: fabrice.ville@insalyon.fr
Received:
20
June
2019
Accepted:
10
September
2019
The study of rolling contact fatigue in spur gears requires a good comprehension of all the phenomena occurring at the material scale. On a numerical point of view, a realistic representation of the material and of the load repartition function of the local microgeometries is needed. However the resulting models are often complex and timeconsuming. So, this work aims at developing a model meeting these specificities. Thus, different sections of the spur gear material granular geometry are simulated first. Secondly, the contact pressure fields are computed accurately relatively to the simulated surface microgeometry. Then, the influence of several parameters on their rolling contact fatigue life is highlighted. Among friction, sliding coefficient, load variation and roughness, these individual or combined parameters are taken into account in the model, tested and their impact stressed out. Finally, a fatigue criteria based on rolling contact fatigue microcracks nucleation at grain boundaries is proposed in order to compare simulations and influencing parameters to the reference.
Key words: Numerical study / rolling contact fatigue / spur gear / rough microgeometry
© G. Vouaillat et al., published by EDP Sciences 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Spur gears are major components of aeronautical transmission boxes. The correct functioning of the aircraft systems is directly dependent of its transmission operational state. However, conception or manufacturing approximations of these parts may lead to component failures, micropitting being the most recurrent one. In the literature, micropitting occurs for rollingsliding contacts and is characterized by typical 10 μm length and depth pits [1–5].
Besides, OILA has pointed out in [6] through twindiscbased experimental investigations, seven factors influencing this specific failure. Among them, pressure, roughness, friction and slidetoroll ratio are predominant. More recently, RYCERZ [7] has also identified slidetoroll ratio as a major influencer in micropitting occurrence and particularly on the rough microcontact stress cycles. This asperity stress history is highlighted as a potential basis for a micropitting criterion. Thus, the present study actually focuses on the influencing factors introduced above first, and on a shear stressesbased micropitting criterion then.
In all the previously quoted studies, micropitting is more generally described as the combined result of a lubrication film decrease and an overstressed material surface from a rough surface finish. Since the lubrication is ignored in the present paper to limit the influencing parameters, an accurate description of both the contact microgeometry and the material microstructure is needed so as to correctly simulate the failure phenomena.
Regarding the material, SADEGHI's working team focuses on the granular description of some specific bearing or gear steels in [8–10]. The austenitic grain boundaries from typical bearing or gear steels are represented. The same concern for the microstructure representation is brought in this study.
Concerning the contact, MORALES has reviewed in [11] the different methods in literature to take into account rough surfaces in computations. His own contribution was based on GREENWOOD and JOHNSON's work [12] and consists of a semianalytical method to process complex roughness profiles. Through a Fourier analysis of these profiles coupled to a transient solution, MORALES proposed in [13–16] a pressure field calculation for complex rough contacting surfaces. Applied to rolling element bearings or gears, these computed pressure fields are timevarying with the contact moving along the specific geometries. Such a model is used hereafter so as to completely simulate rough contacting gears. Thus, rolling contact fatigue estimations are made possible by criteria or stress analyses as presented in [17].
So, the final aim of this study is to simulate and understand the individual and combined impact of several influencing factors (among pressure, roughness, friction and slideroll ratio) on the nucleation of micropitting for contacting spur gears. The literaturebased model is presented in details for a smooth contact in a first approach. The main contribution consisting of the material microstructure introduction coupled to the timeevolving pressure field computation representing industrial gear specificities is notably introduced. Then, the rough contact is considered in regards with the influencing factors so as to finally simulate a real contacting spur gear.
2 Model description
2.1 Introduction
A numerical approach is built to give an answer to the real influence of different parameters on the life duration of gear surface. The model is divided into three parts identified in Figure 1. First, all the input parameters to be used inside the analytical analysis are given. The pressure field computation is then performed function of both the macrogeometry and microone of contacting bodies. These fields are exploited as boundary conditions of the finite element analysis. Finally, stresses and fatigue criterion analyses are operated through a model initially developed by NOYEL [18,19]. The details of the model construction and computations are explained hereafter.
Fig. 1 Organigram of the model. 
2.2 Basics of the model
A granular geometry is set as presented by NOYEL [19] so as to represent roughly typical austenitic grain boundaries observed in basic gear materials. This particular geometry is presented in Figure 2. Grains are generated numerically in this model (grain orientation is randomly set) and their size (25 μm) is chosen to fit with usual industrial gear core grain length [20]. The equivalent ASTM grain size [21] is almost G_{ASTM} = 7.
The geometry is meshed classically to get the Finite Element Model (FEM). Grain Boundaries (GB) are also meshed with cohesive elements to allow the stress computation along these GB. The plane strain hypothesis is applied for the FEM. The model is bounded on three faces by frictionless supports.
Fig. 2 Illustration of the static Hertzian pressure field of simulation #0 and the dimensionless Tresca stress field associated. 
2.3 Study of a centered and static pressure field
When a basic contact between two surfaces is considered without any roughness, microgeometry nor tangential effect, HERTZ's theory is applied [22]. Thus, equations (1) and (2) are used for a dry line contact.$${p}_{0}=\frac{2\cdot {F}_{N}}{\pi \cdot l\cdot a}$$(1) $$a=\sqrt{\frac{8\cdot {F}_{N}\cdot R}{\pi \cdot l\cdot {E}^{\text{'}}}}.$$(2)
The normal pressure field is then computed along the surface as follows in equation (3) and applied above the meshed geometry in Figure 2.$${P}_{N,\mathrm{t}\mathrm{o}\mathrm{t}}(X)=\sqrt{1{(X)}^{2}}.$$(3)
The Tresca equivalent shear stress for this plane strain issue is obtained in Figure 2 after a Finite Element Analysis (FEA) with the numerical data from simulation #0 as presented in Table 1. The equivalent stress is dimensionless (by the value of the Hertzian pressure p_{0}). Besides, the literature classical 0.3 value for the dimensionless Tresca stress is found for the line contact simulated in #0 (as shown on Tresca stress field colorbar in Fig. 2).
From these static and centered pressure and stress fields, only basic fatigue criteria can be applied. For example, a comparison with the micro yield shear stress with a plastic strain of 20μdef from LAMAGNERE [23] (915 MPa for a 100Cr6 bearing steel) can be commonly used. Thus, Tresca maximum shear stress is computed and compared to the yield stress limit in the whole material. A higher value corresponds physically to a dislocation creation which is linked to a plastic strain.
Table of simulations.
3 Study of a smooth contact
Whether a more precise description of the damaging phenomena occurring inside the material is wanted, a computation of the load passing over the surface is needed. So, the previously mentioned model inspired from NOYEL is used. In a first approach a smooth contact is considered.
3.1 Smooth contact passing over a surface
For a smooth Hertzian contact, the pressure field at each step of time is again computed thanks to equation (3). But comparatively to the previous simulation, the pressure field is applied at different successive positions above the meshed geometry to simulate the load passage. At this step, no friction nor sliding is taken into account in the model.
Then along each grain boundary (GB), the mean shear stress (or Intergranular Shear Stress ISS) is computed. A GB is composed of several cohesive elements. Thus, the ISS is estimated as the average value of the shear stress for each cohesive elements along this GB. This way, a complete simulation of the ISS evolution with time is made possible for every GB of the granular geometry. All the computations details are depicted in [19].
An amplitude ΔΤ of the ISS variations is then calculated. This computation is the basis of the fatigue criterion estimation. This criterion is obtained for each GB, the most severe one is isolated and its value becomes the general criterion for the whole simulation. ΔΤ is estimated with the minimummaximum method i.e. the detection of the two extrema of the shear stress. So, the GB presenting the largest value of ΔΤ is isolated and its position and evaluation of ISS are given as results.
ΔΤ_{1} for simulation #1 is given in Figure 3. It corresponds to a smooth contact without any surface shear stress nor sliding effect. All the simulations results will be presented in this work relatively to the reference simulation #1.
The fatigue criterion is based on experimental results found in the literature [24] and represented as a S–N curve in Figure 4. Modeling this S–N curve by a straight line gives directly a relation between the shear stress amplitude and the number of cycles to fatigue in (4).$$\mathrm{log}\left(\Delta \tau \right)=A\cdot \mathrm{log}\left({n}_{i}\right)+\mathrm{log}\left(B\right).$$(4)
These constants are expressed function of the variables and the material properties as given in (5) [24].$$\{\begin{array}{l}A=m\hfill \\ B={\sigma}_{r}\cdot {\left(1+m\right)}^{m}\hfill \end{array}.$$(5)
Then, formulation (6) comes directly from calculations so as to isolate the number of cycles to fatigue n_{i}. The dimensionless version is given in (7).$${n}_{i}={\left[\left(m+1\right)\cdot {\left(\Delta {\tau}_{i}/{\sigma}_{r}\right)}^{m}\right]}^{1}$$(6) $${N}_{i}=\frac{{n}_{i}}{{n}_{\mathrm{r}\mathrm{e}\mathrm{f}}}={\left(\frac{\Delta {\tau}_{i}}{\Delta {\tau}_{\mathrm{r}\mathrm{e}\mathrm{f}}}\right)}^{m}={\left(\Delta T\right)}^{m}.$$(7)
All the material parameters (m and σ_{r}) are found in RAJE's work [1]. With an empiric point of view, this GB corresponds to the location of the first microcrack in the material.
Two other simulations have been conducted to validate the correct evolution of the fatigue criterion through the variations of the Hertzian pressure p_{0}. Thus, simulation #1′ has a maximum pressure P_{max} equals to 0.84 (lower than #1). Simulation #1′ has a P_{max} of 3.63 (higher than #1). It also appears that those values have been chosen as the maximum (for 3.63) and minimum (for 0.84) values of P_{max} that will be used in all the other simulations. So, #1′ and #1′ are defined as the extreme boundaries in ΔΤ and N for the simulations with the same input parameters i.e. from #0 to #10 (cf. Tab. 1).
Consequently to the P_{max} diminution in #1′ [respectively augmentation in #1′], ΔΤ_{1'} computed in Figure 3 is smaller than ΔΤ_{1} [and ΔΤ_{1′ } is larger than ΔΤ_{1}]. This results in an increase for #1′ [respectively a decrease for #1′] in life duration to fatigue initiation. The microcracks initiate each time at the Hertzian depth in the material (cf. Tab. 1 and Fig. 7). So concerning the maximum pressure evolution, the fatigue criterion evolution is consistent with the expectations.
Fig. 3 Influence of the contact pressure and the traction coefficient variations for a smooth contact on ISS, ΔΤ and N results. 
Fig. 7 Illustration of microcracks nucleation location for smooth contacts. 
3.2 Insertion of a friction coefficient
The insertion of a tangential component influence is investigated hereafter. Sliding and friction are two main parameters to take into account when a complete simulation of the contacting gears is needed.
For a smooth contact and among friction and sliding, only friction will impact the fatigue results. Sliding will actually be of no influence on the displacement of local microgeometries above the two contacting surfaces since none is present. Simulation #3 presents the lone influence of friction on a smooth contact. The details of those simulations parameters are shown in Table 1.
The surface shear stress is computed as the product of the normal one for each step and a traction coefficient μ. A maximum value of 0.1 is chosen for the traction coefficient. Such a value is high for a lubricated contact and rather corresponds to a typical boundary regime. It is also assumed that no microslip is present during sliding and only full slip occurs. A positive or a negative sign is attributed to the traction coefficient depending on the sliding direction. This formulation is detailed in equation (8) $${P}_{T,\mathrm{t}\mathrm{o}\mathrm{t}}(X)=\pm \mu \cdot {P}_{N,\mathrm{t}\mathrm{o}\mathrm{t}}(X).$$(8)
The ISS and number of cycles to microcrack nucleation are shown in Figure 3 and compared to the reference one. The most stressed grain boundary (GB) is located both for simulations #1 and #3 in the depth of the material (Z = 0.5 in Fig. 7 and Tab. 1).
In terms of results, the comparison of the curves of #1 and #3 gives a really limited influence of the surface shear stress on the value of ISS: the global value of ΔT is almost not modified. The slight left shift of the ISS of #3 is linked to the position of its damaged GB. This one is located on the left hand side compared to the damaged GB of #1. Thus, GB from #3 is stressed earlier in time than GB from #1. Moreover, the slight shift to the bottom is linked directly to the insertion of the surface shear stress. This component introduces a constant shear contribution to the global stress field.
Finally, the number of cycles to fatigue initiation is almost not modified. So, no major change can be observed from the insertion of a surface shear stress with a constant traction coefficient on a smooth contact.
3.3 Load variation above the surface
The introduction of the theory of spur gears implies a line contact between the two contacting surfaces. Since gear macrogeometry is realized from a circle involute, the equivalent curvature radius will evolve along the gear action line. The same way, the number of meshing teeth is evolving function of the gear parameters. Thus, the normal load is directly influenced. The traction coefficient, speeds and several other parameters are impacted by the evolution of the gear geometrical data.
An example of those different variations for a smooth contact is illustrated in Figure 5 for data given in Table 2.
It is important to note that all the simulation gear parameters (except for simulation #13) are chosen to simulate a symmetrical behavior of resulting parameters around the pitch point I.
Three different cases are isolated on the previous figure. For this specific study, they correspond to three different areas that will be analyzed. The previously introduced parameters will evolve function of the studied area. On a gear action line, these area have been selected to correspond to specific points of the meshing and the damage phenomena or kinetics might evolve from a zone to another one. These areas are defines as:

The area around A: the approach contact point. This area analysis will be indexed in simulations as #i_A.

The area around I: the rolling without sliding pitch point indexed as #i_I.

The area around B: the recess contact point indexed as #i_B.
These three sections are systematically kept to compare the results between the locations on the gear teeth. The global aim of these section identifications is to potentially identify for a complete gear simulation, a fatigue initiation difference between the tooth addendum, dedendum or pitch point.
The modifications of the input parameters depending on the analyzed area will have a direct influence on the contact and the damaging parameters previously mentioned as the sliding rate, the load or the surface shear stress for example. Thus, the Hertzian pressure computed particularly from the normal load and the equivalent curvature radius is evolving along the gear action line.
The influence of the normal load variation is studied in simulation #6. The contact is again a smooth contact and only the normal load evolution is studied here. The ISS variations are presented in Figure 6. Both results for the reference #1 simulation and those for the three areas are shown here. The stress and the number of cycles to fatigue initiation are also available on the figure.
Concerning the area around the pitch point, the ISS evolution and number of cycles are equivalent between #6_I and #1 as parameters from #1 have been chosen deliberately equal to those at the pitch of #6.
The input normal load being lower around approach and recess contact points, the maximum Hertz pressure are also smaller in these zones than at the pitch point. Thus the lower ISS variations and the higher number of cycles to microcrack nucleation are explained directly. The same way as presented above, the grain boundary (GB) stressed around A [respectively around B] is located more on the left [on the right] but at the same dimensionless depth than the one at the pitch. So the slight shifts in ISS are explained.
So, the variation of the normal load leads to an increase of life duration around A and B because of the lower normal pressure imposed. An equivalent life duration around I is obtained compared to the lone Hertzian contact in #1. It is also shown that the load variation analysis to get an influence on the fatigue criteria is relevant but in a limited scale for this specific gear geometry. All the microcracks locations are represented in Figure 7.
Fig. 5 Evolution of the maximum Hertzian pressure evolution, the traction coefficient, the speeds and identification of the three studied areas along the action line with their corresponding pressure field passage. 
Symmetrical gear characteristics.
Fig. 6 Influence of the number of teeth meshing and the friction coefficient evolution for a smooth contact on ISS, ΔΤ and N results. 
3.4 Complete smooth contact
In order to simulate a complete geared smooth contact, simulation #7 is set. The combined influence a normal load variation and a nonconstant traction coefficient is studied. This typicalshaped traction coefficient has been computed from DIAB's work [25] and is illustrated in Figure 5. It is included on the chosen sections in the following ranges:± [0.083; 0.1] (around the addendum and dedendum) and [–0.05; 0.05] (around pitch point). At the exact position of the pitch point, the value is null but increases rapidly around (cf. Fig. 5). It results in quite high traction coefficients close to point I.
All the ISS and number of cycles to fatigue results are presented for #7 in Figure 6. Figure 7 is illustrating the microcracks nucleation locations in the material relatively to the reference one.
Around the approach and recess point, the respective ISS are shifted to the top [and the bottom] of the graph. This is the direct result of the signed traction coefficient introduction. However, the global ΔΤ and N_{i} values are unchanged between simulations #6 and #7 around A and B. So, the traction coefficient and the normal load evolution combined aren't influencing the fatigue life of this simulated smooth contact. Compared to #1 results, only the load evolution is influencing the global stresses inside the material and its resulting fatigue life. In terms of initiation location, microcracks nucleate closer to the surface (cf. Fig. 7) because of the smallest values of normal load developed in the corresponding simulations.
Around the pitch point, the very low difference between results of simulations #1, #6 and #7 confirms the conclusions depicted above.
3.5 Preliminary conclusions
From all the previous simulations run for a smooth contact, both the influence of a surface shear stress and the variation of the normal load are tested. The amplitude of the load imposed is clearly shown to directly influence the stresses in the material and the number of cycles to fatigue initiation. The conclusion is not the same for the friction coefficient. It is indeed demonstrated above that fatigue results are just slightly or even not modified by the introduction of a surface shear stress. Figure 8 sums up all the simulations run for a smooth contact and the subsequent results (number of cycles to fatigue and ISS) obtained. So, the simulation group around reference #1 confirms the previously depicted conclusions.
Fig. 8 Results for several parameters influence on ΔΤ and N values with smooth contacts. 
4 Study of a rough contact
When a more realistic representation of classical surfaces in contact is wanted, the introduction of a rough component is needed. Surfaces microgeometry is indeed highlighted in literature [7,26,27] as a major parameter to take into account so as to correctly forecast micropitting. So, a sinusoidal roughness geometry (one frequency 1/λ) is set [28]. The approximation made on the roughness profile is serving in this first approach the proper analysis of the results. Thus, equation (9) describes this geometry with Amp and Λ respectively the signal dimensionless amplitude and wavelength.$$\stackrel{}{\mathcal{}}}(X)=\text{Amp}\cdot \mathrm{sin}\left(\frac{2\cdot \pi \cdot X}{\Lambda}\right).$$(9)
Inspired from Johnson [29], the computation of the resulting pressure field is developed in [28] thanks to Johnson's parameter.$$\chi =\frac{{\pi}^{2}\cdot \text{Amp}}{2\cdot \Lambda}.$$(10)Two different cases are then identified.
1. χ < 1: Continuous pressure field
The pressure field is continuous and the rough component of pressure can be computed directly from equation (11).$$\Delta P=\chi .$$(11)
2. χ > 1: Discontinuous pressure field (cf. Fig. 9)
The pressure field is discontinuous and the rough component of pressure is also computed directly from equation (12).$$\Delta P=2\cdot \sqrt{\chi}1.$$(12)
The global pressure field is finally obtained as the direct superposition of the Hertzian pressure field for each time step and the rough component computed.
4.1 Basic rough contact
As for a smooth contact, the global load passage over the surface is divided into several steps of time. For each step on a rough contact, the rough component is considered as fixed (for pure rolling conditions) and only the Hertzian component is moving (with respect to the reference coordinate system). The resulting global pressure field is then computed for each step individually.
Simulation #2 has been set to represent a basic rough contact without sliding nor tangential effect. The corresponding pressure field passage is represented in Figure 9. On the figure, the roughness is indeed fixed relatively to the granular geometry and added to the Hertzian pressure to represent the global normal pressure.
In terms of stresses, the equivalent Tresca stress is also shown in Figure 9. Compared to stresses in #0 and #1, the introduction of roughness to the geometry makes maximum stresses to move closer to the surface. The values of Z presented for both simulations in Table 1 and in Figure 16 go along with this analysis. The ISS variations and the number of cycles to fatigue nucleation is notably modified with this simulation as shown in Figure 10.
Thus, the introduction of roughness produces more damage and at a lower number of cycles than the lone Hertzian contact and in the close subsurface of the material. The ISS shape is explained with the position of the area of maximum stress in the material. The stress elevation corresponds to the moment when the global pressure passes just over the GB. The value stays quite constant since the roughness is not moving (no sliding occurs here) until the rough overpressure goes further than the GB.
Fig. 9 Passage of rough contact of simulation #2 over the finite element model with illustration of the Tresca dimensionless stress field associated. 
Fig. 10 Influence of the traction coefficient for a rough contact on ISS, ΔΤ and N results. 
4.2 Friction in a rough contact
The same way as presented for a smooth contact, friction and sliding are two main and linked parameters to take into account in order to fully simulate a contacting gear. A constant traction coefficient is introduced first in simulation #4. Then, simulation #8 shows the influence of a variable traction coefficient on the different gear sections. The influence of sliding will be studied after.
Once again, the way the tangential component is computed is similar for a rough contact to the smooth one. Equation (8) shows the formulation used. The ISS curves obtained for #4 and their number of cycles to fatigue initiation are presented in Figure 10.
A small modification of ΔT is observed between #2 and #4. Because the damaged GB is almost the same in both simulations (position and orientation), the introduction of a surface shear stress with a constant traction coefficient creates a slight increase of shear stress for a rough contact. The number of cycles to fatigue initiation is directly influenced and the same diminution of life duration is obtained from #2 to #4. However the modifications are too small to conclude that a surface shear stress has an important influence on the initiation and the life duration.
When a contacting gear is considered through simulation #8, it has already been shown that the contact parameters evolve function on the position on the gear action line. Thus, the three previously introduced gear sections (addendum, pitch point and dedendum) are used to describe the different fatigue behaviors function of the influencing parameters.
Results appear to be similar on the three gear sections (Fig. 10). First, the traction coefficient values are equal between the sections around the approach (point A) and recess points (point B) with opposite signs. This explains the similar results observed. The traction coefficient around the pitch point is small. This explains that the results in terms of ISS and number of cycles to fatigue is similar to the one obtained in #2 (without any traction coefficient). Moreover, the decreasing value of the traction coefficient from 0.1 around A and B tends to explain the difference with results of simulation #4.
Finally, concerning the initiation location, the microcracks nucleate at the near surface of the material as shown in Figure 16. The stresses increase due to roughness explains this result and the unmodified initiation depth. So, from these simulations, nearly no influence is observed of the traction coefficient evolution on the initiation of fatigue microcracks for a rough contact.
4.3 Influence of sliding in a rough contact
Sliding is then introduced in simulation #5 (without any tangential effect here). Between two time steps, sliding is set in the pressure fields as a shift in the position of the sinusoidal roughness (Δs defined in Fig. 11). The Hertzian component is still moving the same way as presented before.
The ISS computed by the numerical model presents a very different shape as the one observed for the previous simulations in Figure 12. Some stresses oscillations are observed. They result from the passage of the successive roughness peaks over the initiated GB. This change in the ISS shape will influence directly the fatigue criterion and should be taken into account for its estimation as RYCERZ suggested in [8]. Thus, the way ΔΤ is calculated is modified in order to take into account all the damaging phenomena in the material.
So, a rainflow algorithm is set to fit these new ISS shapes. This algorithm is described in [30]. The principle allows the detection of half and full cycles in a random signal. All the local extrema in cycleshaped signal are isolated this way.
The number of oscillations detected corresponds to the number of rough peaks seen by the damaged GB and can be computed analytically. To do so, the time Δt for a point at the surface of the material to see the whole contacting passing at the sliding speed V_{slide} is considered as described in equation (13). The corresponding distance is the contact width$$\Delta t=\frac{2\cdot a}{{V}_{\text{slide}}}=\frac{2\cdot a}{{u}_{1}{u}_{2}}.$$(13)
Both contacting bodies are respectively moving during Δt by a distance Δs_{1} and Δs_{2} described in equation (14) at their respective speed u_{1} and u_{2}.$$\{\begin{array}{c}\hfill \Delta {s}_{1}=\Delta t\cdot {u}_{1}\hfill \\ \hfill \Delta {s}_{2}=\Delta t\cdot {u}_{2}\hfill \end{array}.$$(14)
The relative dimensionless sliding shift ΔS is then computed in (15) and the number of oscillations seen by the considered point is deduced from (16).$$\Delta S=\frac{\Delta {s}_{1}\Delta {s}_{2}}{a}$$(15) $${N}_{\text{oscillations}}=\frac{\Delta S}{\Lambda}.$$(16)
The speeds values directly depend on the slidetoroll ratio. It is fixed in the first approach of simulation #5. So, all the linked parameters are constant and their variations will be studied later. The analytical data from equation (17) (with numerical data presented in Tab. 1) confirm the number of rough peaks in simulation #5 found in the numerical simulation of six oscillations.$$\begin{array}{c}\hfill \Delta S=2\hfill \\ \hfill \Lambda =0.316\hfill \end{array}\Rightarrow {N}_{\text{oscillations}}\approx \mathrm{6.3.}$$(17)
Two ΔΤ results are presented for simulation #5 in Figure 12. They correspond to the stress variation obtained from the two different extrema detection methods. ΔΤ_{5,th} is given by the maximumminimum detection method described previously. This value is compared to ΔΤ_{5,re} the equivalent stress variation obtained for the oscillating ISS. ΔΤ_{5,re} is computed from the following method:

The total number of oscillations for each GB is estimated with the rainflow algorithm (each extrema is identified in Fig. 12 by a full circle)

For each oscillation in each GB, a local ΔΤ_{i,j } (i being the index of simulation and j the oscillation one) is estimated with the minimummaximum method

The global number of cycles N_{i,}_{re} is then calculated with equation (18)

Finally, the equivalent ΔΤ_{i,}_{re} is obtained from N_{i,}_{re} and equation (7).
These results show that initiation is importantly modified with the rainflow method. The ISS obtained when all the oscillations are taken into account is higher but the difference is more important when the number of cycles to fatigue is considered. The life duration is indeed significantly reduced (around ten times).
The sliding effect on the rough contact pressure is introduced through the shift explained before. On simulation #5, a constant slidetoroll ratio is imposed. However, a contacting gear simulation requires the evolution of this parameter and of the resulting calculated shifts. Depending on the respective speed of each body and the sliding ratio, this shift can be positive or negative from a gear section to another. Thanks to the geometric relations in gears (cf. Fig. 13) [31], the shift can be calculated for each teeth in contact as:$$\{\begin{array}{c}\hfill \Delta {s}_{1}=M\text{'}N\text{'}MN\hfill \\ \hfill \Delta {s}_{2}=M\text{'}L\text{'}ML\hfill \end{array}\Rightarrow \{\begin{array}{c}\hfill \Delta {s}_{1}=MM\text{'}\cdot \frac{2\cdot {T}_{1}M+MM\text{'}}{2\cdot R{b}_{1}}\hfill \\ \hfill \Delta {s}_{2}=MM\text{'}\cdot \frac{2\cdot {T}_{2}MMM\text{'}}{2\cdot R{b}_{2}}\hfill \end{array}.$$(19)
With MM ′ the distance on the gear action line between two successive time steps, T_{1} M and T _{2}M the distances between the points of tangency (with base circles) respectively on gear 1 and gear 2 and the point at the considered step M. R_{b1 } and R_{b2 } are the base circle radii of each gear. Δs_{1} and Δs_{2} are the shifts computed for a fixed value of MM′ for each gear. Those values are inserted inside the pressure computations to get the final sliding fields as described in Figure 11.
So, simulation #9 is set to study the influence of the combined variations of the traction coefficient and the slidetoroll ratio on the material fatigue life. Figure 12 presents the results in terms of ISS and number of cycles to fatigue.
Concerning the results around the pitch point, ΔΤ_{9_I} is slightly lower than ΔΤ_{2}, the simulation for a unique rough contact. The ISS variation is directly linked to the presence of sliding even at a small rate. The sign change in the sliding rate at the pitch point is indeed directly responsible for this stress variation decrease. So, the number of cycles to fatigue is influenced and a bit smaller than for a lone rough contact.
Around the approach and recess points A and B, ΔΤ_{9_A} and ΔΤ_{9_B} are almost similar to ΔΤ_{5}. Thus, the results in terms of fatigue life initiation are very close for those three simulations. Given that, friction has the same limited impact as in a smooth contact and sliding rate is clearly designated as the most influencing one parameter on the fatigue life. Concerning the microcracks location, Figure 16 shows that they also initiate at the near surface of the material.
Fig. 11 Details about the introduction of sliding inside the computed pressure fields. 
Fig. 12 Influence of the sliding rate for a rough contact on ISS, ΔΤ and N results. 
Fig. 13 Illustration of the geometries of meshing gears. 
4.4 Complete symmetrical rough contact
Lastly with the symmetrical gear data, simulation #10 is set with all the influencing parameters (roughness, sliding rate, traction coefficient and normal load evolution).
The results for stresses evolutions and life duration are given in Figure 14. Around the pitch point I, the results are similar to those obtained just before without any load variation. Since this variation is very limited around this point, no difference in the stresses is observed and the fatigue life is strictly the same. The microcrack location remains unchanged.
Around points A and B, a difference is observed with the introduction of the normal load evolution. Less oscillations are indeed observed than on the lone rough sliding simulation #5. This phenomenon is explained by the number of teeth meshing. While only one tooth meshes in the pitch point section, two teeth mesh around the approach and the recess points. So, the normal load decreases and the contact semiwidth is reduced function of the load and the equivalent curvature radius.
The global number of steps seen by the damaged GB is directly reduced and the number of oscillations influenced. With numerical data of equation (20) and theoretical equations (13)–(16), 3.4 oscillations are found for simulations #10_A and #10_B (mean values are taken for variating parameters resulting from the gear geometry).$$\begin{array}{c}\hfill \text{mean}(\Lambda )=0.59\hfill \\ \hfill \left\text{mean}(\Delta S)\right=2\hfill \end{array}\Rightarrow {N}_{\text{oscillations}}\approx \mathrm{3.4.}$$(20)
So the ISS for those two simulations includes between three and four oscillations as shown in Figure 14. The exact value depends on the GB location in the material relatively to the rough overstressed peaks position inside the contact.
The difference observed between N_{10_A} and N_{10_B} comes from additional oscillation noticed on the ISS around the recess point B. Thus, the GB is more stressed at point B than at point A and the number of cycles influenced. This phenomenon is the result of the opposite direction of the combined sliding and friction coefficient around A and B. So, fatigue will potentially initiates earlier around B (which corresponds to the tooth dedendum) than at other sections. However the difference between the two sections results is small and is not representative of the slideroll direction influence. This point is discussed in the conclusions.
Fig. 14 Influence of the sliding rate and the traction coefficient evolutions for a rough contact on ISS, ΔΤ and N results. 
4.5 Complete asymmetrical rough gear geometry
Finally, an asymmetrical gear is set in simulation #13. The associated parameters are presented in Table 1 and are inspired from a FZGC gear [32]. The numerical results are to be compared to the experimental ones. So, all the sections of the driving gear are simulated in order to investigate and potentially highlight a fatigue life difference. The associated ISS and number of cycles to fatigue are presented in Figure 15. The numbers of ISS oscillations are computed again for the two extreme sections in equations (21) and (22). Because of the local variations of the speeds resulting from the gear geometry mean results are presented in equations.$$\text{ for}\text{\#}13\_\text{A}:\text{\u2593}\begin{array}{c}\hfill \text{mean}(\Lambda )=0.382\hfill \\ \hfill \left\text{mean}(\Delta S)\right=2\hfill \end{array}\Rightarrow {N}_{\text{oscillations}}\approx 5.2$$(21) $$\text{ for}\#13\_\text{B}:\text{\u2593}\begin{array}{c}\hfill \text{mean}(\Lambda )=0.21\hfill \\ \hfill \left\text{mean}(\Delta S)\right=2\hfill \end{array}\Rightarrow {N}_{\text{oscillations}}\approx \mathrm{9.5.}$$(22)
Around the pitch point, ISS shapes are similar to those obtained for a symmetrical gear. The limited sliding rate in this section only creates an increasingdecreasing behavior. The first respective microcracks initiate at the near surface and latter (in terms of number of cycles to fatigue) than at the approach and recess points.
At point A, the ISS presents a bit more than five oscillations when nine ones are observed around point B. This actually corresponds to the analytical results. This difference is linked to the velocity variation between the sections resulting from the gear geometry. Thus, the considered driving gear is at the dedendum around point A and at the addendum around point B. The velocity along the tooth profile being higher at the addendum, the shifting distance resulting from sliding is higher around the recess point B than around the approach point A for the same duration. So, according to equation (16), more oscillations are expected in simulation #13_B.
Moreover, the global ISS oscillations have a greater amplitude for #13_A than for #13_B. So, in spite of a greater number of oscillations around point B, the equivalent stresses variation ΔΤ_{13_A} is higher than ΔΤ_{13_B}. The resulting fatigue lives to initiation give an earlier microcrack nucleation at the gear tooth dedendum. Figure 16 also shows that the microcracks location is still close to the surface because the roughness introduced.
This numerical conclusion is comparable to the experimental ones described in [32]. The authors have indeed investigated the micropitting initiation on FZGC type gears and the failure appears to initiate preferentially at the tooth dedendum. So, the earlier initiation simulated above around the driving gear dedendum is the first step of the damage empirically observed after microcrack propagation
Fig. 15 Study of a complete asymmetrical gear: influence of the number of teeth evolution, the sliding rate and the traction coefficient for a rough contact on ISS, ΔΤ and N results. 
Fig. 16 Illustration of microcracks nucleation location for rough contacts. 
5 Conclusions
So, a finite element model oriented for fatigue analyses is presented in this study. It is based on the coupled simulation of a standard granular steel microstructure and on a timeevolving contact pressure field representing the surface microgeometries. The individual and combined influences of several factors (among the surface roughness, the load variation, the traction coefficient and the sliding rate) are also considered. Materials fatigue life is deduced from intergranular shear stress calculation at grain boundaries.
All the rough simulations presented and their results (ISS and number of cycles to fatigue) are summarized in Figure 17.
The most preponderant conclusions resulting from this study are:

All the results around the pitch point section demonstrate similar stress and fatigue behaviors as those obtained for the simplest references. So, a classical calculation with a fixed Hertzian pressure field is appropriate for a smooth contact so as to study the fatigue analysis in this section in any case of influencing parameters (except roughness). The conclusion is the same for the rough contact with only a fixed rough pressure field.

The traction coefficient even at a high rate is of no major influence in the nucleation of fatigue microcracks neither in smooth nor rough contacts.

The most influencing factor on materials fatigue life is contact pressure. The load variations simulated in the study have shown to directly modify the calculated pressure fields, the linked stresses and fatigue criterion. A control in the gear geometry design is needed to take into account this effect at a macro scale. The overpressures resulting from the microgeometrical parameters have to be regulated with surface finish improvement. OILA has given the same conclusions through empirical investigations [7].

A representative gear damage simulation is made possible with this model. Different sections of the meshing gears are represented through their evolving parameters. The study results clearly demonstrate a fatigue life difference (more than a magnitude order) between the microcrack nucleation at the tooth dedendum and addendum resulting from gear geometrical specificities in both sections. These results are correlated to MARTINS's ones [32].

The slidetoroll ratio is also a major influencing factor in fatigue microcracks nucleation. When associated to a rough contact, the relative movement between the two contacting surfaces creates stresses oscillations. These additional local stress cycles worsen the material damage and significantly reduce its fatigue life. This study shows for example a decrease in the number of cycles to microcracks nucleation of two magnitude orders for a simple rough contact with and without a slidetoroll ratio. It confirms RYCERZ's conclusions in [8].
The difference highlighted by RYCERZ in [8] which consists of a shorter fatigue life for contacts subjected to negative slidetoroll ratios rather than positive ones is also observed here on a small scale. However, this phenomenon is actually described in the literature [33,34] as part of the microcrack propagation instead of their initiation. So the results obtained here cannot be an explanation to the previous conclusion. This point is a potential prospect of this study.
Fig. 17 Results for several parameters influence on ΔΤ and N values with rough contacts. 
Nomenclature
amp: Amplitude of roughness [m]
Amp: Dimensionless amplitude of roughness
E_{i } : Young modulus of body i [Pa]
E^{ ′} : Equivalent young modulus [Pa]
L_{model} : Length of the load passage over the Finite Element model (FEM) [m]
M_{0} : Gear normal module [m]
m, σ_{r} : Material parameters form S–N curve [Ø] and [Pa]
n_{i} : Number of cycles to fatigue: first microcrack nucleation
n_{ref} : Number of cycles to fatigue for the reference simulation
N_{i } : Number of cycles for simulation #i reported to the reference
N_{oscillations} : Number of local oscillations in Intergranular Shear Stress graphs
p_{0} : Hertzian pressure [Pa]
p_{max} : Maximum pressure [Pa]
P_{max} : Dimensionless maximum pressure
P_{N,tot} : Dimensionless total normal pressure
P_{T,tot} : Dimensionless total surface shear stress
Δp : Rough peak of pressure [Pa]
ΔP : Dimensionless rough peak of pressure
R_{xi } : Curvature radius of body i in the rolling direction x [m]
R_{yi } : Curvature radius of body i in the gear width direction y [m]
R_{x } : Equivalent curvature radius in x direction [m]
R_{y } : Equivalent curvature radius in y direction [m]
R : Equivalent curvature radius [m]
Z_{i} : Number of teeth of pinion i
Δs_{i } : Microscopic displacement resulting from sliding of body i
ΔS : Global microscopic displacement resulting from sliding
Δt : Time for the contact to pass over a fixed point of the material [s]
Δτ_{i } : Intergranular Shear Stress (ISS) variation of simulation #i [Pa]
Δτ_{ref} : Reference ISS variation [Pa]
ΔT_{i } : Dimensionless ISS variation of simulation #i
u_{i } : Velocity of body i [m/s]
V_{slide} : Gear sliding velocity [m/s]
x : Position of the contact in (O,x,y) [m]
X : Dimensionless position in (O,x,y)
z_{max} : Depth of the maximum τ_{xy} [m]
Z_{max} : Dimensionless depth of τ_{xy} maximum
α_{0} : Gear pressure angle [°]
λ : Sinusoidal roughness wavelength [m]
Λ: Dimensionless roughness wavelength
ω : Gear input rotation velocity [rpm]
ν_{i } : Poisson's coefficient of body i
References
 R.S. DwyerJoyce, J.C. Hamer, J.M. Hutchinson, E. Ioannides, R.S. Sayles, A pitting fatigue life model for gear tooth contacts, Tribol. Ser. 18, 391–400 (1991) [CrossRef] [Google Scholar]
 J.A. Brandão, R. Martins, J.H.O. Seabra, J. Castro, An approach to the simulation of concurrent gear micropitting and mild wear, Wear 324, 64–73 (2015) [Google Scholar]
 A.V. Olver, The mechanism of rolling contact fatigue: an update, Proc. Inst. Mech. Eng. Part J. J. Eng. Tribol. 219, 313–330 (2005) [CrossRef] [Google Scholar]
 R.W. Snidle, H.P. Evans, M.P. Alanou, M.J.A. Holmes, Understanding Scuffing and Micropitting of Gears, Williamsburg, USA, 2004 [Google Scholar]
 R.L. Errichello, Morphology of Micropitting, 2012 [Google Scholar]
 A. Oila, S.J. Bull, Assessment of the factors influencing micropitting in rolling/sliding contacts, Wear 258, 1510–1524 (2005) [Google Scholar]
 P. Rycerz, A. Kadiric, The influence of slideroll ratio on the extent of micropitting damage in rolling sliding contacts pertinent to gear applications, Tribol. Lett. 2018, 1 (2018) [Google Scholar]
 N.R. Paulson, J.A.R. Bomidi, F. Sadeghi, R.D. Evans, Effects of crystal elasticity on rolling contact fatigue, Int. J. Fatigue 61, 67–75 (2014) [Google Scholar]
 N.R. Paulson, F. Sadeghi, W. Habchi, A coupled finite element EHL and continuum damage mechanics model for rolling contact fatigue, Tribol. Int. 107, 173–183 (2017) [Google Scholar]
 N. Raje, F. Sadeghi, J. Rateick, G. Richard, M.R. Hoeprich, A numerical model for life scatter in rolling element bearings, J. Tribol. 130, 011011 (2007) [Google Scholar]
 G.E. MoralesEspejel, Surface roughness effects in elastohydrodynamic lubrication: a review with contributions, Proc. Inst. Mech. Eng. Part J. J. Eng. Tribol. 228, 1217–1242 (2013) [CrossRef] [Google Scholar]
 J.A. Greenwood, K.L. Johnson, The behaviour of transverse roughness in sliding elastohydrodynamically lubricated contacts, Wear 153, 107–117 (1992) [Google Scholar]
 J.A. Greenwood, G.E. MoralesEspejel, The behaviour of transverse roughness in EHL contacts, Proc. Inst. Mech. Eng. Part J. J. Eng. Tribol. 208, 121–132 (1994) [CrossRef] [Google Scholar]
 G.E. MoralesEspejel, C.H. Venner, J.A. Greenwood, Kinematics of transverse real roughness in elastohydrodynamically lubricated line contacts using Fourier analysis, Proc. Inst. Mech. Eng. Part J. J. Eng. Tribol. 214, 523–534 (2000) [CrossRef] [Google Scholar]
 G.E. MoralesEspejel, A.W. Wemekamp, A. FélixQuiñonez, Microgeometry effects on the sliding friction transition in elastohydrodynamic lubrication, Proc. Inst. Mech. Eng. Part J. J. Eng. Tribol. 224, 621–637 (2010) [CrossRef] [Google Scholar]
 G.E. MoralesEspejel, A. Gabelli, The progression of surface rolling contact fatigue damage of rolling bearings with artificial dents, Tribol. Trans. 58, 418–431 (2015) [CrossRef] [Google Scholar]
 G.E. MoralesEspejel, P. Rycerz, A. Kadiric, Prediction of micropitting damage in gear teeth contacts considering the concurrent effects of surface fatigue and mild wear, Wear 398–399, 99–115 (2018) [Google Scholar]
 J.P. Noyel, Analyse Des Mécanismes d'initiation de Fissures En Fatigue de Contact: Approche Mésoscopique, Thèse, Université de Lyon, 2015 [Google Scholar]
 J.P. Noyel, F. Ville, P. Jacquet, A. Gravouil, C. Changenet, Development of a granular cohesive model for rolling contact fatigue analysis: crystal anisotropy modeling, Tribol. Trans. 59, 469–479 (2016) [CrossRef] [Google Scholar]
 M. Le, Influence Des Liserés de Carbures Induits Par La Nitruration Gazeuse Sur Les Mécanismes de Fissuration de Fatigue de Contacts Roulants, 2015 [Google Scholar]
 ASTM E11213, Standard Test Methods for Determining Average Grain Size, 2013 [Google Scholar]
 H. Hertz, Über Die Berührung Fester Elastische Körper Und Über Die Harte, Verhandlungen des Vereins zur Beförderung des, 1882 [Google Scholar]
 P. Lamagnere, R. Fougeres, G. Lormand, A. Vincent, D. Girodin, G. Dudragne, F. Vergne, A physically based model for endurance limit of bearing steels, ASME 120, 421–426 (1998) [Google Scholar]
 N. Raje, F. Sadeghi, R.G. Rateick, A statistical damage mechanics model for subsurface initiated spalling in rolling contacts, J. Tribol. 130, 042201 (2008) [Google Scholar]
 Y. Diab, F. Ville, P. Velex, Prediction of power losses due to tooth friction in gears, Tribol. Trans. 49, 260–270 (2006) [CrossRef] [Google Scholar]
 A. Fabre, L. Barrallier, M. Desvignes, H.P. Evans, M.P. Alanou, Microgeometrical influences on micropitting fatigue damage: multiscale analysis, Proc. Inst. Mech. Eng. Part J. J. Eng. Tribol. 225, 419–427 (2011) [CrossRef] [Google Scholar]
 G.E. MoralesEspejel, V. Brizmer, Micropitting modelling in rollingsliding contacts: application to rolling bearings, Tribol. Trans. 54, 625–643 (2011) [CrossRef] [Google Scholar]
 A. Labiau, F. Ville, P. Sainsot, E. Querlioz, A.A. Lubrecht, Effect of sinusoidal surface roughness under starved conditions on rolling contact fatigue, Proc. Inst. Mech. Eng. Part J. J. Eng. Tribol. 222, 193–200 (2008) [CrossRef] [Google Scholar]
 K.L. Johnson, Contact Mechanics, Cambridge university Press, 1987 [Google Scholar]
 I. Rychlik, Simulation of load sequences from rainflow matrices: Markov Method, Int. J. Fatigue 18, 429–438 (1996) [Google Scholar]
 J. Dufailly, Etude Géométrique Des Engrenages Cylindriques de Transmission de Puissance, 1997 [Google Scholar]
 R. Martins, C. Locatelli, J.H.O. Seabra, Evolution of tooth flank roughness during gear micropitting tests, Ind. Lubr. Tribol. 63, 34–45 (2011) [CrossRef] [Google Scholar]
 M. Kaneta, H. Yatsuzuka, Y. Murakami, Mechanism of crack growth in lubricated rolling/sliding contact, ASLE Trans. 28, 407–414 (1985) [CrossRef] [Google Scholar]
 A.F. Bower, The influence of crack face friction and trapped fluid on surface initiated rolling contact fatigue cracks, J. Tribol. 110, 704 (1988) [Google Scholar]
Cite this article as: G. Vouaillat, J.P. Noyel, F. Ville, X. Kleber, S. Rathery, From Hertzian contact to spur gears: analyses of stresses and rolling contact fatigue, Mechanics & Industry 20, 626 (2019)
All Tables
All Figures
Fig. 1 Organigram of the model. 

In the text 
Fig. 2 Illustration of the static Hertzian pressure field of simulation #0 and the dimensionless Tresca stress field associated. 

In the text 
Fig. 3 Influence of the contact pressure and the traction coefficient variations for a smooth contact on ISS, ΔΤ and N results. 

In the text 
Fig. 4 S–N curve in complete reverse torsion [1]. 

In the text 
Fig. 7 Illustration of microcracks nucleation location for smooth contacts. 

In the text 
Fig. 5 Evolution of the maximum Hertzian pressure evolution, the traction coefficient, the speeds and identification of the three studied areas along the action line with their corresponding pressure field passage. 

In the text 
Fig. 6 Influence of the number of teeth meshing and the friction coefficient evolution for a smooth contact on ISS, ΔΤ and N results. 

In the text 
Fig. 8 Results for several parameters influence on ΔΤ and N values with smooth contacts. 

In the text 
Fig. 9 Passage of rough contact of simulation #2 over the finite element model with illustration of the Tresca dimensionless stress field associated. 

In the text 
Fig. 10 Influence of the traction coefficient for a rough contact on ISS, ΔΤ and N results. 

In the text 
Fig. 11 Details about the introduction of sliding inside the computed pressure fields. 

In the text 
Fig. 12 Influence of the sliding rate for a rough contact on ISS, ΔΤ and N results. 

In the text 
Fig. 13 Illustration of the geometries of meshing gears. 

In the text 
Fig. 14 Influence of the sliding rate and the traction coefficient evolutions for a rough contact on ISS, ΔΤ and N results. 

In the text 
Fig. 15 Study of a complete asymmetrical gear: influence of the number of teeth evolution, the sliding rate and the traction coefficient for a rough contact on ISS, ΔΤ and N results. 

In the text 
Fig. 16 Illustration of microcracks nucleation location for rough contacts. 

In the text 
Fig. 17 Results for several parameters influence on ΔΤ and N values with rough contacts. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.