Proposal of a Probabilistic Model on Rotating Bending Fatigue Property of a Bearing Steel in a Very High Cycle Regime

: In S-N diagrams for high strength steels, the duplex S-N curves for surface-initiated failure and interior inclusion-initiated failure were usually conﬁrmed in the very high cycle regime. This trend is more distinct in the loading type of rotating bending, due to the stress distribution across the section. In the case of interior failure mode, the ﬁsh-eye is usually observed on the fracture surface and an inclusion is also observed at the center of the ﬁsh-eye. In the present work, the authors attempted to construct a probabilistic model on the statistical fatigue property in the interior failure mode, based on the distribution characteristics of the location and the size of the interior inclusion at the crack initiation site. Thus, the P-S-N characteristics of the bearing steel (SUJ2) in the very high cycle regime were successfully explained.


Introduction
From both viewpoints of economical design of mechanical structures and reduction of the environmental load, elongation of the design life for every product in the wide range of industries, is an important subject. In such a circumstance, for train axles and wheels, load bearing parts of automobiles, turbine rotors of the power plants and so on, the number of loading cycles reaches over 10 9 cycles, during the service period of the respective mechanical components [1][2][3][4]. One of the typical aspects of fatigue property in the very high cycle regime for structural steels is the fact that the duplex S-N characteristics consisting of the respective S-N curves for the surface-initiated and the interior-initiated failures are observed, as reported by many researchers [3,[5][6][7]. The term "S-N curve" refers to the relationship between the stress level [S] and the number of stress cycles to failure [N]. Expressions of S-N property and S-N characteristics are also used in this paper, based on a similar sense. In the case of interior-initiated failure, an inclusion is usually found at the core portion of the fish-eye formed on the fracture surface [3,5,8]. Thus, the interior inclusion at the crack initiation site plays a dominant role to govern the fatigue strength and fatigue life of the specimen.
It is well-known that the fatigue limit of the metallic material can be evaluated by combining the material hardness and the concept of √ area for the defect size proposed by Murakami et al. [9]. Another important factor to control the fatigue property is the location of the inclusion, since the stress distribution across the section has a steep slope in rotating bending [10]. Accordingly, the fatigue property of those metallic materials should be analyzed as a function of the size and the location of the interior inclusions, together with the hardness. However, both the size and the location of the inclusions inside the material have particular distributions, depending on the fabrication process of the actual metallic materials [11][12][13].
From this point of view, the authors attempted to construct a probabilistic model on the very high cycle fatigue property for the bearing steel of SUJ2 in rotating bending, by combining the distribution characteristics of the size and the location of the interior inclusion. In this model, a two-parameter Weibull distribution was accepted for the distribution pattern of the inclusion size, whereas the random distribution was accepted for the location of the inclusion. Since the scatter of the hardness is sufficiently small, the hardness of this steel was supposed to be constant for the sake of simplicity. Although some number of papers analyzing the statistical aspects of S-N property were published [14][15][16][17][18][19], the main target of those research is the statistical property within the usual life region, which is less than 10 7 cycles. Considering this point, the main aim of the present study was set to the physical interpretation of the statistical aspects of the S-N characteristics of a bearing steel, in the very high cycle regime.

Experimental Fatigue Test Data Referred in This Study
Sufficient number of fatigue test data are required to analyze the statistical fatigue property of metallic materials. However, the fatigue test in the very high cycle regime occupies a long time. For example, at the usual loading frequency of 50-60 Hz, it takes around 200 days to reach the loading cycles of N = 10 9 . In such a circumstance, the authors extracted a series of fatigue test results obtained as round robin tests by the Research Group on Statistical Aspects of Materials Strength [RGSAMS] [3]. These fatigue tests were performed in rotating bending in very high cycle regime over N = 10 9 cycles, for a bearing steel of JIS:SUJ2, by using the same type of testing machine and common specimens. Since such fatigue test results give useful data for discussion on the very high cycle fatigue property of metallic materials, they are often referred to by many researchers in this area [3,7,20,21].
All experimental results are plotted as an S-N diagram in Figure 1, in which fatigue test data obtained at different laboratories in seven universities in Japan were indicated altogether. Open symbols represent the S-N data for the surface-induced fracture, whereas solid symbols indicate the results for the interior-induced fracture. Several data points attached arrows around N = 10 9 , which indicate runout specimens without fracture. The solid and dotted curves indicate the S-N curves fitted to the respective fracture modes, using the JSMS Standard of "Standard Regression Method of S-N Curves", accepting the semi-logarithmic bilinear model [22]. The fatigue limit for the surface-induced fracture was σ w = 1278 MPa, as given by the horizontal portion of the solid S-N curve. This conventional fatigue limit cannot provide the true fatigue limit, since fatigue fractures take place in the very high cycle region, at stress levels lower than the horizontal portion of the solid S-N curve. For the experimental data in the interior-induced fracture mode, the linear regression line was indicated by a dotted line. Although the horizontal line was also indicated from the point of N = 10 9 along the dotted S-N curve, further experimental reconfirmation would be required in order to make clear the appearance of the horizontal portion in the fracture mode by continuing fatigue tests, until further long life region. Typical examples of fracture surfaces in the respective fracture modes are shown in Figure 2, where (a) gives the SEM observation of the fracture surface in the surface-induced fracture mode and (b) is the fracture surface in the interior-induced fracture mode [23]. The fatigue crack occurs at the specimen surface in the short life region at a relatively Typical examples of fracture surfaces in the respective fracture modes are shown in Figure 2, where (a) gives the SEM observation of the fracture surface in the surface-induced fracture mode and (b) is the fracture surface in the interior-induced fracture mode [23]. The fatigue crack occurs at the specimen surface in the short life region at a relatively high stress level, as shown in Figure 2a. However, a clear fish-eye is found on the fracture surface in the very high cycle region, at low stress levels, as shown in Figure 2b. In such a case, an inclusion is also found at the core site of the fish-eye. This fact means that the fatigue crack takes place around the inclusion at the crack initiation site, in a very high cycle regime. The main aim of this study was to analyze the statistical distribution characteristics of the S-N property in the interior inclusion-induced fracture mode, by combining distribution characteristics of the inclusion size and the inclusion depth. Typical examples of fracture surfaces in the respective fracture modes are shown in Figure 2, where (a) gives the SEM observation of the fracture surface in the surface-induced fracture mode and (b) is the fracture surface in the interior-induced fracture mode [23]. The fatigue crack occurs at the specimen surface in the short life region at a relatively high stress level, as shown in Figure 2a. However, a clear fish-eye is found on the fracture surface in the very high cycle region, at low stress levels, as shown in Figure 2b. In such a case, an inclusion is also found at the core site of the fish-eye. This fact means that the fatigue crack takes place around the inclusion at the crack initiation site, in a very high cycle regime. The main aim of this study was to analyze the statistical distribution characteristics of the S-N property in the interior inclusion-induced fracture mode, by combining distribution characteristics of the inclusion size and the inclusion depth. For the sake of important reference, static mechanical properties such as the tensile strength and the Vickers' hardness were examined, prior to performing fatigue tests of this bearing steel. Experimental results of these properties are indicated in Table 1. For the sake of important reference, static mechanical properties such as the tensile strength and the Vickers' hardness were examined, prior to performing fatigue tests of this bearing steel. Experimental results of these properties are indicated in Table 1.

Effects of Size and Depth of Inclusions on the Fatigue Strength
According to the √ area model proposed by Y. Murakami [9], the fatigue limit of a metallic material, σ w , was well provided by the defect size ρ and the Vickers' hardness HV (kgf/mm 2 ), as follows; where area indicates the area of the defect projected perpendicular to the longitudinal direction of the specimen, and the factor of α was given by α = 1.43 for surface defect and α = 1.56 for interior defect, respectively. When this concept of √ area was applied to the fatigue strength at N = 10 9 cycles as an attempt, assuming a spherical interior inclusion with a radius of ρ, the fatigue strength at N = 10 9 , σ w9 , was given by the following expression, σ w9 = 1.56(HV + 120)/ πρ 2 Since the Vickers' hardness of this bearing steel (JIS:SUJ2) was reported as HV = 778 [3], Equation (2) was reduced as follows; Equation (3) indicates the relationship between the fatigue strength at N = 10 9 and the inclusion size (radius) for this steel.
Such a relationship is depicted in Figure 3. When the inclusion size is small, those inclusions have no harmful effect on the fatigue strength, as reported by many researchers [24,25]. From this point of view, if the critical small size is assumed as ρ = 3 µm in this study, the upper bound of the fatigue strength given by Equation (3) becomes σ w9 = 1060 MPa, as indicated in Figure 3. On the other hand, an inclusion larger than 40 µm (ρ > 40 µm) is seldom found in the present steel. In such an assumption, the lower bound of the fatigue strength given by Equation (3) becomes σ w9 = 689 MPa, as indicated by the horizontal dotted line in Figure 3. Thus, the distribution range of the fatigue strength at N = 10 9 due to the distribution of the inclusion size would be 689 MPa~1060 MPa.

Distribution Characteristics of Inclusion Size and Inclusion Depth
The inclusion size is not constant, and therefore, the size has its own dist characteristics peculiar to the fabrication process of the metallic material [11course, since the inclusion size is always positive, we have > 0. In other wo lower bound of the inclusion size is assumed "0". In addition, a size larger than 4 seldom found in the bearing steel. Distribution pattern of the inclusion size with th feature is well-represented by a two-parameter Weibull distribution [27], withou tion parameter. The probability density function and the cumulative distribution f are given as follows: where a is the shape parameter and b is the scale parameter, respectively. When the number of inclusions included within the critical volume of a spec denoted by n, the inclusion with the maximum size would be the crack starter un cyclic loadings. Accordingly, the distribution characteristics of the maximum size n of inclusions in the critical volume become the most important factor to contro tigue strength and fatigue life of the metallic material. Based on the concept of ex distribution [28], the probability density function of the maximum size among In the case of the rotating bending, the stress distribution across the specimen section has a distinct gradient, such that the stress on the surface gives the maximum and the stress becomes zero at the center of the section. Accordingly, the inclusion depth is another important factor to govern the fatigue strength of σ w9 , even if the inclusion size is fixed to a certain value. When the fatigue strength at the depth ξ is denoted as σ w9 , the fatigue strength evaluated by nominal bending stress on the specimen surface, σ * w9 , is given as follows [26]; where r indicates the radius of the critical portion of the specimen. Here, substituting Equation (3) into Equation (4), we have the following expression; Consequently, the fatigue strength of σ * w9 for the specimen with an inclusion size of ρ and a depth of ξ could be calculated by Equation (5).

Distribution Characteristics of Inclusion Size and Inclusion Depth
The inclusion size is not constant, and therefore, the size has its own distribution characteristics peculiar to the fabrication process of the metallic material [11][12][13]. Of course, since the inclusion size ρ is always positive, we have ρ > 0. In other words, the lower bound of the inclusion size is assumed "0". In addition, a size larger than 40 µm is seldom found in the bearing steel. Distribution pattern of the inclusion size with the above feature is well-represented by a two-parameter Weibull distribution [27], without a location parameter. The probability density function and the cumulative distribution function are given as follows: and where a is the shape parameter and b is the scale parameter, respectively. When the number of inclusions included within the critical volume of a specimen is denoted by n, the inclusion with the maximum size would be the crack starter under the cyclic loadings. Accordingly, the distribution characteristics of the maximum size among n of inclusions in the critical volume become the most important factor to control the fatigue strength and fatigue life of the metallic material. Based on the concept of extremes' distribution [28], the probability density function of the maximum size among n inclusions, f (ρ n ), is provided by Although we have to know the original distribution of the inclusion size ρ itself, it is difficult to directly observe such a distribution pattern. In such a circumstance, the authors determined the respective parameters of the original size distribution, so as to satisfy the following conditions: This condition for the inclusion size was introduced here as a possible assumption supposed from the observation results by Toriyama et al. [12]. Thus, we obtain the results for the Weibull parameters as a = 1.139 and b = 7.214, respectively. In addition, since the number of inclusions inside the critical volume, n, is also unknown, the authors analyzed the distribution of the extremes by supposing various values as the number of inclusions, like n = 2, 3, 5, 6, 10, 20, and 40, respectively. Probability density functions of the maximum inclusion size thus calculated by Equation (8) are depicted in Figure 4, where f (ρ 1 ) plotted by a fine dotted curve in the most left hand side indicates the original probability density function of the inclusion size. The mode of the largest inclusion size, i.e., the inclusion size at each peak of the probability density function, tends to shift into the right hand side, with an increase in the number of inclusions in the critical volume of the specimen.
As described in Section 3.1, the inclusion depth at the crack initiation site, ξ, plays an important role to control both the fatigue strength and fatigue life of the metallic material, in the case of the rotating bending, due to the steep slope of the stress distribution across the section. Figure 5 indicates the distribution feature of the inclusions projected to the specimen section, assuming random distribution (uniform distribution) [23]. Here, if the inclusion depth from the surface is denoted as ξ 1 , ξ 2 , . . . , ξ n , from the most shallow inclusion, one can obtain the one-dimensional distribution of the inclusion depth. As reported in another paper [26], the probability density function of the inclusion depth ξ can be represented by the following expression: ( 1 ) plotted by a fine dotted curve in the most left hand side indicates the original probability density function of the inclusion size. The mode of the largest inclusion size, i.e., the inclusion size at each peak of the probability density function, tends to shift into the right hand side, with an increase in the number of inclusions in the critical volume of the specimen. As described in Section 3.1, the inclusion depth at the crack initiation site, ξ, plays an important role to control both the fatigue strength and fatigue life of the metallic material, in the case of the rotating bending, due to the steep slope of the stress distribution across the section. Figure 5 indicates the distribution feature of the inclusions projected to the specimen section, assuming random distribution (uniform distribution) [23]. Here, if the inclusion depth from the surface is denoted as 1 , 2 , · ·, , from the most shallow inclusion, one can obtain the one-dimensional distribution of the inclusion depth. As reported in another paper [26], the probability density function of the inclusion depth ξ can be represented by the following expression: In Equation (9), r indicates the radius of the specimen. The probability density of Equation (9) becomes a maximum at the specimen surface (ξ = 0), whereas the density becomes zero at the center of the specimen section. As reported by several researchers [5,29], the depth of the crack initiation site is restricted within the thin surface layer of ξ < 250 μm in rotating bending. Therefore, the probability density function of Equation (9) should be normalized as the conditional distribution, under the condition of ξ < 250 μm. Thus, we have the actual probability density function of ξ, as follows: where indicates the probability giving the condition of 0 < < 250 μm and its probability becomes = 0.3055 in the condition of the present work. The probability density function thus normalized by Equation (10) is depicted in Figure 6.

Joint Distribution of Inclusion Size and Inclusion Depth and Analysis of Fatigue Strength Distribution
In the previous section, distribution characteristics on the size and the depth of the inclusion are discussed, and the probability density functions of ( ) and 0 ( ) are derived as Equations (8) and (10) In Equation (9), r indicates the radius of the specimen. The probability density of Equation (9) becomes a maximum at the specimen surface (ξ = 0), whereas the density becomes zero at the center of the specimen section. As reported by several researchers [5,29], the depth of the crack initiation site is restricted within the thin surface layer of ξ < 250 µm in rotating bending. Therefore, the probability density function of Equation (9) should be normalized as the conditional distribution, under the condition of ξ < 250 µm. Thus, we have the actual probability density function of ξ, as follows: where F c indicates the probability giving the condition of 0 < ξ < 250 µm and its probability becomes F c = 0.3055 in the condition of the present work. The probability density function thus normalized by Equation (10) is depicted in Figure 6.

Joint Distribution of Inclusion Size and Inclusion Depth and Analysis of Fatig Distribution
In the previous section, distribution characteristics on the size and th inclusion are discussed, and the probability density functions of ( ) and rived as Equations (8)

Joint Distribution of Inclusion Size and Inclusion Depth and Analysis of Fatigue Strength Distribution
In the previous section, distribution characteristics on the size and the depth of the inclusion are discussed, and the probability density functions of f (ρ n ) and f 0 (ξ) are derived as Equations (8) and (10). Since both ρ n and ξ are statistically independent random variables, the joint probability density function, h(ρ n , ξ), is given by By putting n = 5, a = 1.139, b = 7.214, r = 1.5, and F c = 0.3055 along the previous example, the joint probability density function of h(ρ n , ξ) was numerically calculated. The result of h(ρ n , ξ) thus obtained is depicted in Figure 7, where the joint probability function of h(ρ n , ξ) gives the curved surface like a mountain range. In the figure, ∆H(ρ n , ξ) corresponding to the volume of the vertical column indicates the probability that ρ n and ξ are within the square region of ∆ρ n × ∆ξ. The dashed line appearing on the ξ − h(ρ n , ξ) plane corresponds to the marginal distribution for the inclusion depth ξ, as given by Equation (10). On the other hand, another dashed curve appearing on the ρ n − h(ρ n , ξ) plane corresponds to the marginal distribution for the inclusion size of ρ n , given by Equation (12).
In this place, let us divide the ρ n − ξ plane into fine meshes of ∆ρ n × ∆ξ and consider the mesh at the arbitrary point of (ρ n , ξ), as shown in Figure 7. Then, the probability that the inclusion size ρ n and the inclusion depth ξ yield within the area of ∆ρ n × ∆ξ is provided by the volume of the vertical column raising at this mesh of ∆ρ n × ∆ξ. In other words, the volume of this column corresponds to the probability that the fatigue strength at N = 10 9 is given by Equation (5). Thus, we have the following equation, Based on the repetition of numerical calculations by Equation (13), the probability density function f (σ * w9 ) and the cumulative distribution function F(σ * w9 ) can be analyzed numerically.

Distribution Characteristics of Fatigue Strength 9 *
Since the number of inclusions included within the critical volume of a specimen is unknown, the authors analyzed the distribution characteristics of the fatigue strength, 9 * , assuming several numbers of inclusions like n = 2, 3, 5, 6, 10, 20, and 40, respectively. Among them, the analytical results of the fatigue strength distribution obtained for n = 2, 5, and 20, by the method in the previous section are indicated in Figure 8 as typical examples. Figure 8a shows the probability density functions under these three cases of n = 2, 5, and 20, whereas Figure 8b indicates the cumulative distribution functions corresponding to the respective cases. Weibull parameters a and b in Equation (8) are provided to satisfy the condition of [1 μm < < 15 μm] = 0.8. It is found that the peak (mode) and the standard deviation of the fatigue strength distribution tends to decrease with an increase in the inclusion number.

Distribution Characteristics of Fatigue Strength σ * w9
Since the number of inclusions included within the critical volume of a specimen is unknown, the authors analyzed the distribution characteristics of the fatigue strength, σ * w9 , assuming several numbers of inclusions like n = 2, 3, 5, 6, 10, 20, and 40, respectively. Among them, the analytical results of the fatigue strength distribution obtained for n = 2, 5, and 20, by the method in the previous section are indicated in Figure 8 as typical examples. Figure 8a shows the probability density functions under these three cases of n = 2, 5, and 20, whereas Figure 8b indicates the cumulative distribution functions corresponding to the respective cases. Weibull parameters a and b in Equation (8) are provided to satisfy the condition of P[1 µm < ρ < 15 µm] = 0.8. It is found that the peak (mode) and the standard deviation of the fatigue strength distribution tends to decrease with an increase in the inclusion number.  In order to investigate the effect of the condition for the dispersion of the inclusion size of , similar analyses were performed, giving some other conditions of [1 μm < < 10 μm] = 0.8 and [1 μm < < 5 μm] = 0.8. Among these analyses, only the results of the cumulative distribution functions are indicated in Figure 9, due to article page limit. As seen in Figures 8b and 9, median of the fatigue strength distribution tends to increase with a decrease of the dispersion of the inclusion size, while standard deviation of the fatigue strength distribution tends to decrease a little, depending on the decrease of the dispersion of the inclusion size.
Comparing these analytical distribution characteristics of the fatigue strength 9 * in Figures 8 and 9 with the experimental result in Figure 1, the distribution feature of the fatigue strength In order to investigate the effect of the condition for the dispersion of the inclusion size of ρ, similar analyses were performed, giving some other conditions of P[1 µm < ρ < 10 µm] = 0.8 and P[1 µm < ρ < 5 µm] = 0.8. Among these analyses, only the results of the cumulative distribution functions are indicated in Figure 9, due to article page limit. As seen in Figures 8b and 9, median of the fatigue strength distribution tends to increase with a decrease of the dispersion of the inclusion size, while standard deviation of the fatigue strength distribution tends to decrease a little, depending on the decrease of the dispersion of the inclusion size.

Expansion of the Probabilistic Model to Analyze P-S-N Characteristics in Interior-Induced Fracture
Target of the probabilistic model described in Section 3 is the distribution characteristics of the fatigue strength only at stress cycles of N = 10 9 . Accordingly, the analytical Comparing these analytical distribution characteristics of the fatigue strength σ * w9 in Figures 8 and 9 with the experimental result in Figure 1, the distribution feature of the fatigue strength σ * w9 under the conditions of P[1 µm < ρ < 15 µm] = 0.8 and n = 5 was roughly in agreement with the experimental distribution aspect of the fatigue strength at N = 10 9 . The most preferable number of inclusions, n, is further discussed again in Section 4.4.

Expansion of the Probabilistic Model to Analyze P-S-N Characteristics in Interior-Induced Fracture
Target of the probabilistic model described in Section 3 is the distribution characteristics of the fatigue strength only at stress cycles of N = 10 9 . Accordingly, the analytical model should be conceptionally expanded to interpret the whole statistical aspect of the S-N property, in the interior inclusion-induced fracture mode. From this point of view, denoting the fatigue strength at any stress cycle, N f ix , by σ w , Equation (3) is rewritten as follows: Based on the dotted S-N curve in Figure 1, parameters λ and γ in Equation (15) were determined from two points N = 10 6 , σ * w9 = 1424 MPa and N = 10 9 , σ * w9 = 920 MPa as λ = −168 and γ = 1512, respectively. Substituting these values into Equations (14) and (15), Equation (14) is rewritten as Here, σ w − ρ relationships under the fixed numbers of the fatigue life of N = 10 6 , N = 10 7 , N = 10 8 , and N = 10 9 are depicted in Figure 10. As suggested from Equations (14)- (16), analytical curves in Figure 10 tend to shift in parallel along the ordinate. Of course, the curve tends to shift downwards with an increase in the fixed number of stress cycles. Applying the concept of Equation (4), one could calculate the fatigue strength σ * w by Equation (17) includes arbitrarily given number of stress cycles , together with the inclusion size and the inclusion depth ξ. Therefore, combining the joint probability function in Figure 7 with Equation (17), one could analyze the distribution characteristics of the fatigue strength at any number of stress cycles by repeating the numerical calculations. Probability density functions of the fatigue strength thus obtained at N = 10 6 , 10 7 , 10 8 , and 10 9 are indicated at the respective numbers of the stress cycles in Figure 11. These analytical results are roughly in agreement with the overall feature of the experimental aspect on the statistical fatigue characteristics. The scale of the axis for the probability density is adjusted to correspond to each value listed in the table attached in the righthand side in Figure 11. The percentile points of = 1%, 10%, 50%, 90%, and 99% for the Inclusion radius (μm) Figure 10. σ wρ relationships at several given numbers of stress cycles.
Equation (17) includes arbitrarily given number of stress cycles N f ix , together with the inclusion size ρ and the inclusion depth ξ. Therefore, combining the joint probability func-tion in Figure 7 with Equation (17), one could analyze the distribution characteristics of the fatigue strength at any number of stress cycles N f ix by repeating the numerical calculations.
Probability density functions of the fatigue strength thus obtained at N = 10 6 , 10 7 , 10 8 , and 10 9 are indicated at the respective numbers of the stress cycles in Figure 11. These analytical results are roughly in agreement with the overall feature of the experimental aspect on the statistical fatigue characteristics. The scale of the axis for the probability density is adjusted to correspond to each value listed in the table attached in the right-hand side in Figure 11. The percentile points of F = 1%, 10%, 50%, 90%, and 99% for the fatigue strength are indicated by marks of 's along the vertical axis of the probability density functions. The respective thin dotted lines passing through the percentile points corresponding to the same probability give the P-S-N curves representing the statistical fatigue characteristics of this steel. All data points yield within the range of 1-99% in Figure 11. This fact suggests that the probabilistic model developed in this study has an availability to interpret the physical meaning of the statistical fatigue property in the very high cycle regime.

Analysis of the Fatigue Life Distributions in Interior-Induced Fracture Mode
Although distribution characteristics of the fatigue strength * at the given number of stress cycles were analyzed in the previous section, the fatigue life distributions could also be analyzed based on Equation (17), by giving any value of the stress amplitude * . Applying this method, fatigue life distributions were calculated at * = 1400 MPa, 1200 MPa, 1000 MPa, and 900 MPa, respectively. The probability density functions thus obtained were indicated at the respective stress levels in Figure 12. Since the fatigue life increased logarithmically with decrease of the stress level, the probability density at the peak of the density function tended to decrease such that the integrated value of the density function maintained unity. Values of the peak density at the respective stress levels are indicated on the right-hand side, in Figure 12. As shown in this table, the peak density varied drastically, depending on the stress level, due to the above described reason, but the distribution curves were drowned when they had the same height for the sake of convenience.

Analysis of the Fatigue Life Distributions in Interior-Induced Fracture Mode
Although distribution characteristics of the fatigue strength σ * w at the given number of stress cycles were analyzed in the previous section, the fatigue life distributions could also be analyzed based on Equation (17), by giving any value of the stress amplitude σ * w . Applying this method, fatigue life distributions were calculated at σ * w = 1400 MPa, 1200 MPa, 1000 MPa, and 900 MPa, respectively. The probability density functions thus obtained were indicated at the respective stress levels in Figure 12. Since the fatigue life increased logarithmically with decrease of the stress level, the probability density at the peak of the density function tended to decrease such that the integrated value of the density function maintained unity. Values of the peak density at the respective stress levels are indicated on the right-hand side, in Figure 12. As shown in this table, the peak density varied drastically, depending on the stress level, due to the above described reason, but the distribution curves were drowned when they had the same height for the sake of convenience.
In order to distinguish the fatigue strength distribution and the fatigue life distribution, the notations of f (σ w ) and F(σ w ) are used for the fatigue strength distribution, whereas other notations of g(σ w ) and G(σ w ) are accepted for the fatigue life distribution. Similar to Figure 11, the percentile points of G = 1%, 10%, 50%, 90%, and 99% for the fatigue life are indicated by marks of 's along the horizontal lines, at the respective stress levels. The respective thin dotted lines passing through the percentile points corresponding to the same probability give the P-S-N curves of this steel. All experimental data in the interior-induced fracture appear within the range of 1-99% in Figure 12. Thus, it was finally noted that the statistical fatigue property in the interior-induced fracture of this steel could be well explained from either one of the fatigue strength distribution and the fatigue life distribution, through the probabilistic model proposed in this study.
of stress cycles were analyzed in the previous section, the fatigue life distributions could also be analyzed based on Equation (17), by giving any value of the stress amplitude * . Applying this method, fatigue life distributions were calculated at * = 1400 MPa, 1200 MPa, 1000 MPa, and 900 MPa, respectively. The probability density functions thus obtained were indicated at the respective stress levels in Figure 12. Since the fatigue life increased logarithmically with decrease of the stress level, the probability density at the peak of the density function tended to decrease such that the integrated value of the density function maintained unity. Values of the peak density at the respective stress levels are indicated on the right-hand side, in Figure 12. As shown in this table, the peak density varied drastically, depending on the stress level, due to the above described reason, but the distribution curves were drowned when they had the same height for the sake of convenience.

Reconfirmation of the Number of Inclusions in the Critical Volume
As described in Section 4.1, since the number of inclusions included in the critical volume is unknown, the analysis was carried out by setting three conditions of n = 2, 5, and 20, as the inclusion number. Then, the condition of n = 5 was tentatively selected as a preferable number of inclusions. However, in order to confirm the most reasonable number of inclusions, the fatigue life distributions were additionally analyzed, under conditions of n = 6, 8, and 10, respectively. Based on the series of analytical results, percentile points of the fatigue life distribution corresponding to G = 1%, 10%, 50%, 90%, and 99% were determined; the relationship between those percentile points and the number of inclusions are depicted in Figure 13. Comparing five solid curves carefully with the dispersion aspect of the experimental results in an interior-induced fracture in Figure 12, it was finally noted that the analytical result at n = 6 was well fitted compared to the result at n = 5.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 16 In order to distinguish the fatigue strength distribution and the fatigue life distribution, the notations of ( ) and ( ) are used for the fatigue strength distribution, whereas other notations of ( ) and ( ) are accepted for the fatigue life distribution. Similar to Figure 11, the percentile points of = 1%, 10%, 50%, 90%, and 99% for the fatigue life are indicated by marks of ◇'s along the horizontal lines, at the respective stress levels. The respective thin dotted lines passing through the percentile points corresponding to the same probability give the P-S-N curves of this steel. All experimental data in the interior-induced fracture appear within the range of 1-99% in Figure 12. Thus, it was finally noted that the statistical fatigue property in the interior-induced fracture of this steel could be well explained from either one of the fatigue strength distribution and the fatigue life distribution, through the probabilistic model proposed in this study.

Reconfirmation of the Number of Inclusions in the Critical Volume
As described in Section 4.1, since the number of inclusions included in the critical volume is unknown, the analysis was carried out by setting three conditions of n = 2, 5, and 20, as the inclusion number. Then, the condition of n = 5 was tentatively selected as a preferable number of inclusions. However, in order to confirm the most reasonable number of inclusions, the fatigue life distributions were additionally analyzed, under conditions of n = 6, 8, and 10, respectively. Based on the series of analytical results, percentile points of the fatigue life distribution corresponding to = 1%, 10%, 50%, 90%, and 99% were determined; the relationship between those percentile points and the number of inclusions are depicted in Figure 13. Comparing five solid curves carefully with the dispersion aspect of the experimental results in an interior-induced fracture in Figure 12, it was finally noted that the analytical result at n = 6 was well fitted compared to the result at n = 5.
Based on this evidence, the statistical fatigue properties (P-S-N properties) in Figures  11 and 12 were numerically analyzed by putting n = 6 for the number of inclusions included in the critical volume of the specimen.

Mutual Relationship between Both Distributions of Fatigue Strength and Fatigue Life
In this study, the authors showed that both fatigue strength distribution and fatigue life distribution could be derived from the joint probability density function of ℎ( , ), given by Equation (11). It was also confirmed that the analytical results were in good Based on this evidence, the statistical fatigue properties (P-S-N properties) in Figures 11 and 12 were numerically analyzed by putting n = 6 for the number of inclusions included in the critical volume of the specimen.

Mutual Relationship between Both Distributions of Fatigue Strength and Fatigue Life
In this study, the authors showed that both fatigue strength distribution and fatigue life distribution could be derived from the joint probability density function of h(ρ n , ξ), given by Equation (11). It was also confirmed that the analytical results were in good agreement with the statistical aspect of the experimental fatigue data in the very high cycle regime. In this place, mutual relationship between the distribution aspects of the fatigue strength and the fatigue life is further discussed to reconfirm the reasonability of the present probabilistic model for the P-S-N characteristics. The mutual relationship between both distributions is illustrated in Figure 14, where the ordinate is the fatigue strength σ * w and the abscissa is the fatigue life N. When any point P(N fix , σ * w ) is taken in this diagram, the probability that the fatigue life at the stress level σ * w is less than N fix always corresponds to the probability that the fatigue strength at the definite stress cycles N fix is lower than σ * w . In other words, the cumulative probability of N fix , G(N fix ,σ * w ), always corresponds to the cumulative probability of σ * w , F(N fix ,σ * w ) [30]. G(N fix ,σ * w ) gives the distribution function of the fatigue life at σ * w and F(N fix ,σ * w ) gives the distribution function of the fatigue strength at N fix . Thus, we have G(N fix , σ * w ) = F(N fix , σ * w )  Based on this equality of Equation (18), one could obtain the S-N curve corresponding to any level of the fracture probability. Thus, the probabilistic S-N curves at the respective fracture probabilities could be analyzed as P-S-N curves. As described in Sections 4.2 and 4.3, the distribution characteristics of the experimental S-N data were well explained from either viewpoint of the fatigue strength distribution or the fatigue life distribution. The equality of Equation (18) is the reason why the reasonable P-S-N curves could be equally obtained from different viewpoints of the fatigue strength distribution and the fatigue life distribution.

Conclusions
Main conclusions obtained in this study are summarized as follows.
1. The probability density functions of the inclusion size at the crack initiation site, ( ), was successfully derived by combining the Weibull distribution and the concept of extreme distribution. In addition, the probability density function of the inclusion depth, 0 ( ), was also derived from the uniform distribution of the location of the inclusion in the material space. 2. Since the inclusion size and the crack depth ξ are statistically independent, the joint probability density function of these random variables, ℎ( , ), is given by the Equation (18) implies that the dashed area G(N fix ,σ * w ) is always equal to the other dashed area F(N fix ,σ * w ) in Figure 14. Based on this equality of Equation (18), one could obtain the S-N curve corresponding to any level of the fracture probability. Thus, the probabilistic S-N curves at the respective fracture probabilities could be analyzed as P-S-N curves. As described in Sections 4.2 and 4.3, the distribution characteristics of the experimental S-N data were well explained from either viewpoint of the fatigue strength distribution or the fatigue life distribution. The equality of Equation (18) is the reason why the reasonable P-S-N curves could be equally obtained from different viewpoints of the fatigue strength distribution and the fatigue life distribution.

Conclusions
Main conclusions obtained in this study are summarized as follows.

1.
The probability density functions of the inclusion size at the crack initiation site, f (ρ n ), was successfully derived by combining the Weibull distribution and the concept of