Effect of characteristics of non-metallic inclusions on fatigue life of high-speed rail bearings based on multiple linear regression

In the process of smelting and manufacturing, bearing steel inevitably have defects such as non-metallic inclusions. Moreover, inclusion geometry (type, quantity, size, distribution, and shape) has a significant effect on the fatigue failure of bearing steel. In order to study effect extent of the size and distribution of inclusion on the fatigue life of bearings, this paper uses the test high-speed rail bearings LY-Z028 to research and analysis. The micro-model of contact between single inclusion and the matrix is established, and then performs ANSYS simulation on stress concentration which caused by size of single inclusion and relative surface depth. Further, multiple linear regression method is used to compare the effect extent of these two factors on the fatigue life. The results show that size of single inclusion has a greater impact on the fatigue life than relative surface depth, which can provide an effective reference for smelting and production requirements.


Introduction
As one of the key components in the operation of high-speed railway trains, it is always inevitable that there are defects such as non-metallic inclusions in the process of smelting and manufacturing bearing steel. Due to evolution of these micro-defects, equipment shutdowns and structural failures caused great economic losses and safety accidents. Then people have studied effect of non-metallic inclusions on bearing fatigue durability. It is believed that existence of such defects destroy continuity of the matrix material, and easily occur local stress concentration. Because of stress concentration, the material lattice glides, and slip bands are generated inside the matrix, which accelerates the generation of initial fatigue cracks, and ultimately leads to a significant reduction in the fatigue life of bearings [1].
In addition, more and more researchers are turning their attention to the effect of various characteristics of inclusions on bearing fatigue life. Meng Bo and Dai Jingjun [2] studied the local micromechanical behavior of typical inclusions by using three-dimensional finite element method, then analyzed effect of the distribution state, location and size of inclusions on the local micromechanical behavior. Koumi [3] proposed a new numerical method to calculate stress distribution around an ellipsoid inclusion. The results show that when the ellipsoid's long axis symmetry is not parallel to the contact surface, stress distribution is asymmetric. Feng Lei, Xuan Fuzhen [4] analyzed effect of inclusions on local stress concentration and fatigue crack growth  [5] also studied effect of the type, shape, size and distance of inclusions on stress distribution.
In summary, series of studies have found that researchers have made many simulation models for stress concentration which generated by non-metallic inclusions, and have analyzed stress magnitude brought by various characteristics of inclusions to bearing subsurface. But effect extent of different characteristics of inclusions on fatigue life have not been deeply studied and analyzed, and it is difficult to provide a direct and effective reference for actual manufacturing. In view of the above limitations, this paper uses multiple linear regression analysis. A multiple linear regression equation which includes diameter of inclusions and relative surface depth as explanatory variables and the fatigue life of bearings as analysis object is established to compare their effect extent on bearing fatigue life.

Characteristics of non-metallic inclusions
Non-metallic inclusions in steel refer to oxides, sulfides, silicates or nitrides that do not have metallic properties in steel. In most cases, the existence of non-metallic inclusions has certain harmfulness, which can destroy the continuity of metal materials and thereby reducing the plasticity, toughness and fatigue life of metal materials.
At present, the classification of inclusions in bearing steel standards is mainly grouped into Class A (sulfide), Class B (oxide, including nitride), Class C (silicate), Class D (spherical oxide) and Class Ds (large single-particle Class D inclusions). Class A: single gray inclusion with high ductility and a wide range of aspect ratio (length/width), generally, the ends are rounded; Class B: most of them have no deformation , Angular, small form ratio (usually less than 3), black or bluish particles, arranged in a row along the rolling direction (at least 3 particles); Class C: a single black or dark gray inclusion with a wide range of aspect ratios (generally greater than 3) and high ductility, and generally an sharp angle at the ends; Class D(or Class Ds): non-deformed, angular or rounded, relatively small in shape (generally less than 3), black or bluish, irregularly distributed particles [6].
Commonly, non-metallic inclusions in steel are divided into hard inclusions and soft inclusions. Hard inclusions mainly include spherical oxides and nitrides, and soft inclusions mainly contain sulfides and silicates. Figure 1 shows the spherical inclusions in the bearing spalling pits which observed by SEM after the test roller bearing LY-Z028 has fatigued.

Effect of inclusions on subsurface stress of bearings
The existence and distribution of inclusions in bearing steel are complicated. In order to better probe effect of inclusions on subsurface stress distribution of the raceway, this paper simplifies the inclusion model and studies its effect on the subsurface stress of bearings from the perspective of single spherical inclusion.
When non-coordinated elastic contact objects (objects with dissimilar shapes) are in contact, their contact area is usually very narrow. At this time, contact stress decreases rapidly as the distance away 3 from the area, but it is highly concentrated near this area. Therefore, when studying contact characteristics of objects, the researches often focus on the contact area and its vicinity. Based on this conclusion, this paper processes each object as a semi-infinite elastic solid (also called elastic halfspace) bounded by a flat surface, so that the stress and strain can be calculated very approximately.
Consequently, based on the above theory, this paper simplifies the bearing model into an elastic half space, and considers stress distribution under the subsurface of bearing raceway when the most stressed roller is in contact with the inner and outer rings of bearing during the operation. According to Murakami's research, if the equivalent length S of double inclusions is equal (S refers to the shadow area of the defects on a vertical plane to direction of the maximum tensile stress), then its shape has little effect on fatigue work of the material. Hence, the equivalent spherical defect is used here instead of complex-shaped inclusion. According to the above assumptions, a single inclusion simulation model is established, as shown in Figure 2. Among them, Z inc means depth of relative contact surface, D inc means diameter of inclusion, b is the contact half-width, P(x) means surface contact stress between the most stressed roller and the raceway of ring, and the direction is from left to right. Calculated according to Hertz line contact theory, the expression of P(x) is (1)  M is maximum equivalent contact stress value, which is 742.128MPa by static simulation. The specific parameters of simulation model of inclusion in the elastic semi-infinite objects can be seen in Table 1. Table 1. Specific parameters of the simulation model of inclusion.

Effect of diameter of inclusion on the subsurface stress
According to related literature, current steel smelting technology can control the average diameter of inclusions below 15-20μm, but it is far from enough to improve fatigue life of bearings. Yang Zhenguo deduced and verified rationality of the "critical inclusion" theory of high-strength steel based on the inclusion shadow area model which proposed by Murakami. The "critical inclusion" theory considers that high-strength steel has a critical inclusion size. If actual inclusion size is larger than the critical inclusion, fatigue cracks are likely to initiate from the inclusion; on the contrary, actual inclusion size is less than the critical inclusion size, it is very likely that fatigue cracks will not initiate at the inclusion. The calculation formula of critical inclusion size derived from this theory is: Among them, means the critical inclusion size; C is the modified inclusion position coefficient, when inclusion is on the subsurface, the value is 0.528; HV means the Vickers-hardness of bearing steel.
At present, according to relevant research, the Vickers-hardness of GCr15 is above 720, so under Formula (2), diameter of the critical inclusion can be calculated to be less than 1μm. Based on the above data, this paper mainly considers diameter of inclusion in the range of 1-20μm, and changes the diameter to reflect impact of subsurface stress.
Assuming that relative surface depth of inclusions remains unchanged (the depth is 0.6b), the Poisson's ratio is equal to 0.29 and the friction coefficient between inclusion and steel matrix is 0.25, changing diameter of inclusion (from 2μm to 16μm), this paper considers the subsurface dynamic shear stress and Von Mises stress distribution in two cases of soft inclusion (105GPa) and hard inclusion (315GPa). Figure 3 shows dynamic shear stress near the inclusion when the size is 6μm.  It can be observed from the Figure 3 that whether it is a soft or hard inclusion, the maximum point of dynamic shear stress is distributed in the surrounding contact area, stress concentration is obvious to generate, and the maximum dynamic shear stress is basically located at about ±30 o and ±150 o of inclusion. On the basis of the simulation results, when relative surface depth of inclusion is constant, and the diameter of inclusion gradually changes from 2μm to 16μm, stress distribution structure diagram is basically the same as Figure 3.
For Von Mises stress, following Figure 4 shows corresponding stress distribution diagram when the size of inclusion is 6μm.  Next, drawing a line chart by MATLAB, as shown in Figure 5, which visually shows effect of different inclusion diameters (varying from 2μm to 16μm) on the subsurface dynamic shear stress and Von Mises stress distribution. Under the above figure, when keeping relative surface depth of inclusion unchanged and the diameter changed, the maximum stress on the subsurface of contact area increases if diameter of inclusion increases. In general, stress concentration caused by hard inclusion is greater than that of soft inclusion. This is because soft inclusions have certain plastic toughness, and counteract a less stress.

Effect of relative surface depth of inclusion on the subsurface stress
Related references record that when relative surface depth exceeds a certain size, impact of inclusions on the bearing subsurface stress is inessential [7]. Therefore, based on the model shown in Figure 2, this paper maintains the elastic modulus of inclusion at 315GPa, 6μm in diameter and the friction coefficient between inclusion and steel matrix is 0.25, while changing relative surface depth of inclusion (increasing from 0.1b to 1.5b) to analyze stress change. And Figure 6 shows dynamic shear stress and Von Mises stress distribution if relative surface depth is 0.6b.  It can be seen from the above figure that near the inclusion and their contact area, dynamic shear stress and Von Mises stress reach the maximum, and stress concentration is obvious to generate. And according to the simulation results, as relative surface depth of inclusion increases, stress distribution is similar.  Figure 7 shows the maximum dynamic shear stress and Von Mises stress varying with changes of relative surface depth. It can be seen from the above figure that if relative surface depth of inclusion is in the range of 0.4b to 0.9b, extent of effect on the stress is not very large. When relative surface depth is large enough, the stress gradually decreases.
Prior to this, Bryan Allison and Anup Pandkar [8] also did related researches. By comparing and referring to both of the maximum Von Mises stress, the maximum stress values obtained by simulation in this paper are relatively reliable.

Prediction and calculation of fatigue life of bearings
Traditional calculation method of bearing life is no longer applicable to the non-uniformity model containing inclusions. For the non-uniformity model, there are two main types of methods to determine the fatigue life: test method and analysis method. Although the test method has high reliability, it is still difficult to grasp some local or micro variables, and very costly. The analysis method is to establish a certain model to calculate fatigue life based on fatigue performance of the bearing material and load on the material or structure [9][10][11][12].
There are many kinds of fatigue life analysis methods, but the mainstream analysis methods include stress-strain field intensity method, local stress-strain method and so on. According to the advantages and disadvantages of these methods, after comprehensive consideration, this paper uses local stress-strain method to analyze effect of diameter of inclusion and relative surface depth on bearing life.
According to local stress and strain theory, in the constant fatigue test, the Basquin formula based on equivalent stress is usually used for life prediction. The Basquin formula is (3) In this formula, N f means the fatigue life, and  ae is equivalent stress amplitude (Von Mises stress is used in this paper). Considering the average stress,  ae can be calculated by the following formula: Where  ad is equivalent stress amplitude when average stress is zero,  me is average stress, and  b is ultimate tensile strength of the material.
In Formula (3), means fatigue strength coefficient. If tensile fatigue strength limit of the material is unknown, it can be estimated by = b +350. For bearing steel GCr15,  b =1617MPa, where can be taken as 1967MPa.  means fatigue strength coefficient, under slope method, and  = -0.12.
The Basquin formula is an idealized formula based on material properties. Taking into account non-uniformity of materials, it is necessary to modify Formula (3) using fatigue strength reduction factor K f . Modified formula for this: In the formula, a means the Neuber constant, and its value is related to  b , as shown in Figure 8.
When the material is GCr15, a is 0.0025; ρ is root radius of the notch, where it can be replaced by inclusion radius R inc ; K T is the microscopic stress concentration factor, which can be obtained from the formula: Among them,  max means the maximum local elastic stress and  0 is nominal stress, where K T can be obtained from the finite element simulation results above. According to the above calculation method of fatigue life prediction, the graph is drawn by using MATLAB, which shows the correlation more intuitively. The result is shown in Figure 9: (a) diameter (b) relative surface depth Figure 9. Graph of effect of inclusion diameter and relative surface depth.
From the above graph, it can be seen that if keeping relative surface depth of inclusion unchanged, the predicted fatigue life decreases significantly as diameter of inclusion increases; If diameter of inclusion unchanged, when relative surface depth is between 0.4b and 0.9b, the predicted fatigue life fluctuates around 3210 7 cycles, and the extent of effect is small.

Model construction
According to Figure 1, the inclusions of LY-Z028 bearing are generally spherical oxides, which are hard inclusions by microscopic observation after fatigue failure. In the following, this paper takes hard spherical oxide inclusions as the research object, and predicted fatigue life of high-speed rail bearings N f is interpreted object, and the diameter of inclusion x 1 and relative surface depth x 2 are used as explanatory variables. Then, using multiple linear regression analysis to compare effect extent of two explanatory variables on fatigue life of bearings, a regression equation is established. Suppose N f has the following relationship with x 1 and x 2 : In the above formula,  is random disturbance term.  1 and  2 are regression coefficients, and their absolute values indicate extent of interpretation of the variable to predicted fatigue life N f , and the positive and negative indicate direction of the variable's influence on N f . In order to accurately study the more important factors that affect N f , effect of dimensions needs to be eliminated. A standardized regression model is used here [13][14][15][16].

Model checking
Before performing regression analysis, this paper first examines a multicollinearity test on variables. As shown in Figure 10, the VIF values are all far less than 5, and there is no multicollinearity, that is, there is no mutual interference relationship between all explanatory variables . Then, the explanatory variables are tested for heteroscedasticity. As shown in Figure 11, the P value obtained by BP test is bigger than 0.05, indicating that the hypothesis is accepted at a 95% confidence level (the hypothesis is that disturbance term does not have heteroscedasticity), that is, disturbance term is considered to have no heteroscedasticity. Through the above test and analysis, the established regression equation is statistically significant.

Result of regression analysis
This paper uses Stata15.0 to perform multiple linear standardized regression after inputting the relevant data. The standardized regression coefficients are used to compare strength of the explanatory variable's influence. The result is shown in Figure 12. It can be drawn from the figure: at the 95% confidence level, Prob=0.0022 is less than 0.05, indicating that the hypothesis is rejected (the hypothesis is  1 = 2 =0), so the equation has significant regression coefficients  1 and  2 . The final standardized linear regression equation is as follows: In this model, the absolute value of  1 is bigger than  2 , and it can be analyzed that the diameter of inclusion x 1 has a greater impact on the fatigue life of bearing than relative surface depth x 2 . Therefore, in the smelting and manufacturing process of bearing steel, it is necessary to control and reduce size of non-metallic inclusions as much as possible to improve fatigue life of bearing effectively.

Conclusions
Based on the results of finite element analysis, inclusions size or distribution cause local stress to rise, and the stress concentration occurs on the subsurface of bearing steel, which eventually leads to bearings failure. For the study of effect extent of single spherical inclusion geometry (size and distribution) on bearings fatigue life, this paper uses multiple linear regression method. And through the above research, the following points are drawn: (1) Within a certain range of single spherical inclusion (diameters between 8-16m), stress concentration caused by hard inclusion is greater than soft inclusion.
(2) If relative surface depth of single spherical inclusion is in the range of 0.4b to 0.9b, extent of effect on the stress is not very large. If relative surface depth is deep enough (much bigger than 1.6b), stress concentration caused by it becomes less and less.
(3) Based on the linear regression equation obtained above, it is concluded that under certain conditions and ranges, the size of single spherical inclusion has greater effect on fatigue life than relative surface depth.
However, effect of single spherical inclusion geometry (size, distribution) on the overall fatigue limit is relatively small. There are other influence factors such as initial scratches on the work surface [17][18], inclusion orientation, heterogeneity of the bearing steel, temperature effects, etc. This means that if only the contact model used above is considered, then a much simpler analytical method is sufficient to estimate the fatigue limit of bearings.