Microstructure-Based Fatigue Modeling with Residual Stresses: Effect of Inclusion Shape on Very High Cycle Fatigue Life

: When considering the effect of inclusions on fatigue life, the size effect of inclusions is well recognized. However, many of these studies overlooked or decoupled the size effect from the shape features. Therefore, in this study, the inﬂuence of the shape characteristics of inclusions with 3 equivalent sizes of 26.6 µ m, 13.3 µ m, and 4.2 µ m on the very high cycle fatigue life of high-strength steels is investigated based on a microstructure-sensitive modeling approach, considering residual stresses. A shape parameter, unifying the aspect ratio and tilting angle of inclusion, is introduced. Based on this parameter, a new formulation of fatigue life with respect to inclusions is also proposed, extending the former one to consider the shape effect of inclusions. It is concluded that the general trend that the fatigue life increases with the decrease in inclusion size is still valid, while the shape features in terms of aspect ratio and tilting angle complicate the quantitative inﬂuence of inclusions size signiﬁcantly. Even for a constant inclusion size, the combination of shape factor and tilting angle could change the fatigue life with one order of magnitude compared with the commonly assumed round shape. These ﬁndings would enhance the precision for the fatigue life estimation based on pre-inclusion analysis and also eventually provide new dimensions for inclusion engineering to improve fatigue resistance, as size will not be the only design parameter for fatigue life.


Introduction
Fatigue property is an important mechanical property for many kinds of materials. With the improvement of material quality, more focuses are moving from low cycle fatigue to high cycle fatigue (HCF) and very high cycle fatigue (VHCF). With the increase in the circle number before fatigue failure, the main influencing factors on fatigue behavior tend to change. When the circle number is not that large, the most evident factor influencing the fatigue property is the surface quality. Rough surfaces always lead to bad fatigue behavior; meanwhile, for VHCF, the most evident factor turns to internal defects, such as inclusions. Li et al. [1,2] claimed that the inclusion-induced fatigue cracks are in the majority in their research on VHCF.
To reveal the effect of inclusions on fatigue properties, many studies were conducted via experiments and simulations [3][4][5][6]. A majority of them focused on the effect of inclusion size on fatigue life and it is proven that inclusion size indeed plays a significant role in determining the fatigue life of metals [7][8][9][10]. However, recent research also shows that the effect of other features of the inclusions also cannot be neglected, such as the type and shape of inclusions. Fu et al. [11] studied the fatigue behavior under the influence of different types of inclusions and suggested that early fatigue of bearings is governed by SiO 2 fragmentation and late fatigue by TiN inclusions. Gu et al. [12] found that the fatigue properties of bearing steels with different amounts of inclusions are dramatically different and the critical sizes of inclusions with different types that trigger fatigue cracks are distinct. Ulrike et al. [13] claimed that different 18Ni maraging steels containing different inclusions also show different fatigue properties and the effect of inclusions on fatigue behavior is related to the combined states between the inclusions and the matrix. Some advanced experimental methods were also adopted in the analysis of fatigue failure and inclusions. Naragani et al. [14] studied the crack initiation mechanism in a Ni-based superalloy specimen with high energy synchrotron X-rays and in situ mechanical testing and found that crack nucleation was easily created due to the inclusion, specifically the residual stress state and local bonding state at the inclusion-matrix interface.
Although the fatigue mechanism, concerning inclusions, has been studied rather thoroughly, it is still challenging to analyze the quantitative effect of inclusions on fatigue life due to the complex features of inclusions and the uncertainty of the inclusion distribution in the matrix materials. To solve this problem, microstructure-based models with inclusions were developed. This type of study for fatigue property was initially introduced by Dunne et al. [15] and McDowell and Dunne [16], and most of them were focusing on the microstructural features of the matrix materials, as they mainly studied the HCF properties. Concerning the VHCF simulation, the investigation team has been continuously developing microstructure-based models by including various inclusions as well as residual stresses between the inclusions and matrix [5,6]. The predicted fatigue life shows a good agreement with the experiment data in terms of the calcium aluminate inclusions and also reflects a clear effect of the inclusion size on the fatigue life.
As mentioned before, in addition to size, inclusions in steels have multiple features, including shape and mechanical/physical properties (Young's modulus, Poisson's ratio, thermal expansion coefficient, etc.). The mechanical properties of inclusions are governed by their compositions, which are connected to the generation process of inclusions [17][18][19][20]. The size and shape features belong to geometrical characteristics of inclusions. Except for the size effect, the quantitative correlation between the other features of inclusions and fatigue life is still missing. Therefore, the aim of the current study is to develop a modeling approach, correlating the shape characteristics of inclusions with fatigue life, and eventually establishing a quantitative relation between them for a comprehensive and accurate evaluation of the fatigue life of engineering materials.
The effect of shape features of inclusions on fatigue life is actually not decoupled from the size analysis. There are a few methods of these two features. The simplest one is to adopt the projected area of the inclusion perpendicular to the fatigue loading direction or calculate the equivalent diameter with the projected area. Concerning the inclusion shape, some researchers [21] use eccentricity, defined in Equation (1), to describe the features.
where a 0 is the distance of the intersection of the radius of the two centers of the ellipse where both have the same length; b 0 is the trace that is perpendicular to the intersection and the trace between the two centers. The range of ε is from 0 to 1. When a 0 equals b 0 , the ellipse has the shape of a circle. With an eccentricity of ε = 1, the ellipse would transform into a straight line. In this case, the inclusion size can be calculated with Equation (2), as follows: With these methods, some of the inclusion features can be described well when the shape of inclusions is not complicated. When the local shapes of the inclusion are distinctly different, these methods will lose details during the description. Gu et al. [22] offered a method to recreate the inclusion shape based on the coordinates of inclusion edges. The method can reproduce details of the local inclusion shape well. However, this method is designed for digitalization of existing inclusions in the microstructure model and shows obvious weakness on the analytical assistant of the total effect of a single inclusion on fatigue life due to its complexity. To solve this problem, a new method to describe the inclusion shape is offered in the present study. Using this method, the shape factors, such as aspect ratio and rotation angle, are unified in one single parameter. With the assistance of the new shape parameter, the quantitative relationship between inclusion and the corresponding fatigue life is clarified based on the simulation results of microstructure-sensitive modeling.

Modeling Strategy
The modeling strategy of the microstructure-based fatigue life prediction model with residual stresses can be divided into five steps, as follow: (a) generation of inclusions with various features in representative volume element model; (b) simulation of residual stress and fatigue process; (c) extraction of fatigue indicator parameter (FIP) data to derive the fatigue life for each configuration; (d) conduction of sensitivity study about the shape influence on fatigue life; (e) formulating a general equation to correlate the fatigue life with inclusions shapes. The flow chart of the modeling strategy is shown in Figure 1. is designed for digitalization of existing inclusions in the microstructure model and shows obvious weakness on the analytical assistant of the total effect of a single inclusion on fatigue life due to its complexity. To solve this problem, a new method to describe the inclusion shape is offered in the present study. Using this method, the shape factors, such as aspect ratio and rotation angle, are unified in one single parameter. With the assistance of the new shape parameter, the quantitative relationship between inclusion and the corresponding fatigue life is clarified based on the simulation results of microstructure-sensitive modeling.

Modeling Strategy
The modeling strategy of the microstructure-based fatigue life prediction model with residual stresses can be divided into five steps, as follow: (a) generation of inclusions with various features in representative volume element model; (b) simulation of residual stress and fatigue process; (c) extraction of fatigue indicator parameter (FIP) data to derive the fatigue life for each configuration; (d) conduction of sensitivity study about the shape influence on fatigue life; (e) formulating a general equation to correlate the fatigue life with inclusions shapes. The flow chart of the modeling strategy is shown in Figure 1.

Reconstruction of Microstructure
The reconstruction of the geometry model of the steel matrix was based on the results from electron backscatter diffraction (EBSD) tests with the assistance of representative volume element (RVE). Instead of direct digitalization of the microstructure graphs, the statistical method was applied during the reconstruction [23]. The characteristic parameters of the virtual microstructure were obtained through the grain size distribution. Further information on grain distribution and a detailed method of RVE generation can be referred to in our previous study [5]. Based on this method, 10 different RVEs were generated to avoid the unreasonable results led by the occasionality of the microstructure in the previous study.

Reconstruction of Microstructure
The reconstruction of the geometry model of the steel matrix was based on the results from electron backscatter diffraction (EBSD) tests with the assistance of representative volume element (RVE). Instead of direct digitalization of the microstructure graphs, the statistical method was applied during the reconstruction [23]. The characteristic parameters of the virtual microstructure were obtained through the grain size distribution. Further information on grain distribution and a detailed method of RVE generation can be referred to in our previous study [5]. Based on this method, 10 different RVEs were generated to avoid the unreasonable results led by the occasionality of the microstructure in the previous study.

Digitalization of Inclusions
To describe the inclusion shape accurately, the inclusions were concreted with coordinates of their edges. With the assistance of coordinates, the inclusions were all inserted in the center of the RVEs. In this study, three features of the inclusions were considered: size, aspect ratio, and the rotation angle relative to the loading direction. The inclusion sizes, which were expressed with the square roots of the inclusion area, √ area, were determined as 26.6 µm, 13.3 µm, and 4.2 µm. For each inclusion shape, the variations of two features, aspect ratio A and the tilting angle θ, are shown in Table 1, as well as their geometry sketches. All the changes of inclusion geometries in Table 1 were based on circle and ellipse. None of the inclusions with sharp edges were considered in this study. In Table 1, aspect ratios are expressed with A (A ≥ 1). When A = 1, the shape of the inclusion is a perfect circle. The tilting angle, θ, as a representation of inclusion orientation, is the angle between the major axis and the loading direction (0 • ≤ θ ≤ 180 • ). Due to the symmetry of 0 • ≤ θ ≤ 90 • and 90 • ≤ θ ≤ 180 • , only the case with 0 • ≤ θ ≤ 90 • are discussed in the present study. For a better understanding of these two parameters, a diagram of an example inclusion is shown in Figure 2. All the 12 cases of inclusions with different √ area, A, and θ are inserted in different RVEs constructed in Section 2.1. Figure 3 gives an example of one RVE with different inclusions in the center (the dark grey part represents the inclusion). All the edges of the inclusions were assumed to be tightly bonded with the steel matrix in this study. Table 1. Geometry sketches of inclusions with different shapes in this study and the corresponding values of aspect ratio A and the rotation angle θ.

Digitalization of Inclusions
To describe the inclusion shape accurately, the inclusions were concreted with coordinates of their edges. With the assistance of coordinates, the inclusions were all inserted in the center of the RVEs. In this study, three features of the inclusions were considered: size, aspect ratio, and the rotation angle relative to the loading direction. The inclusion sizes, which were expressed with the square roots of the inclusion area, √ , were determined as 26.6 μm, 13.3 μm, and 4.2 μm. For each inclusion shape, the variations of two features, aspect ratio A and the tilting angle θ, are shown in Table 1, as well as their geometry sketches. All the changes of inclusion geometries in Table 1 were based on circle and ellipse. None of the inclusions with sharp edges were considered in this study. In Table 1, aspect ratios are expressed with A (A ≥ 1). When A = 1, the shape of the inclusion is a perfect circle. The tilting angle, θ, as a representation of inclusion orientation, is the angle between the major axis and the loading direction (0° ≤ θ ≤ 180°). Due to the symmetry of 0° ≤ θ ≤ 90° and 90° ≤ θ ≤ 180°, only the case with 0° ≤ θ ≤ 90° are discussed in the present study. For a better understanding of these two parameters, a diagram of an example inclusion is shown in Figure 2. All the 12 cases of inclusions with different √ , A, and θ are inserted in different RVEs constructed in Section 2.1. Figure 3 gives an example of one RVE with different inclusions in the center (the dark grey part represents the inclusion). All the edges of the inclusions were assumed to be tightly bonded with the steel matrix in this study. Table 1. Geometry sketches of inclusions with different shapes in this study and the corresponding values of aspect ratio A and the rotation angle θ.

Geometry sketches of inclusions
Parameters A = 1 θ: unconstrained

Digitalization of Inclusions
To describe the inclusion shape accurately, the inclusions were concreted with coordinates of their edges. With the assistance of coordinates, the inclusions were all inserted in the center of the RVEs. In this study, three features of the inclusions were considered: size, aspect ratio, and the rotation angle relative to the loading direction. The inclusion sizes, which were expressed with the square roots of the inclusion area, √ , were determined as 26.6 μm, 13.3 μm, and 4.2 μm. For each inclusion shape, the variations of two features, aspect ratio A and the tilting angle θ, are shown in Table 1, as well as their geometry sketches. All the changes of inclusion geometries in Table 1 were based on circle and ellipse. None of the inclusions with sharp edges were considered in this study. In Table 1, aspect ratios are expressed with A (A ≥ 1). When A = 1, the shape of the inclusion is a perfect circle. The tilting angle, θ, as a representation of inclusion orientation, is the angle between the major axis and the loading direction (0° ≤ θ ≤ 180°). Due to the symmetry of 0° ≤ θ ≤ 90° and 90° ≤ θ ≤ 180°, only the case with 0° ≤ θ ≤ 90° are discussed in the present study. For a better understanding of these two parameters, a diagram of an example inclusion is shown in Figure 2. All the 12 cases of inclusions with different √ , A, and θ are inserted in different RVEs constructed in Section 2.1. Figure 3 gives an example of one RVE with different inclusions in the center (the dark grey part represents the inclusion). All the edges of the inclusions were assumed to be tightly bonded with the steel matrix in this study. Table 1. Geometry sketches of inclusions with different shapes in this study and the corresponding values of aspect ratio A and the rotation angle θ.

Geometry sketches of inclusions
Parameters A = 1 θ: unconstrained

Digitalization of Inclusions
To describe the inclusion shape accurately, the inclusions were concreted with coordinates of their edges. With the assistance of coordinates, the inclusions were all inserted in the center of the RVEs. In this study, three features of the inclusions were considered: size, aspect ratio, and the rotation angle relative to the loading direction. The inclusion sizes, which were expressed with the square roots of the inclusion area, √ , were determined as 26.6 μm, 13.3 μm, and 4.2 μm. For each inclusion shape, the variations of two features, aspect ratio A and the tilting angle θ, are shown in Table 1, as well as their geometry sketches. All the changes of inclusion geometries in Table 1 were based on circle and ellipse. None of the inclusions with sharp edges were considered in this study. In Table 1, aspect ratios are expressed with A (A ≥ 1). When A = 1, the shape of the inclusion is a perfect circle. The tilting angle, θ, as a representation of inclusion orientation, is the angle between the major axis and the loading direction (0° ≤ θ ≤ 180°). Due to the symmetry of 0° ≤ θ ≤ 90° and 90° ≤ θ ≤ 180°, only the case with 0° ≤ θ ≤ 90° are discussed in the present study. For a better understanding of these two parameters, a diagram of an example inclusion is shown in Figure 2. All the 12 cases of inclusions with different √ , A, and θ are inserted in different RVEs constructed in Section 2.1. Figure 3 gives an example of one RVE with different inclusions in the center (the dark grey part represents the inclusion). All the edges of the inclusions were assumed to be tightly bonded with the steel matrix in this study. Table 1. Geometry sketches of inclusions with different shapes in this study and the corresponding values of aspect ratio A and the rotation angle θ.

Geometry sketches of inclusions
Parameters A = 1 θ: unconstrained

Digitalization of Inclusions
To describe the inclusion shape accurately, the inclusions were concreted with coordinates of their edges. With the assistance of coordinates, the inclusions were all inserted in the center of the RVEs. In this study, three features of the inclusions were considered: size, aspect ratio, and the rotation angle relative to the loading direction. The inclusion sizes, which were expressed with the square roots of the inclusion area, √ , were determined as 26.6 μm, 13.3 μm, and 4.2 μm. For each inclusion shape, the variations of two features, aspect ratio A and the tilting angle θ, are shown in Table 1, as well as their geometry sketches. All the changes of inclusion geometries in Table 1 were based on circle and ellipse. None of the inclusions with sharp edges were considered in this study. In Table 1, aspect ratios are expressed with A (A ≥ 1). When A = 1, the shape of the inclusion is a perfect circle. The tilting angle, θ, as a representation of inclusion orientation, is the angle between the major axis and the loading direction (0° ≤ θ ≤ 180°). Due to the symmetry of 0° ≤ θ ≤ 90° and 90° ≤ θ ≤ 180°, only the case with 0° ≤ θ ≤ 90° are discussed in the present study. For a better understanding of these two parameters, a diagram of an example inclusion is shown in Figure 2. All the 12 cases of inclusions with different √ , A, and θ are inserted in different RVEs constructed in Section 2.1. Figure 3 gives an example of one RVE with different inclusions in the center (the dark grey part represents the inclusion). All the edges of the inclusions were assumed to be tightly bonded with the steel matrix in this study. Table 1. Geometry sketches of inclusions with different shapes in this study and the corresponding values of aspect ratio A and the rotation angle θ.

Geometry sketches of inclusions
Parameters and ellipse. None of the inclusions with sharp edges we Table 1, aspect ratios are expressed with A (A ≥ 1). When A is a perfect circle. The tilting angle, θ, as a representation angle between the major axis and the loading direction (0 metry of 0° ≤ θ ≤ 90° and 90° ≤ θ ≤ 180°, only the case with present study. For a better understanding of these two par ple inclusion is shown in Figure 2. All the 12 cases of inclu and θ are inserted in different RVEs constructed in Section of one RVE with different inclusions in the center (the dark sion). All the edges of the inclusions were assumed to be matrix in this study. Table 1. Geometry sketches of inclusions with different shapes in values of aspect ratio A and the rotation angle θ.

Geometry sketches of inclusions
Parameters

Residual Stress Simulation
During the heat treatment process of quenching, the temperature of the steel decreases sharply, which introduces residual stresses between the inclusion and steel matrix due to different thermal expansion properties and mechanical properties. The residual stress has been proven to affect the fatigue process significantly [24,25]. In this study, the residual stress distribution was simulated and set as the initial condition for fatigue loading, element by element. The simulation was conducted with Abaqus/standard. The inclusions type was chosen as calcium aluminate, which is one of the most common inclusion types existing in steels. The properties of calcium aluminate inclusions are obtained from references and shown in Table 2. For the inclusions, only the elastic properties were considered; while for the steel matrix, elastic and plastic properties were both considered. The mechanical properties of the steel matrix and detailed numerical study on the effect of inclusions with different shapes and types on residual stress distribution were investigated in the previous study [22]. The temperature variation in this study is from 835 °C to 20 °C, which is consistent with the quenching process.

Crystal Plasticity Model and Parameter Calibration
In the microstructure-based fatigue simulation, the response behavior of the steel matrix to fatigue loading is described by the CP model. The formulation of the model is based on the well-established phenomenological description of slip, which can be found

Residual Stress Simulation
During the heat treatment process of quenching, the temperature of the steel decreases sharply, which introduces residual stresses between the inclusion and steel matrix due to different thermal expansion properties and mechanical properties. The residual stress has been proven to affect the fatigue process significantly [24,25]. In this study, the residual stress distribution was simulated and set as the initial condition for fatigue loading, element by element. The simulation was conducted with Abaqus/standard. The inclusions type was chosen as calcium aluminate, which is one of the most common inclusion types existing in steels. The properties of calcium aluminate inclusions are obtained from references and shown in Table 2. For the inclusions, only the elastic properties were considered; while for the steel matrix, elastic and plastic properties were both considered. The mechanical properties of the steel matrix and detailed numerical study on the effect of inclusions with different shapes and types on residual stress distribution were investigated in the previous study [22]. The temperature variation in this study is from 835 • C to 20 • C, which is consistent with the quenching process.

Crystal Plasticity Model and Parameter Calibration
In the microstructure-based fatigue simulation, the response behavior of the steel matrix to fatigue loading is described by the CP model. The formulation of the model is based on the well-established phenomenological description of slip, which can be found in detail in the review paper by Roters et al. [27]. The model has shown great predictive capabilities in correlating the deformation anisotropy [28], damage [29], and residual stresses [30] with microstructure features. The model was extended to include the kinematic hardening in our previous study [31]. The main equations in this model are shown in Table 3.

Calculation of Shear Stress
n a -normal to the slip plane; m a -slip direction; S-second Piola-Kirchhoff stress tensor.
In Table 3, the calibration parameters are τ i , τ s , h 0 , a, . γ 0 , m, G 1 , and G 2 , and the elastic parameters are C 11 , C 12 , and C 44 . The calibration was based on the hysteresis loops from low cycle fatigue tests with a fixed strain. The calibration process was iterative with the trial-and-error method. The calibration results for the steel matrix are reported in our previous studies [5,6], which is also shown in Table 4.

Fatigue Indicator Parameter and Fatigue Simulation
A fully reversed tension-compression (R = −1) fatigue loading in the y-direction ( Figure 2) was applied in this simulation. The stress amplitude is predefined, which is 1200 MPa in this study. Six cycles of fatigue loading were acted on the model since the saturated plastic slip per cycle can be generally reached after two cycles, according to multiple studies [32,33]. The local accumulated dislocation slip was recognized as fatigue indicator parameter (FIP) according to the studies of Mughrabi et al. [34][35][36], Dunne et al. [15], Manonukul and Dunne [37], and Cheong and Busso [38]. To obtain the information of FIP accurately and clear the effect of meshing, the maximum grain-level-averaged value was calculated. The local accumulated dislocation slip p acc and the maximum grain-level-averaged value P max are calculated with Equations (7) and (8) [31].
where the parameter i is the identifier of the finite element of the involved grain, N Gr E is the number of elements within the grain, and V Gr i is the volume of the element i of the grain. The position of P max is identified as the fatigue crack initiation site.
Researchers also studied the relation between P max distribution and the cycles for fatigue crack. Gillner et al. [21,31] proposed a method to calculate fatigue life with P max and other fitting parameters (Equations (9) and (10)), as follow: where A m and B m are the distribution parameters of P max under stress amplitude indicating number m; k is the number of RVEs; N c is the number of cycles to calculate P max ; d gr is the mean grain size; α p and u are fitting parameters calibrated with at least two fatigue life data; σ m is the stress amplitude value under m.
To adopt the method offered by Gillner et al. [21,31], over three stress amplitudes are necessary for the simulation. Besides, at least two experimental fatigue life data under different stress amplitudes are needed for curve calibration. These facts narrow the application scope of this method.
Compared with Gillner et al. [21,31], Boeff et al. [39] proposed a relatively simple method, claiming that the number of cycles to initiate a microcrack is inversely proportional to P max (Equation (11)), as follows: where p c is the critical accumulated plastic slip. ∆P is the accumulated plastic slip per cycle in the stable regime. Since only one stress amplitude was adopted in this simulation, and the method of Boeff et al. [39] was chosen to analyze the relationship between P max and cycle N f .

Residual Stress Distribution with Different Inclusions
Since the residual stress distribution will be the initial state before the fatigue stress loading and the residual stress will significantly affect the fatigue behavior, the residual stress state concerning inclusions with different sizes and shapes will be analyzed in this section. Figure 4 shows typical residual stress distributions around different inclusions. The residual stress distributions are calculated with the same RVE even though the inclusion features are different. As shown in the figure, the maximum residual stress exists in the boundary between the inclusion and steel matrix. The value of residual stress decreases gradually when the position turns far away from the inclusion. It is also obvious that the residual stress is larger for inclusions with larger sizes when the shape keeps consistent. While for inclusions with the same size, the maximum residual stress around ellipse inclusions appears larger than that around circular inclusion. For elliptical inclusions, the When the maximum residual stress with different inclusions in different RVEs are calculated and analyzed, another interesting fact is noticed. Figure 5 shows the avera maximum residual stress with different RVEs and the same inclusion and the standa deviation of these data. The standard deviations for inclusions with √ = 4.3 μm distinctly larger than the larger inclusions no matter what the inclusion shapes are. T fact indicates that for smaller inclusions, especially when the inclusion size is close to grain size, the value of the maximum residual stress is affected by the microstruct around the inclusion more seriously, compared with the inclusion features themselv While for larger inclusions, the fluctuation of the maximum residual stress turns sma under different local features of the microstructure. In this case, the inclusion feature more important in residual stress analysis.  When the maximum residual stress with different inclusions in different RVEs are all calculated and analyzed, another interesting fact is noticed. Figure 5 shows the average maximum residual stress with different RVEs and the same inclusion and the standard deviation of these data. The standard deviations for inclusions with √ area = 4.3 µm are distinctly larger than the larger inclusions no matter what the inclusion shapes are. This fact indicates that for smaller inclusions, especially when the inclusion size is close to the grain size, the value of the maximum residual stress is affected by the microstructure around the inclusion more seriously, compared with the inclusion features themselves. While for larger inclusions, the fluctuation of the maximum residual stress turns smaller under different local features of the microstructure. In this case, the inclusion feature is more important in residual stress analysis. When the maximum residual stress with different inclusions in different RVEs are all calculated and analyzed, another interesting fact is noticed. Figure 5 shows the average maximum residual stress with different RVEs and the same inclusion and the standard deviation of these data. The standard deviations for inclusions with √ = 4.3 μm are distinctly larger than the larger inclusions no matter what the inclusion shapes are. This fact indicates that for smaller inclusions, especially when the inclusion size is close to the grain size, the value of the maximum residual stress is affected by the microstructure around the inclusion more seriously, compared with the inclusion features themselves. While for larger inclusions, the fluctuation of the maximum residual stress turns smaller under different local features of the microstructure. In this case, the inclusion feature is more important in residual stress analysis.

Fatigue Crack Initiation Site Prediction
To analyze the crack initiation site with inclusions of different characters, Figure 6 shows the p acc contour distributions of steel matrix in typical RVEs with different inclusions. The white arrows point at the position of P max . As shown in the figure, the positions of P max are all around the inclusions; meanwhile, the specific positions are not exactly the same as those of the maximum residual stress, which are more connected with the loading directions. This conclusion can be proven in the cases of the Ellipse-2 and Ellipse-3 inclusions. In Figure 6c,d,g,h,k,l, the positions of P max are not in the sites with the largest curvature of the ellipse. The positions of P max are all at the end of the ellipse perpendicular to the loading direction, where the stress concentration is caused by the fatigue loading peaks. In the cases of the present study, the residual stress introduced by heat treatment around inclusions can just affect the value of fatigue life, but cannot constitute a key factor of the fatigue crack initiation site.

Fatigue Crack Initiation Site Prediction
To analyze the crack initiation site with inclusions of different characters, Figure 6 shows the pacc contour distributions of steel matrix in typical RVEs with different inclusions. The white arrows point at the position of Pmax. As shown in the figure, the positions of Pmax are all around the inclusions; meanwhile, the specific positions are not exactly the same as those of the maximum residual stress, which are more connected with the loading directions. This conclusion can be proven in the cases of the Ellipse-2 and Ellipse-3 inclusions. In Figure 6c,d,g,h,k,l, the positions of Pmax are not in the sites with the largest curvature of the ellipse. The positions of Pmax are all at the end of the ellipse perpendicular to the loading direction, where the stress concentration is caused by the fatigue loading peaks. In the cases of the present study, the residual stress introduced by heat treatment around inclusions can just affect the value of fatigue life, but cannot constitute a key factor of the fatigue crack initiation site.

Results of Fatigue Life Prediction with Microstructure Modeling
To illustrate the effect of the shape features (aspect ratio and tilting angle) on the fatigue life, the material parameter of c shall be firstly calibrated. This calibration requires one data point in the S-N curve from experimental data. In this study, the experimental results for a circular inclusion with the size of √ = 13.3 μm under the stress amplitude of 1200 MPa are used. The morphology of the fatigue fracture and the composition of the inclusion inducing this fracture is shown in Figure 7. It is emphasized here that c is considered to be a material intrinsic parameter. Therefore, with one special case of circular inclusion and one stress amplitude, it is plausible to extend it for other extrinsic applications with different inclusion staples and stress amplitudes.

Results of Fatigue Life Prediction with Microstructure Modeling
To illustrate the effect of the shape features (aspect ratio and tilting angle) on the fatigue life, the material parameter of p c shall be firstly calibrated. This calibration requires one data point in the S-N curve from experimental data. In this study, the experimental results for a circular inclusion with the size of √ area = 13.3 µm under the stress amplitude of 1200 MPa are used. The morphology of the fatigue fracture and the composition of the inclusion inducing this fracture is shown in Figure 7. It is emphasized here that p c is considered to be a material intrinsic parameter. Therefore, with one special case of circular inclusion and one stress amplitude, it is plausible to extend it for other extrinsic applications with different inclusion staples and stress amplitudes.   During the calculation of the values of ∆ , the stable cycles in this study are decided from cycle 4 to cycle 6, according to previous studies [5,6]. The average fatigue life of different inclusions is shown in Figure 8a

Inclusion Parameters Unifying and the Formulation of Fatigue Life concerning Inclusion Parameters
In previous studies, the shapes of inclusions were considered as total areas or the projection areas along the loading direction when their effects on fatigue behavior were analyzed, which are quite reasonable and effective under some conditions. However, these assumptions and simplifications cannot fully describe the features of the inclusions, especially for the inclusions with the same projection area but different local shapes. In this section, a new method to describe the size and shape of inclusions is offered with three parameters-√ area, A, and θ-based on the four kinds of shapes shown in Section 2.2. In the end, an equation S shape = f ( √ area, A, θ)-is proposed, where S shape is the equivalent area when the inclusion shape is taken into consideration. Firstly, this method is designed from circular inclusions. For the circular inclusion, the change of θ is meaningless, as long as A = 1. To meet this characteristic, the parameter A and θ is represented with θ × (A − 1). Therefore, the entirety of θ × (A − 1) remains to be 0 no matter how θ changes. However, θ × (A − 1) cannot meet the demands from elliptical inclusions. For elliptical inclusions with θ = θ 0 and θ = 180 • − θ 0 , the effects on fatigue behavior are the same due to the symmetry effect. Therefore, the trigonometric function sinθ is introduced accordingly. θ × (A − 1) is improved to be sinθ × (A − 1). Under this circumstance, the values of sinθ × (A − 1) for circle inclusion and Ellipse-3 inclusion are both 0. However, the effects of these two inclusions on fatigue are obviously different. To solve this problem, the parameter A is recalled. Additionally, the value 0 of sinθ × (A − 1) will also cause trouble when it shows in the denominator in the function, hence the function to describe inclusion shape becomes A sinθ × (A − 1) . The values of these expressions are shown in Table 5 for better understanding. Table 5. Values of expressions in the function to describe inclusion shape.

Inclusion Parameters Unifying and the Formulation of Fatigue Life Concerning Inclusion Parameters
In previous studies, the shapes of inclusions were considered as total areas or the projection areas along the loading direction when their effects on fatigue behavior were analyzed, which are quite reasonable and effective under some conditions. However, these assumptions and simplifications cannot fully describe the features of the inclusions, especially for the inclusions with the same projection area but different local shapes. In this section, a new method to describe the size and shape of inclusions is offered with three parameters-√ , A, and θ-based on the four kinds of shapes shown in Section 2.2. In the end, an equation- where Sshape is the equivalent area when the inclusion shape is taken into consideration.
Firstly, this method is designed from circular inclusions. For the circular inclusion, the change of θ is meaningless, as long as A = 1. To meet this characteristic, the parameter A and θ is represented with θ × (A − 1). Therefore, the entirety of θ × (A − 1) remains to be 0 no matter how θ changes. However, θ × (A − 1) cannot meet the demands from elliptical inclusions. For elliptical inclusions with θ = θ0 and θ = 180° − θ0, the effects on fatigue behavior are the same due to the symmetry effect. Therefore, the trigonometric function sinθ is introduced accordingly. θ × (A − 1) is improved to be sinθ × (A − 1). Under this circumstance, the values of sinθ × (A − 1) for circle inclusion and Ellipse-3 inclusion are both 0. However, the effects of these two inclusions on fatigue are obviously different. To solve this problem, the parameter A is recalled. Additionally, the value 0 of sinθ × (A − 1) will also cause trouble when it shows in the denominator in the function, hence the function to describe inclusion shape becomes A sinθ × (A − 1) . The values of these expressions are shown in Table 5 for better understanding. Table 5. Values of expressions in the function to describe inclusion shape.

Geometry sketches of inclusions
According to the data regulations, when the shapes of inclusions were considered as total areas or the projection areas along the loading direction, the increase in A sinθ × (A − 1) will lead to the increase in the inclusion effect on fatigue life. Therefore, √ shape = f(√ , A, θ) is expressed with Equations (12) and (13), as follow: where is the inclusion shape parameter. It is noted that an extra parameter h is also inserted in this equation. This parameter is designed for adjusting the influence weight of

Inclusion Parameters Unifying and the Formulation of Fatigue Life Concerning Inclusion Parameters
In previous studies, the shapes of inclusions were considered as total areas or the projection areas along the loading direction when their effects on fatigue behavior were analyzed, which are quite reasonable and effective under some conditions. However, these assumptions and simplifications cannot fully describe the features of the inclusions, especially for the inclusions with the same projection area but different local shapes. In this section, a new method to describe the size and shape of inclusions is offered with three parameters-√ , A, and θ-based on the four kinds of shapes shown in Section 2.2. In the end, an equation- where Sshape is the equivalent area when the inclusion shape is taken into consideration.
Firstly, this method is designed from circular inclusions. For the circular inclusion, the change of θ is meaningless, as long as A = 1. To meet this characteristic, the parameter A and θ is represented with θ × (A − 1). Therefore, the entirety of θ × (A − 1) remains to be 0 no matter how θ changes. However, θ × (A − 1) cannot meet the demands from elliptical inclusions. For elliptical inclusions with θ = θ0 and θ = 180° − θ0, the effects on fatigue behavior are the same due to the symmetry effect. Therefore, the trigonometric function sinθ is introduced accordingly. θ × (A − 1) is improved to be sinθ × (A − 1). Under this circumstance, the values of sinθ × (A − 1) for circle inclusion and Ellipse-3 inclusion are both 0. However, the effects of these two inclusions on fatigue are obviously different. To solve this problem, the parameter A is recalled. Additionally, the value 0 of sinθ × (A − 1) will also cause trouble when it shows in the denominator in the function, hence the function to describe inclusion shape becomes A sinθ × (A − 1) . The values of these expressions are shown in Table 5 for better understanding. Table 5. Values of expressions in the function to describe inclusion shape.

Geometry sketches of inclusions
According to the data regulations, when the shapes of inclusions were considered as total areas or the projection areas along the loading direction, the increase in A sinθ × (A − 1) will lead to the increase in the inclusion effect on fatigue life. Therefore, √ shape = f(√ , A, θ) is expressed with Equations (12) and (13), as follow: where is the inclusion shape parameter. It is noted that an extra parameter h is also inserted in this equation. This parameter is designed for adjusting the influence weight of

Inclusion Parameters Unifying and the Formulation of Fatigue Life Concerning Inclusion Parameters
In previous studies, the shapes of inclusions were considered as total areas or the projection areas along the loading direction when their effects on fatigue behavior were analyzed, which are quite reasonable and effective under some conditions. However, these assumptions and simplifications cannot fully describe the features of the inclusions, especially for the inclusions with the same projection area but different local shapes. In this section, a new method to describe the size and shape of inclusions is offered with three parameters-√ , A, and θ-based on the four kinds of shapes shown in Section 2.2. In the end, an equation- where Sshape is the equivalent area when the inclusion shape is taken into consideration.
Firstly, this method is designed from circular inclusions. For the circular inclusion, the change of θ is meaningless, as long as A = 1. To meet this characteristic, the parameter A and θ is represented with θ × (A − 1). Therefore, the entirety of θ × (A − 1) remains to be 0 no matter how θ changes. However, θ × (A − 1) cannot meet the demands from elliptical inclusions. For elliptical inclusions with θ = θ0 and θ = 180° − θ0, the effects on fatigue behavior are the same due to the symmetry effect. Therefore, the trigonometric function sinθ is introduced accordingly. θ × (A − 1) is improved to be sinθ × (A − 1). Under this circumstance, the values of sinθ × (A − 1) for circle inclusion and Ellipse-3 inclusion are both 0. However, the effects of these two inclusions on fatigue are obviously different. To solve this problem, the parameter A is recalled. Additionally, the value 0 of sinθ × (A − 1) will also cause trouble when it shows in the denominator in the function, hence the function to describe inclusion shape becomes A sinθ × (A − 1) . The values of these expressions are shown in Table 5 for better understanding. Table 5. Values of expressions in the function to describe inclusion shape.

Geometry sketches of inclusions
According to the data regulations, when the shapes of inclusions were considered as total areas or the projection areas along the loading direction, the increase in A sinθ × (A − 1) will lead to the increase in the inclusion effect on fatigue life. Therefore, √ shape = f(√ , A, θ) is expressed with Equations (12) and (13), as follow: where is the inclusion shape parameter. It is noted that an extra parameter h is also inserted in this equation. This parameter is designed for adjusting the influence weight of

Inclusion Parameters Unifying and the Formulation of Fatigue Life Concerning Inclusion Parameters
In previous studies, the shapes of inclusions were considered as total areas or the projection areas along the loading direction when their effects on fatigue behavior were analyzed, which are quite reasonable and effective under some conditions. However, these assumptions and simplifications cannot fully describe the features of the inclusions, especially for the inclusions with the same projection area but different local shapes. In this section, a new method to describe the size and shape of inclusions is offered with three parameters-√ , A, and θ-based on the four kinds of shapes shown in Section 2.2. In the end, an equation- where Sshape is the equivalent area when the inclusion shape is taken into consideration.
Firstly, this method is designed from circular inclusions. For the circular inclusion, the change of θ is meaningless, as long as A = 1. To meet this characteristic, the parameter A and θ is represented with θ × (A − 1). Therefore, the entirety of θ × (A − 1) remains to be 0 no matter how θ changes. However, θ × (A − 1) cannot meet the demands from elliptical inclusions. For elliptical inclusions with θ = θ0 and θ = 180° − θ0, the effects on fatigue behavior are the same due to the symmetry effect. Therefore, the trigonometric function sinθ is introduced accordingly. θ × (A − 1) is improved to be sinθ × (A − 1). Under this circumstance, the values of sinθ × (A − 1) for circle inclusion and Ellipse-3 inclusion are both 0. However, the effects of these two inclusions on fatigue are obviously different. To solve this problem, the parameter A is recalled. Additionally, the value 0 of sinθ × (A − 1) will also cause trouble when it shows in the denominator in the function, hence the function to describe inclusion shape becomes A sinθ × (A − 1) . The values of these expressions are shown in Table 5 for better understanding. Table 5. Values of expressions in the function to describe inclusion shape.

Geometry sketches of inclusions
According to the data regulations, when the shapes of inclusions were considered as total areas or the projection areas along the loading direction, the increase in A sinθ × (A − 1) will lead to the increase in the inclusion effect on fatigue life. Therefore, √ shape = f(√ , A, θ) is expressed with Equations (12) and (13), as follow: where is the inclusion shape parameter. It is noted that an extra parameter h is also inserted in this equation. This parameter is designed for adjusting the influence weight of According to the data regulations, when the shapes of inclusions were considered as total areas or the projection areas along the loading direction, the increase in A sinθ × (A − 1) will lead to the increase in the inclusion effect on fatigue life. Therefore, is expressed with Equations (12) and (13), as follow: where ξ is the inclusion shape parameter. It is noted that an extra parameter h is also inserted in this equation. This parameter is designed for adjusting the influence weight of A and θ and balancing the synergistic effects of these two parameters, which will be fitted later.
To decide the value of parameter h, multiple values of h are applied to the calibration trials. The variations of shape indicator ξ under different shapes are shown in Figure 9. For circle inclusion, the value of ξ keeps constant. While for Ellipse-1, Ellipse-2, and Ellipse-3 inclusion, the value of ξ decreases, and the decreasing rate is slightly flattened with the increase in h. A and θ and balancing the synergistic effects of these two parameters, which will be fitted later.
To decide the value of parameter h, multiple values of h are applied to the calibration trials. The variations of shape indicator ξ under different shapes are shown in Figure 9. For circle inclusion, the value of ξ keeps constant. While for Ellipse-1, Ellipse-2, and Ellipse-3 inclusion, the value of ξ decreases, and the decreasing rate is slightly flattened with the increase in h. As presented in previous studies, the logarithm of fatigue life, Nf, can be expressed with a power function of stress amplitude and √ . In this study, the stress amplitude is not considered as a variable, thus Nf will be fitted with a power function of a single variable √ shape , where √ shape is a replacement of √ , as shown in Equation (14).
where a and b are the fitting parameters for power function. Combined with Equation (11), Equation (14) becomes the following: Based on this variation trend, Equation (15) is calibrated with data in Figure 8a and lg( c ) = 3.2 under different values of h. Several representative calibration results are shown in Figure 10. With the increase in h, the quality of calibration improves sharply at first, which is also proved by the fitting degrees R 2 ( Figure 11). When the value of h exceeds 0.8, R 2 begins decreasing. Therefore, the formulation of fatigue life concerning inclusion parameters can be expressed with Equation (16).
To further express the representative meaning of inclusion shape parameter, , Figure 12 shows the values of under different aspect ratios, A, tilting angle, θ, and h. Based on the fact that the increase in will lead to the increase in the inclusion effect on fatigue property and the decrease in fatigue life, the geometrical and physical representativeness of parameter and h can be described as follows. When A equals one, the value of and the inclusion effect keeps constant no matter how θ changes. When A remains constant and larger than one, the value of the inclusion effect keeps growing with the increase in θ. However, when θ keeps constant, the changing trend of the value of and the inclusion effect with the increase in A becomes different, which is distinguished with light As presented in previous studies, the logarithm of fatigue life, N f , can be expressed with a power function of stress amplitude and √ area. In this study, the stress amplitude is not considered as a variable, thus N f will be fitted with a power function of a single variable S shape , where S shape is a replacement of √ area, as shown in Equation (14).
where a and b are the fitting parameters for power function. Combined with Equation (11), Equation (14) becomes the following: Based on this variation trend, Equation (15) is calibrated with data in Figure 8a and lg(p c ) = 3.2 under different values of h. Several representative calibration results are shown in Figure 10. With the increase in h, the quality of calibration improves sharply at first, which is also proved by the fitting degrees R 2 ( Figure 11). When the value of h exceeds 0.8, R 2 begins decreasing. Therefore, the formulation of fatigue life concerning inclusion parameters can be expressed with Equation (16).
To further express the representative meaning of inclusion shape parameter, ξ, Figure 12 shows the values of ξ under different aspect ratios, A, tilting angle, θ, and h. Based on the fact that the increase in ξ will lead to the increase in the inclusion effect on fatigue property and the decrease in fatigue life, the geometrical and physical representativeness of parameter ξ and h can be described as follows. When A equals one, the value of ξ and the inclusion effect keeps constant no matter how θ changes. When A remains constant and larger than one, the value of ξ the inclusion effect keeps growing with the increase in θ. However, when θ keeps constant, the changing trend of the value of ξ and the inclusion effect with the increase in A becomes different, which is distinguished with light grey color in Figure 12. When sinθ is larger than h, the increase in A will lead to the increase in the value of ξ and the inclusion effect; when sinθ is smaller than h, the increase in A will lead to the decrease in the value of ξ and the inclusion effect. According to the change of h, the different values of critical θ can be calculated. In the above text, the value of h has already been determined based on the calibration quality between fatigue life and √ area ξ , which indicates that the critical value of h = arcsin(0.8) = 53 • . When θ > 53 • , the effect of inclusion on the fatigue property will increase when A increases; when θ < 53 • , the effect of inclusion on the fatigue property will decrease when A increases.
Crystals 2022, 12, x FOR PEER REVIEW 13 of 17 grey color in Figure 12. When sinθ is larger than h, the increase in A will lead to the increase in the value of and the inclusion effect; when sinθ is smaller than h, the increase in A will lead to the decrease in the value of and the inclusion effect. According to the change of h, the different values of critical θ can be calculated. In the above text, the value of h has already been determined based on the calibration quality between fatigue life and √ , which indicates that the critical value of h = arcsin(0.8) = 53°. When θ > 53°, the effect of inclusion on the fatigue property will increase when A increases; when θ < 53°, the effect of inclusion on the fatigue property will decrease when A increases.

Validation of the Formulation of Fatigue Life concerning Inclusion Parameters
According to the research of Murakami [10], the 10 7 fatigue limit for stress ratio, R = −1, depends on the square root of inclusion area to the −1/6th power, shown in Equation (17): where HV is the Vickers hardness of the steel; √ area is the square root of the inclusion area, µm; a is the position coefficient of the inclusion, a = 1.43 for surface inclusions, a= 1.41 for inclusions in contact with the sample surface, and a = 1.56 for internal inclusions. For high cycle fatigue and very high cycle fatigue properties, Furuya et al. [40] reported that the 10 9 fatigue limit also obeyed this −1/6th power rule for inclusions larger than 15 µm.
In our previous study, a formulation of fatigue life with respect to stress amplitude and inclusions size was also proposed based on the experimental data and microstructure-based fatigue modeling. The formulation is shown as Equation (18), as follows: where σ represents the stress amplitude in MPa; √ area represents the square root of the inclusion area in µm; and F 1 , F 2 , and F 3 are the fitting constants, where F 1 = 15,402 MPa 0.98 ·µm 1/4.7 , F 2 = −0.98, and F 3 = −1/4.7.
The inclusions in the previous study were calcium aluminate circle inclusion. To compare the previous formulation and Equation (16). The value of ξ is calculated first with the data of circle inclusion (A = 1, (sinθ − h) × (A − 1) = 0). Therefore, for circle inclusion, Equation (16) becomes Equation (19).
It is noted that the exponent changed to −1/26. With this exponent, the slope of the curve turns flatter, which means the effect of inclusion size change on fatigue life decreases in this formulation. Since the inclusions considered in Equations (18) and (19) are actual inclusions observed in the fatigue failure. The shapes of these inclusions are not that regular. The local irregularity directly leads to a stronger impact on fatigue life; meanwhile, in the present study, artificial inclusions with perfect regular shapes are inserted in the model for numerical study. The absence of the roughness in the inclusion shapes is the major cause of the smaller value of the exponent. In the following work, the surface roughness of inclusions will also be included in the simulation to extend the applicability of the formulation.

Conclusions
(1) In this study, the effect of the shape features of inclusions on fatigue life is systematically investigated via a microstructure-based modeling approach. Based on the findings, an analytical formulation to correlate the fatigue life with the size and shape features of inclusions is proposed. This formulation extends the former fatigue equation to shape parameters (aspect ratio, A, and tilting angle, θ). (2) To describe the inclusion shape, a new parameter, ξ, which unifies the aspect ratio and tilting angle of the inclusion, is introduced. This parameter, ξ, shows its representativeness in both the geometrical and physical sides. This new parameter and inclusion area jointly determine the equivalent inclusion area, which was adopted in the formulation of fatigue life as a single variable. (3) A critical θ is also offered. When θ > 53 • , the effect of inclusion on the fatigue property will increase when A increases; when θ < 53 • , the effect of inclusion on the fatigue property will decrease when A increases. (4) For smaller inclusions, the value of the maximum residual stress around the inclusion is affected by the microstructure of the steel matrix around the inclusion more seriously compared with the inclusion features itself; meanwhile, for larger inclusions, the inclusion shape is more important in residual stress analysis. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.