*Article* **Empirical Solution of Stress Intensity Factors for the Inclined Inner Surface Crack of Pipe under External Pressure and Axial Compression**

**Xi-Ming Yao 1,2, Yu-Chen Zhang 1,2, Qi Pei 1,2, Li-Zhu Jin 1,2, Tian-Hao Ma 1,2, Xiao-Hua He 1,2 and Chang-Yu Zhou 1,2,\***


**Abstract:** Based on fracture mechanics theory, a finite element method was used to determine the stress intensity factors of the inclined crack on the inner surface of the pipe under axial compression load and external pressure. The effects of different influencing factors on the stress intensity factor along the crack front considering crack closure were systematically explored, which were different to those under internal pressure. The effects of high aspect ratio on *KII*, the crack inclination asymmetry caused by curvature and the effects of the friction coefficient on the stress intensity factors of the pipe with an inclined inner surface crack under axial compression load and external pressure were explored in this paper. To be fit for defect assessment, the solutions for stress intensity factors *KII* and *KIII* were derived, and new correction factors *f<sup>θ</sup>* and *f<sup>μ</sup>* were proposed in the empirical solutions to accommodate the crack inclination asymmetry and the friction coefficient, respectively.

**Keywords:** inclined surface crack; stress intensity factor; pipe; crack closure; external and axial pressure; finite element analysis

**Citation:** Yao, X.-M.; Zhang, Y.-C.; Pei, Q.; Jin, L.-Z.; Ma, T.-H.; He, X.-H.; Zhou, C.-Y. Empirical Solution of Stress Intensity Factors for the Inclined Inner Surface Crack of Pipe under External Pressure and Axial Compression. *Materials* **2023**, *16*, 364.

https://doi.org/10.3390/ma16010364 Academic Editor: Roberto Citarella

Received: 25 November 2022 Revised: 19 December 2022 Accepted: 27 December 2022 Published: 30 December 2022

**Copyright:** © 2022 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/).

#### **1. Introduction**

Pipes are widely used in various fields of petroleum, chemical industry, and natural gas. In engineering practice, there are a large number of shells subjected to external pressure loads such as submarines, aerospace simulators, deep-sea pipelines, and buried pipelines. Due to the limitations of the manufacturing process, the structures may have tiny defects such as cracks. Lin and Smith [1] found that arbitrarily shaped cracks become semi-elliptical after a certain period of expansion. The stress intensity factor is an important fracture parameter for judging crack propagation and failure in defect assessment. Therefore, since the last century, scholars have conducted many studies on this subject.

In the early years, Raju and Newman [2] established a stress intensity factor empirical solution for mode I internal and external surface cracks of cylindrical vessels through FEA. The formula gave the influence coefficients based on different shell and crack sizes. Navid and Fenner [3] calculated the stress intensity factors of a thick-walled cylinder with internal pressure for different crack sizes based on the boundary integral equation method. Afterward, various numerical techniques were developed to calculate the stress intensity factor of pipe surface cracks such as the line spring model [4] and the weight function method [5–7]. Kamaya [8] combined the finite element alternating method and the finite element analysis to calculate the stress intensity factors of shell surface cracks. Wallbrink [9] proposed a semi-analytical method of conformal transformation to predict the stress intensity factors of circumferential surface cracks. However, most of these studies were limited to simple loading and the crack direction was axial or circumferential. In fact, in most cases, the loading conditions of the pipelines are complex and the crack direction is inclined.

The research on mixed mode cracks first started from a flat plate. Murakami [10] analyzed *KI*, *KII*, *KIII* by the body force method, and provided a formula for the maximum value *Kθmax* at the crack front according to the projected area of the crack in the direction of the maximum principal stress. Zeng and Dai [11] proposed a simplified analytical model of an inclined surface crack under biaxial stress, and proposed a closed solution for *KI* and *KIII* at the deepest point of the crack front. Some scholars have also used FEA [12,13], the numerical influence function method [14], and the experimental method [15] to study the mixed mode cracks of flat plates. The study of mixed mode cracks in shell can be divided into two types. First, the crack direction is determined, usually axial or circumferential, and the mixed mode crack is formed by applying a complex far-field load on the structure. Shahani and Habibi [16] considered the stress intensity factor of cylinders with a circumferential crack under the action of axial force, bending moment, and torque. Predan et al. [17] studied the stress intensity factor for circumferential semi-elliptical surface cracks in a hollow cylinder subjected to pure torsion. Ramezani et al. [18] demonstrated the empirical solution of the stress intensity factor of cylinder surface crack under pure torsion. Second, the mixed mode cracks were formed by changing the inclination angle of the surface cracks. For example, Li et al. [19,20] studied the stress intensity factor of the inclined crack in the pipeline under the far field tension and tension-bending, and provided an empirical formula to calculate the complete value of the stress intensity factor at the crack front through the influence coefficient under different aspect ratios. Li and Mao [21] calculated the stress intensity factor of the inclined surface crack of the outer wall of the heat exchanger by numerical simulation.

In the research on cracks, it is inevitable to encounter the problem of complete or partial closure of the cracks. For simplicity, more papers have chosen to ignore the surface contact of the cracks. However, the contact closure of the crack face significantly changes the stress–strain distribution at the crack front, thereby affecting the fracture behaviors of the structure. In fact, since the contact area is time-varying, the contact problem is a typical nonlinear problem with boundary conditions. In static fracture mechanics, the contact friction of the crack face has an obvious influence on the stress intensity factor. At present, there have been few studies on mixed mode cracks on the friction surfaces under compressive loading, which have mainly focused on flat plates [22–24]. Bowie and Freese [25] proposed a crack closure technique to correct the solution of overlapping crack surfaces, but their method does not consider the sliding between crack surface, and is only suitable for the case of a large friction coefficient. Liu and Tan [26] used the boundary element method to study the effect of the interaction between the friction and crack surface on the stress intensity factor. Hammouda et al. [27,28] used finite element analysis to study the effect of crack surface friction and crack inclination on the stress intensity factor of the central cracked plate under unidirectional compressive load. In addition, Dorogoy and Banks-Sills [29] investigated the effects of loading angle and friction coefficient on the stress intensity factor and the crack length of Brazilian disc cracks under concentrated loading using a finite difference solution.

To the best of the authors' knowledge, researchers have not focused on the surface crack stress intensity factor of shell structures subjected to compressive loads, however, in engineering practice, fracture damage occurs in shell structures subjected to simultaneous external and axial pressures. In this paper, the stress intensity factors of the shell subjected to simultaneous external and axial pressures were evaluated by using the three-dimensional finite element analysis method, and the effects of the relative depth of the surface cracks, aspect ratio, crack inclination, and friction coefficient on the stress intensity factor of the crack fronts were investigated, which were different from those under internal pressure. The empirical solutions of the mode II and mode III stress intensity factors along the crack front are given using numerical methods, respectively.

#### **2. Finite Element Model of Pipe Crack**

#### *2.1. Meshing of the Pipe with Inclined Crack on the Inner Surface*

For the study of mixed mode cracks, the most widely used and accurate method is FEA, which involves the extended finite element method (XFEM) and contour integral. Based on the static loading model in this paper, the fluctuation of results caused by inaccurate meshing can be better avoided by using contour integrals. The commercial code ABAQUS [30] was used for FEA. The material used in this paper was TA2 with a modulus of elasticity of 101,901 MPa and a Poisson's ratio of 0.348. The pipe size was fixed at *Ri* = 100 mm, *R*<sup>0</sup> = 110 mm, *t* = 10 mm. In order to avoid the influence of the boundary effect on the stress intensity factor, the pipe length should be six times that of *Ri* and 20 times that of the crack half-length c, *l* = 600 mm [31,32]. The crack size is described by dimensionless parameters: the relative crack depth (*a/t*) and the ratio of long and minor axes of semi-elliptical cracks (*a/c*). The crack inclination angle is *θ*. The pipe model is shown in Figure 1a. The external pressure of the pipe *P*<sup>0</sup> = 100 MPa and the axial pressure *p* = 525 MPa. The boundary conditions were that one end of the pipe was completely fixed and the other end restricted its circumferential constraint, as shown in Figure 1b.

**Figure 1.** (**a**) Pressure pipe with an inclined crack on the inner surface. (**b**) Loads and boundary conditions.

The current general semi-elliptical crack meshing method was used by the authors [31,33]. The six-node triangular prism element C3D6 was used to sweep along the path around the crack tip to simulate the stress–strain singularity at the crack front. C3D6 belongs to a fully integrated element with linear interpolation in all directions, the number of nodes is 6, and

the number of integration points is 6. The eight-node hexahedral element C3D8R was used in the nearby area. C3D8R belongs to a linear reduced integral element with eight nodes and only one integral point exists in the center of the cell, which is equivalent to a constant stress cell. The reduced integration element can effectively avoid stress discontinuity problems. The meshes at the crack were subdivided, the mesh size was controlled within 1 mm, and the mesh size of the rest zone was controlled within 10 mm to reduce the calculation time. The mesh division is shown in Figure 2.

**Figure 2.** Meshing of the pipe with the inclined internal surface crack:(**a**) global view of the cracked pipe; (**b**) equal divisions at the crack front; (**c**) hexahedron elements around the crack front; (**d**) wedgeshaped elements at the crack front.

#### *2.2. Penalty Function Method for Contact Problem*

Since the pipe is not only subjected to external pressure but also to axial compressive load, contact should be set between the crack surfaces. Generally speaking, the crack surfaces in the contact state have three characteristics:


In ABAQUS, the penalty function method, Lagrange multiplier method, and static friction–kinetic friction index decay method can be used to solve the contact problem. In this paper, a penalty function was used to characterize the tangential friction force between the crack surfaces. The penalty function method requires the normal and tangential friction coefficients, which is similar to setting a "spring" between the contact surfaces, and the "spring" works only when the contact surfaces are closed. In this way, the normal contact pressure can be expressed as:

$$P\_n = \begin{cases} 0 & \mu\_n > 0\\ \
k\_n \mu\_n & \mu\_n \le 0 \end{cases} \tag{1}$$

According to Coulomb's law of friction, the frictional stress at the contact surface can be written as:

$$\begin{aligned} \pi\_{\mathfrak{s}} &= k\_{\mathfrak{s}} \mu\_{\mathfrak{s}} - \mu P\_{\mathfrak{n}} \\ \pi\_{\mathfrak{t}} &= k\_{\mathfrak{t}} \mu\_{\mathfrak{t}} - \mu P\_{\mathfrak{n}} \end{aligned} \tag{2}$$

where *t, n, s* represent the tangential, normal and sub-normal directions, respectively, as shown in Figure 3. The advantage of the penalty function is that it does not increase the degree of freedom of the problem. The disadvantage is that when the friction coefficient is too large, the convergence will become more difficult. The values of the friction coefficient *μ* in this paper were 0, 0.2, and 0.4, respectively.

**Figure 3.** Local coordinate system of the quarter-elliptic crack front.

#### *2.3. Mesh Independence Verification*

Raju and Newman [2] first proposed a formula to calculate the stress intensity factor of a mode I crack in a pipe, and Chun-Qing Li [19] extended this formula to calculate the mixed-mode stress intensity factor of an inclined crack, as follows:

$$\mathbf{K} = \sigma\_0 \sqrt{\pi \frac{a}{Q}} \mathbf{F} \left( \frac{a}{t}, \frac{a}{c}, \frac{t}{R\_i}, \xi, \theta \right) \tag{3}$$

where *σ*<sup>0</sup> = *PR/t* is the average hoop stress of the pipe; *K* and *F* are the stress intensity factor and influence coefficient function under different cracking modes, respectively; *Q* is the shape factor obtained from the second type of elliptic integral; and the empirical formula is given by Shiratori et al. [34].

$$\mathbf{K} = \left\{ \begin{array}{c} \mathbf{K}\_{I} \\ \mathbf{K}\_{II} \\ \mathbf{K}\_{III} \end{array} \right\} \tag{4}$$

$$F\left(\frac{a}{t}, \frac{a}{c}, \frac{t}{R\_i}, \xi, \theta\right) = \left\{ \begin{array}{c} F\_I \\ F\_{II} \\ F\_{III} \end{array} \right\} \tag{5}$$

$$\begin{array}{ll} Q = 1 + 1.464 \left( \frac{a}{c} \right)^{1.65} & \frac{a}{c} \leqslant 1\\ Q = \left[ 1 + 1.464 \left( \frac{c}{d} \right)^{1.65} \right] \left( \frac{a}{c} \right)^2 & \frac{a}{c} > 1 \end{array} \tag{6}$$

In order to unify the evaluation indicators, it is necessary to normalize the stress intensity factor and the position along the crack front. In this paper, the normalization factor *K*<sup>0</sup> was defined as [18]:

$$K\_0 = \sigma \sqrt{\pi \cdot a\_0} \quad a\_0 = 1 \tag{7}$$

where *σ* = *PR*/2*t* is the axial compressive stress. The position along the crack front can be normalized according to the number of equal meshes. The advantage of this normalization method is that *K0* becomes a constant, which can truly reflect the variation trend of *K*.

Figure 4 shows the normalized stress intensity factors at the crack front of three kinds of mesh sizes at different angles. The mesh size is shown in Table 1. It can be seen that the results of the first two mesh sizes were basically the same. A comparison with the normalized *KII* and *KIII* results with Chun-Qing Li [19] after changing the load using mesh2 is given in Figure 5. The average differences were 8.98% and 2.55%, respectively, which demonstrates the accuracy of the results.

**Figure 4.** Effect of the mesh density on the normalized SIF at different angles, *a/t* = 0.2, *a/c* = 0.2. (**a**) Normalized *KII*, (**b**) normalized *KIII*.

**Table 1.** Crack front meshing.


**Figure 5.** Comparison of the stress intensity factors with mixed modes, *a/t* = 0.2, *a/c* = 0.4, *θ* = 45◦. (**a**) Normalized *KII*, (**b**) normalized *KIII* [19].

#### *2.4. Comparison of SIFs between Crack Opening and Crack Closing*

In order to compare the change in the stress intensity factor for crack opening and crack closing, two sets of cracks with different angles were selected in this paper. The contact surface friction coefficient *μ* = 0. The pipe load for crack opening was changed to the internal pressure and axial tension, and the results of the comparison are shown in Figure 6. To facilitate comparison of the differences, the normalized stress intensity factors under external pressure and axial compression were expressed symmetrically about the transverse coordinate. It can be seen that the absolute values of the normalized stress intensity factors at crack closure were obviously higher for both the mode II crack and mode III crack. For this phenomenon, the authors used the stress field calculation method for the plastic zone of the crack tip mentioned by Shlyannikov and Tumanov [35] to convert the three stress states obtained from the finite element calculations to give the sub-normal stresses and tangential stresses at the crack front in these two cases, corresponding to mode II and mode III cracks, respectively, as shown in Figure 7a,b, and it can be seen that *τ<sup>s</sup>* and *τ<sup>t</sup>* at the crack front of the pipe under axial compression and external pressure were higher than those under axial tension and internal pressure, which is the fundamental reason for the difference in SIFs. This difference also proves that it is necessary to study the stress intensity factor of cracks under axial compression load and external pressure.

**Figure 6.** Comparison of SIFs under tension and compression loads with different angles, *μ* = 0, *a/t* = 0.2, *a/c* = 0.8. (**a**) Normalized *KII*, (**b**) normalized *KIII*.

**Figure 7.** Comparison of shear stress along the crack fronts under tension and compression loads. (**a**) Shear stress along the crack front of mode II. (**b**) Shear stress along the crack front of mode III.

#### **3. Results and Discussion**

The stress intensity factor for mode I crack is 0 under the compression load [36,37]. However, due to the crack inclination and external pressure, the *KII* and *KIII* generated by the interaction of tangential and frictional forces at the crack front cannot be neglected. In this paper, the SIFs of pipe surface crack fronts were obtained by FEA under different influencing factors, which included the relative depth (*a/t*), the aspect ratio (*a/c*), the crack

inclination angle (*θ*), and the friction coefficient (*μ*). The specific calculation scheme is shown in Table 2. A total of 336 finite element models were calculated in this paper.


**Table 2.** Finite element calculation scheme.

#### *3.1. Effects of Crack Geometry and Inclination Angle*

Figure 8 shows the variation in the mode II and mode III SIFs with relative depth. The long axes of the semi-elliptical cracks were fixed as *c* = 5 mm and *c* = 10 mm, respectively, the friction coefficient *μ* was 0. For the convenience of presentation, only the comparison results of the maximum value along the crack front are given, and the rules were consistent for the remaining positions. It can be seen that the absolute value of the mode II and mode III stress intensity factors increased with the deepening of the crack, and the increase in the absolute values was more obvious when the crack size was larger. The surface crack was most dangerous at *θ* = 45◦ and became safer with the offset of crack inclination to the axial and circumferential directions.

**Figure 8.** Effect of the relative depth on the normalized SIFs, *μ* = 0. (**a**) Normalized *KII* (*c* = 5 mm, *a/c* = 0.4, 0.8, 1.6); (**b**) normalized *KIII* (*c* = 5 mm, *a/c* = 0.4, 0.8, 1.6); (**c**) normalized *KII* (*c* = 10 mm, *a/c* = 0.2, 0.4, 0.8); (**d**) normalized *KIII* (*c* = 10 mm, *a/c* = 0.2, 0.4, 0.8).

The variation in the mode II and mode III SIFs with the aspect ratio is given in Figure 9, and the minor axis of the semi-elliptical crack was fixed as *a* = 4 mm and *a* = 6 mm, respectively, the friction coefficient *μ* was 0. For the mode II cracks, when *a/c* < 1, the absolute value of the normalized *KII* increased with the aspect ratio, and when *a/c* > 1, the absolute value of the normalized *KII* decreased with the increase in the aspect ratio. The *KII* values of the crack fronts for different aspect ratios at two depths are given in Figure 9e,f. For semi-elliptical cracks with lower aspect ratios, the maximum value of *KII* always appeared at the surface point of the crack front. However, for semi-elliptical cracks with higher aspect ratios, the maximum value of *KII* appeared near the surface point, which led to a smaller value at the surface point. These variations are consistent with the results of Yang et al. [38] for high aspect ratio cracks. For mode III cracks, the absolute value of *KIII* decreased monotonically with the increase in aspect ratio, with the most pronounced decrease at *θ* = 45◦.

**Figure 9.** Effect of the aspect ratio on the normalized SIFs (*μ* = 0). (**a**) Normalized *KII* (*a* = 4 mm, *a*/*t* = 0.4); (**b**) normalized *KIII* (*a* = 4 mm, *a/t* = 0.4); (**c**) normalized *KII* (*a* = 6 mm, *a*/*t* = 0.6); (**d**) normalized *KIII* (*a* = 6 mm, *a/t* = 0.6); (**e**) normalized *KII* near the crack surface (*a* = 4 mm, *a*/*t* = 0.4); (**f**) normalized *KII* near the crack surface (*a* = 6 mm, *a*/*t* = 0.6).

The variation in the mode II and mode III SIFs with the crack inclination angle is given in Figure 10, the friction coefficient *μ* was 0. It can be seen that in all cases, the absolute values of *KII* and *KIII* were symmetrical about the deepest point. The maximum value of the stress intensity factor always appeared at the surface point or the deepest point at *θ* = 45◦, and the SIF value was smaller when the crack inclination was closer to the axial and circumferential directions, which was due to the more obvious closure effect under compression load, and the tangential component parallel to the crack surface became smaller. It should be noted that for small-size cracks, *KII* and *KIII* were almost identical when *θ* was symmetric about 45◦ (*θ* = 30◦, 60◦ and *θ* = 15◦, 75◦), and the pipe surface cracks at this time could be approximated as flat plate cracks, which has been verified in the study of mixed mode cracks in flat plates and pipes [13,20]. However, for large size pipe surface cracks, the difference caused by the angle was more obvious. The crack surface had a different curvature on the inner wall surface of the pipe at different inclination angles, which made the stress intensity factors no longer consistent for *θ* = 30◦, 60◦ and *θ* = 15◦, 75◦. As shown in Figure 10e,f, this difference due to curvature was more pronounced in spherical shells [39].

**Figure 10.** Effect of the inclination angle on the normalized SIFs, *μ* = 0. (**a**) Normalized *KII* (*a/t* = 0.2, *a/c* = 1.6); (**b**) normalized *KIII* (*a/t* = 0.2, *a/c* = 1.6); (**c**) normalized *KII* (*a/t* = 0.4, *a/c* = 0.8); (**d**) normalized *KIII* (*a/t* = 0.4, *a/c* = 0.8); (**e**) normalized *KII* (*a/t* = 0.6, *a/c* = 0.4); (**f**) normalized *KIII* (*a/t* = 0.6, *a/c* = 0.4).

#### *3.2. Effect of Friction Coefficient on Contact Surface*

The variation in the mode II and mode III SIFs with the friction coefficient is given in Figure 11. Due to the symmetry of the curves, only the results for the half-length along the crack front are shown in the figure. The absolute values of *KII* and *KIII* decreased with an increasing friction coefficient for different crack sizes and crack inclination angles. At *θ* = 45◦, the friction coefficient had the greatest effect on the SIFs, and this effect also increased as the crack size increased. For the mode II cracks, the friction coefficient had the greatest effect on the SIF at the surface point, and this effect decreased as the location of the crack front moved forward, reaching a minimum at the deepest point. In contrast to mode II cracks, the friction coefficient significantly changed the SIF at the deepest point of mode III, but had a dropping effect on the surface point.

**Figure 11.** Effect of the coefficient of friction on normalized SIFs. (**a**) Normalized *KII* (*a/t* = 0.2, *a/c* = 1.6); (**b**) normalized *KIII* (*a/t* = 0.2, *a/c* = 1.6); (**c**) normalized *KII* (*a/t* = 0.6, *a/c* = 0.4); (**d**) normalized *KIII* (*a/t* = 0.6, *a/c* = 0.4).

#### **4. Empirical Solution of SIFs for the Pipe with Inclined Inner Surface Cracks under External Pressure and Axial Compression**

Although FEA is one of the most effective methods for calculating SIFs, it is difficult to apply in engineering practice due to the complexity of modeling and time-consuming calculations. In this paper, a new form of empirical solution was proposed by least-squares fitting based on the FEA results to give the influence coefficients at different *a/t* with the normalized crack front position *ξ* as the basis function. A new correction factor for the inclination angle, *fθ*, was proposed based on the effect caused by the curvature of the pipe above-mentioned in the effect of the crack inclination angle. The friction coefficient influence coefficient, *fμ*, was proposed because the pipe is subjected to external pressure and axial pressure. The influence coefficients are shown in Tables 3 and 4, and the fitted equations are shown in Equations (8)–(16).

$$\mathbf{K} = \sigma \sqrt{\pi \frac{a}{Q}} \mathbf{F} \left( \frac{a}{t}, \frac{a}{c}, \tilde{\xi}, f\_{\theta'} f\_{\mu} \right) \tag{8}$$

where *σ* is the far-field compressive stress, *K* = *KI I KIII* ,*F a t* , *a <sup>c</sup>* , *ξ*, *fθ*, *f<sup>μ</sup>* = *FI I FIII* .


**Table 3.** The sub-curve fitting constants in Equation (10) *hi*.

**Table 4.** The sub-curve fitting constants in Equation (14) *hi*.


Influence coefficient functions for *KII*:

$$F\_{II} = \left(H\_1 + H\_2 \mathfrak{J} + H\_3 \mathfrak{J}^2 + H\_4 \mathfrak{J}^3\right) f\_\theta f\_\mu \tag{9}$$

$$H\_{i(i=1,2,3,4)} = h\_1 + h\_2(\frac{a}{c}) + h\_3(\frac{a}{c})^2 + h\_4(\frac{a}{c})^3 \tag{10}$$

$$f\_{\theta} = \left\{ 1 + \frac{ac}{600} [\sin(\theta - 45)] \right\} \sin 2\theta \tag{11}$$

$$f\_{\mu} = (1 + 0.575\mu)(ac)^{-0.505\mu} \tag{12}$$

Influence coefficient functions for *KIII*:

$$F\_{III} = \left(H\_1 + H\_2 \mathfrak{J} + H\_3 \mathfrak{J}^2\right) f\_\theta f\_\mu \tag{13}$$

$$H\_{i(i=1,2,3)} = h\_1 + h\_2(\frac{a}{c}) + h\_3(\frac{a}{c})^2 + h\_4(\frac{a}{c})^3 \tag{14}$$

$$f\_{\theta} = \left\{ 1 + \frac{ac}{1500} [\sin(45 - \theta)] \right\} \sin 2\theta \tag{15}$$

$$f\_{\mu} = (1 + 0.281\mu)(ac)^{-0.511\mu} \tag{16}$$

The fitting results for *θ* = 45◦ and *μ* = 0 at the three sizes are given in Figure 12a,b. The results of the fit at different angles for the same crack size are given in Figure 12c,d, where *μ* = 0. The fitting results for the same crack size with different angles and different friction coefficients are given in Figure 12e,f; due to the symmetry of the curves, only the results for the half-length of the crack front are shown in the figure. It can be seen that the fitting results achieved good agreement with the FEA results, which can be applied for the defect assessment.

**Figure 12.** Comparison of SIFs fitting results with (**a**,**b**) different crack sizes (*a/t* = 0.2, 0.4, 0.6, *a/c* = 1.6, 0.8, 0.4, *θ* = 45◦, *μ* = 0); (**c**,**d**) different inclination angles (*a/t* = 0.6, *a/c* = 0.4, μ = 0); (**e**,**f**) different friction coefficients (*a/t* = 0.4, *a/c* = 0.8, *μ* = 0, 0.2, 0.4).

#### **5. Conclusions**

In this paper, the SIFs of the inclined crack on the inner surface of the pipe under combined axial pressure load and external pressure were investigated by the finite element method, and the main conclusions are as follows.

(1) By comparing the SIF values along the crack fronts, it can be found that *KII* and *KIII* for the closed surface cracks under external and axial pressure were higher than *KII* and *KIII* for the open surface cracks under internal pressure and axial tension.

(2) The effects of the relative depth and aspect ratio on the SIFs of the inclined inner crack were analyzed. The results showed that the larger the *a* and *c*, the easier it is for the crack to expand, and when *a/c* >1, *KIImax* does not appear at the surface point of the crack, but near the surface point.

(3) The effects of the crack inclination angle on the SIFs of the inclined inner crack were analyzed. The results showed that the SIFs at the crack front were the largest at *θ* = 45◦, and their values decreased with an inclination angle toward the axial and circumferential directions. The larger the crack size, the more obvious the asymmetry of the SIFs about *θ* = 45◦.

(4) The effects of the friction coefficient on the SIFs of surface crack were analyzed. The larger the friction factor, the smaller the SIFs along the crack front. The friction coefficient had the greatest effect on the surface point of mode II cracks and the deepest point of mode III cracks, and the greatest effect on cracks with *θ* = 45◦.

(5) Based on the above results, a new solution for stress intensity factors *KII* and *KIII* were proposed, and the corresponding coefficients for different crack sizes, an angle correction factor *fθ*, and a friction correction factor *f<sup>μ</sup>* are given.

**Author Contributions:** X.-M.Y.: Conception, methodology, software, writing—manuscript, and designed the research. Y.-C.Z., Q.P. and L.-Z.J.: Software. T.-H.M.: Investigation. X.-H.H.: Funding acquisition. C.-Y.Z.: Formal analysis, conception and funding acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper was funded by the National Natural Science Foundation of China (51975271).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (51975271).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**



#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Analysis of the Vehicle Chassis Axle Fractures**

**Živile Decker ˙ 1, Vitalijus Rudzinskas 1, Kazimierz Drozd 2, Jacek Caban 2,\*, Jurijus Tretjakovas 3, Aleksander Nieoczym <sup>2</sup> and Jonas Matijošius <sup>4</sup>**


**Abstract:** With the rapid development of the road transport industry, trucks with semi-trailers have become the main means of transporting goods by road. High quality, durability and reliability of the construction are the main requirements for the production of trailers. Trailer and semi-trailer axles are one of the main and most important components of a truck. Due to the fact that semi-trailer axles are subjected to additional static and dynamic loads during operation, their proper construction is extremely important, therefore they should be carefully designed and tested. The durability of the suspension components refers to the duration of the onset of fatigue. This article presents an analysis of damage to the rear axle of the semi-trailer using macroscopic observations of the damage site and dynamic FEA of stress distribution in the axle material. In order to identify the probable cause of the damage, eight cases of loading the semi-trailer axle were considered. Analytical solutions have shown that in various cases the yield point is exceeded and the strength limit of the modeled semi-trailer axle is reached. The risk of damage to the vehicle's suspension system components increases on poor roads (bumps and winding road sections).

**Keywords:** failure analysis; FEA; fracture mechanics; macroscopic research; semi-trailers

#### **1. Introduction**

With the rapid development of the road transport industry, trucks with semi-trailers have become the main means of transporting goods by road. The economic goal of all transport companies is to transport as many goods as possible and with as few journeys as possible. High quality, durability, and structural reliability are the main requirements for the production of trailers. The axles of trailers and semi-trailers are one of the main and most important elements, which must be carefully designed and tested experimentally under static and dynamic loads, as the axle is subjected to additional loads in the event of road roughness or off-road [1]. The durability of the suspension parts refers to the duration of the onset of fatigue, defined as the number of cycles up to a certain component cracking length under cyclic loads [2,3].

The reliability of individual safety systems of a given means of transport translates directly into road traffic safety [4–8]. Much attention is also paid to the diagnostics of individual vehicle systems that affect safety. It is worth mentioning here the research presented by Savchenko et al. [9] and Gnap et al. [10], who demonstrated the possibility of using MEMS sensors in vehicles. Hudec et al. [11], Gajek [12] and Tucki et al. [13] paid attention to the monitoring of the technical condition of vehicles. On the other hand, research on the impact of heavy goods vehicle load during braking is presented in detail in [14–17].

**Citation:** Decker, Ž.; Rudzinskas, V.; Drozd, K.; Caban, J.; Tretjakovas, J.; Nieoczym, A.; Matijošius, J. Analysis of the Vehicle Chassis Axle Fractures. *Materials* **2023**, *16*, 806. https:// doi.org/10.3390/ma16020806

Academic Editors: Grzegorz Lesiuk and Dariusz Rozumek

Received: 9 December 2022 Revised: 5 January 2023 Accepted: 6 January 2023 Published: 13 January 2023

**Copyright:** © 2023 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/).

The topic of scientific information is the determination of the causes of vehicle chassis axle defects and fractures, which can be found in a wide variety of ways. The literature begins to address this problem along with the first defects in the chassis axles. A review of recent research shows that researchers are looking to find the most optimal design solutions based on durability, lightness of construction, strength, and low cost. It also examines the causes of defects by evaluating the chemical composition of the metal, microstructure, mechanical properties, acting on loads, the location of the object under investigation in the whole system, and so on.

One of the most popular test methods used to detect cracks, fractures, and other defects is fractographic testing. Such investigations determine the complex causes of decay: the influence of the shape of the part, the load conditions, the interaction of the elements, the peculiarities of the structure of the material, and other factors are determined from the fractures.

In addition to traditional qualitative fractography, quantitative fractography aims to measure the topographic features of the fracture surface, revealing significant fracture surface characteristics. Modern fracture image analysis systems play an important role in advancing and successfully achieving these goals, not only to speed up measurement procedures but also to perform operations that would not be possible in other ways [18].

Fractography is a method in failure analysis for studying the fracture surface of materials. The nature of decay and fracture is characterized by specific characteristics of the mechanisms of decomposition (brittleness, plasticity, etc.). Each substance has the property of decomposing. As the material is exposed to various loads and environmental influences, the nature of decomposition is usually multifaceted. The nature of fracture of a part and the loading conditions has an interface that is characterized by: the nature of the fracture, the mechanism of crack formation, and the relief of the fracture surface [19].

Due to the high load on the rear axle, especially on the tractor, its service life is shortened. Fractographic studies show that the main cause of axle failure is fatigue. Fatigue cracks have been observed near welds. The results show that the axles break due to bending fatigue caused by improper welding [20,21]. Improper welding in the heat-affected zone (HAZ) reduces the plasticity of the material, resulting in structural stress concentration points and inclusions that subsequently affect the cracks. Pre-treatment of pre-weld and post-weld heat treatment of medium carbon steel is necessary to control the hardness level of HAZ and to reduce residual stress [20].

In the literature, it is often recommended to start the engineering calculations of axes from an analytical theoretical model, which can be used to estimate displacements, specific deformations, stresses, and the self-frequency spectrum [1,22]. Theoretical calculations of axle strength are one of the most important tasks in vehicle design. The axles are exposed to different external loads with a certain frequency, which depends on the speed, the actual load, the road conditions, and many other factors. At the same time, resonant phenomena are possible, which can lead to higher than nominal stresses and many other adverse phenomena. Variable external loads cause periodic changes in stresses, which contribute to the growth of fatigue cracks leading to fatigue fractures [1,23].

During the operation of the axle, the load acting on the axle housing in the vertical direction has a significant effect on the fatigue life of the components [24]. Cracks are caused by a constant stress concentration in the axle housing, resulting in fatigue at the concentration point. If the load exceeds a certain threshold, microscopic cracks will begin to form at the stress concentration point. Later, the cracks grow due to the cyclic load of fatigue. Eventually, the cracks reach a critical limit and then the structure suddenly breaks [21,25–28]. Axle housing failures are also affected by factors such as uneven load effect, housing slope, and eccentricity.

The paper presents an analysis of damage to the rear axle of a truck semi-trailer. For this purpose, FEA preceded by a theoretical introduction was used. Several load cases were considered in order to determine the most likely point of damage initiation to the semitrailer axle. The structure of the work is as follows: the part covering materials and methods

is presented in Section 2. Then, in Section 3, the results of analytical calculations and FEM numerical simulations for axle loads with various excitation (forces or displacement) are presented. Macroscopic analysis of the damaged element and material properties tests are also presented in Section 3. As a result of the tests and analyses carried out, conclusions were developed, which are included in Section 4.

#### **2. Materials and Methods**

#### *2.1. Analytical Calculations*

The first step in the analysis of the damaged semi-trailer axle was analytical calculations, which were then used as input data for dynamic stress analysis, which was then performed using the finite element method. The calculated values of individual parameters were obtained using appropriate mathematical relationships, taking into account the assumptions of engineering knowledge and the experience contained in the available professional literature was used. The results of the analysis are presented in Section 3.1.

#### *2.2. FEM Simulations*

Simulations of the operation of the semi-trailer's non-driven axle were carried out using the Abaqus Explicit module. The model (Figure 1), apart from the axle itself, included two hinges, two bushings, two spring-damping elements and two points reflecting the contact points of the wheels with the road. The distribution of forces adopted for the simulation assumed an even position of the nominal load over the entire space of the platform.

**Figure 1.** Complete model of the axle for simulations (front in Z direction).

The prepared model contains simplifications compared to the actual structure [29,30], i.e., the possibility of deformation of the wishbones, which were modeled with Rigid Body elements, was omitted. All degrees of freedom were fixed but with rotation around the Y axis. Each of the wishbones is coupled with a bushing that is mounted on the axle. This coupling also does not take into account the possibility of deformation of these two parts (rocker arm and bush) relative to each other. The bushings are connected to the side members by means of integrated spring/dashpot elements, while in reality the shock absorber works between the swingarm and the front side member, and the air bellows between the swing arm and the side member behind the axle. The model does not take into account the stiffness and damping of the road wheels. The contact points of the wheels with the road have been connected with Coupling elements to the ends of the axles.

In the presented model, the connection of the bushing with the axle tube was modeled using welds with a shape and location similar to the real structure. In addition, contact was modeled between the inner surfaces of the bushing and the outer cylindrical surface of the axle with a friction coefficient of 1.0. It was assumed that the spring elements of the suspension have a stiffness of 4.5 kN/mm and the damping coefficient of each shock absorber is 10 N s/mm. These parameters were selected in such a way that the wheel does not lose its grip with the ground and that the deflection of the suspension after its load corresponding to the permissible axle load is in the range of 15–20 mm.

The mesh for the axles and bushing was created using linear elements of the C3D8R type, which are bricks with an integration point reduced to its center. In the area of contact between these two parts, the mesh of elements was densified (Figure 2). The bushing has 11,496 elements and 17,943 nodes, while the axle has 10,120 elements and 20,424 nodes.

**Figure 2.** Axle meshed for simulations (front in Z direction).

Simulations were carried out for eight axle load cases. For the load with the force of the wheels (wheels) in each case, the value of 44.1 kN was used, resulting from the nominal load capacity of the axle, which increases in 8 ms. A vehicle traveling at a speed of 50 km/h covers about 110 mm of road in this time.

One of the ways of loading concerned only the left wheel, i.e., as if the right wheel fell into a hole in the road. In the second case, both wheels were loaded simultaneously. Four consecutive cases consisted of vertical loading of both wheels, followed by: braking of the left wheel, braking of both wheels simultaneously, side loading of the left wheel, side loading of both wheels. The adopted values of forces (44.1 kN for each wheel) reflect the situation as if the coefficient of friction between the tires and the road surface was 1.0 both when braking and when cornering. One-wheel braking and lateral force on only one wheel can occur when the other wheel loses grip with the road surface.

Simulations were also carried out of the case when the left wheel overruns a small triangular mogul with a height of 30 mm, and then both wheels simultaneously overcome such an obstacle. Each of these cases was simulated as a kinematic excitation occurring in 8 ms.

The simulation results were analyzed in terms of the distribution of stresses and deformations of the axle material and compared with analytical calculations. Next, the results of the simulations (HMH max stress value, point/ node of max stress generation, number of stress cycles in the simulation time for max stressed nodes) were validated with the initiation point of the fracture observed on the real object.

#### *2.3. Fracture Analysis and Fractography*

The damaged semi-trailer axle was inspected. For the fractographic analysis, a fragment of the material from the damaged place was taken, then this element was cleaned of dirt and rust using the ultrasonic method. Fractographic analysis of the crack surface was performed using the MBS-2 stereoscopic microscope. Photographic documentation of the breakthrough was prepared using the Optikam Microscopy (Ponteranica, Italy) digital USB

camera. The test sample was cut using a liquid-cooled Optimum (Long Island City, NY, USA) band saw. The sample was then subjected to further analyses.

#### *2.4. Strength Tests*

The obtained material samples from the damaged semi-trailer axle were subjected to a static tensile test. The tests of the mechanical properties of the steel were carried out in the Laboratory of Strength Mechanics of the Vilnius University of Technology in accordance with the procedure contained in the ISO 6892-1 standard. Young modulus was calculated as a regression coefficient. Tests for obtaining the curve were conducted according to ISO 6892-1:2019 annex G. The Instron 8801 servohydraulic fatigue testing system (with dynamic extensometer) was used for the mechanical characterization of the sample material.

#### **3. Results and Discussion**

#### *3.1. Static Analyses of the Strength of the Axle Material*

Bending around the horizontal axis was described in [31] and the results concerning bending around horizontal axis stresses are used here.

Additional bending due to potholes around the horizontal axis is investigated.

If the truck goes to the left (or to the right), the wheel could be in way of the pothole. The most dangerous case in Figure 3 was found to be when the trajectory of the truck has less curvature. Therefore, the speed of the track is minimal. In this case there could be an additional bending moment acting around the horizontal axis, which could cause the bending moment to increase or decrease.

**Figure 3.** The scheme of the trailer turning.

Pull force is the force that a truck or prime mover can exert onto a transporter, or any type of trailer for that matter.

To go from pull force, many terms and conversion factors are thrown on the table, including the number of driven axles, gearbox ratio, rear end ratio, tire size, truck weight, and fifth wheel capacity.

The track pool force according to [32] is 11.5 Mg for truck axles and for truck trailer axles it is no more than 9.0 Mg.

$$F\_R = 9.0 \cdot 1000 \cdot 9.81 \cdot \frac{1.0}{2} = 44.1 \text{ kN},\tag{1}$$

The reaction force of the wheel (Figure 3) is derivable from the equation of equilibrium. The additional bending moment acting around the horizontal axis is shown in (Figure 4).

**Figure 4.** Reaction force.

$$M\_H = 44.1 \cdot \frac{0.995}{2} = 22.0 \text{ kNm},\tag{2}$$

Load bending for one side of the axle:

$$M\_L = 9.0 \cdot 1000 \cdot 9.81 \cdot 262.5 = 23.2 \text{ kNm},\tag{3}$$

$$
\sigma\_L = \frac{M\_L}{W} = \frac{23.2 \cdot 10^3}{125 \cdot 10^{-6}} = 185.0 \text{ MPa},\tag{4}
$$

The full stresses due to bending moment acting around the horizontal axis:

$$
\sigma\_A = \frac{M\_H + M\_L}{W} = \frac{(22.0 + 23.1) \cdot 10^3}{125 \cdot 10^{-6}} = 360.7 \cdot 10^6 \text{ Pa} = 360.7 \text{ MPa},\tag{5}
$$

#### **Bending around the vertical axis**

The bending moment acting around the vertical axle is the result of the sudden braking of the truck. The important fact is that the most significant movement occurred just after the first intensive braking [31].

The inertia force of the truck for longitudinally forward, when braking, according to standard EN 12195-1 [32], is *Fi = m*·*g*·0.8.

For a semi-trailer, it is approximately 55% load [31] and one axle inertia load is equal:

$$F\_{\text{Iside}} = \frac{m \cdot g}{2} = \frac{9000 \cdot 9.81}{2} = 44.1 \cdot 10^3 \text{ N} = 44.1 \text{ kN},\tag{6}$$

Especially, this internal force is decreasing in dynamic shock–image a wheel in a pothole during the braking of the truck. The authors accepted the dynamics coefficient according to [33] as being approximately equal to 1.2.

Bending moment acting around a vertical axis during braking:

$$M\_V = 44.1 \cdot 0.515 \cdot 1.2 = 27.3 \text{ kNm},\tag{7}$$

where: 0.515 m from the center of the wheel to the center of the rocker arm, (49.2 kNm taking into account the braking force).

Stresses due to bending moment around the vertical axis during braking:

$$
\sigma\_V = \frac{M\_V}{W} = \frac{27.3 \cdot 10^3}{125 \cdot 10^{-6}} = 218.2 \cdot 10^6 \text{ Pa} = 218.2 \text{ MPa} \tag{8}
$$

This value will be used for the calculation of compound stresses of the axle for finding the critical region. Torque from load and braking and stress

$$M\_{LV} = \sqrt{M\_L^2 + M\_V^2} = \sqrt{(23.1 \cdot 10^3)^2 + (27.3 \cdot 10^3)^2} = 35.8 \text{ kNm}$$

$$\sigma\_V = \frac{M\_{LV}}{W} = \frac{35.8 \cdot 10^3}{125 \cdot 10^{-6}} = 286.1 \cdot 10^6 \text{ Pa} = 286.1 \text{ MPa}.\tag{9}$$

**Bending stresses**

The bending moments for tubular section may be summed up superimposed.

$$\sigma\_{B} = \frac{\sqrt{M\_{H} + M\_{V}}}{W} = \frac{\sqrt{\left(22.0 \cdot 10^{3}\right)^{2} + \left(27.3 \cdot 10^{3}\right)^{2}}}{125 \cdot 10^{-6}} = 280.1 \cdot 10^{6} \text{ Pa} = 280.1 \text{ MPa} \tag{10}$$

The bending stresses are maximal in the front down quarter (Figure 3) of the axle. **Torsion of the axle**

In a scenario with potholes, the wheels of the truck are on different levels. According to [34], air cushion throw is up to 410 mm, while in [35], the nominal ride height is, on average, 300 mm.

When one wheel is in the background and another is in the pothole, this results in torsion of the axle (Figure 5). The angle of twist is:

**Figure 5.** Air cushion throw, (L ≤ 300 mm [35]).

$$
\phi = \frac{200}{500 + 385} = 0.339 \text{ rad},
\tag{11}
$$

The angle of twist, internal torque, polar moment, length, and shear modulus are related by the formula:

$$
\Delta\phi = \frac{T \cdot l}{G \cdot I\_p \prime} \tag{12}
$$

where: *Ip* = <sup>2</sup> · *<sup>I</sup>* = <sup>2</sup> · <sup>912</sup> = 1824 cm4.

Shear modulus is obtained from a relationship that exists between *G*, *E* and *υ* (Poisson's ratio).

$$G = \frac{E}{2 \cdot (1 + v)} = \frac{204}{2 \cdot (1 + 0.30)} = 78.5 \text{ GPa} \tag{13}$$

where *E*—modulus of elasticity (*E* = 204 GPa); *υ*—Poisson's ratio (*υ* = 0.30).

Internal torque:

$$T = F\_R \cdot 0.500 = 44.1 \cdot 0.5 = 22.1 \cdot 10^3 \text{Nm}.\tag{14}$$

Shearing stresses in slow motion:

$$
\tau = \frac{T}{W\_p} = \frac{T}{\frac{l\_p}{4\sqrt{2}}} = \frac{22.1 \cdot 10^3}{250058} = 88.3 \cdot 10^6 \text{ Pa} = 88.3 \text{ MPa}.\tag{15}
$$

Shearing stresses in fast motion:

$$
\pi = 1.2 \cdot 88.3 = 105.9 \,\text{MPa},
\tag{16}
$$

Stress in axle when trailer crosses potholes:

$$
\sigma\_{\mathcal{P}} = \sqrt{\sigma\_{\mathcal{L}}^2 + 3 \cdot \tau\_{\mathcal{F}}^2} = \sqrt{185.0^2 + 3 \cdot 105.9^2} = 260.6 \text{ MPa},
\tag{17}
$$

#### **Compound stresses**

The influence of bending and torsion with a truck's turning in slow motion is of an equation:

$$
\sigma\_d = \sqrt{\sigma\_B^2 + 3 \cdot \pi^2} = \sqrt{280.1^2 + 3 \cdot 105.9^2} = 334.9 \text{ MPa},
\tag{18}
$$

A compound dynamic stress of 335 MPa is calculated for driving at a speed 90 km/h on roads with potholes and more turns.

#### *3.2. Effects of Simulations*

As a result of the simulation, it was found that most often the area of maximum stress (HMH) in the axle material occurs in the place of its cooperation with the sleeve, where it is connected by a weld. This maximum lies closer to the wheel, on the outside of the axle, not between the wishbones. In addition, it usually occurs in the rear part of the axle, which, as it is further away from the axis of rotation of the wishbones, is exposed to higher torsional stresses when the deflection of the axle ends is different.

One of the ways of loading (Figure 6) concerned only the left wheel, i.e., as if the right wheel fell into a hole in the road. Four subsequent cases concerned the vertical load on both wheels (Figure 7), followed by braking the left wheel (Figure 8), braking both wheels (Figure 9), side load on the left wheel (Figure 10), and side load on both wheels (Figure 11). The adopted values of forces (*F*1side = 44.1 kN for each wheel) reflect the situation as if the coefficient of friction between the tires and the road surface was 1.0, both when braking and when cornering. One-wheel braking and lateral force on only one wheel can occur when the other wheel loses grip with the road surface.

**Figure 6.** HMH stress distribution in the axle while only left wheel loaded with *F*1side force. Trail front in Z direction.

**Figure 7.** Stress distribution in the axle while both wheels loaded with *F*1side force.

**Figure 8.** Stress distribution in the axle while only left wheel break is active.

**Figure 9.** Stress distribution in the axle while only both wheels break.

**Figure 10.** Stress distribution in the axle and max strengthen node 4384 while side force during turning acts to the left wheel.

**Figure 11.** Stress distribution in the axle while side force during turning acts to both wheels.

The analysis of the distribution of stress in Figures 6–11 demonstrates that the max values in bushing were always near the welding. For the max strengthened node, the lowest value 136 MPa was observed for the symmetrical loading of two wheels (Figure 7). A slightly larger value (151 MPa) was found for the load on one (left) axle wheel. After that, when side (Figures 8 and 9) or break (Figures 10 and 11) force was added to the simulations, for the most loaded nodes, HMH stress was found to be at least twice as big (312 MPa).

The last simulations concerned the case where the left wheel overruns a small trough (Figure 12) of a triangle shape with a height of 30 mm, and then both wheels simultaneously overcome such an obstacle (Figure 13). Each of these cases was simulated as a kinematic excitation occurring in 8 ms.

**Figure 12.** Stress distribution in the axle and max strengthened node 1275 while left wheel passes the triangular road deformation.

**Figure 13.** Stress distribution in the axle and max strengthened node 12,707 while both wheels pass the triangular road deformation.

When forcing the vertical displacement of the left wheel, the maximum stresses occur at the point of contact between the sleeve and the axle in the lower part of the axle (Figure 12). It can be seen from the graph that apart from the average stresses in the range of 100–200 MPa, there are sudden jumps in their values (peaks) by about 200 units to over 300 MPa and many load cycles.

The value of stresses in the axle (436 MPa in Figure 13) for the case of symmetrical excitation of both wheels is higher than for the displacement of only one wheel (392 MPa in Figure 12). However, their course over time is less abrupt for node 12,707 (Figure 13) than for node 1275 (Figure 12).

For the excitation carried out as displacements of both wheels, the highest stresses of 436 MPa were recorded. They occurred in the central part of the axis, with its deformation corresponding to the first form of natural vibrations.

For analytical calculations, each case can be considered separately: load, braking, side forces.

#### *3.3. Visual Fractographic Analysis*

The test object is a broken axle of a semi-trailer chassis. The axle broke in February 2016. The total mileage of the axle is 326,516 km. The axle breaks near the weld (near the mounting location). The broken axle is analyzed in the condition in which it was delivered for service under warranty. After inspecting the axle, a fracture of the axle is visible near the coupling joint in which the axle is fastened with the bracket. Non-destructive methods show only part of the fracture (Figure 14). Fractographic analysis of the fracture surface was performed using the MBS-2 stereo microscope and an Optikam Microscopy digital USB camera.

**Figure 14.** Axis fracture location.

For fractographic analysis, a specimen is excised from the fracture site. The ultrasonic bath removes dirt and rust from the fracture surface before inspecting the specimen. The specimen was cut using an Optimum liquid-cooled band saw. Liquid cooling during cutting was used to prevent the structural transformation of the steel.

Mechanical properties of steel tests were performed in the Laboratory of Strength Mechanics of Vilnius Gediminas Technical University according to standard ISO 6892-1. For further analysis using FEM (Finite Element Method), a numerical model of the considered trailer axle was prepared and simulations were carried out using the Abaqus Explicit module.

Fractures of a complex nature predominate in the fracture as the structural material is subjected to environmental influences and various deformations (Figure 15). In the manufacturing process, for example in the welded joints, areas of different mechanical properties are formed. Where there is a decrease in the plasticity of the material or in the area of higher stresses, individual voids may appear and a process of long-term plastic deformation may begin. By inspecting the fracture surface, three areas of fracture can be distinguished: fatigue before fracture, the onset of fracture, and major fracture.

**Figure 15.** Fracture surface: Areas 1 to 5 are used for fractographic analysis.

The first is the area of fatigue before fracture or decay; it has a wavy relief (Figures 16 and 17).

**Figure 16.** Fracture area: the area of the focal point of decay.

**Figure 17.** Fracture area: wider view of the area of the focal point of decay.

In the area of fatigue before the fracture, in the base and at the weld metal, near the weld where the highest stresses are applied, cracks are observed and a fragile fine-grained structure is visible. The formation and propagation of a crack depending on the type of deformation, the structure of the material, the level of the load, the shape of the part, and many other factors. The fracture is thought to have been initiated by long-term loading, as cracks and ribbed fracture morphology are characteristic of the effects of long-term loading.

Figures 16–20 show the results of fractographic studies of individual fracture areas (Figure 15). Fracture surfaces: Figure 16—an area no. 1; Figure 17—an area no. 2; Figure 18—an area no. 3; Figure 19—an area no. 4; Figure 20—an area no. 5.

**Figure 18.** Fracture area: the area of gradual crack growth.

**Figure 19.** Fracture area: The area of the final fracture.

**Figure 20.** Fracture area: area of plastic deformation.

The accumulation of plastic deformation occurs during long-term operation. It is known that during plastic deformation, the density of dislocations in the metal increases. New dislocations are due to new internal sources, the best known of which is the Franco-Reed source. The increase in the density of dislocations affects the constant increase in the

hardness of the metal. To achieve this, the process of increasing the hardness of the metal is accompanied by an increase in the brittleness of the metal. As a result, micro-cracks appear on the surface in the area of maximum hardness. During subsequent operations, these cracks turn into macro-cracks.

Fatigue fracture occurs when the surface of the fracture is perpendicular to the direction of the maximum stresses and has characteristic areas: the first—the decay foci, the second—the gradual increase in the crack, the third—the final fracture. In the area of the decay foci, a fine crystalline relief structure is visible, in the second area, the relief is crystalline in structure, and the third area is also crystalline in structure but fragile. Cracks are noticeable on the surface of the part.

The second is the area at the beginning of the fracture. Grooves and fatigue thresholds are observed, which are characteristic of low-cycle fatigue fractures (Figure 18). Fatigue thresholds and grooves are among the signs of macroscopic decay. These signs indicate the direction of crack propagation, which is associated with plastic deformations, a decrease or increase in the rate of crack propagation, depending on the effects of the environment.

The relief of the main or final fracture area is crystalline in structure, brittle, and the surface is rough and porous (Figure 19).

If all materials were absolutely plastic or absolutely brittle, plastic, or brittle fracture would occur during the tensile tests. Since there are no absolutely plastic and absolutely brittle materials, structural plastic materials are fragile when they decompose.

At the edges of the fracture, the shear characteristic of plastic deformation is seen, as well as grooves and fatigue thresholds characteristic of fatigue fracture (Figure 20).

Fractographic analysis shows that it is a fatigue fracture characterized by three fracture areas and the fracture surface oriented perpendicular to the direction of maximum stresses. It can be assumed that the fracture may have been initiated by a long-term load.

#### *3.4. Strength Properties of the Axle Material*

The mechanical properties of the steel of the axle were obtained experimentally according to standard ISO 6892-1. The middle part of the axle was used for testing (Figure 21).

**Figure 21.** Tested middle part of axel.

The blanks for specimens were cut from the middle. It was important to know the numbers of blanks and their position in the axle (Figure 22a,b). The specimen for the tension test was produced by the milling process (Figure 22c).

**Figure 22.** The blanks for specimens: (**a**) number of blanks; (**b**) orientations of blanks (**c**) specimen.

To identify the cross sections area of the specimens (Table 1), the cross sections parameters such as length *a* and thickness *b* were measured. The Young's modulus *E =* 204 GPa, yield stress *σ<sup>y</sup>* = 581 MPa and ultimate stress *σ<sup>u</sup>* = 663 MPa (Table 1) were obtained as average from the experimental curves.

**Table 1.** Mechanical characteristics of axle.


The axis in this place puts the difficult working conditions of the surface in contact with the sleeve. This can cause corrosion of the axle surface (compare Figures 14 and 15) and fretting. In addition, the end of the axle is located near the wheel, from which sand and gravel can be thrown onto the surface of the material and damage the paint coating. Corrosion pits cause additional weakening of the material and surface notches, which are dangerous, especially when the material is subjected to fatigue, i.e., in the conditions of axle operation. The safety factor related to the min yield stress value obtained from strength tests is 567/392 = 1.45. It should be noted that the result of 567 MPa was obtained for a standardly prepared tensile sample (Figure 22c). On the surface of the axis, near the fracture (Figure 14), corrosion pits are clearly visible, reducing the strength of the material, and causing geometrical and structural stress concentrations.

Poisson's ratio *ν* = 0.30, which is standard for steel was taken.

After tensions, tests of the specimens were measured and machine diagrams force F- displacement L were transformed to stress σ–strain ε diagrams. Finally, the behavior of the steel used for the axle is described by the non-linear stress σ–strain ε diagram (Figure 23), which is an average of results.

**Figure 23.** Experimental stress–strain diagrams.

In any case, this result has shown that about 60% (335/567) of elastic behavior of steel is used. Note that the analysis Equation (18) did not evaluate possible stress concentration factors due to welding impact on an axle and damage to its surface by corrosion (compare Figure 14).

In addition, the surface of the central part of the pipe is well protected with an anticorrosion coating (Figure 21) and there is no danger of its abrasion in contact with the sleeve (compare Figure 1) or the impact of gravel thrown by a wheel rolling nearby. For

these reasons, a safety factor of 567/436 = 1.30 (Figure 13) may be sufficient to achieve permanent fatigue strength.

Table 2 shows that only in the case of the simplest analytical calculations (bending stress from permissible load for one side of the axle *σL*) were the calculated stresses greater than those resulting from the simulation. In all other cases, it was the alternative: the simulation resulted in stresses higher than those calculated analytically. The reason for this may be that the simulation takes into account all types of axle loads, while the analytical calculations (except *σp*) did not take into account possible torsion (different deflection of the wishbones).


**Table 2.** Comparison of maximum stresses calculated analytically (A) and at simulations (S).

\* load (wheel displacement) of left side only.

Apart from the case of overcoming a hole in the road, the absolute difference between the stress values calculated and those resulting from the simulation was about 30 MPa, i.e., about 5% of the yield stress value.

#### **4. Conclusions**

On the basis of the non-destructive tests as well as analyses and simulations, the following conclusions were made.


**Author Contributions:** Conceptualization, Ž.D., K.D. and J.C.; methodology, Ž.D., V.R., K.D. and A.N.; software, K.D.; validation, V.R., A.N. and J.M.; formal analysis, K.D., J.T. and A.N.; investigation, Ž.D., V.R., K.D. and J.T.; resources, Ž.D., K.D. and J.C.; data curation, K.D. and J.T.; writing—original draft preparation, Ž.D., K.D., J.C., J.T. and A.N.; writing—review and editing, Ž.D., K.D., J.C., A.N. and J.M.; visualization, J.C. and J.M.; supervision, K.D., J.C. and J.M.; project administration, J.C. and J.M.; funding acquisition, J.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work was partly prepared as part of Jacek Caban and Kazimierz Drozd scientific internship at the Institute of Mechanical Science of Vilnius Gediminas Technical University which took place from 26 July to 6 August 2021.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Takayuki Shiraiwa \*, Fabien Briffod and Manabu Enoki**

Department of Materials Engineering, School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

**\*** Correspondence: shiraiwa@rme.mm.t.u-tokyo.ac.jp

**Abstract:** The 7075 aluminum alloy is a promising material for the aerospace industry due to its combination of light weight and high strength. This study proposed a method for predicting fatigue crack initiation of the 7075 aluminum alloy by crystal plasticity finite element analysis considering microstructures. In order to accurately predict the total fatigue life, it is necessary to calculate the number of cycles for fatigue crack initiation, small crack growth, and long crack growth. The long crack growth life can be estimated by the Paris law, but fatigue crack initiation and small crack growth are sensitive to the microstructures and have been difficult to predict. In this work, the microstructure of 7075 aluminum alloy was reconstructed based on experimental observations in the literature and crystal plasticity simulations were performed to calculate the elasto-plastic deformation behavior in the reconstructed polycrystalline model under cyclic deformation. The calculated local plastic strain was introduced into the crack initiation criterion (Tanaka and Mura, 1981) to predict fatigue crack initiation life. The predicted crack initiation life and crack morphology were in good agreement with the experimental results, indicating that the proposed method is effective in predicting fatigue crack initiation in aluminum alloys. From the obtained results, future issues regarding the prediction of fatigue crack initiation were discussed.

**Keywords:** aluminum alloy; fatigue; crystal plasticity; finite element method; crack initiation

**Citation:** Shiraiwa, T.; Briffod, F.; Enoki, M. Prediction of Fatigue Crack Initiation of 7075 Aluminum Alloy by Crystal Plasticity Simulation. *Materials* **2023**, *16*, 1595. https:// doi.org/10.3390/ma16041595

Academic Editors: Grzegorz Lesiuk and Dariusz Rozumek

Received: 23 November 2022 Revised: 12 January 2023 Accepted: 13 February 2023 Published: 14 February 2023

**Copyright:** © 2023 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/).

#### **1. Introduction**

The 7000 (Al-Zn-Mg) and 2000 (Al-Cu-Mg) series aluminum alloys are widely used as structural materials in aircraft and space applications. The stable precipitate phases in Al-Zn-Mg alloys are MgZn2 (η phase) and Mg3Zn3Al2 (T phase) [1]. The metastable phases (η', T' phases) that appear during the transition from the spherical GP zone to the stable phase are fine and produce substantial lattice strains, resulting in excellent precipitation strengthening. The 7075 alloy (Al-5.6Zn-2.5Mg-1.6Cu-0.23Cr in wt%) is a representative 7000 series alloys, further enhanced with the addition of Cu and Cr, and is commonly utilized in a fully precipitated state by T6 heat treatment. While the 7075 alloys have the highest tensile strength among aluminum alloys, they are inferior in fatigue strength and corrosion resistance, making them difficult to use in parts subjected to tensile fatigue loading. Their application is limited to parts subjected to compressive fatigue loading, such as the upper wing surface. Since aircraft are subjected to cyclic loading due to takeoffs, landings, and changes in atmospheric pressure, further improvement in fatigue performance is required.

The fatigue failure of metallic materials occurs when the material is subjected to repeated loads, resulting in irreversible cyclic plastic deformation, crack initiation, the small crack growth, the long crack growth, and final failure. The behavior of long crack growth can be predicted based on linear fracture mechanics. However, the prediction of fatigue crack initiation and the microstructurally small fatigue crack growth remains challenging. Since the work by Forthyth [2], many observations have been made regarding fatigue crack initiation in aluminum alloys. Under the positive stress ratio (tensile fatigue loading), fatigue crack initiations from inclusions [3] and voids [4] have been reported in 7000 series alloys, but the effect of the type and size of the inclusions on crack initiation is not clear. On the other hand, it has been observed that under conditions of stress ratios below zero, as used in practical applications, fatigue cracks are generated by the formation of transgranular facets at the surface of the material [5,6]. In order to predict such fatigue crack initiation within grains, it is necessary to analyze the cyclic plastic deformation behavior of polycrystalline materials considering the grain morphologies and crystallographic orientations.

The crystal plasticity finite element method (CPFEM) is a suitable numerical method for such purposes and has been used to predict fatigue properties in steels [7], titanium alloys [8], and aluminum alloys [9–11]. In these analyses, fatigue crack initiation life is often predicted by quantifying the driving force for fatigue crack initiation as a fatigue indicator parameter (FIP) [12], which is defined by the Tanaka-Mura model based on dislocation theory [13] and/or the Fatemi-Socie model based on critical plane theory [14]. However, methods for quantifying the FIP from the stress-strain field computed by CPFEM are still not well established. Previous papers on fatigue prediction of aluminum alloys have several drawbacks, such as the inability to predict crack shape considering the plastic anisotropy of each grain and the inability to relate the slip band length to grain morphology in the Tanaka-Mura model. Our research group has proposed a method to solve these problems by introducing potential crack paths parallel to the slip bands, and applied the method to pure iron [15], titanium alloys [16], and magnesium alloys [17].

The purpose of this paper is to investigate the applicability of fatigue crack initiation prediction by the crystal plasticity simulation to aluminum alloys. Among the 7000 series alloys, the 7075 alloy was selected as the target of this study because of the abundance of experimental fatigue data. The fatigue prediction results by CPFEM and potential crack paths are compared with the experimental data in the literature. The physical interpretation of the fatigue crack initiation criterion will be also discussed.

#### **2. Materials and Methods**

#### *2.1. Materials*

The material investigated in this study is 7075 aluminum alloy. Numerical models in this paper were constructed and calibrated with reference to experimental data of the 7075 aluminum alloy from the literature [5,18]. In the literature, a rolled 7075 aluminum alloy was subjected to T6 heat treatment, i.e., solution heat treatment at 480 ◦C for 1 h followed by aging at 120 ◦C for 24 h. The grains were elongated in the rolling direction, with average dimensions of 600 μm in the rolling direction (RD) and 108 μm in the transverse direction (TD). The tensile strength was 623 MPa, and elongation was 16%. The inverse polar figure map obtained by electron back-scatter diffraction (EBSD) experiments [18] are shown in Figure 1a.

**Figure 1.** (**a**) Inverse pole figure map of 7075 aluminum alloy in the TD-RD plane (reprinted with permission from Ref. [18]), (**b**) finite element model of polycrystalline aggregate and boundary conditions.

#### *2.2. Generation of Synthetic Microstructures*

As mentioned in the previous section, the material has an anisotropic grain morphology. To reproduce such grain morphology, a polycrystalline aggregate model was constructed by anisotropic weighted Voronoi tessellation [15]. In this tessellation method, ellipses are sampled from assumed statistical distributions of grain shapes and positioned in model space by a relaxed random sequential addition (RSA) algorithm [19]. Then, an anisotropic tessellation is performed using the center positions of the ellipses as seeds, and finally the crystal orientation is assigned. In this study, the statistical distributions of grain shape parameters were assumed based on the experimental data. The major axis length was assumed to be normally distributed with a mean of 600 μm and a standard deviation of 50 μm, the aspect ratio of the major axis to the minor axis was normally distributed with a mean of 0.2 and a standard deviation of 0.05, and the angle between the major axis and RD was normally distributed with a mean of 0◦ and a standard deviation of 5◦. The crystal orientation was assumed uniformly random for simplicity. A synthetic polycrystalline microstructure generated using these statistical distributions is shown in Figure 1b. It is noted that the polycrystalline structure is periodic in the x- and y-directions (RD and TD). Considering the computational cost, the model size and mesh size were set to 1500 μm × 3000 μm and 5 μm, respectively. The number of elements was 180,000. Three-dimensional eight-node hexahedral elements (C3D8) were used for meshing.

#### *2.3. Constitutive Law and Parameter Calibration*

To analyze the elasto-plastic deformation behavior of the polycrystalline model under cyclic loading, a constitutive law based on crystal plasticity theory was used in the finite element analysis. The total deformation gradient **F** can be decomposed into an elastic part **F**e and a plastic part **F**p:

$$\mathbf{F} = \mathbf{F\_{e}F\_{p}} \tag{1}$$

where **F**e represents the elastic deformation and the rigid body rotation, and **F**p represents the plastic deformation of the crystal lattice. The plastic velocity gradient **L** is defined as:

$$\mathbf{L}\_{\mathbf{P}} = \dot{\mathbf{F}}\_{\mathbf{P}} \mathbf{F}\_{\mathbf{P}}^{-1} = \sum\_{a=1}^{N} \dot{\gamma}^{a} \mathbf{m}^{a} \odot \mathbf{n}^{a} \tag{2}$$

where . *γα* , **m***<sup>α</sup>* and **n***<sup>α</sup>* are the plastic shear rate, slip direction vector and normal to the slip plane of the slip system *α*, respectively. A total of 12 slip systems of {111} + 110, was considered. The shear rate of each slip system is modeled by a phenomenological power law:

$$\dot{\gamma}\_s^a = \dot{\gamma}\_0 \left| \frac{\mathbf{r}^a}{\mathbf{r}\_c^a} \right|^n \text{sgn}(\mathbf{r}^a - \chi^a) \tag{3}$$

where . *<sup>γ</sup>*<sup>0</sup> is the reference shear rate ( . *γ*<sup>0</sup> = 10−<sup>3</sup> s), *τ<sup>α</sup>* is the resolved shear stress of the slip system *α*, *τ<sup>a</sup> <sup>c</sup>* is the critical resolved shear stress (CRSS), and *n* is the strain rate exponent. The work hardening behavior due to slip-slip interaction was formulated as follows:

$$
\tau\_c^\alpha = \tau\_0^\alpha + \int\_t \dot{\tau}\_s^\alpha \, \_{s \to s} \mathbf{d}t \tag{4}
$$

$$\dot{\tau}\_{\text{s}\to\text{s}}^{a} = \frac{\mathbf{d}\,\tau\_{\text{s}\to\text{s}}^{a}}{\mathbf{d}\,\Gamma\_{\text{s}}} \sum\_{\beta=1}^{N\_{\text{s}}} h^{a\beta} \left| \dot{\gamma}\_{\text{s}}^{a} \right| \tag{5}$$

$$\tau\_{\mathbf{s}\to\mathbf{s}}^{a} = \tau\_1^a \left\{ 1 - \exp\left( - \left| \frac{b\_1^a}{\tau\_1^a} \right| \Gamma\_\mathbf{s} \right) \right\} \tag{6}$$

where *τ<sup>α</sup>* <sup>0</sup> is the initial CRSS, *<sup>τ</sup><sup>α</sup>* <sup>1</sup> is the saturated CRSS, *<sup>b</sup><sup>α</sup>* <sup>1</sup> is the initial hardening rate, Γ<sup>s</sup> is the total shear strain, and *hαβ* is the interaction coefficient of the slip systems *α* and *β*. To simulate cyclic deformation behavior, the back stress evolution was defined by the Armstrong-Frederick model as follows:

$$
\dot{\chi}^a = A \dot{\gamma}^a - B \left| \dot{\gamma}^a \right| \chi^a \tag{7}
$$

where *A* and *B* are constants, each associated with a dynamic hardening and a dynamic recovery.

To calibrate the crystal plasticity parameters, strain-controlled low-cycle fatigue simulations were performed on the polycrystalline model generated in the previous section. The boundary conditions are depicted together with the model in Figure 1b. Periodic boundary conditions were applied in the x- and y-directions (RD and TD), respectively. To provide the periodic boundary conditions, equation constraints in the finite element code Abaqus were used to constrain the displacement of each node on the left end to match the displacement of the corresponding node on the right end. The top and bottom ends were similarly constrained. While maintaining these constraints, the average value of the x-displacement at the left end was fixed, and displacement (or load) was applied on the right end to simulate strain-controlled (or stress-controlled) cyclic fatigue tests. In the strain-controlled tests, the polycrystalline model was subjected to cyclic displacement in the RD with a maximum strain of 1.5%, a strain ratio of *R*<sup>ε</sup> = −1, and the number of cycles at 10. The stress-strain hysteresis loop was compared with experimental data [18] and the crystal plasticity parameters were calibrated to minimize the discrepancy. The calibrated crystal plasticity parameters are shown in Table 1.

**Table 1.** Elastic constants and crystal plasticity parameters for aluminum alloy 7075.


#### *2.4. High-Cycle Fatigue Simulation and Crack Initiation Prediction*

Stress-controlled fatigue simulations were performed to predict crack initiation in high-cycle fatigue tests. A periodic boundary condition was applied to the polycrystalline model as in the previous section, and cyclic loading was applied in the RD direction with a maximum stress of 190 to 340 MPa and a stress ratio of R = −1. The number of cycles was set to 10.

To assess the crack initiation life, an analysis based on the Tanaka-Mura model [20] was applied to the results of the crystal plasticity finite element simulations. A schematic of the crack initiation analysis is shown in Figure 2. The Tanaka-Mura model considers that a crack occurs on a slip band when the accumulated strain energy on the slip band exceeds a critical value. To specify the slip bands in the finite element model, the intersecting lines of the 2D finite element model and the slip plane were drawn at regular intervals for each slip system of each grain, as shown in Figure 2a. These lines were defined as potential crack paths. As an example, Figure 2b shows potential crack paths in a certain grain. For each potential crack path, the fatigue indicator parameter (FIP) was defined as a function of the plastic shear strain amplitude at the last cycle averaged over the path, Δ*γ*α, and the length of the path, *d*.

$$FIP = d\left(\frac{\Delta\gamma^a}{2}\right)^2\tag{8}$$

**Figure 2.** (**a**) Schematic of a potential crack path defined as the intersection of the slip plane and the 2D finite element model. (**b**) An example of the potential crack paths for a single crystalline grain. The path with the highest fatigue indicator parameter (FIP) is highlighted with a bold line.

Following the Tanaka-Mura model, the crack initiation life, *Ni*, is given by:

$$N\_i = \frac{8GW\_c}{\pi (1 - \nu)d(\Delta \tau - 2\tau\_c)^2} = \frac{A\_{\rm TM}}{FIP} \tag{9}$$

where *G* is the shear modulus, *Wc* is the fracture surface energy per unit area, *ν* is the Poisson's ratio. Using the FIP defined above, the crack initiation life can be predicted by the inverse FIP and a single material constant, *A*TM. The material constant was determined by fitting with experimental data on crack initiation life at maximum stress of 270 MPa [5]. For all potential crack paths in the polycrystalline model, *Ni* was calculated and it was assumed that the initial crack occurs in the path with the smallest value of *Ni*.

#### **3. Results and Discussion**

#### *3.1. Low-Cycle Fatigue Behavior*

Low-cycle fatigue simulations were repeated with varying crystal plasticity parameters until the difference in the stress-strain hysteresis loop from the experimental data was sufficiently small. The first ten hysteresis loops obtained from low-cycle fatigue experiments performed in the rolling direction are shown in Figure 3a. Basically, it is difficult to uniquely determine the solution for multiple crystal plasticity parameters based only on macroscopic stress-strain behavior. To avoid this non-uniqueness problem, among the crystal plasticity parameters, *<sup>n</sup>*, . *γ*0, *A*, and *B* were fixed using literature values [9], while the other parameters *τ*0 <sup>α</sup>, *τ*<sup>1</sup> <sup>α</sup>, and *b*<sup>1</sup> <sup>α</sup> were varied. The stress-strain hysteresis loops calculated with these crystal plasticity parameters are shown in Figure 3b. The obtained stress-strain behavior agreed well with the experimental data [18], including the yield stress, the maximum stress, the loop shape, the slight cyclic hardening at the beginning, and almost saturation by 10 cycles. Other experimental investigations have also reported that the hysteresis loop reaches a saturated state after about 10 cycles at most strain levels [21–23].

#### *3.2. High-Cycle Fatigue Behavior*

Stress-controlled fatigue simulations were performed using the crystal plasticity parameters calibrated in the previous section. Based on the results of the low-cycle fatigue simulation, the number of cycles was set to 10, and the stress-strain state in the last cycle was assumed to be the same as the steady state in the high-cycle fatigue test (e.g., half the number of cycles to failure). As an example, the Mises stress distribution at maximum load in the last cycle under the conditions of a maximum stress of 340 MPa and a stress ratio of *R* = −1 is shown in Figure 4a. When the applied stress was 340 MPa, the local stress exceeded 400 MPa due to the anisotropy of the elasto-plastic behavior of the crystal grains. The macroscopic stress-strain curve is shown in Figure 4b. The highest strain appeared in the tensile loading of the first cycle, and the loops after the second cycle overlapped each other. On the compression side, almost the same loops showed from the first cycle. To examine the local deformation behavior in more detail, stress-strain curves at higher and lower stress positions were displayed in Figure 4c. At the higher stress position, the maximum tensile stress in the loading direction was 421 MPa, and at the lower stress position, the stress was 269 MPa, indicating that the stress varied about ±20% from the applied stress. Additionally, the local stress ratio showed slightly different values for each position. For example, the two elements used in Figure 4c had a stress ratio of −0.95 at the higher stress position and a stress ratio of −1.07 at the lower stress position.

#### *3.3. Prediction of Fatigue Crack Initiation*

The plastic shear strain amplitudes for each slip system were extracted from the high-cycle fatigue simulations in the previous section to calculate the FIP values for all potential crack paths. The cumulative plastic shear strain is shown in Figure 5 (the slip system corresponding to each contour plot is listed in Figure 2b). As can be seen from the figure, the plastic strain showed a quite heterogeneous distribution due to the plastic activities on slip systems. The FIP values of all paths were computed and the path with the highest FIP was assumed to be the crack initiation site. The predicted initial crack is indicated by an arrow in the figure. It was oriented close to 45◦ to the loading direction and occurred in a grain with a significant accumulation of plastic strain. This inclination of the crack to the loading direction is consistent with experimental observations [5]. Fatigue prediction studies using FIP often use volume averaging, but such methods can eliminate the intragranular heterogeneity of plastic strain. Especially for coarse grains

with anisotropic geometry, the slip band size is highly dependent on the slip plane and grain shape. The concept of potential crack paths used in this paper can take into account the effect of grain shape anisotropy on crack initiation. The same crack initiation analysis was performed by changing the maximum stress in the range from 190 Mpa to 340 Mpa. The material constant *A*TM in the Tanaka-Mura model was determined to be 0.66 mm by fitting based on the experimental data of the crack initiation life at an intermediate stress level, the maximum stress of 270 Mpa [5]. This value was used to predict fatigue crack initiation life at other stress levels. The predicted crack initiation life is shown in Figure 6 together with the experimental S-N curve. The predicted crack initiation life at the maximum stress of 340 Mpa was in good agreement with the experimental value of *Ni*. Under loading conditions with no experimental crack initiation data, the predicted results of crack initiations were not contradicted with experimental failure life, *N*<sup>f</sup> (the predicted crack initiation life was smaller than the *N*f). These results suggest that the crystal plasticity simulation proposed in this study can predict the fatigue crack initiation life in 7075 aluminum alloy. It should be noted that experimental results may be affected not only by the microstructures, but also by various factors such as specimen geometry, surface roughness, and alignment of the fatigue testing machine.

**Figure 4.** (**a**) Mises stress distribution at maximum stress applied (*σ*max = 340 MPa, R = −1). (**b**) Macroscopic stress-strain curve. (**c**) Stress-strain curves at the higher and lower stress locations indicated by arrows in (**a**).

**Figure 5.** Distribution of cumulative plastic shear strain for 12 slip systems in the last calculation step (*σ*max = 340 MPa, R = −1). The predicted crack initiation site is indicated by an arrow and a red solid line.

**Figure 6.** Predicted fatigue crack initiation life and experimental fatigue crack initiation and failure life taken from the literature [5].

The proposed method has the advantage of accurately predicting crack initiation with only a single fitting parameter (*A*TM) and is useful for predicting fatigue crack initiation under various loading conditions. On the other hand, several issues remain to be addressed in order to apply the proposed method more universally to various situations. The first issue is the effect of free surfaces. In the experiments, extrusion and intrusion on an unconstrained surface can lead to fatigue crack initiation [2]. Therefore, it is necessary to include surface topology in the numerical simulations. A study that applied a crystal plasticity simulation similar to this paper to a polycrystalline model with a free surface reported that crack initiation was predicted on the surface [24]. The second issue is the crystallographic texture. It has been shown experimentally in aluminum alloys that the crystallographic texture has a significant effect on fatigue propagation and total fatigue life [25,26]. Virtual experiments based on numerical simulations also showed that the effect of the texture on fatigue crack initiation is not negligible [27]. The third issue is to model the various sources of crack initiation. It has been reported that cracks initiate from non-metallic inclusions under conditions where the stress ratio *R* = 0 or higher [3,23]. There are hard inclusions such as Al2Cu and soft inclusions such as Mg2Si, both of which can be sources of fatigue crack initiation. Another issue is the number of grains included in the polycrystalline model. Since fatigue crack initiation is a phenomenon involving scattering, fatigue crack initiation analysis should be performed on models with sufficient volume and the highest value of the FIP should be adopted. Since computational cost often limits the model size, applying an extreme value analysis [28] can be helpful. To predict total fatigue life, it is also necessary to accurately predict microstructural small fatigue crack propagation [10].

#### *3.4. Physical Meaning of Tanaka-Mura Parameter*

Let us consider the physical meaning of the fitted parameter, *A*TM. This consideration is important to make the proposed method more generally applicable to various materials. In the fatigue analysis proposed in this paper, the plastic shear strain induced by external forces and the microstructure-dependent slip length are quantified in the FIP. Hence, the *A*TM is considered to involve the nature of the material. Basically, a higher *A*TM means higher resistance to fatigue crack initiation. Pioneering work by Feltner et al. [29] more than 50 years ago pointed out that the stacking fault energy (SFE) affects dislocation structure changes and cyclic stress-strain responses, which are precursors to fatigue crack initiation. Then, many studies have reported that in materials with high SFE, cross slip is promoted, and subgrain structures of dislocations form under cyclic loading, leading to fatigue crack initiation along subgrain boundaries [30–33]. In their experiments, it has been difficult to separate the effect of SFE from the effect of yield stress on fatigue properties, making quantitative evaluation challenging. In the analysis performed in this paper, the effects of yield stress are separated to the FIP and the effect of SFE is considered to be included in the *A*TM. The *A*TM values derived for various metallic materials in previous studies are listed in Table 2 together with the SFE in the literature [34,35]. Overall, the *A*TM tended to increase with decreasing the SFE. This trend is roughly consistent with previous studies, as a lower SFE implies superior fatigue properties. The *A*TM value of 7075 aluminum alloy obtained in this study appear to deviate from the trend. One of the reasons is that the experimental values for crack initiation life were taken from literature 30 years ago, and the spatial resolution of the observations was lower than for the other three alloys obtained more recently, resulting in an overestimation of *Ni*. Another reason is that the 7075 aluminum alloy is precipitation-hardened, and the finely dispersed particles on the atomic order promote planar slips, as examined in detail by Gerold et al. [36]. As Mughrabi points out [37], explaining fatigue behavior with SFE is not inherently appropriate and should be compared to the ease of cross slip and the nature of dislocation slip (plane or wavy slip) because short-range order and short-range clustering in solid solution have a more significant effect on promoting planar slip than SFE. To this end, it is important to quantitatively evaluate the dislocation behavior of materials subjected

to cyclic deformation, which should be discussed with numerical simulations based on discrete dislocation dynamics (DDD) [38]. In DDD, the challenge is how to formulate the dislocation mobility law, which would require MD calculations [39] and inverse analysis of dislocation motion by acoustic emission (AE) [40].

**Table 2.** Tanaka-Mura parameter and stacking fault energy (SFE). The SFE is the literature value of the base metal of each material.


It should be noted that the *A*TM also depends on the distance between the slip bands assumed in the model. Recently, high-resolution digital image correlation (HR-DIC) has made it possible to observe slip bands [41]. This technique would be useful for introducing reasonable slip-band spacing into the numerical model. Additionally, it should be kept in mind that the Tanaka-Mura model (Equation (9)) is derived from strong hypotheses, such as the irreversibility of dislocation movement, uniform shear stress on the slip band, and no work hardening. The assumption regarding slip irreversibility is inconsistent with the recent observations made by atomic force microscopy [42,43]. The observation suggested that the slip irreversibility depends on the plastic strain amplitude. One possible approach to incorporate this phenomenon is to define the fraction of the slip irreversibility as a function of plastic strain amplitude and introduce it into the Tanaka-Mura model, but there is a lack of sufficient experimental knowledge to formulate such a function. Experimental methods to measure slip irreversibility more easily are required. The DIC technique and AE method may be useful for such purposes. Previous studies have suggested that in cyclic testing of pure aluminum, the generated AE signals were associated with dislocations that move during stress reversal [44]. Further quantitative discussion of dislocation behavior during cyclic deformation will be necessary to better understand the physics of the fatigue crack initiation.

#### **4. Conclusions**

This work proposed a method to predict fatigue crack initiation of 7075 aluminum alloy by calculating the cyclic deformation of a polycrystalline model based on crystal plasticity simulations and substituting the results into a crack initiation criterion of the Tanaka-Mura model. The proposed method requires the crystal plasticity parameters of the material and the critical value of the Tanaka-Mura model. The following conclusions can be drawn:


parameter on fatigue crack initiation. The importance of analyzing the dynamic dislocation behavior under cyclic loading was pointed out.

**Author Contributions:** Conceptualization, T.S., F.B. and M.E.; methodology, T.S. and F.B.; software, T.S. and F.B.; validation, T.S.; formal analysis, T.S. and F.B.; investigation, T.S.; resources, M.E.; data curation, T.S. and F.B.; writing—original draft preparation, T.S.; writing—review and editing, F.B. and M.E.; visualization, T.S. and F.B.; supervision, M.E.; project administration, T.S., F.B. and M.E.; funding acquisition, T.S. and M.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by grant-in-aid from Japan Aluminium Association (JAA), Light Metal Educational Foundation of Japan, JSPS KAKENHI for Scientific Research (B) (Grant Number 21H01648), and JST SICORP (Grant Number JPMJSC21E1).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The raw/processed data are available from the corresponding author upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Michał Böhm 1,\*, Adam Niesłony 1, Szymon Derda 1, Robert Owsi ´nski 1, Miloslav Kepka, Jr. 2, Ivana Zetkova 2, Miroslav Zetek 2, Šárka Houdková <sup>3</sup> and Mariusz Prazmowski ˙ <sup>1</sup>**


**Abstract:** At present, due to advanced fatigue calculation models, it is becoming more crucial to find a reliable source for design S–N curves, especially in the case of new 3D-printed materials. Such obtained steel components are becoming very popular and are often used for important parts of dynamically loaded structures. One of the commonly used printing steels is EN 1.2709 tool steel, which has good strength properties and high abrasion resistance, and can be hardened. The research shows, however, that its fatigue strength may differ depending on the printing method, and may be characterized by a wide scatter of the fatigue life. This paper presents selected S–N curves for EN 1.2709 steel after printing with the selective laser melting method. The characteristics are compared, and conclusions are presented regarding the resistance of this material to fatigue loading, especially in the tension–compression state. A combined general mean reference and design fatigue curve is presented, which incorporates our own experimental results as well as those from the literature for the tension–compression loading state. The design curve may be implemented in the finite element method by engineers and scientists in order to calculate the fatigue life.

**Keywords:** fatigue of materials; S–N curves; 1.2709 steel; 3D-printed materials; SLM 3D printing

#### **1. Introduction**

The proper use of S–N curve data is one of the most important steps in the process of fatigue life assessment. The reliability of the results depends on the engineer's experience, but mostly on the quality of the S–N data. Therefore, it is important to work on reliable data sets when implementing S–N curves in the fatigue life estimation process, usually together with the finite element method (FEM). The process is especially interesting when dealing with 3D-printed materials. EN 1.2709 tool steel can be found in the literature under different descriptions and names such as EOS MS1: US classification 18% Ni Maraging 300, European 1.2709 and German X3NiCoMoTi 18-9-5. It is used in the additive manufacturing process of 3D elements. EOS MS1 is delivered only as a powder. The material is praised for its good strength properties and high abrasion resistance as well as its potential to be hardened. Due to the fact that the material exists under different names, it is sometimes difficult to find the required experimental data in order to predict the fatigue life of designed components printed using it. Due to the 3D printing technology, an increasing number of mechanical elements are being designed and produced with the selective laser melting method. This method is regarded as the most promising in terms of additive manufacturing. The analysis performed within this paper consists of fatigue S–N data that were obtained for specimens produced with the use of the selective laser melting (SLM) method. The influence of different parameters was taken into account during the data analysis. Parameters such as different heat treatments and defect sizes and locations as well as the specimen smoothness

**Citation:** Böhm, M.; Niesłony, A.; Derda, S.; Owsi ´nski, R.; Kepka, M., Jr.; Zetkova, I.; Zetek, M.; Houdková, Š.; Prazmowski, M. General Reference ˙ and Design S–N Curves Obtained for 1.2709 Tool Steel. *Materials* **2023**, *16*, 1823. https://doi.org/10.3390/ ma16051823

Academic Editor: Feng Qiu

Received: 25 January 2023 Revised: 17 February 2023 Accepted: 20 February 2023 Published: 23 February 2023

**Copyright:** © 2023 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/).

and stress concentrators in the specimen shape and geometry were also described for the literature data. All data analyses were performed for the case of cyclic constant amplitude stress loading conditions. This served in the process of the comparison between the S–N curves taken from the literature and the obtained experimental curve for the tension– compression tests.

As we sorted through the literature, we noticed that papers about the fatigue life of EN 1.2709 tool steel are difficult to find, especially due to the fact that the steel can be found under different names. After noticing this, we observed that this topic is being pursued by some researchers, especially for 3D printing with the SLM technology. In this section, a number of related papers are presented, which take into account different parameters that influence the obtained experimental fatigue results; some of them also refer to the powder bed fusion (PBF) process or direct metal laser sintering (DMLS), which are equal to the SLM technology. Interesting papers that discuss the maraging steel fatigue properties for casted specimens are also mentioned. The authors mostly focus on papers that published an S–N curve. We can divide the papers into subcategories that present parameters such as:


As for different heat treatments, we found a large number of papers, such as the paper by Bouzakis et al. [1], where the authors additionally discussed the corrosion fatigue problem. Another group of authors led by Kostic et al. [2] analyzed the effect on heattreated and untreated specimens under bending loading. Elangeswaran et al. [3] also focused on the heat treatment effects as well as the surface post-treatment processes that might increase the fatigue life. It can be concluded that heat treatment increases the fatigue life of specimens.

As for the defect size and location, we found the paper by Bai et al. [4], where the authors discussed the influence of hole-type defects inside the printed material. No specific S–N curve was presented, which could have been widely discussed or used in the analysis. It may be concluded that laser power has a positive effect on the surface quality, which serves in prolonging the fatigue life of specimens.

As for powder contamination, we found the paper by Gatto et al. [5], where the authors analyzed this influence on fatigue life in comparison with forged specimens. They stated that powder contamination has a great influence on the fatigue life of 3D-printed specimens in terms of its decrease. The presented experimental results were found for the case of non-zero mean stress under tension–compression.

As for the specimen shape, we found the papers by Branco et al., where the authors discussed the effect of notches on the fatigue life [6,7]. This effect of a stress concentrator in the specimen affects the fatigue life performance of SLM-printed materials. Dörfert et al. [8] discussed the fatigue of conventionally obtained and 3D-printed materials. The S–N curves showed that the conventionally obtained material had higher fatigue strength characteristics than the 3D-printed material. Fitzka et al. [9] analyzed the size effect in terms of the fatigue life for two different sizes of sheet specimens that were tested under ultrasonic frequencies obtained from a coil material. Guo et al. [10] presented a review in which they compared the fatigue results of 3D-printed specimens for notched and unnotched cases, as well as for cases using PBF printing technology. They presented an S–N curve for a variety of experimental results under tension–compression taken from the literature. Meneghetti et al. [11] presented a paper focused mostly on the effects of the build orientation of specimens and heat treatments, along with the fatigue results for the tension– compression loading state with the use of DMLS. Solberg et al. [12] discussed the effect of the directional build orientation under tension–compression for plate specimens. The comparison was performed for different build angles, and surfaces were as-built during the

printing process or machined. The machined specimens had superior fatigue properties in comparison to the as-built ones. Tshabalala et al. [13] also performed a tension–compression test for cylindrical specimens with different build orientations. They observed that tool steels develop compressive mean stress during axial tests. It can be concluded that the shape, dimensions and distortion in terms of stress concentrators in the form of notches can cause a decrease in the fatigue life.

As for the loading type, we found the paper by Branco et al. [14], where the authors presented fatigue results under constant and variable amplitude loading. The same authors discussed the problem of multiaxial loading conditions concerning the studied material [15]. Crocollo et al. [16] presented experimental results for specimens under bending loading conditions. They focused on the sensitivity of the technology in regard to the build orientation. The paper by Damon et al. [17] presented tension–compression fatigue results under different temperatures and discussed the effect of the vertical and horizontal orientation of the printing procedure. Vilchez et al. [18] presented the fatigue results for maraging steel 300 in the tension–compression state under ultrasonic frequencies. The material was not printed but obtained from solid bars, but the obtained results showed the enormous possibilities in terms of the obtained fatigue life. It can be concluded that variable loads can cause a decrease in the fatigue life of specimens.

As for non-zero mean stress loading conditions, we found the paper by Chang [19] et al., where the authors presented multiple S–N curves for four different cycle asymmetry ratios, R. Their paper also strongly focused on the description of the defects that may influence the fatigue life of the components, as well as the temperature effects. Nevertheless, these tests were performed on a casted material, where the specimens were obtained from casted rods. Schuller et al. [20] presented fatigue S–N curves for two different R values under ultrasonic loading tests. It may be concluded that positive mean stresses influence the fatigue life of printed specimens by decreasing their lifetime.

Summarizing the literature review, it is found that the material has been tested by various scientists under different conditions. The shortcoming is the fact that it is difficult to formulate a "safe" design curve on the basis of individual literature results. It is important to note that the authors analyzed the fatigue results for 3D-printed specimens of this steel under different loading states inter alia for bending and tension–compression. The study presents never-before-published experimental fatigue results for 3D-printed specimens of 1.2709 tool steel under tension–compression on a fatigue test stand under the stress ratio R = −1. The most important part of this paper is the presentation of new experimental fatigue results as well as the presentation of a mean reference curve for the obtained experimental data and literature data as well as a reference design curve that was prepared according to the British Standard: BS 7608:1993 for the same data. There are currently no reference curves that represent a data set as large as the one presented within this paper. The design curve may be implemented in FEM by engineers and scientists in order to calculate the fatigue life in a reliable way, especially for simple objects that have been printed, e.g., plates with the dimensions reported by Van Vihn et al. [21].

#### **2. Materials and Methods**

The material used in the experimental investigation was EN 1.2709 tool steel. The mechanical properties of this steel are presented in Table 1. The powder of this material was used in order to print specimens, whose geometry and dimensions are presented in Figure 1. The specimens were not subject to additional machining, but they were in the state after solution annealing.


**Table 1.** Important mechanical properties of 1.2709 tool steel.

**Figure 1.** Geometry and dimensions of the specimens used in this experimental research for the tension–compression tests.

The machine used during the printing process was the EOS M290 with EOSTATE Exposure OT, MeltPool Monitoring systems and is presented in Figure 2. The default EOS parameter set was used for the printing of the samples in the Z direction (laser power of 258 W, scanning rate of 960 m/s, hatching of 0.11 mm, layer thickness of 40 micrometers, protective atmosphere of nitrogen gas). The samples were tested after solution annealing (temperature of 820 ◦C/1 h, slow cooling in the furnace).

**Figure 2.** Selective laser sintering 3D printing machine EOSM290 used in the process of specimen preparation.

The specimens underwent fatigue tests under tension–compression performed on the Instron 8852 test stand presented in Figure 3. The tests were performed with a zero mean stress value. Therefore, no initial compression or tension was added. In Figure 4, a photograph of the specimens before and after the fatigue tests is shown.

**Figure 3.** Fatigue test stand for tension–compression axial loading, Instron 8852, at Opole University of Technology.

**Figure 4.** Photograph of test specimens before and after axial fatigue testing, OUTech Opole.

The basic formula used for the description of the S–N curve is

$$N = B(\sigma\_a)^{-m},\tag{1}$$

where *N* is the number of cycles until failure, *B* is the S–N curve constant and *m* is the slope of the S–N curve.

An examination of the microstructure was also performed. The metallographic specimens were ground with abrasive papers or diamond pastes with decreasing gradation. Then, they were polished using an aqueous suspension of Al2O3 and etched using an etching reagent (10 mL HF, 30 mL HNO3, 50 mL H2O).

#### **3. Results**

In this section, the obtained experimental results and S–N curves are presented, as well as the S–N curves taken from the literature in order to compare other states of loading conditions. The results of the experimental fatigue tests are presented in Table 2. The microstructure of the material is presented in Figure 5. The microstructure resembles a Widmanstätten austenite structure. The Widmanstätten structure is formed in the weld after remelting and cooling. This might be the case for this material as the SLM method works on the basis of this principle. In Figure 6, the fatigue fracture surfaces of three specimens under different loading amplitude values are presented. It can be noticed that the fractures obtained using the low-cycle (MS1\_AN\_PLATF\_5\_Pr012 and MS1\_AN\_PLATF\_5\_Pr015) and high-cycle regimes (MS1\_AN\_PLATF\_5\_Pr004) in the fatigue curve have an expected fracture area of propagation and initiation on the surface of the specimens, which is the expected case for the tension–compression loading state. In Figure 7, the authors present the S–N curves obtained for the tension–compression loading state under R = −1 for smooth cylindrical specimens.


**Table 2.** Fatigue test results under tension–compression loading and R = −1.

Beside the experimental results of the fatigue tests, this section also presents an analysis of different S–N curves obtained for additively manufactured specimens reported by other authors. A Wöhler diagram is used to describe the fatigue life. The diagram presents the fatigue life in terms of cycles until failure *Nf* on the x axis. The constant stress amplitude *σ<sup>a</sup>* on the y axis represents the value at which the material is loaded during the constant tests. The literature sometimes presents the S–N diagrams in relation to the maximum stress *σmax*. In Figure 8, the experimental results of the tension–compression tests under R = −1 performed by Branco et al. [14] for specimens created with the use of the SLM technology are shown.

**Figure 5.** Microstructure of the 3D-printed specimens under different magnifications.

**Figure 6.** Fatigue fracture surfaces of three specimens under different loading amplitudes: (**a**) MS1\_AN\_PLATF\_5\_Pr004; (**b**) MS1\_AN\_PLATF\_5\_Pr012; (**c**) MS1\_AN\_PLATF\_5\_Pr015.

**Figure 7.** S–N curve for the experimental tension–compression fatigue test results of 1.2709 steel obtained by the authors under R = −1.

**Figure 8.** S–N curve for the experimental tension–compression fatigue test results (R = −1) reported by Branco et al. [14].

In Figure 9, the experimental results obtained in the work of Croccolo et al. for specimens created with the direct metal laser sintering method and tested under a rotating bending load (R = −1) are shown [16]. The authors tested the sensitivity of the material in terms of the build direction. In Figure 10, the experimental results obtained by Cruces et al. for the tension–compression (R = −1) of smooth specimens with hollow cylinders are shown [7]. The authors analyzed smooth specimens with two differently sized circular stress concentrators. It can be noticed that their results apply to notched specimens with 0.4 and 1 mm notches.

**Figure 9.** S–N curve for the experimental rotating bending fatigue test results reported by Croccolo et al. [16].

**Figure 10.** S–N curve for the experimental tension–compression fatigue test results reported by Cruces et al. [7].

In Figures 11–14, the results of the bending fatigue tests performed by Kostic et al. are presented [2]. The figures show the results for four different material conditions. The first is the condition of no heat treatment and no machining of the surface; the second is the condition of heat treatment but no machining; the third is the condition of no heat treatment and machining; and the fourth is the condition of heat treatment and machining. The experimental tension–compression S–N curves for the tests performed by Damon et al. [17] under R = −1 with two different temperature conditions, namely, room temperature (25 ◦C) and 400 ◦C, are shown in Figures 15–18. The specimens were created with the use of the SLM technology in the horizontal and vertical directions. In Figure 19, the case of an S–N curve with non-zero mean stress under a stress ratio of R = 0 is presented. In Figure 20, the case of non-zero mean stress with a value of R = 0.55 is presented. Both these cases apply to tension–compression loading.

**Figure 11.** S–N curve for the experimental bending fatigue test results reported by Kostic et al. [2] for the case of no heat treatment and no machining.

**Figure 12.** S–N curve for the experimental bending fatigue test results reported by Kostic et al. [2] for the case of heat treatment and no machining.

**Figure 13.** S–N curve for the experimental bending fatigue test results reported by Kostic et al. [2] for the case of no heat treatment and machining.

**Figure 14.** S–N curve for the experimental bending fatigue test results reported by Kostic et al. [2] for the case of heat treatment and machining.

**Figure 15.** S–N curve for the experimental tension–compression (R = −1) fatigue test results reported by Damon et al. [17] under the conditions of room temperature and a horizontal direction of the printing orientation.

**Figure 16.** S–N curve for the experimental tension–compression (R = −1) fatigue test results reported by Damon et al. [17] under the conditions of room temperature and a vertical direction of the printing orientation.

**Figure 17.** S–N curve for the experimental tension–compression (R = −1) fatigue test results reported by Damon et al. [17] under the conditions of a high temperature of 400 ◦C and a horizontal direction of the printing orientation.

**Figure 18.** S–N curve for the experimental tension–compression (R = −1) fatigue test results reported by Damon et al. [17] under the conditions of a high temperature of 400 ◦C and a vertical direction of the printing orientation.

**Figure 19.** S–N curve for the experimental tension–compression (R = 0) fatigue test results reported by Gatto et al. [5] for the case of non-zero mean stress.

**Figure 20.** S–N curve for the experimental tension–compression (R = 0.55) fatigue test results reported by Tshabalala et al. [13] for the case of non-zero mean stress.

#### **4. Discussion**

In terms of practical use and the determination of parameters that might be used by engineers and designers, we can notice that there are some specific types of S–N curves that might be applicable during their praxis. Most of them incorporate simple constant amplitude tension–compression or bending states. These types of S–N curves should be divided in terms of heat treatment, loading type, specimen smoothness and stress concentrators. From the variety of S–N curves that can be found in the literature, some specific discussion points may arise. We can notice that simple fatigue tests performed under rotating bending or pure bending allow the reader to gain more information about the behavior of printed maraging steel. However, due to the fact that these results are very rare, it is currently difficult to create a reference fatigue curve for a large amount of data. On the other hand, it is possible to form such a curve for tension–compression data. On the basis of this research, including the authors' own experimental data and the literature data, we can try to create a reference fatigue curve for the stress asymmetry ratio R = −1 with no temperature effect. This type of curve can be used in the FEM calculation procedure in order to simulate simple as well as more sophisticated loading states. Such a reference curve is presented in Figure 21. The curve was calculated on the basis of the 55 experimental points taken from our own results and the literature results. For the bending or other specific loading cases, we can notice that the results may still be insufficient to present a reference fatigue curve, which would be based on the data from only one or two papers. What can be noticed while comparing the stress amplitudes in terms of the loading state is that the bending and other states have lower curves than the tension–compression S–N curves in terms of the stress amplitude values. Table 3 presents the fatigue curve values obtained in the experimental and literature results. It also presents the reference tension–compression fatigue curve data.

**Table 3.** Fatigue test results under tension–compression loading and R = −1.



As an addition to the mean curve based on the above data, there is a calculated design S–N curve with a certainty of survival equal to 97.7%. A lognormal distribution and therefore a relationship between the standard deviation and probability are assumed. The design fatigue curve is shifted by two standard deviations below the generic mean curve. This is carried out using a similar method to that of the British Standard [22].

#### **5. Conclusions and Observations**

On the basis of the obtained results, we can formulate the following conclusions and observations:


**Author Contributions:** Conceptualization, A.N. and M.B.; methodology, M.B. and A.N.; software, M.B.; validation, M.B. and M.K.J.; formal analysis, A.N. and M.B.; investigation, S.D.; resources, I.Z. and M.Z.; data curation, M.B. and M.K.J.; writing—original draft, M.B. and A.N.; writing—review and editing, M.B., R.O., M.K.J., S.D., I.Z., M.Z. and M.P.; visualization, M.B. and R.O.; supervision, A.N.; project administration, Š.H.; funding acquisition, Š.H.; microstructure, M.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was co-financed by: (1) National Science Centre, Poland, 2020/02/Y/ST8/00093 within the M-ERA NET 2 program, project DePriSS; (2) Technology Agency of the Czech Republic (TA CR) under the EPSILON Programme within the M-Era.Net 2 co-fund call. ˇ

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
