**1. Introduction**

In mechanical systems, gears are widely used components to transmit torque and motion (i.e., mechanical power) between non-coaxial shafts [1]. Due to their working principle (i.e., meshing of conjugate profiles), teeth are subject to various damage mechanisms that can lead to the failure of the entire mechanical system [2,3]. Wear, scuffing, and (micro) pitting in the teeth flank are just a few examples of failure modes that, in turn, can be attributable to high contact pressures and/or insufficient lubrication [4–8]. However, Tooth (root) Bending Fatigue (TBF) is the most dangerous one [9,10].

TBF leads to the nucleation and propagation of a crack in the Tooth Root Radius (*ρf P*) due to the varying stress induced by tooth bending during meshing [11,12]. Therefore, a fundamental aspect to be considered while designing gears is the capability to withstand cyclic bending loads [13]. With this respect, different standards support the gear design to avoid TBF failures, e.g., ISO 6336-3 [14,15] and ANSI/AGMA [16]. The above-mentioned standards support gear design through the determination of Tooth Bending Strength (TBS).

**Citation:** Concli, F.; Maccioni, L.; Fraccaroli, L.; Bonaiti, L. Early Crack Propagation in Single Tooth Bending Fatigue: Combination of Finite Element Analysis and Critical-Planes Fatigue Criteria. *Metals* **2021**, *11*, 1871. https://doi.org/10.3390/ met11111871

Academic Editors: Filippo Berto, Ricardo Branco and Shengchuan Wu

Received: 28 October 2021 Accepted: 17 November 2021 Published: 21 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

According to the Method B of ISO 6336-3 [15], the maximum stress *σF* at the *ρf P* due to pure bending has to not exceed the permissible bending stress *σFP* that, in turn, is proportional to the material strength *σFlim*, i.e., a material property usually determined through Single Tooth Bending Fatigue (STBF) tests [17–20].

In STBF tests, teeth belonging to a gear sample made of the material to be characterized are loaded with two pulsating, competing, parallel, and discordant forces applied through two anvils (having parallel faces) mounted on a universal (fatigue) testing machine [21,22]. Exploiting the Wildhaber property [1], these forces are applied perpendicularly to two tooth flanks, resulting in tangent to the base circle [23]. STBF tests are interrupted if a tooth fails or if it withstands the run-out condition (106–108) cycles [24–28]. The statistical elaboration of failures and run-outs (the load levels are defined with the staircase approach [29,30]) leads to the determination of the load-carrying capacity [20]. The methods for translating this load into stresses acting at the *ρf P* (i.e., *<sup>σ</sup>Flim*) can be different: (1) through the reverse application of the standard (e.g., [31–34]), (2) by means of experimental measures (e.g., exploiting strain gauges [35]), (3) via Finite Element (FE) simulations (e.g., [36–39]).

With respect to the FE simulations, on the one hand, they allow us to obtain relevant information on the principal stresses in the *ρf P*, for each loading condition. On the other hand, results of FE simulation have to be further elaborated to estimate the fatigue behavior of a specific gear design [40,41]. In other words, numerical simulations of STBF tests and further data elaboration based on fatigue criteria using material data obtained via standard tests (i.e., torsion, bending and traction quasi-static and fatigue tests on standard specimens) seems to be a valuable alternative to long experimental campaigns.

In recent studies, the authors pointed out that the FE simulation results of the STBF configuration can be analyzed and elaborated via different fatigue criteria based on the critical plane approach [41–43]. These allow for evaluating the criticality of each point along the *ρf P*. Moreover, it permits to individuate the potential propagation direction of the crack after nucleation. Nevertheless, it has been observed that different fatigue criteria could lead to different results in terms of TBS and/or crack propagation direction [42].

The goal of the present paper is to evaluate the most appropriate fatigue criteria for characterizing the fatigue behavior in terms of the individuation of the nucleation point and the determination of the direction of early propagation of the crack in real mechanical components characterized by non-proportional multiaxial states of stress. This stress state, i.e., any state of time varying stress where the orientation of the principal axes changes with respect to a reference system integral with the studied component, can be found in gears [38]. In this respect, STBF tests described in [35,38–40] have been numerically reproduced and the FE results have been analyzed through different fatigue criteria based on the critical plane, i.e., Findley [44], Matake [45], McDiarmid [46], Papadopoulos [47], and Susmel et al. [48]. The outcomes of the elaboration have been compared with the cracks observed in the above-mentioned experimental campaigns [35,38–40].

The present paper is organized as follows. In Section 2, the mathematical elaboration of a generic time-dependent stress tensor *σ*(*t*) according to different fatigue criteria is presented. In Section 3, FE modeling, numerical data processing, and experimental data acquisition are shown. Comparison of numerical and experimental results are presented in Section 4. Discussions and conclusion can be found in Section 5.

### **2. Background: Mathematical Modeling of Fatigue Criteria Based on the Critical Plane**

In the present section, the mathematical modeling is presented of the main fatigue criteria based on critical plane (i.e., Findley [44], Matake [45], McDiarmid [46], Papadopoulos [47], and Susmel et al. [48]). Each fatigue criterion starts from the time-dependent stress tensor *σ*(*t*) (Equation (1)) referred to the point whose fatigue behavior has to be evaluated.

$$
\overline{\sigma}(t) = \begin{bmatrix}
\sigma\_{xx}(t) & \tau\_{xy}(t) & \tau\_{xz}(t) \\
\tau\_{yx}(t) & \sigma\_{yy}(t) & \tau\_{yz}(t) \\
\tau\_{zx}(t) & \tau\_{zy}(t) & \sigma\_{zz}(t)
\end{bmatrix} \tag{1}
$$

More specifically, a generic plane (including the point to be evaluated) can be defined by means of its normal vector *n* that, in turn, can be defined according with its spherical coordinates (*ϕ<sup>n</sup>*, *<sup>ϑ</sup>n*) with respect to a generic reference system (Figure 1).

**Figure 1.** Components of *Pn*(*ϕ<sup>n</sup>*, *ϑ<sup>n</sup>*, *t*) on the plane *<sup>n</sup>*(*ϕ<sup>n</sup>*, *<sup>ϑ</sup>n*).

According to Equation (2), it is possible to calculate the stress vector *Pn* acting on the afore-mentioned plane; the vector can, in turn, be decomposed into a normal component *σn* and into a tangential component *τn* (Figure 1).

$$P\_n(\varphi\_{n,\prime}\vartheta\_n, t) = \overline{\sigma}(t) \,\,\mathfrak{n}(\varphi\_{n,\prime}\vartheta\_n) \tag{2}$$

On the one hand, *σn* (which can be calculated through Equation (3)) presents a fixed direction and a time-dependent modulus. On the other hand, *τn* has a time-dependent modulus and direction. Therefore, *τn* has to be further decomposed into its components along the *u* and *v* directions (Figure 2). The unitary vectors *n*, *u*, *v* are defined as in Equation (4). In Equation (5), *τn* is defined.

$$
\sigma\_{\mathfrak{n}}(\varphi\_{\mathfrak{n}}, \theta\_{\mathfrak{n}}, t) = \mathfrak{n}^{\mathsf{T}}(\varphi\_{\mathfrak{n}}, \theta\_{\mathfrak{n}}) \overline{\overline{\sigma}}(t) \,\mathfrak{n}(\varphi\_{\mathfrak{n}}, \theta\_{\mathfrak{n}}) \tag{3}
$$

$$\mathfrak{u}(\boldsymbol{\varrho}\_{\boldsymbol{\theta}\boldsymbol{\eta}},\boldsymbol{\theta}\_{\boldsymbol{\eta}}) = \begin{bmatrix} -\sin\theta\_{\boldsymbol{\theta}\boldsymbol{\eta}} \\ \cos\boldsymbol{\varrho}\_{\boldsymbol{\theta}\boldsymbol{\eta}} \\ 0 \end{bmatrix} \boldsymbol{\mathfrak{v}}(\boldsymbol{\varrho}\_{\boldsymbol{\theta}\boldsymbol{\eta}},\boldsymbol{\theta}\_{\boldsymbol{\eta}}) = \begin{bmatrix} -\cos\theta\_{\boldsymbol{\theta}\boldsymbol{\eta}}\cos\boldsymbol{\varrho}\_{\boldsymbol{\theta}} \\ -\cos\theta\_{\boldsymbol{\theta}\boldsymbol{\eta}}\sin\boldsymbol{\varrho}\_{\boldsymbol{\theta}\boldsymbol{\eta}} \\ \sin\boldsymbol{\theta}\_{\boldsymbol{\eta}} \end{bmatrix}; \boldsymbol{\mathfrak{u}}(\boldsymbol{\varrho}\_{\boldsymbol{\theta}\boldsymbol{\eta}},\boldsymbol{\theta}\_{\boldsymbol{\theta}}) = \begin{bmatrix} \sin\theta\_{\boldsymbol{\theta}\boldsymbol{\eta}}\cos\boldsymbol{\varrho}\_{\boldsymbol{\theta}\boldsymbol{\eta}} \\ \sin\boldsymbol{\theta}\_{\boldsymbol{\theta}}\sin\boldsymbol{\varrho}\_{\boldsymbol{\theta}\boldsymbol{\eta}} \\ \cos\boldsymbol{\theta}\_{\boldsymbol{\theta}} \end{bmatrix} \tag{4}$$

$$\pi\_{\boldsymbol{\mathfrak{u}}}(\boldsymbol{\varrho}\_{\boldsymbol{\mathfrak{u}}\_{\boldsymbol{\mathfrak{t}}}}\boldsymbol{\vartheta}\_{\boldsymbol{\mathfrak{t}}},t) = \boldsymbol{\mathfrak{u}}^{\boldsymbol{\mathfrak{t}}}(\boldsymbol{\varrho}\_{\boldsymbol{\mathfrak{u}}\_{\boldsymbol{\mathfrak{t}}}}\boldsymbol{\mathfrak{t}}\_{\boldsymbol{\mathfrak{t}}}) \overline{\overline{\boldsymbol{\mathfrak{t}}}}(\boldsymbol{\mathfrak{t}}) \boldsymbol{\mathfrak{u}}(\boldsymbol{\varrho}\_{\boldsymbol{\mathfrak{t}}},\boldsymbol{\mathfrak{t}}\_{\boldsymbol{\mathfrak{t}}}) + \boldsymbol{\mathfrak{v}}^{\boldsymbol{\mathfrak{t}}}(\boldsymbol{\varrho}\_{\boldsymbol{\mathfrak{t}}},\boldsymbol{\mathfrak{t}}\_{\boldsymbol{\mathfrak{t}}}) \overline{\overline{\boldsymbol{\mathfrak{t}}}}(\boldsymbol{\mathfrak{t}}) \boldsymbol{\mathfrak{v}}(\boldsymbol{\varrho}\_{\boldsymbol{\mathfrak{t}}},\boldsymbol{\mathfrak{t}}\_{\boldsymbol{\mathfrak{t}}}) \tag{5}$$

**Figure 2.** *u* and *v* on the plane *<sup>n</sup>*(*ϕ<sup>n</sup>*, *<sup>ϑ</sup>n*) and definition of the curve Γ*<sup>n</sup>*.

It is worth noting that, for periodic stresses (having a period *T*), the point of the arrow of the vector *Pn* describes a closed tridimensional curve. Consequently, *τn* describes a closed curve in the plane (Figure 2). In the Figure, this curve is indicated as Γ*<sup>n</sup>*. On the one hand, the normal stress *σn* ranges from a minimum *<sup>σ</sup>n*,*min* to a maximum *<sup>σ</sup>n*,*max* value (Figure 2). Therefore, it is possible to define the value of the alternating stress (acting on the plane having normal *n*) as *<sup>σ</sup>n*,*<sup>a</sup>* defined according to Equation (6). On the other hand, to define the value of alternate tangential stress *<sup>τ</sup>n*,*<sup>a</sup>* (acting on the plane having normal *n*), literature reports different methods. The most diffused one is the Minimum Circumscribed

Circle (MCC) (Equation (7)) [49]. Considering that the curve Γ*n* is representative of the tangential stresses acting on the studied plane during the entire loading cycle, the MCC method suggests determining *<sup>τ</sup>n*,*<sup>a</sup>* as the radius of the smallest circle that can entirely contain the curve Γ*n* (Figure 3).

$$
\sigma\_{\mathfrak{n},\mathfrak{a}} = \max\_{\widetilde{T}} \{ \sigma\_{\mathfrak{n}}(t) \} - \min\_{\widetilde{T}} \{ \sigma\_{\mathfrak{n}}(t) \} = \sigma\_{\mathfrak{n},\max} - \sigma\_{\mathfrak{n},\min} \tag{6}
$$

$$\tau\_{\mathfrak{n},\mathfrak{a}} = M\_{\widetilde{T}} C \{ \tau\_{\mathfrak{n}}(t) \} \tag{7}$$

**Figure 3.** Minimum Circumscribed Circle (MCC) method.

By varying the spherical coordinates (*ϕ<sup>n</sup>*, *<sup>ϑ</sup>n*) systematically, it is possible to define a series of different planes passing through the point to be evaluated. For each plane it is possible to calculate the related stress parameters, i.e., *<sup>τ</sup>n*,*a*, *<sup>σ</sup>n*,*min*, *<sup>σ</sup>n*,*max*, and *<sup>σ</sup>n*,*a*. Based on these stress parameters, it is possible to individuate the critical plane having specific spherical coordinates (*ϕ<sup>c</sup>*, *<sup>ϑ</sup>c*). According to the Matake, the Susmel et al., the Papadopoulos, and the McDiarmid criteria, the critical plane is defined as the plane that displays the maximum value of *<sup>τ</sup>n*,*<sup>a</sup>* (Equation (8)).

$$(\mathcal{q}\_{\mathbb{C}'}, \theta\_{\mathbb{C}}) \to \max\_{\boldsymbol{\varrho}, \boldsymbol{\theta}} \{ \tau\_{\boldsymbol{n}, \boldsymbol{a}}(\boldsymbol{\varrho}, \boldsymbol{\theta}) \} \tag{8}$$

Conversely, according to the Findley criterion, the critical plane is defined as the plane that presents the maximum value of the damage parameter (*DP*) defined as in Equation (9). This is a function of the alternating tangential stress (*<sup>τ</sup>n*,*<sup>a</sup>*) and the maximum stress reached in a cycle (*<sup>σ</sup>n*,*max*). Therefore, the Findley criterion could lead to a critical plane having a different orientation with respect to the critical plane according to the other fatigue criteria.

$$\tau\_{\mathbf{t}}(\varphi\_{\mathbf{C},\mathbf{t}},\theta\_{\mathbf{C}}) \to \max\_{\boldsymbol{\uprho},\boldsymbol{\uprho}} \left\{ \tau\_{\mathbf{t},\mathbf{d}}(\boldsymbol{\uprho},\boldsymbol{\uptheta}) + \frac{2r\_{\pi/\sigma}-1}{2\left(\sqrt{r\_{\pi/\sigma}-r\_{\pi/\sigma}^{2}}\right)} \sigma\_{\mathbf{t},\max}(\boldsymbol{\uprho},\boldsymbol{\uptheta}) \right\} \tag{9}$$

where *<sup>r</sup>τ*/*σ* is the ration between the material fatigue limit at symmetrical alternating torsional loading (*<sup>τ</sup>f*) and material fatigue limit at symmetrical alternating bending loading (*<sup>σ</sup>f*) as in Equation (10). It is worth noticing that these material properties can be estimated through simple fatigue tests.

$$
\tau\_{\tau/\sigma} = \tau\_f / \sigma\_f \tag{10}
$$

Once the critical plane (*ϕ<sup>c</sup>*, *<sup>ϑ</sup>c*) has been identified, the various criteria require that the damage parameter on this plane be calculated. In the present paper, the stress parameters related with the critical plane are labeled with the subscript *c*, i.e., *<sup>τ</sup>c*,*a*, *<sup>σ</sup>c*,*max*, *<sup>σ</sup>c*,*a*.

The various criteria differ on how the damage parameter (*DP*) is calculated. According to Findley, the damage parameter (*DPFindley*) is defined as in Equation (11). According to the Matake criteria, the *DPMatake* is affected by the alternating (tangential) stress *<sup>σ</sup>c*,*<sup>a</sup>* (*<sup>τ</sup>c*,*<sup>a</sup>*) (Equation (12)). The Susmel et al. criteria requires calculating the damage parameter (*DPSusmel et al*.) according to Equation (13). With respect to the criteria proposed by McDiarmid, it is necessary to consider the ultimate tensile stress *σR* (Equation (14)).

$$DP\_{\text{Finding}} = \tau\_{c,\mu} + \frac{2r\_{\tau/\sigma} - 1}{2\left(\sqrt{r\_{\tau/\sigma} - r\_{\tau/\sigma}^2}\right)} \sigma\_{c,\text{max}} \tag{11}$$

$$DP\_{\text{Matake}} = \pi\_{\text{c,a}} + (2r\_{\pi/\sigma} - 1)\sigma\_{\text{c,a}} \tag{12}$$

$$DP\_{\text{Susnel et al.}} = \tau\_{\text{c,a}} + \left(\tau\_f - \frac{\sigma\_f}{2}\right) \frac{\sigma\_{\text{c,max}}}{\tau\_{\text{c,a}}} \tag{13}$$

$$DP\_{\text{McDiarmid}} = \pi\_{c,a} + \frac{\pi\_f}{2\sigma\_R} \sigma\_{c,\text{max}} \tag{14}$$

To implement the Papadopoulos' criteria, it is important to define the maximum octahedral stress *<sup>σ</sup>h*,*max* (in the time window *T*). It can be calculated through Equation (15), where *σO* is a vector with the principal stresses, i.e., the stresses that, for the same time instant *t*, satisfies Equation (16). *I* is the identity matrix. The *DPPapadopoulos* can be calculated according to Equation (17).

$$
\sigma\_{h,max} = \max\_{T} \left\{ \frac{1}{3} \sum\_{i=1,2,3} \sigma\_{O\_i} \right\} \tag{15}
$$

$$\det \left| \overline{\overline{\sigma}}(\mathbf{t}) - \sigma\_O \overline{I} \right| = 0 \tag{16}$$

$$DP\_{Popadopulos} = \pi\_{c,a} + \left(\frac{3}{2}(2r\_{\pi/\sigma} - 1)\right)\sigma\_{b,\text{max}}\tag{17}$$

Eventually, each fatigue criteria (based on critical plane) state that the component works safely as long as the value of the damage parameter, in each point, is below a specific threshold. Therefore, it is possible to calculate a safety factor (*SF*) (for each criterion and in each position) which formulation depends on the implemented criterion. For example, *SFFindlay* is defined in Equation (18). In Equations (19)–(22) the *SF* for the Matake, Susmel et al., McDiarmid, and Papadopoulos criteria can be found, respectively. *SF* > 1 means that the analyzed stress state has not reached the critical value according to the studied criterion (and vice versa for *SF* < 1).

$$S\_{Finding} = \frac{\frac{\tau\_f}{2\left(\sqrt{r\_{\pi/\sigma} - r\_{\pi/\sigma}^2}\right)}}{DP\_{Finding}}\tag{18}$$

$$S\_{F\_{Matake}} = \frac{\tau\_f}{DP\_{Matake}}\tag{19}$$

$$S\_{F\_{\text{Susnel et al.}}} = \frac{\mathbf{r}\_f}{DP\_{\text{Susnel et al.}}}\tag{20}$$

$$S\_{F\_{McDiramid}} = \frac{\tau\_f}{DP\_{McDiramid}}\tag{21}$$

$$S\_{F\_{\text{Papadopvolos}}} = \frac{\mathfrak{r}\_f}{DP\_{\text{Papadopvolos}}} \tag{22}$$
