© G. Vouaillat et al., published by EDP Sciences 2019

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 rolling-sliding contacts and is characterized by typical 10 μm length and depth pits [1–5].

Besides, OILA has pointed out in [6] through twin-disc-based experimental investigations, seven factors influencing this specific failure. Among them, pressure, roughness, friction and slide-to-roll ratio are predominant. More recently, RYCERZ [7] has also identified slide-to-roll ratio as a major influencer in micropitting occurrence and particularly on the rough micro-contact 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 stresses-based 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 semi-analytical 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 time-varying 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 slide-roll ratio) on the nucleation of micropitting for contacting spur gears. The literature-based 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 time-evolving 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 macro-geometry and micro-one 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, micro-geometry 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₀). 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.

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 minimum-maximum 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.

ΔΤ₁ 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).$$\log \left(\Delta \tau \right)=A\cdot \log \left(n_i\right)+\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 micro-crack 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₀. 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 ΔΤ₁ [and ΔΤ_1′ is larger than ΔΤ₁]. This results in an increase for #1′ [respectively a decrease for #1′] in life duration to fatigue initiation. The micro-cracks 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. 4 S–N curve in complete reverse torsion [1].

Fig. 7 Illustration of micro-cracks 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 micro-slip 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 micro-crack 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 macro-geometry 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 micro-crack 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 micro-cracks 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.

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 non-constant traction coefficient is studied. This typical-shaped 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 micro-cracks 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, micro-cracks 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.$${\displaystyle \stackrel{}{\mathcal{}}}(X)=\text{Amp}\cdot \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 micro-cracks 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 micro-cracks 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 cycle-shaped 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₁ and Δs₂ described in equation (14) at their respective speed u₁ and u₂.$$\{\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 slide-to-roll 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 maximum-minimum 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 minimum-maximum 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).

$$N_{i,\mathrm{r}\mathrm{e}}=\frac{{\left(\Delta {\tau }_{\mathrm{r}\mathrm{e}\mathrm{f}}/{\sigma }_r\right)}^m}{{\displaystyle \sum _j}{\left(\Delta {\tau }_{i,j}/{\sigma }_r\right)}^m}.$$(18)

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 slide-to-roll 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}M-MM\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₁ M and T ₂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₁ and Δs₂ 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 slide-to-roll 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 ΔΤ₂, 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 ΔΤ₅. 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 micro-cracks 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 micro-crack 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 semi-width 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 slide-roll 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 FZG-C 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{▓}\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{▓}\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 increasing-decreasing behavior. The first respective micro-cracks 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 micro-crack nucleation at the gear tooth dedendum. Figure 16 also shows that the micro-cracks 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 micro-pitting initiation on FZG-C 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 micro-crack 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 micro-cracks 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 time-evolving 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 micro-cracks 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 micro-geometrical 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 micro-crack 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 slide-to-roll ratio is also a major influencing factor in fatigue micro-cracks 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 micro-cracks nucleation of two magnitude orders for a simple rough contact with and without a slide-to-roll 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 slide-to-roll 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 micro-crack 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

a : Contact semi-width [m]

amp: Amplitude of roughness [m]

Amp: Dimensionless amplitude of roughness

b : Gear width [m]

δ : Maximum deformation [m]

E_i : Young modulus of body i [Pa]

E^′ : Equivalent young modulus [Pa]

F_N : Normal load [N]

L_model : Length of the load passage over the Finite Element model (FEM) [m]

M₀ : Gear normal module [m]

m, σ_r : Material parameters form S–N curve [Ø] and [Pa]

n_i : Number of cycles to fatigue: first micro-crack 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₀ : 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]

T₀ : Input torque [N 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

α₀ : Gear pressure angle [°]

χ : Johnson's parameter

λ : Sinusoidal roughness wavelength [m]

Λ: Dimensionless roughness wavelength

ω : Gear input rotation velocity [rpm]

ν_i : Poisson's coefficient of body i

References

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 micro-cracks 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 micro-cracks 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