Residual Stress Distribution Design for Gear Surfaces Based on Genetic Algorithm Optimization

The rolling contact fatigue of gear surfaces in a heavy loader gearbox is investigated under various working conditions using the critical plane-based multiaxial Fatemi–Socie criterion. The mechanism for residual stress to increase the fatigue initiation life is that the compressive residual stress has a negative normal component on the critical plane. Based on this mechanism, the genetic algorithm is used to search the optimum residual stress distribution that can maximize the fatigue initiation life for a wide range of working conditions. The optimum residual stress distribution is more effective in increasing the fatigue initiation life when the friction coefficient is larger than its critical value, above which the fatigue initiation moves from the subsurface to the surface. Finally, the effect on the fatigue initiation life when the residual stress distribution deviates from the optimum distribution is analyzed. A sound physical explanation for this effect is provided. This yields a useful guideline to design the residual stress distribution.


Introduction
The gearbox is the most important transmission system in many applications such as heavy equipment, automobiles, helicopters, wind turbines, etc. [1,2]. A key aspect in designing a gearbox is to prolong the gear tooth longevity. Rolling contact fatigue (RCF) that results from the cyclic rolling-sliding contact in the gear meshing pair is a common failure mode of the gear tooth surface [3]. It cannot be avoided and must occur as long as gears are operated long enough. RCF first leads to the initiation of fatigue microcracks. These microcracks then propagate and merge into multiple macrocracks, which can form pieces of material removed from the bulk, i.e., pitting. Pitting can accelerate tooth surface deterioration [4], cause transmission error [5], and increase vibration and noise [6], etc. Therefore, it is imperative to investigate the behavior of RCF.
Sharif et al. [7] investigated the influence of surface roughness on RCF under elastohydrodynamic (EHL) contact conditions. They pointed out that the ratio of the film thickness over the surface roughness is the key index of the fatigue performance. Li and Kahraman [8] proposed a novel model for the progression of pitting behavior in spur gear contacts. The novelty lies in the model's capability of simultaneously considering various time-varying variables such as contact load, contact radii, and surface velocities. Paulson et al. [9] found that Hertzian pressure distribution, compared with EHL under the same contact loading, reduces fatigue life. Moreover, the material degradation affected the contact pressure distribution and somehow increased the fatigue life. Wang et al. [10,11] performed extensive finite element simulations of RCF in carburized gear surfaces and revealed that proper carburization process can improve the fatigue performance. Golmohammadi and Sadeghi [12] studied the effect of surface dents on RCF and found that a dent of sharp edge and high pile-up can significantly deteriorate the fatigue performance.
The work in [7][8][9][10][11][12] provided useful insights into RCF. Unfortunately, they ignored the residual stress (RS) induced by shot peening, which is an essential treatment to enhance the fatigue performance in gear manufacturing (see Figure 1, where the typical RS distribution induced by shot peening is shown [13]). Liu et al. [14] revealed that compressive residual stress (CRS) can increase fatigue initiation life in contrast to tensile residual stress (TRS). Paladugu and Hyde [15] additionally considered the material microstructure under dry and lubricated conditions and confirmed the conclusion in [14]. Wang et al. [16] compared TRS and CRS of the same magnitude and found that the decrease in fatigue initiation life induced by TRS is four times the increase by CRS. Xu et al. [17] presented a model to predict the initiation and propagation of cracks due to RCF and showed that the prediction considering residual stresses is in better agreement with the experiments. Ooi et al. [18] theoretically and experimentally demonstrated that the effect of residual stress on the RCF life of carburized AISI 8620 steel varies with the percentage of restrained austenite. vealed that proper carburization process can improve the fatigue performance. Golmohammadi and Sadeghi [12] studied the effect of surface dents on RCF and found that a dent of sharp edge and high pile-up can significantly deteriorate the fatigue performance. The work in [7][8][9][10][11][12] provided useful insights into RCF. Unfortunately, they ignored the residual stress (RS) induced by shot peening, which is an essential treatment to enhance the fatigue performance in gear manufacturing (see Figure 1, where the typical RS distribution induced by shot peening is shown [13]). Liu et al. [14] revealed that compressive residual stress (CRS) can increase fatigue initiation life in contrast to tensile residual stress (TRS). Paladugu and Hyde [15] additionally considered the material microstructure under dry and lubricated conditions and confirmed the conclusion in [14]. Wang et al. [16] compared TRS and CRS of the same magnitude and found that the decrease in fatigue initiation life induced by TRS is four times the increase by CRS. Xu et al. [17] presented a model to predict the initiation and propagation of cracks due to RCF and showed that the prediction considering residual stresses is in better agreement with the experiments. Ooi et al. [18] theoretically and experimentally demonstrated that the effect of residual stress on the RCF life of carburized AISI 8620 steel varies with the percentage of restrained austenite. Although RS was included in the investigation of RCF in [14][15][16][17][18], the effect of the distribution of residual stress was not addressed. Guan et al. [19] and Zhang et al. [20] examined the effect of CRS distribution resulting from different shot peening processes and pointed out the potential of tailoring CRS distribution to improve the fatigue performance. However, no further study on this issue was carried out. Walvekar and Sadeghi [21] studied the effect of CRS distributions through a parametric study and demonstrated the dependence of fatigue performance on variables associated with the CRS distribution profile, such as the depth and magnitude of the peak CRS and the CRS layer thickness (see Figure 1). Morris et al. [22] included the effect of material microstructure and reported that higher CRS near the surface results in deeper initiation of fatigue cracks, which therefore enhances the fatigue performance.
As can be seen from the above literature review, relevant works addressing the effect of RS distribution [19][20][21][22] have provided useful information for the RS distribution design. However, they only considered the effect of a single variable of the distribution profile. No efforts were made to optimize the RS distribution to obtain the best fatigue performance. Therefore, the current study adopts a genetic-algorithm-based optimization scheme to search for the optimum RS distribution. Based on the optimization results, some recommendations for the RS distribution design are provided. As an example, this work is performed for the surfaces of gears used in the gearbox of a heavy loader. Although RS was included in the investigation of RCF in [14][15][16][17][18], the effect of the distribution of residual stress was not addressed. Guan et al. [19] and Zhang et al. [20] examined the effect of CRS distribution resulting from different shot peening processes and pointed out the potential of tailoring CRS distribution to improve the fatigue performance. However, no further study on this issue was carried out. Walvekar and Sadeghi [21] studied the effect of CRS distributions through a parametric study and demonstrated the dependence of fatigue performance on variables associated with the CRS distribution profile, such as the depth and magnitude of the peak CRS and the CRS layer thickness (see Figure 1). Morris et al. [22] included the effect of material microstructure and reported that higher CRS near the surface results in deeper initiation of fatigue cracks, which therefore enhances the fatigue performance.
As can be seen from the above literature review, relevant works addressing the effect of RS distribution [19][20][21][22] have provided useful information for the RS distribution design. However, they only considered the effect of a single variable of the distribution profile. No efforts were made to optimize the RS distribution to obtain the best fatigue performance. Therefore, the current study adopts a genetic-algorithm-based optimization scheme to search for the optimum RS distribution. Based on the optimization results, some recommendations for the RS distribution design are provided. As an example, this work is performed for the surfaces of gears used in the gearbox of a heavy loader.

Residual Stress Distribution
The components of residual stress are depicted in Figure 2a. Both normal components in the x and z directions exhibit identical distributions, while the component in the y direction is small enough to be ignored [23]. The typical RS distribution along the depth (see, e.g., [13]) induced by shot peening is shown in Figure 1. However, the distribution is too complex and requires too many variables to define its profile. This disables efficient optimization of RS distribution (introduced in Section 5). Thus, as was done in [21], the distribution profile is simplified by ignoring the detrimental TRS and assuming the profile to consist of straight lines, as shown in Figure 2b. This simplification reduces the number of variables of the distribution profile without losing the basic features of the RS distribution. The design variables for the optimization of the RS distribution profile are as follows: (1) σ surface -the RS at the surface; (2) σ max -the peak RS; (3) y max -the depth of σ max ; (4) y core -the depth where RS vanishes. nents in the x and z directions exhibit identical distributions, while the component in the y direction is small enough to be ignored [23]. The typical RS distribution along the depth (see, e.g., [13]) induced by shot peening is shown in Figure 1. However, the distribution is too complex and requires too many variables to define its profile. This disables efficient optimization of RS distribution (introduced in Section 5). Thus, as was done in [21], the distribution profile is simplified by ignoring the detrimental TRS and assuming the profile to consist of straight lines, as shown in Figure 2b. This simplification reduces the number of variables of the distribution profile without losing the basic features of the RS distribution. The design variables for the optimization of the RS distribution profile are as follows: (1) σsurface-the RS at the surface; (2) σmax-the peak RS; (3) ymax-the depth of σmax; (4) ycore-the depth where RS vanishes.
The values of these variables can be controlled to formulate different RS distributions by adjusting shot peening process parameters such as the shot velocity, shot size, coverage, etc. [24].  Figure 3a presents the contact problem in a pair of meshing gears used in the 157 kW gearbox of a heavy loader. The geometry and material properties of the gears are presented in Table 1. The contact of a helical gear pair is equivalent to that of a tapered roller and a deformable half-space, to which the solution can be found in [25]. The present study focuses on the methodology for the design of RS distribution. The methodology itself is universal and valid for any type of gear. Therefore, for the convince of demonstrating this methodology, the helix angle is not considered and both the pinion and gear are regarded as spur gears, the simplest type of gear. The contact of meshing teeth is equivalent to that of two cylinders of varying radii (indicated by two dashed arcs), which depends on the position of the meshing point on the line of action (LOA). Since the gear contact is elastic, the contact of two cylinders can be simplified as that of a rigid cylinder of equivalent radius R and a deformable half-space of equivalent Young's modulus E * in plane strain condition [23], as shown in Figure 3b. The values of these variables can be controlled to formulate different RS distributions by adjusting shot peening process parameters such as the shot velocity, shot size, coverage, etc. [24]. Figure 3a presents the contact problem in a pair of meshing gears used in the 157 kW gearbox of a heavy loader. The geometry and material properties of the gears are presented in Table 1. The contact of a helical gear pair is equivalent to that of a tapered roller and a deformable half-space, to which the solution can be found in [25]. The present study focuses on the methodology for the design of RS distribution. The methodology itself is universal and valid for any type of gear. Therefore, for the convince of demonstrating this methodology, the helix angle is not considered and both the pinion and gear are regarded as spur gears, the simplest type of gear. The contact of meshing teeth is equivalent to that of two cylinders of varying radii (indicated by two dashed arcs), which depends on the position of the meshing point on the line of action (LOA). Since the gear contact is elastic, the contact of two cylinders can be simplified as that of a rigid cylinder of equivalent radius R and a deformable half-space of equivalent Young's modulus E * in plane strain condition [23], as shown in Figure 3b.   The contact model in Figure 3b aims to simulate the contact in the region close to the pitch point, where it is more susceptible to RCF since there is only one pair of teeth in contact [26]. Thus, the rated meshing force per unit length Prated = 423 N/mm, which can be easily obtained using the parameters in Table 1 [27], is the normal contact load. To gain basic physical insights, the variation of R is ignored and R = 15.231 mm at the pitch point is used for the rigid cylinder. Due to the involute geometry, the gear surface is subjected to rolling-sliding contact. This type of contact loading can be represented by Hertzian pressure with shear traction:  The contact model in Figure 3b aims to simulate the contact in the region close to the pitch point, where it is more susceptible to RCF since there is only one pair of teeth in contact [26]. Thus, the rated meshing force per unit length P rated = 423 N/mm, which can be easily obtained using the parameters in Table 1 [27], is the normal contact load. To gain basic physical insights, the variation of R is ignored and R = 15.231 mm at the pitch point is used for the rigid cylinder. Due to the involute geometry, the gear surface is subjected to rolling-sliding contact. This type of contact loading can be represented by Hertzian pressure with shear traction:

Gear Contact Model
where p is the normal pressure, P is the applied normal load, x c is the x-coordinate of the normal pressure center, b = (4PR/(πE * )) 0.5 is the half contact width, t is the shear traction, and µ is the friction coefficient ranging from 0 to 0.3, representing ideally lubricated to dry friction conditions. The analytical solution for the stress/strain field in the half-space was provided in [28] (see Appendix A). When the RS distribution is incorporated in the half-space, it can be demonstrated that the deformation in the half-space is still elastic and the RS distribution is not altered as the number of loading cycles increases. Therefore, the RS distribution can be directly added to the stress solution provided in [28] to obtain the stress field with the presence of RS. Since material points at the same depth are subjected to the same loading cycle, the fatigue analysis will be conducted only for target points at x = 0 (see the red line in Figure 3b). x c is varied from −30b to 30b so that the pressure profile moves over a sufficiently long distance to ensure that the material points at x = 0 undergo a complete loading cycle.

Fatemi-Socie Multiaxial Fatigue Criterion
The complex multiaxial stress state exists in the half-space subjected to rolling-sliding contact. Moreover, this kind of contact loading is a typical non-proportional loading that results in stress components varying non-proportionally during the loading process, i.e., the principal directions constantly change. To analyze the RCF, a proper multiaxial fatigue criterion that takes the effect of the loading non-proportionality into account is needed. The critical plane-based Fatemi-Socie fatigue criterion [29] showed good agreement between the predicted and experimental results in RCF in gears [30] and was adopted extensively in many theoretical studies on gear RCF (see, e.g., [14]). Therefore, this criterion is also used in the present study.
The Fatemi-Socie fatigue criterion states that a fatigue microcrack at a material point is initiated along the critical plane that has the maximum shear strain range over a loading cycle. However, each material point usually has multiple critical planes. In this case, the fatigue microcrack occurs along the critical plane with the maximum fatigue damage. The fatigue damage of one material point on an arbitrary plane can be quantified by the Fatemi-Socie damage parameter (FSDP) as follows: where ∆γ and σ n,max are the maximum shear strain range and maximum normal stress on this plane over one loading cycle, respectively; Y is the material yield strength, and k is a fitting parameter ranging from 0 to 1 [31]. For high cycle fatigue, which is the case in the present study, k approaches unity. Thus, k = 1 is adopted. The normal stress component is included to consider the fact that compressive stress tends to retard fatigue microcrack initiation while tensile stress accelerates the initiation. The maximum FSDP on all the critical planes is denoted by D FS . The fatigue initiation life N f of a material point is related to D FS by: where N f is the number of loading cycles when the fatigue microcrack is initiated; G = E/(2 + 2ν); τ f and γ f are the shear fatigue strength coefficient and the shear fatigue ductility coefficient, respectively; m and n are the fatigue strength exponent and the fatigue ductility exponent, respectively. τ f = 1296 MPa, γ f = 0.437, m = −0.087, and n = −0.58. Note that temperature may have an impact on the above parameters. However, this is not considered since this paper only focuses on demonstrating the methodology for the design of RS distribution. One can easily incorporate a temperature-dependent fatigue citation in the proposed methodology to obtain a temperature-dependent optimum RS distribution. Figure 4a shows a flow chart that illustrates the procedure to predict D FS and N f for one material point. The critical plane needs to be first identified. To this end, the stressstrain history over a loading cycle for a material point should be obtained. This can be done by taking the coordinates of this point into Equations (A1)-(A7) and varying x c from −30b to 30b. Then, all the candidate planes that pass this point are examined to find the critical plane. Each candidate plane is characterized by the angle α between its normal direction and the x-axis (see Figure 4b), ranging from 0 • to 180 • with intervals of 0.2 • . The stressstrain components on each candidate plane can be obtained by coordinate transformation.
Each candidate plane's Δγ and σn,max can thus be obtained. Then, the critical planes can be determined. DFS is the largest FSDP among all the critical planes. Finally, the fatigue initiation life of this point, Nf, can be calculated using Equation (6).  Each candidate plane's ∆γ and σ n,max can thus be obtained. Then, the critical planes can be determined. D FS is the largest FSDP among all the critical planes. Finally, the fatigue initiation life of this point, N f , can be calculated using Equation (6).

Optimization Scheme
The genetic algorithm (GA) is an optimization methodology inspired by the "survival of the fittest" rule in the theory of evolution. It is an efficient and robust method to solve complex optimization problems [32][33][34][35]. In gear-related problems, it has been used to optimize the tooth profile [35], maintain the reliability of gear transmission [36], and improve the load-sharing performance of planetary gears [37]. In this paper, the GA is used to search for the optimum RS distribution that maximizes the fatigue initiation life of the gear surface, which is (N f ) min , i.e., the minimum N f among all the points. Since (N f ) min increases with decreasing (D FS ) max according to Equation (6), the objective of this optimization problem is to minimize (D FS ) max . As can be seen from Figure 2, σ surface , σ max , y max , and y core are the four design variables for the optimization, and thus, the constraints are applied on these four variables.
The mathematical description of this optimization problem is as follows: where x = (σ surface , σ max , y max , y core ) is called an individual in GA. Each component of the individual is called a "gene". The constraints in Equation (8b-g) are summarized from the existing experimental [38][39][40][41] and theoretical [13,42] results for RS distribution profiles produced using the shot peening process. The basic steps in the GA [43] are illustrated by the flow chart in Figure 5.

Mechanism for RS to Increase (Nf)min
As mentioned, since the fatigue damage along the x-axis direction exhibits a uniform distribution, fatigue analysis was performed for material points at x = 0. Figure 6 shows the shear strain range Δγ, the maximum normal stress over a loading cycle σn,max, and the FSDP on all the candidate planes (0° ≤ α ≤ 180°) of material points at x = 0 for three working conditions where P/Prated = 1 with μ = 0, 0.1, and 0.3 without the presence of RS. These results are normalized by their maximum absolute values. The critical planes of each material point along the depth are denoted by the solid cyan lines in the Δγ and FSDP distribution.
As would be expected from the theory of elasticity [44], each point has two critical planes with equal maximum Δγ. These two planes have an α with a difference of 90°. For Step 1: Population initialization. Generate a random initial population containing N individuals that satisfy the constraints in Equation (8b-g).
Step 2: Fitness evaluation. Calculate (D FS ) max for each individual and use 1/(D FS ) max to measure its fitness. Individuals having higher fitness are more likely to survive.
Step 3: Population evolution. Three genetic operators are needed for population evolution to select and generate individuals with high fitness and eliminate those with low fitness. The selection operator is first applied to specify individuals of higher fitness as parents for the next generation. Then, two parents exchange one or multiple genes and generate two new individuals using the crossover operator. The crossover fraction, P c , decides how many individuals are generated in the next generation. Finally, genes of individuals are altered randomly using the mutation operator.
Step 4: Termination. The optimization process is terminated once one of the following two criteria is satisfied: (1) The generation number reaches the prescribed upper limit G max ; (2) The relative change in the highest fitness over G s generations is less than the function tolerance value e. The optimum solution is the individual with the highest fitness in the last generation. The input parameters for the GA are listed in Table 2. As mentioned, since the fatigue damage along the x-axis direction exhibits a uniform distribution, fatigue analysis was performed for material points at x = 0. Figure 6 shows the shear strain range ∆γ, the maximum normal stress over a loading cycle σ n,max , and the FSDP on all the candidate planes (0 • ≤ α ≤ 180 • ) of material points at x = 0 for three working conditions where P/P rated = 1 with µ = 0, 0. since the σn,max distribution over α is not uniform, the FSDP on the two critical planes is different according to Equation (5). Thus, there is only one critical plane having the largest FSDP (i.e., DFS) and the orientation of the critical plane corresponding to DFS is denoted by αc.  Figure 7 shows a typical RS distribution and the normal component σn,r introduced by this RS distribution on all candidate planes. The σn,r distribution in Figure 7b can be added to the σn,max distribution without RS in Figure 6 to obtain the new σn,max distribution with RS. In this way, the FSDP distribution can be changed by RS. Nonetheless, there is an exception that when αc is close to 90 • (see Figure 6, μ = 0), Nf will remain unchanged since RS will have zero normal components on the αc plane and will not be able to decrease DFS in such a case. As would be expected from the theory of elasticity [44], each point has two critical planes with equal maximum ∆γ. These two planes have an α with a difference of 90 • . For µ = 0, the critical plane orientations are independent of the depth. For µ = 0, they depend on the depth due to the shear traction induced by friction. The effect of shear traction decays with increasing depth. The effect of σ n,max on the fatigue damage is demonstrated by the difference between the distributions of ∆γ and FSDP. σ n,max is compressive for µ = 0 and tensile for µ = 0. For µ = 0 and 0.1, the effect of σ n,max is negligible. Due to the mild friction on the surface, the maximum damage occurs below the surface. On the other hand, for µ = 0.3, the effect of σ n,max is appreciable and the high friction gives rise to the maximum damage appearing on the surface. For all the three cases (P/P rated = 1; µ = 0, 0.1, and 0.3), since the σ n,max distribution over α is not uniform, the FSDP on the two critical planes is different according to Equation (5). Thus, there is only one critical plane having the largest FSDP (i.e., D FS ) and the orientation of the critical plane corresponding to D FS is denoted by α c . Figure 7 shows a typical RS distribution and the normal component σ n,r introduced by this RS distribution on all candidate planes. The σ n,r distribution in Figure 7b can be added to the σ n,max distribution without RS in Figure 6 to obtain the new σ n,max distribution with RS. In this way, the FSDP distribution can be changed by RS. Nonetheless, there is an exception that when α c is close to 90 • (see Figure 6, µ = 0), N f will remain unchanged since RS will have zero normal components on the α c plane and will not be able to decrease D FS in such a case.  Figure 7 shows a typical RS distribution and the normal component σn,r introduced by this RS distribution on all candidate planes. The σn,r distribution in Figure 7b can be added to the σn,max distribution without RS in Figure 6 to obtain the new σn,max distribution with RS. In this way, the FSDP distribution can be changed by RS. Nonetheless, there is an exception that when αc is close to 90 • (see Figure 6, μ = 0), Nf will remain unchanged since RS will have zero normal components on the αc plane and will not be able to decrease DFS in such a case.    Figure 8a, as µ increases, the location of (D FS ) max moves from the depth of around 0.13 mm to the surface. µ c can be defined as the critical friction coefficient. For µ < µ c , the location of (D FS ) max is below the surface. For µ > µ c , the location of (D FS ) max is on the surface. It should be noted that the D FS distributions for different normal loadings ranging from 0.5P rated to 1.5P rated share similar profiles. For this reason, Figure 8a is only shown for a fixed normal load with different µ and µ c = 0.23 is independent of normal loading.  Figure 8a, as μ increases, the location of (DFS)max moves from the depth of around 0.13 mm to the surface. μc can be defined as the critical friction coefficient. For μ < μc, the location of (DFS)max is below the surface. For μ > μc, the location of (DFS)max is on the surface. It should be noted that the DFS distributions for different normal loadings ranging from 0.5Prated to 1.5Prated share similar profiles. For this reason, Figure 8a is only shown for a fixed normal load with different μ and μc = 0.23 is independent of normal loading. Figure 8b presents αc associated with (DFS)max as a function of μ under the rated normal load. For 0 < μ ≤ 0.3, αc associated with (DFS)max is not close to 90 • . This provides the solid theoretical basis that (DFS)max can be decreased or (Nf)min can be increased by introducing a proper RS distribution. Since there are four variables to define the RS distribution profile, the optimum RS distribution that can maximize (Nf)min was sought using the optimization scheme introduced in the previous section.  Figure 9a presents the percentage increase in (Nf)min induced by the optimum RS distribution compared to the case without RS. From the discussion on Figure 8, the DFS dis-  Figure 8b presents α c associated with (D FS ) max as a function of µ under the rated normal load. For 0 < µ ≤ 0.3, α c associated with (D FS ) max is not close to 90 • . This provides the solid theoretical basis that (D FS ) max can be decreased or (N f ) min can be increased by introducing a proper RS distribution. Since there are four variables to define the RS distribution profile, the optimum RS distribution that can maximize (N f ) min was sought using the optimization scheme introduced in the previous section. Figure 9a presents the percentage increase in (N f ) min induced by the optimum RS distribution compared to the case without RS. From the discussion on Figure 8, the D FS distribution depends on the working condition and so does the optimum RS distribution. Thus, the percentage increase in (N f ) min varies with the working condition. As the normal loading and µ increase, the percentage increase in (N f ) min increases. Compared with the increasing normal loading, µ has a higher impact on the percentage increase in (N f ) min . Moreover, when µ > µ c , the percentage increase in (N f ) min is dramatically higher than that when µ < µ c . engineering practice, it is almost impossible to precisely achieve the optimum RS distribution in this way. Thus, it is desired that the optimum RS distribution can be varied to some extent while (Nf)min is not sacrificed too much. By this means, it is more apt to adjust the shot peening process parameters to obtain an acceptable RS distribution and achieve a satisfying (Nf)min. To this end, it is necessary to investigate how the variables of the RS distribution profile affect (Nf)min when they deviate from those of the optimum distribution.

Effect on (Nf)min When the RS Distribution Deviates from the Optimum
Since the optimum RS distribution depends on the DFS distribution, which has a weak dependence on the normal load P, the effect on (Nf)min when the RS distribution deviates from the optimum is also weakly dependent on P. In light of this, the optimum RS distributions under investigation are those for fixed P with different μ. Figure 10 shows the optimum RS stress distribution for P/Prated = 1 and μ = 0.1, 0.2, and 0.3. The optimum distributions exhibit two common features. They all have nearly the same (ymax)opt of 0.5b, which is the depth of the maximum Δγ [21], and the same (σmax)opt = −1000 MPa, which is the lower limit of σmax. (σsurface)opt and (ycore)opt do not exhibit a simple dependence on μ since they do not monotonically increase or decrease with μ. The effect on (Nf)min when the RS distribution deviates from the optimum shall be investigated by deviating a single component of (σsurface, σmax, ymax, and ycore) from its corresponding optimum value. Figure 11 shows the percentage decrease in (Nf)min when ymax and ycore deviate from their optimum values. As can be seen from Figure 11a, the deviation from (ymax)opt causes a slight decrease in (Nf)min. As μ increases, (Nf)min becomes more sensitive to the deviation of ymax. From Figure 11b, it can be seen that the deviation of ycore almost has no effect on (Nf)min. Figure 12 shows and explains the effect on (Nf)min when σsurface deviates from (σsurface)opt. As can be seen from Figure 12a, for relatively low μ, deviation of σsurface does not affect (Nf)min since (DFS)max occurs in the subsurface. However, for μ = 0.3, (Nf)min is very sensitive to the deviation of σsurface unless σsurface/(σsurface)opt becomes larger than 0.96. This can be explained by Figure 12b. As |σsurface| increases, (DFS)max, which appears on the surface, decreases and the location of (DFS)max moves from the surface to the subsurface when σsurface/(σsurface)opt reaches 0.96. Thereafter, (DFS)max in the subsurface is not affected by σsurface  Figure 9b presents the percentage increase in the depth of (D FS ) max induced by the optimum RS distribution for µ < µ c . For µ > µ c , (D FS ) max moves from the surface to the subsurface. In this case, the depth of (D FS ) max when the optimum RS distribution is applied is presented. The increase in the depth of (D FS ) max is beneficial in increasing the propagation life since more loading cycles are needed for the fatigue microcrack to propagate to the surface before a pitting particle is formed [14]. Provided that the relation between the depth of (D FS ) max and the propagation life is quantified, the RS distribution can be optimized to maximize the sum of the initiation life, (N f ) min , and propagation life. Nonetheless, this is out of the scope of the present study. Figure 9a demonstrates that the optimum RS distribution, which is obtained using the GA, can indeed substantially increase (N f ) min . Theoretically speaking, the shot peening process parameters can be adjusted to achieve the optimum RS distribution. However, in engineering practice, it is almost impossible to precisely achieve the optimum RS distribution in this way. Thus, it is desired that the optimum RS distribution can be varied to some extent while (N f ) min is not sacrificed too much. By this means, it is more apt to adjust the shot peening process parameters to obtain an acceptable RS distribution and achieve a satisfying (N f ) min . To this end, it is necessary to investigate how the variables of the RS distribution profile affect (N f ) min when they deviate from those of the optimum distribution.

Effect on (N f ) min When the RS Distribution Deviates from the Optimum
Since the optimum RS distribution depends on the D FS distribution, which has a weak dependence on the normal load P, the effect on (N f ) min when the RS distribution deviates from the optimum is also weakly dependent on P. In light of this, the optimum RS distributions under investigation are those for fixed P with different µ. Figure 10 shows the optimum RS stress distribution for P/P rated = 1 and µ = 0.1, 0.2, and 0.3. The optimum distributions exhibit two common features. They all have nearly the same (y max ) opt of 0.5b, which is the depth of the maximum ∆γ [21], and the same (σ max ) opt = −1000 MPa, which is the lower limit of σ max . (σ surface ) opt and (y core ) opt do not exhibit a simple dependence on µ since they do not monotonically increase or decrease with µ. The effect on (N f ) min when the RS distribution deviates from the optimum shall be investigated by deviating a single component of (σ surface , σ max , y max , and y core ) from its corresponding optimum value.    Figure 11 shows the percentage decrease in (N f ) min when y max and y core deviate from their optimum values. As can be seen from Figure 11a, the deviation from (y max ) opt causes a slight decrease in (N f ) min . As µ increases, (N f ) min becomes more sensitive to the deviation of y max . From Figure 11b, it can be seen that the deviation of y core almost has no effect on (N f ) min . Figure 12 shows and explains the effect on (N f ) min when σ surface deviates from (σ surface ) opt . As can be seen from Figure 12a, for relatively low µ, deviation of σ surface does not affect (N f ) min since (D FS ) max occurs in the subsurface. However, for µ = 0.3, (N f ) min is very sensitive to the deviation of σ surface unless σ surface /(σ surface ) opt becomes larger than 0.96. This can be explained by Figure 12b. As |σ surface | increases, (D FS ) max , which appears on the surface, decreases and the location of (D FS ) max moves from the surface to the subsurface when σ surface /(σ surface ) opt reaches 0.96. Thereafter, (D FS ) max in the subsurface is not affected by σ surface and neither is (N f ) min .
As shown in Figure 10, (σ max ) opt for all the three cases reached its lower limit, −1000 MPa. To investigate the effect on (N f ) min when σ max deviates from (σ max ) opt , σ max has to be increased. However, since the constraint in Equation (8b) must be satisfied, increasing σ max should be coupled with increasing σ surface . To this end, we simply increased both σ max and σ surface by the same amount from (σ max ) opt and (σ surface ) opt , respectively. In this way, the effect of deviation of σ max was studied. Figure 13 shows (N f ) min and the orientation of the critical plane corresponding to (D FS ) max , α c , as functions of |σ max | for four values of µ. The effect of deviation of σ max on (N f ) min is salient in agreement with the experimental results reported in [45][46][47] that the peak RS plays a major role in the fatigue performance.
For µ < µ c , the location of (D FS ) max remains in the subsurface. As shown in Figure 6 (µ = 0.1), the material point at this location has two critical planes of maximum ∆γ, which are denoted as planes 1 and 2 with orientations of α 1 and α 2 close to 0 • and 90 • , respectively. Initially, without RS, i.e., |σ max | = 0, plane 1 has a higher σ n,max than plane 2 (see Figure 6). According to Equations (5) and (6), the FSDP on plane 1 is larger and is thus equal to (D FS ) max . Consequently, α c = α. Since the magnitude of the negative normal component of RS on plane 1 is larger than on plane 2 (see Figure 7b), as |σ max | increases, the FSDP on plane 1 decreases faster and is exceeded by the FSDP on plane 2. As a result, the FSDP on plane 2 becomes (D FS ) max so that α c switches to α 2 . This results in the two distinct stages of (N f ) min vs. |σ max | in Figure 13a,b.
For µ ≥ µ c , the location of (D FS ) max moves from the surface to the subsurface as |σ max | increases. As shown in Figure 6 (µ = 0.3), when (D FS ) max is on the surface, the two critical planes of maximum ∆γ of the material point at the location of (D FS ) max are denoted as planes 1 and 2. When (D FS ) max moves to the subsurface, the counterparts of planes 1 and 2 are denoted by planes 3 and 4. As can be seen, since α 1 ≈ 45 • and α 2 ≈ 135 • , the negative normal component on planes 1 and 2 are almost the same (see Figure 7b). This results in almost identical FSDP vs. |σ max | for these two planes. Thus, unlike Figure 13a,b, the switch of α c from α 2 to α 1 (see Figure 13d) is not reflected as the deflection on the curve of (N f ) min vs. |σ max |. Therefore, the two distinct stages of (N f ) min vs. |σ max | in Figure 13c,d result from the change in the location of (D FS ) max , which is different from the case of µ < µ c in Figure 13a   As shown in Figure 10, (σmax)opt for all the three cases reached its lower limit, −1000 MPa. To investigate the effect on (Nf)min when σmax deviates from (σmax)opt, σmax has to be   As shown in Figure 10, (σmax)opt for all the three cases reached its lower limit, −1000 MPa. To investigate the effect on (Nf)min when σmax deviates from (σmax)opt, σmax has to be increased. However, since the constraint in Equation (8b) must be satisfied, increasing σmax should be coupled with increasing σsurface. To this end, we simply increased both σmax and σsurface by the same amount from (σmax)opt and (σsurface)opt, respectively. In this way, the effect of deviation of σmax was studied. Figure 13 shows (Nf)min and the orientation of the critical plane corresponding to (DFS)max, αc, as functions of |σmax| for four values of μ. The effect of deviation of σmax on (Nf)min is salient in agreement with the experimental results reported in [45][46][47] that the peak RS plays a major role in the fatigue performance. For μ < μc, the location of (DFS)max remains in the subsurface. As shown in Figure 6 (μ = 0.1), the material point at this location has two critical planes of maximum Δγ, which are denoted as planes 1 and 2 with orientations of α1 and α2 close to 0° and 90°, respectively. Initially, without RS, i.e., |σmax| = 0, plane 1 has a higher σn,max than plane 2 (see Figure 6). According to Equations (5) and (6), the FSDP on plane 1 is larger and is thus equal to (DFS)max. Consequently, αc = α. Since the magnitude of the negative normal component of RS on plane 1 is larger than on plane 2 (see Figure 7b), as |σmax| increases, the FSDP on plane 1 decreases faster and is exceeded by the FSDP on plane 2. As a result, the FSDP on plane 2 becomes (DFS)max so that αc switches to α2. This results in the two distinct stages of (Nf)min vs. |σmax| in Figure 13a,b.
For μ ≥ μc, the location of (DFS)max moves from the surface to the subsurface as |σmax| increases. As shown in Figure 6 (μ = 0.3), when (DFS)max is on the surface, the two critical planes of maximum Δγ of the material point at the location of (DFS)max are denoted as planes 1 and 2. When (DFS)max moves to the subsurface, the counterparts of planes 1 and 2 are denoted by planes 3 and 4. As can be seen, since α1 ≈ 45 • and α2 ≈ 135 • , the negative normal component on planes 1 and 2 are almost the same (see Figure 7b). This results in almost identical FSDP vs. |σmax| for these two planes. Thus, unlike Figure 13a,b, the switch of αc from α2 to α1 (see Figure 13d) is not reflected as the deflection on the curve of (Nf)min Despite the reason leading to the two distinct stages of (N f ) min vs. |σ max | being different for µ < µ c and µ ≥ µ c , it is interesting to investigate |σ max,t |, at which the two stages are divided. As demonstrated in Figure 14, |σ max,t | increases with µ significantly faster when µ > µ c . The increasing rates of (N f ) min with |σ max | before and after |σ max,t | are also presented in this figure. For µ < µ c , (N f ) min vs. |σ max | is bilinear and the increasing rate is thus the slope. For µ ≥ µ c , (N f ) min vs. σ max is nonlinear when |σ max | < |σ max,t |. Thus, the increasing rate is obtained as the average of ∂(N f ) min /∂|σ max | over 0 < |σ max | < |σ max,t | for qualitative analysis. In general, the increasing rate of (N f ) min with |σ max | when |σ max | < |σ max,t | is higher. This indicates that a further increase in |σ max | over |σ max,t | becomes less effective in increasing (N f ) min . This theoretical finding is in accordance with experimental results showing that RS induced by additional shot peening is less effective in preventing fatigue in 17NiVrMo6-4 steel [48,49] and Austempered ductile iron [50].
The results presented in Figures 13 and 14 provide useful information on the effect on (N f ) min when σ max deviates from (σ max ) opt . It is true that the ideal scenario is that |σ max | = |(σ max ) opt | = 1000 MPa. However, a higher |σ max | requires higher shot velocity that leads to higher surface roughness. Moreover, when |σ max | > |σ max,t |, it becomes less effective in increasing (N f ) min by increasing |σ max |. This is especially problematic for the case of low µ, where α c is close to 90 • (see Figure 13a). In this case, the negative component of RS on this critical plane is trivial and the increasing rate of (N f ) min with |σ max | is also very slow. Therefore, when one tries to increase (N f ) min by increasing |σ max |, he may need to comprehensively consider the slowing increasing rate of (N f ) min with |σ max | and the increasing surface roughness. slope. For μ ≥ μc, (Nf)min vs. σmax is nonlinear when |σmax| < |σmax,t|. Thus, the increasing rate is obtained as the average of ∂(Nf)min/∂|σmax| over 0 < |σmax| < |σmax,t| for qualitative analysis. In general, the increasing rate of (Nf)min with |σmax| when |σmax| < |σmax,t| is higher. This indicates that a further increase in |σmax| over |σmax,t| becomes less effective in increasing (Nf)min. This theoretical finding is in accordance with experimental results showing that RS induced by additional shot peening is less effective in preventing fatigue in 17NiVrMo6-4 steel [48,49] and Austempered ductile iron [50]. Figure 14. |σmax, t| and the increasing rate of (Nf)min with |σmax| before and after |σmax,t| as functions of μ.
The results presented in Figures 13 and 14 provide useful information on the effect on (Nf)min when σmax deviates from (σmax)opt. It is true that the ideal scenario is that |σmax| = |(σmax)opt| = 1000 MPa. However, a higher |σmax| requires higher shot velocity that leads to higher surface roughness. Moreover, when |σmax| > |σmax,t|, it becomes less effective in increasing (Nf)min by increasing |σmax|. This is especially problematic for the case of low μ, where αc is close to 90 • (see Figure 13a). In this case, the negative component of RS on this critical plane is trivial and the increasing rate of (Nf)min with |σmax| is also very slow. Therefore, when one tries to increase (Nf)min by increasing |σmax|, he may need to comprehensively consider the slowing increasing rate of (Nf)min with |σmax| and the increasing surface roughness.

Conclusions
An RCF analysis based on the multiaxial Fatemi-Socie fatigue criterion was conducted for gear surfaces in a heavy loader gearbox under various working conditions with different combinations of meshing force and friction coefficient. The effect of RS on the fatigue performance was intensively studied. As μ increases, the fatigue initiation location moves from the subsurface to the surface. The μ at which this transition occurs was defined as μc. The mechanism for RS to increase the fatigue initiation life is that the compressive RS has a negative component on the critical plane and, thus, suppresses the fatigue damage. Figure 14. |σ max,t | and the increasing rate of (N f ) min with |σ max | before and after |σ max,t | as functions of µ.

Conclusions
An RCF analysis based on the multiaxial Fatemi-Socie fatigue criterion was conducted for gear surfaces in a heavy loader gearbox under various working conditions with different combinations of meshing force and friction coefficient. The effect of RS on the fatigue performance was intensively studied. As µ increases, the fatigue initiation location moves from the subsurface to the surface. The µ at which this transition occurs was defined as µ c . The mechanism for RS to increase the fatigue initiation life is that the compressive RS has a negative component on the critical plane and, thus, suppresses the fatigue damage.
The RS distribution, simplified as a piecewise bilinear function of the depth to the surface, is described by four variables, including y surface , y core , σ max , and σ surface . The optimum RS distribution that can increase the fatigue initiation life of the gear surface (N f ) min to the maximum was sought using the GA. The optimum RS distribution is not universal but depends on the working condition. It was demonstrated that the optimum RS distribution can substantially increase (N f ) min under a wide range of working conditions. This improvement in fatigue performance is more significant for higher µ. Another unexpected benefit brought by the optimum RS distribution is that the fatigue initiation occurs at a deeper depth, which increases the fatigue propagation life.
In practice, it is difficult to exactly achieve the optimum RS distribution by adjusting shot peening process parameters. Thus, it is important to explore to what extent modifications can be made on the optimum distribution without too much sacrificing of (N f ) min . It was found that for any µ, be it larger or smaller than µ c , the deviation of y core and y max from their optimum values has a negligible effect on (N f ) min . On the other hand, σ max has a salient effect and |σ max | should be as large as possible. However, in practice, a cautious decision should be made since increasing |σ max | increases surface roughness and can have a trivial effect on increasing (N f ) min when µ is low or |σ max | > |σ max,t |. The effect of the deviation of σ surface depends on µ. For µ < µ c , σ surface has no effect. For µ > µ c , σ surface should be close to (σ surface ) opt , leaving little room for acceptable deviation of σ surface . This information can serve as a general guidance on RS distribution design by properly modifying the optimum one obtained from GA optimization.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. Plane orientation α c

Conflicts of Interest
The orientation of the critical plane corresponding to D FS ∆γ The maximum shear strain range on a plane over a loading cycle σ n,max Maximum normal stress on a plane over a cycle σ surface RS at the surface σ max Peak RS y max Depth of σ max y core Depth where RS vanishes σ max,t The transition value for σ max Subscripts opt Values associated with the optimum RS distribution 1, 2, 3, . . .