Equivalent method for assessment of bearings at operational types of loadings

Abstract

To investigate the effect of various factors on bearing stress response, Huber-Hencky-von Mises stress serves as a bridge, the equivalent interrelation between radial loading, axial loading, and temperature of bearing is studied using finite element method (FEM). Symbolic regression (SR) algorithm is employed to analyze simulation results, establishing a functional expression between independent and dependent variables by optimizing combinations of variables, constants, and functional forms. The results showed that within the specified force and temperature values, the curved surface of the equivalent correlation function, trained using the SR algorithm, demonstrates smoothness. Both training and validation data exhibit a strong correlation with this curved surface. Among the three factors, temperature exerts the greatest influence on bearing stress values, followed by radial loading, and axial loading components had the smallest impact.

Keywords

Bearing FEM loading-temperature equivalence symbolic regression algorithm correlation function

Introduction

The design of bearing reliability for a high-speed railway transmission system necessitates consideration of key factors such as loading, temperature, vibration, and flow field. To verify whether the bearing meets the reliability design standard, stand testing of the high-speed rail transmission systems is imperative. However, effectively simulating equivalent multi-field coupling service conditions poses a challenge. The bearing stand test is an essential experiment link before the actual service of the bearing. Only through the stand test can the safety and reliability of the bearing during operation be guaranteed.

The finite element method (FEM) is extensively utilized in structural engineering strength analysis due to its convenience, practicality, and effectiveness. However, extracting key information and establishing correlation functions from limited data analysis poses a current challenge.^1–7 Addressing the dynamic coupling effect of the radial force, axial force components, and bending moment of the bearing, Zhao et al.⁸ proposed a prediction method for the service life of the hub bearing under multi-axis loading conditions. Hu et al.⁹ introduced an energy method to solve address shock absorber lifespan concerns under different types of loading, with test results correlating with actual service mileage. Lang et al.¹⁰ established a coupling analysis model for the thermal effect of sliding bearing and accurately determined the equivalent temperature of the bearing’s flow field.

It is necessary to obtain linear or nonlinear fitting expressions for calculation and analysis when dealing with multi-factor interaction problem of in engineering. Machine learning offers a valuable tool for deriving more accurate fitting functions from limited samples. Serving as a significant post-processing analysis method, machine learning effectively addresses the challenge of maximum likelihood estimation with small specimens and enhances generalization performance. Leveraging powerful nonlinear analysis and multivariate processing capabilities, machine learning technology has enabled the widespread application of predictive models across various fields.^11–16 Liao et al.¹⁷ utilized a genetic wavelet neural network to investigate the behavior of fatigue short crack propagation, conducting experiments at two loading frequencies. They calculated the effective short crack density and its propagation rate. The findings demonstrated the reasonable and effective nature of the genetic wavelet neural network in studying fatigue short crack evolution behavior. Zheng et al.¹⁸ established prediction models using four machine learning algorithms, including random forest, ridge regression, and support vector regression based on radial basis kernel function and linear kernel function. These models aimed to predict the density, formation energy, band gap, and crystal volume of perovskite material datasets. The results indicated that various machine learning algorithms exhibited strong prediction capabilities for the given data samples. However, most of the current machine learning methods such as neural networks, random forests, and support vector regression models are black box models that are difficult to explain, and they only show good prediction performance in the range of data sets used for training, with poor extrapolation and interpretability. Establishing a machine learning prediction model that can be extended beyond the training data set is the current focus and difficulty. Schmidt and Lipson¹⁹ derived the equations of motion of a simple pendulum and double pendulum by symbolic regression (SR) algorithm, and they found that even if the sin and cos terms in the initial function of the equation are removed, the final results of SR would appear Taylor expansions of sin and cos, respectively. In contrast to traditional machine learning algorithms, the SR algorithm operates under a single assumption: that there exists a functional relationship between the input and output data. Importantly, the symbolic algorithm serves as an explanatory white box model, with its output being an equation represented by variables, constants, and functions.^20–22

In assessing the service reliability of transmission system bearings, numerous factors come into play. The impact of these factors on bearing stress response is investigated using a random forest algorithm, which is commonly employed to assess the importance of features in a dataset. This algorithm creates multiple decision trees by randomly selecting subsets of features from the data and integrates the results from these trees to determine the contribution of each parameter to the overall accuracy of the model.^23,24 The key to accurately evaluating the service performance of bearings lies in determining whether the life of transmission system bearings under operating conditions is equivalent to that under bench test conditions. Simplifying the stand test loading scheme and equivalent multiple factors to a single factor through the equivalent method is crucial, as it provides essential technical support for the stand equivalent test and comprehensive evaluation of the transmission system bearing’s service performance.

In this paper, a multi-factor and multi-level simulation scheme was designed using an orthogonal test method. The Huber-von Mises-Hencky stress of the concerned part was used as the bridge, and a finite element calculation data set was established. The SR algorithm was then applied for data analysis and processing. This approach aimed to explore and establish the equivalent correlation function between radial and axial loading components, as well as their changes versus temperature values. The technical route of the simulation research is illustrated in Figure 1.

Figure 1.

Research technical route.

Principle and algorithm

Loading-temperature equivalence principle

The internal contact state of the high-speed train bearing is extremely complex in the actual service process, and it is difficult to accurately grasp, which directly affects the safety of the train operation. Therefore, the thermal-mechanical coupling model of the bearing is established to analyze the change in the internal contact state of the bearing by setting different loading-temperature fields. Based on the second law of engineering thermodynamics, heat can be transferred from a high-temperature region to a low-temperature region, while a low-temperature region cannot spontaneously transfer heat to a high-temperature region. Heat transfer always occurs when there is a temperature difference. Temperature differences in objects are widespread in nature, and the heat transfer phenomenon is very common in practical life.²⁵ Heat transfer occurs through three primary modes: conduction, convection, and radiation, which collectively constitute all heat transfer processes. Compared with the approximate calorific value model, the Palmgren calorific value model is closer to the experimental results and has higher accuracy. Therefore, the Palmgren calorific value model is selected to calculate the heating power of the axle box bearing of the high-speed train. The equation is as follows:

M = M_{i} + M_{v}

(1)

M_{i} = f_{1} F_{β} d_{m}

(2)

M_{v} = (\begin{matrix} 10^{- 7} f_{0} {(v_{0} n)}^{2 / 3} d_{m}^{3}, v_{0} n \geq 2000 \\ 160 \times 10^{- 7} f_{0} d_{m}^{3}, v_{0} n \leq 2000 \end{matrix})

(3)

where M_i and M_v are the friction torque caused by the external load and the viscous friction of the lubricant respectively. f₁ is the load coefficient, and f₁ = ¢ (Fs/Cs) £, ¢ and £ are the coefficients related to the bearing structure type, F_s is the equivalent static load, C_s is the basic rated static load, d_m is the diameter of bearing pitch circle. For tapered roller bearings, F_β = 0.4F_r + 0.4cosαF_a, F_a is the radial force, F_r is the axial force, α is the contact angle. v₀ is the kinematic viscosity of the lubricant, n is the bearing speed. For the entire bearing, the heat source continuously passes through the contact area. Therefore, the heat flux density can be calculated on average based on the heat generated on the surface of the contact area. The equation is as follows:

q = H / S

(4)

H = π n (M_{i} + M_{v}) / 30

(5)

S = 2 l \sqrt{\frac{4 F (\frac{1 - μ_{1}^{2}}{E_{1}} + \frac{1 - μ_{2}^{2}}{E_{2}})}{π l \frac{1}{ρ \sum 1}}}

(6)

where H is the calorific value, S is the contact area, μ₁ and μ₂ are the Poisson’s ratio of the outer ring (inner ring) and the roller material, respectively, E₁ and E₂ are the Young’s modulus of the outer ring (inner ring) and the roller material respectively, l is the contact length between the roller and the outer ring (inner ring), 1/ρ∑¹ is the comprehensive curvature radius, F is an external load. At present, the bearing convection heat transfer model is mainly Harris convection heat transfer model.²⁵ The heat transfer coefficient between the traditional bearing cavity and the lubricant given by the Harris convection heat transfer model is as follow:

α = 0.332 {kP}_{r}^{1 / 3} {(\frac{u}{vx})}^{0.5}

(7)

where k is the thermal conductivity, p_r is the Planck constant, u is the following velocity of lubricating oil relative to the heat exchange surface, v is the viscosity of lubricating oil, x is the characteristic length. For the inner raceway, u takes the speed of the inner ring, x takes the diameter of the inner raceway. For the outer raceway, u takes the cage speed and o takes the bearing pitch diameter.

In summary, the corresponding stress, strain, and displacement components can be obtained after applying loading or temperature to the bearing. Taking the Huber-von Mises-Hencky (HMH) stress of the concerned part as the equivalent bridge, the equivalence between radial force and axial force can be realized. The effect of the temperature field caused by the existence of grease between the roller and the rolling road inside the bearing varies with speed. The internal temperature of the bearing is set as a fixed temperature when establishing the thermo-mechanical coupling model, and the internal contact characteristics of the bearing are studied by setting different temperatures. Similarly, the HMH equivalent stress of the concerned part is used as the equivalent bridge to realize the equivalence between radial loading and temperature.

SR algorithm

Orthogonal experimental design is an effective approach for multi-factor experiments, several representative points are selected from the comprehensive test according to the orthogonality of the orthogonal table, and these points have the characteristics of uniform dispersion and neat comparability, but the data are limited and it is difficult to determine the trend of the test. The symbolic regression algorithm can establish a symbolic model for nonlinear systems without relying on prior knowledge or models, so the mathematical relationship expression between various factors can be established by using the SR algorithm.^26,27

SR does not require the assumption of a specific functional form between independent variables and dependent variables and only assumes that their relationship can be expressed functionally. SR searches for the optimal solution by combining variables, constants, and various functions.²⁸ The essence of SR is based on gene programming technology, which is similar to Darwin’s natural selection principle. As shown in Figure 2(a), any formula can be expressed in the form of a binary tree. The structure of the SR tree includes terminal nodes and internal nodes. The terminal nodes represent constants and variables, and the internal nodes represent functional relations.

Figure 2.

SR binary tree, crossover, and variation: (a) symbolic regression binary tree, (b) crossover, and (c) mutation.

The initial population of SR is randomly generated according to the defined function, and the individuals are screened according to their fitness. Individuals with higher fitness are more likely to be selected for the next generation, where they evolve through crossover and mutation. As shown in Figure 2(b) and (c), crossover refers to the random exchange of subtrees between the previous generation equations to generate new offspring individuals. Mutation means that one or more nodes in the equation are randomly adjusted to ensure population diversity to explore better data-fitting equations. Fitness evaluation, selection, crossover, and mutation occur in each generation to generate new individuals. In this continuous iterative process, the complexity of the equation increases with the number of iterations. Concurrently, the mean squared error (MSE) gradually decreases. However, the score improvement will gradually diminish, become more variable and lead to instability after a certain complexity level. The equation may also exhibit singular terms. Therefore, the criterion for terminating the algorithm is when the MSE is less than a certain value, the score improvement is stable, and the equation does not exhibit singular terms. By setting the individual complexity and function’s form to control the iterative solution process, the higher the complexity means the more complex the solution function. The function form setting can be set by referring to the function contained in the prior formula, and can also be set artificially.^29,30

Abaqus 2020 commercial finite element analysis software was used for the simulation test, simulation result data and a data set were collected. Randomly select 80% of the data in the data set, use the SR algorithm for training, randomly generate the first generation of individuals, and then use the remaining 20% of the data in the data set as the test set to test the accuracy of the training model. The training model is reasonably evaluated by setting the complexity of the correlation function and judging the MSE value.

To facilitate the practical application of engineering, we defined the function involved in the result expression in the early stage, that is, addition, subtraction, multiplication, and division in this study. The complexity of the model was set to 20. Individuals with higher fitness can generate new individuals through cross-mutation, and then evaluate their fitness. The number of cycles was equal to 1500. After completing the cycle, the equation with the highest fitness was finally output as the correlation function, the SR algorithm process is shown in Figure 3. The characteristic variables are input, the input variables are random combinations of radial and axial force or temperature, and the output variables are the maximum stress value at a certain point or structure.

Figure 3.

SR flow chart.

Numerical process

Finite element model

The tapered roller bearing is a widely used type of bearing in various mechanical equipment. It is known for its strong bearing capacity, high-temperature resistance, and high rotation accuracy. According to the size of the CRH380 axle box bearing provided by Luozhou LYC Co. Ltd., the main geometric parameters of the bearing are shown in Table 1. The three-dimensional model and finite element model were established by commercial software, as shown in Figure 4. Figure 4(a) is the geometric model of axle box bearing, and Figure 4(b) is the finite element model of axle box bearing. The material properties of each part of the finite element model are shown in Table 2. The material of the inner ring and outer ring is G20CrNi2MoA, and the material of the roller is GCr15. To ensure efficiency and accuracy in calculating the HMH stress values under varying operational conditions, a sub-model was used in the numerical process. In the finite element modeling, the outer surface of the outer ring of the bearing was fixed, the inner ring end face was axially constrained. The circumferential motion of all rollers is constrained, the bearing force of different types of loading ARE applied to the inner surface of the bearing inner ring. The temperature boundary is applied to the roller surface and the raceway, and the convective heat transfer coefficient is applied to the remaining surface.

Table 1.

Main geometric parameters of axle box bearings.

Bearing geometric parameter	Symbol	Value
Outside diameter	D (mm)	240
Bore diameter	d (mm)	130
Pitch circle diameter	d _m (mm)	187.2
Number of single row rollers	Z	17
Roller average diameter	D _m (mm)	26.69
Outer ring contact angle	α (°)	12.08

Figure 4.

CRH380 axle box bearing schematic diagram: (a) geometric model, (b) finite element model, (c) bearing sub-model, and (d) focus position.

Table 2.

Finite element model material information.

Component	Elasticity modulus (MPa)	Density (kg/m³)	Poisson ratio
Roller	2.19 × 10⁵	7.88 × 10³	0.29
Cage	2.59 × 10³	1.12 × 10³	0.3
Cone	2.05 × 10⁵	7.83 × 10³	0.3
Cup	2.05 × 10⁵	7.83 × 10³	0.3

Considering the characteristics of the finite element model of tapered roller bearing and the configuration of the computer, the inner and outer rings were divided by the Hex dominant method, and the global mesh size was defined as 1.2 mm. The roller is meshed by Multizone, and the global mesh size is defined as 1.0 mm. The unit type is C3D8R three-dimensional linear hexahedral unit. Most of the numerical models are hexahedral meshes, and also contain a small number of tetrahedral meshes. The order of the finite element is second-order and quadratic. The final bearing finite element model consists of 896,752 mesh elements and 942,815 mesh nodes. Because the maximum stress position on the same contact line can change under different working conditions, two situations are distinguished in data processing, namely: (1) The maximum stress on the contact line. (2) The stress of a fixed point on the contact line, as shown in Figure 4(c) and (d). For the convenience of analysis, each contact line position was named A/B/C/D, and the arrow direction in the figure was the load loading direction. The fixed point is the maximum stress position on the contact line under only radial loading. After the fixed point is determined, the stress value of the fixed point and the maximum value of the equivalent stress on the contact line are extracted for analysis when the radial loading-axial loading-temperature factors are equivalent.

In the process of bearing temperature field analysis, it is necessary to define the appropriate contact mode.³¹ Otherwise, it is difficult to accurately reflect the internal temperature characteristics. According to the actual working conditions of the bearing, two sets of contact pairs need to be set up, which are the contact pairs between the roller surface the outer ring raceway and the inner ring raceway. In the roller-outer ring contact pair, the outer ring raceway is selected as the target surface, and the conical surface of the roller is selected as the contact surface. In the roller-inner ring contact pair, the inner ring raceway is selected as the target surface, and the conical surface of the roller is selected as the contact surface. The contact types are all non-separating contacts, and the contact algorithm uses the augmented Lagrangian algorithm. The contact area settings are shown in Table 3.

Table 3.

Finite element model contact attribute information.

Component	Normal contact properties	Tangential contact properties	Tangential friction coefficient
Cup-roller	Hard contact	Lagrange	0.01
Cone-roller	hard contact	Lagrange	0.01
Cage-roller	Hard contact	Lagrange	0.005

Working condition setting

The orthogonal experimental design method was adopted by many researchers because of its high efficiency, rapidity, economy, and statistical background.^32–34 In the design process of the simulation scheme, the influence of axial, radial loading components and temperature on the HMH stress of each contact line position inside the bearing is focused on.

The orthogonal experiment method selects some representative points from the comprehensive experiment according to the orthogonality, which has the characteristics of uniform dispersion, neat, and comparable, and is an efficient and rapid test method. When considering the equivalence of radial and axial force components or the equivalence of radial loading and temperature, the traditional orthogonal method was used to set the working conditions, and 25 working conditions can be set. When considering the three-factor five-level design, 25 working conditions can be obtained by using the DOE design method, as shown in Table 4. Factors 1–3 in the table represent radial load (30/40/50/60/70 kN), axial load (−15/−10/0/10/15 kN), and temperature (40°/60°/80°/100°/120°) respectively, and 1–5 represent five different levels.^35–37

Table 4.

Orthogonal design scheme.

Number	Factor1	Factor 2	Factor 3	Number	Factor1	Factor 2	Factor 3
1	1	1	1	14	3	4	5
2	1	2	3	15	3	5	2
3	1	3	5	16	4	1	3
4	1	4	2	17	4	2	5
5	1	5	4	18	4	3	2
6	2	1	5	19	4	4	4
7	2	2	2	20	4	5	1
8	2	3	4	21	5	1	2
9	2	4	1	22	5	2	4
10	2	5	3	23	5	3	1
11	3	1	4	24	5	4	3
12	3	2	1	25	5	5	5
13	3	3	3

Results and discussion

Radial loading-axial loading equivalence

The traditional test bench for axial loading-radial loading dynamic actuators is complex and expensive. Establishing a reasonably designed test bench with limited funds is a challenge. Understanding the equivalent relationship between radial loading and axial loading is crucial for establishing a simplified single-bearing test bench.

To facilitate comparison, the complexity and MSE value of all the training models obtained by the SR algorithm are placed in Table 5. The subscripts FP and MP in the table represent fixed points and maximum points, respectively. 80% of the training data and 20% of the verification data were drawn in the correlation function surface map. The drawing results are shown in Figure 5, where (a)–(h) are the A/B/C/D data of the contact line, Figure 5(a) to (d) is the stress of the fixed point on the contact line, and Figure 5(e) to (h) is the maximum stress on the contact line.

Table 5.

Correlation function complexity and mean square error value.

	A_FP	B_FP	C_FP	D_FP	A_MP	B_MP	C_MP	D_MP
Radial loading-axial loading equivalence
Complexity	19	17	17	19	19	19	19	15
MSE	0.025	0.119	0.007	0.001	3.959	102.803	822.649	29.479
Radial loading-temperature equivalence
Complexity	19	19	19	19	19	19	19	17
MSE	46.828	62.601	26.771	19.589	97.868	92.102	5809.194	5.115
Radial loading-axial loading-temperature equivalence
Complexity	19	19	19	19	19	19	17	17
MSE	39.219	79.851	27.458	39.696	93.509	82.239	777.568	9.871

Figure 5.

Radial-axial loading equivalent correlation function surface: (a–d) corresponding to fixed point results and (e–h) corresponding to maximum point results.

The stress value in Figure 5(f) is observed to be 600 MPa, significantly lower than the stress value in Figure 5(g) (1000 MPa). These positions correspond to the contact lines of the roller with the outer and inner rings, respectively, with identical mechanical parameters. The difference in stress values is attributed to the roller’s force. As per the contact line diagram in Figure 4(d), the tapered roller experiences three contact forces from the inner and outer raceways and the guiding large rib, that is, the roller-ring raceway and the roller-inner ring rib load. During actual service, the maximum force occurs at the roller-inner ring contact point, so the observed results can be explained. All the equivalent correlation function surface obtained by the SR algorithm was smooth, and the training data and verification data scatters had a good correlation with the correlation function surface in Figure 5. With the increase in value of axial loading, the stress value of the fixed position of each contact line does not change significantly, and the maximum value of each contact line increases approximately linearly. The roller shape of the axle box bearing is conical. After the axial loading is added, the roller slides slightly in the axial direction, resulting in a linear gradient change of the stress value.

After obtaining the correlation function between the radial and axial loading components, the axial force equivalent to the radial one using the stress value as the intermediary. This is achieved by maintaining the Γ value before and after the equivalent, with modifications made to either the x or y value. The modifications serve the purpose of achieving the equivalent loading. The engineer only needs to clarify the magnitude of the radial and axial force, and then substitute it into equations (8)–(15) to obtain the stress value of the concern point. If the axial force is expected to be equivalent to the radial force, only the axial one (y value) is regarded as 0, the Γ value is unchanged, and the x value is adjusted to be equivalent.

The expression of the correlation function trained by the SR algorithm is shown in equations (8)–(15), where Γ is the contact line stress value (MPa), x and y represent the radial force (kN) and the axial force (kN), respectively. The MSE values of B_MP and C_MP in Table 5 were larger than other points. The reason is that the stress value on the B contact line and the C contact line corresponding to the maximum point is larger than the stress value of the contact line in other working conditions. The stress values on the B contact line and the C contact line are 600 and 1000 MPa, respectively. The MSE value in this paper is defined by the mean square sum of the error between the simulation value and the training value. Taking the C_MP value as an example, C_MP = 822.649, the average error between the simulation results at the maximum point of contact line C position and the training results of the SR algorithm is √822.649 = 28.682 MPa, which is only 2.868% of the maximum stress value of contact line (the maximum stress value of contact line is 1000 MPa), and the error range is acceptable in engineering practice. Therefore, the mathematical relationship fitted by the SR algorithm can approximately describe the equivalent relationship between the radial loading and the axial loading. After the single stand test of the bearing is calibrated, the above expression can be used to equivalent the radial loading-axial loading.

Contact line A_FP:

\begin{matrix} Γ = - 7.205 \times 10^{- 4} \cdot x^{2} + 1.109 \cdot x - 5.894 \times 10^{- 5} \cdot xy \\ + 6.918 \end{matrix}

(8)

Contact line B_FP:

Γ = 5.377 \cdot 10^{- 3} \cdot x^{2} - 2.441 \times 10^{- 5} \cdot xy + 0.7071 \cdot x + 62.81

(9)

Contact line C_FP:

Γ = 1.536 \cdot x + 7.008 \times 10^{- 2} \cdot y + 31.24

(10)

Contact line D_FP:

\begin{matrix} Γ = 2.269 \times 10^{- 5} \cdot x^{3} + 2.837 \times 10^{- 6} \cdot x^{2} y + 1.013 \cdot x \\ - 9.406 \times 10^{- 3} \cdot y - 0.5417 \end{matrix}

(11)

Contact line A_MP:

\begin{matrix} Γ = - 0.2056 \cdot y^{2} - 1.7027566 \times 10^{- 2} \cdot xy + 2.587 \cdot x \\ + 3.118 \cdot y + 39.97 \end{matrix}

(12)

Contact line B_MP:

Γ = 0.3614 \cdot y^{2} - 24.46 \cdot y + 8.155 \cdot x + 8.927

(13)

Contact line C_MP:

Γ = 0.2755 \cdot y^{2} - 9.839 \times 10^{- 2} \cdot xy + 4.729 \cdot x - 7.459 \cdot y

(14)

Contact line D_MP:

Γ = 6.199 \times 10^{- 2} \cdot y^{2} - 2.829 \cdot y + 1.964 \cdot x - 9.405

(15)

Radial loading-temperature equivalence

To simulate the temperature rise effect, the environmental box needs to be set up when establishing a single bench test bench. If the influence of temperature on the contact line stress inside the bearing can be equivalent to the radial force, it can be much easier to design a single-bearing bench. Considering that the internal temperature of the bearing fluctuates significantly at different speeds, the nonlinear influence of the internal temperature rise of the bearing is not considered when establishing the thermal-mechanical coupling model. The internal temperature of the bearing is set as a fixed value, and the influence of nonlinear temperature rise is considered in the later research. The existence of internal defects in the bearing makes the stress value near the contact line position rise sharply. In this study, the abnormal temperature rise caused by the damage to the inner ring raceway, outer ring raceway, or roller is not considered.

The contact stress corresponding to each contact line is different under different damage degrees in the bearing. Therefore, the temperature range is set to 40°C–120°C in this study. In the later stage, the change of the stress value of each contact line under different damage degrees in the bearing under higher temperature conditions can be studied. Based on the above assumptions, a total of 25 groups of radial loading-temperature conditions were set up. Using the method described in Section “Finite element model,” when the radial loading-temperature is equivalent, the complexity and MSE value of the training model can be obtained by the SR algorithm. The model expression obtained by training is shown in equations (16)–(23), and the surface fitting effect and verification results of the equivalent correlation function are shown in Figure 6. In the formula, Γ is the contact line stress value (MPa), and x and z represent the radial force (kN) and temperature (°C) respectively.

Figure 6.

Radial loading-temperature equivalent correlation function surface: (a–d) corresponding to fixed point results and (e–h) corresponding to maximum point results.

The results in Figure 6(a) to (d) depicted the radial loading-temperature HMH stress outcomes at a fixed position along the contact line. The HMH equivalent stress values range around 400 MPa, while the values can reach 2000 MPa at the maximum point. Compared to the axial loading (the HMH equivalent stress values can reach 1000 MPa), temperature exerts a greater influence on the HMH stress values at the maximum point position. This may be because the temperature change will lead to a change in the length or volume of the material, thus causing the change in stress.³⁸ High temperatures causing the material’s internal structure to loosen and grains to grow, which reduces strength and hardness. Additionally, higher temperatures make the material more susceptible to deformation, thus enhancing its toughness and plasticity.³⁹ The red scatter points in Figure 6(e) to (h) represent the HMH equivalent stress values under specific working conditions. Due to the complexity of dynamic simulation modeling of geometric shape, motion, and lubricant characteristics of rolling elements in bearings, this study does not consider the dynamic changes of boundary conditions and the influence of lubrication methods, which is also the focus of future research.^40,41

The surface of each equivalent correlation function was smooth in Figure 6, the scatter points of training data and verification data had a good correlation with the surface of the correlation function, and the surface of the equivalent correlation function was non-uniform and gentle. It is mentioned in the thermodynamic theory that high-temperature conditions can affect the thermal expansion coefficient of the material, causing the volume change of the roller, and the stress value can change accordingly. Compared with the radial force-axial force equivalence results, the MSE value of the radial force-temperature equivalent correlation function was on the rise as a whole, and it can be inferred that the temperature has a more significant effect on the contact line stress value.

Contact line A_FP:

Γ = 1.265 \cdot x + 2.246 \cdot z - 0.01528 \cdot xz - 89.63

(16)

Contact line B_FP:

Γ = - 0.1623 \cdot x + 1.492 \cdot z + 0.02708 \cdot xz - 8.942

(17)

Contact line C_FP:

Γ = 1.193 \cdot x + 1.994 \cdot z - 55.38

(18)

Contact line D_FP:

Γ = 0.4192 \cdot x + 0.00775 \cdot z^{2} + 25.45

(19)

Contact line A_MP:

Γ = - 0.4342 \cdot x + 3.316 \cdot z + 0.05027 \cdot xz - 76.98

(20)

Contact line B_MP:

Γ = - 4.536 \times 10^{- 5} \cdot x z^{2} + 0.02401 \cdot z^{2} + 0.08346 \cdot x + 71.5

(21)

Contact line C_MP:

Γ = 5.739 \cdot x + 17.96 \cdot z - 486.0

(22)

Contact line D_MP:

Γ = 0.4746 \cdot x + 2.486 \cdot z + 0.00524 \cdot xz - 61.97

(23)

Radial loading-axial loading-temperature equivalence

In the analysis of multi-index orthogonal simulation results, intuitive analysis and variance analysis are generally used. The intuitive analysis method is simple and easy to understand, and the variance analysis method has high accuracy. Combining the advantages of the two methods, the simulation results were analyzed to distinguish the primary and secondary factors.^42,43 The intuitive analysis method is used for analysis in this study.

Among the three factors, the temperature had a great influence on the experimental results, and the primary and secondary relationship of the other two factors needs to be determined by combining range analysis.^44–46 The results of the range analysis are shown in Table 6. To facilitate analysis, the range data in the table retains 0 valid digits. At position C or D, whether it was a fixed point or a maximum point, the primary and secondary relationships of the influencing factors were temperature, axial, and radial force components. At position A and B, there was no uniform law between the primary and secondary factors affecting the fixed point and the maximum point. Although the axial force in the primary and secondary relationship of the influencing factors was before the radial force, the range level of the radial and axial types of loading in the table was approximately at the same level.

Table 6.

Range analysis results and Đ value calculation results.

Range analysis results
Fixed location	A-point			B-point			C-point			D-point
Extreme difference (MPa)	4	8	118	81	17	230	48	3	159	17	4	97
Secondary relation	c b a			c a b			c a b			c a b
Maximum location	A-point			B-point			C-point			D-point
Extreme difference (MPa)	141	68	462	4	12	271	226	41	1441	36	4	219
Secondary relation	c a b			c b a			c a b			c a b
Đ value calculation results
Location	Factor a: Radial load			Factor b: Axial load			Factor c: Temperature
A-fixed point	0.9061			0.9077			1.2777
B-fixed point	3.4320			3.4352			3.8544
C-fixed point	1.9495			1.9735			2.1231
D-fixed point	0.9643			0.9624			1.1395
A-maximum point	10.4315			10.5401			12.2671
B-maximum point	5.3052			5.3068			6.9003
C-maximum point	104.1299			105.9621			116.8856
D-maximum point	2.9678			2.9655			3.6517
F distribution results	Đ _0.1(4,12) = 2.48			Đ _0.05(4,12) = 3.26			Đ _0.01(4,12) = 5.412

Taking the factor level as the abscissa and the average deviation of the simulation result as the ordinate, the trend diagram of each working condition factor index is drawn. As shown in Figure 7(a) and (b), the trend diagram of each contact line position factor and index was drawn. The abscissa a₁–a₅, b₁–b₅, and c₁–c₅ in the diagram represent the five horizontal values of radial force, axial force, and temperature respectively. To accurately estimate the importance of the influence of each factor on the simulation results, the above data were analyzed by variance analysis. In F statistics, the value of Đ is defined as the ratio of inter-group fluctuation to intra-group fluctuation. The value obtained from F statistics was compared with the probability graph obtained from the degrees of freedom of factors and random errors, and the influence degree of each factor was obtained.

Figure 7.

Equivalent results: (a) and (b) is the factor and indicator trend chart and (c–j) is the radial loading-axial loadingtemperature equivalent correlation function surface.

The calculated Đ value is listed in Table 6, and the F distribution table is consulted for the significance test. From the data in Table 6, it can be seen that because Đ_0.01(4,12) = 5.412 was greater than the Đ value of the fixed point, the significance of the factors a, b, and c cannot be distinguished by the data. The SR algorithm is used to evaluate and screen the main influencing characteristics in combination with the random forest model. It is found that the axial force was not sensitive to the stress result. Therefore, the radial force-temperature-stress value data was used to draw the correlation surface. When the radial force-temperature is equivalent, the complexity and MSE value of the training model can be obtained by the SR algorithm. As shown in Figure 7(c) to (f), the training data and verification data were all in the radial force-axial force-temperature equivalent correlation function surface, and the fitting effect was good. The correlation functions corresponding to each position are shown in equations (20)–(23). In the formula, Γ is the contact line stress value (MPa), x is the radial force (kN), and z is the temperature (°C).

Combined with the analysis of the maximum data in Table 6, for the A-maximum point, due to Đ_0.01(4,12) = 5.412 < Đ_A = 10.432, Đ_0.01(4,12) = 5.412 < Đ_B = 10.54, Đ_0.01(4,12) = 5.412 < Đ_C = 12.27, the change of factor a, b, and c levels had a highly significant impact on the simulation results. The difference between the Đ values corresponding to the radial force and the axial force in the table was small. The results obtained by combining the random forest feature screening again show that the axial force was insensitive to the stress value. The model trained by the SR algorithm is tested. The results are shown in Figure 7(g) to (j). The training points and verification points in the figure were evenly distributed in the equivalent correlation function surface of the SR algorithm. The contact line stress value changes approximately linearly with the increase of axial load and temperature. The trained model expression is as follows: equations (24)–(31), where Γ is the contact line stress value (MPa), and x and z represent radial force (kN) and temperature (°C), respectively.

The above results show that temperature is the most important factor affecting the stress value of the contact line. The order of temperature terms in the correlation function obtained by the SR algorithm is relatively higher than that of other influencing factors, and the order of influencing factors in the correlation function corresponding to different attention parts is different. Compared with the results in Table 5, the MSE value of C_MP at the contact line position between the roller and bearing inner ring was larger than that at other points, due to the larger maximum stress value extracted at this position. For example, the maximum stress value in Figure 7(i) is 1750 MPa, which is about six times that of the latter (D_MP = 300 MPa) in Figure 7(j), so the MSE value calculated by C_MP is larger.

Contact line A_FP:

Γ = 0.005465 \cdot z^{2} - 0.002224 \cdot xz + 0.7544 \cdot z + 0.1547

(24)

Contact line B_FP:

Γ = 0.02725 \cdot xz - 0.1736 \cdot x + 1.508 \cdot z - 9.605

(25)

Contact line C_FP:

Γ = 1.205 x + 2.0 z - 56.61

(26)

Contact line D_FP:

Γ = 0.4382 \cdot x + 1.242 \cdot z - 19.06

(27)

Contact line A_MP:

Γ = - 8.189 \cdot 10^{- 5} \cdot x z^{2} + 0.05326 \cdot xz + 3.814 \cdot z - 120.8

(28)

Contact line B_MP:

Γ = 0.02593 \cdot z^{2} - 0.01335 \cdot xz + x + 48.83

(29)

Contact line C_MP:

Γ = - 0.09956 \cdot z^{2} + 34.0 \cdot z + 5.629 \cdot x - 1052.01

(30)

Contact line D_MP:

Γ = 0.8832 \cdot x + 2.745 \cdot z - 82.02

(31)

Discussion

Based on the above analysis, the SR algorithm can better process the simulation test data and provide the mathematical equivalent correlation function of each contact line. Combining the SR algorithm with the finite element analysis, this method has certain feasibility from the perspective of simulation. Through linear interpolation, training, and verification, the original limited data is extrapolated to obtain more data, which provides simulation support for the realization of a single-bearing bench design.

From the perspective of engineering, the increase of stress value at any contact line position inside the bearing can affect the reliability analysis results of the bearing. Combined with the simulation results, whether it is a fixed point or a maximum point, the stress value at the contact line between the inner ring of the bearing and the roller is generally higher than the stress value of other contact lines, which is the focus of attention. In the above results, it is pointed out that the temperature has a great influence on the stress value of the contact line position inside the bearing. The grease can reduce the friction between the roller and the raceway to a certain extent, reduce the temperature between the contact surface of the bearing, and the same effect can also be achieved by reducing the bearing speed. The changing trend of stress value after exceeding the normal temperature range remains to be studied. In the bearing stand test, the infrared thermometer is set up on the upper part of the bearing box to monitor the internal temperature of the bearing. The normal operation state monitoring data and the abnormal temperature rise monitoring data are imported into the finite element software for analysis. The symbol regression method is also used for data processing to establish the bearing health state data set and the abnormal temperature rise data set. The training mathematical correlation function can be used for data screening and bearing temperature alarm device trigger limit, which provides a basis for bearing health state monitoring and evaluation.

In summary, the results showed that the axial force has less influence on the calculation results than the radial force in the process of axial force-radial force-temperature HMH stress equivalence. By referring to the expression of the fixed points correlation function at the contact line position in the outer ring of the roller, the axial load can be equivalent to the radial load. The equivalent effect can be verified by the calibration test of the strain gauge arranged in the outer ring of the bearing. It is meaningful for designing or simplifying the scheme of the bearing stand based on the life consistency equivalent model.

Conclusion

By establishing the sub-model of CRH380 axle box bearing, the radial force-axial force-temperature equivalence is studied, and the Huber-von Mises-Hencky stress is used as a bridge to exploring the equivalent relationship between different factors. The Symbolic regression algorithm is applied for data analysis and processing. The main conclusions of this study are summarized below.

(1) The equivalent correlation function trained by the SR algorithm can well express the equivalent relationship between the equivalent parameters, and 20% of the data of the database is randomly selected to verify the training model, the effect of data randomness is not included in the study.

(2) The axial loading can be equivalent to the radial load through the equivalent correlation function, and the temperature in the three factors has a great influence on the stress value of the contact line position inside the bearing.

(3) It is feasible and generalizable to analyze and process the finite element results with an SR algorithm from the perspective of simulation.

Footnotes

Handling Editor: Sharmili Pandian

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors sincerely thank the support from the National Key Research and Development Program of China (2021YFB3400703), the National Natural Science Foundation of China (Grant No. 52375159), the National Railway Administration of China (KF2023-025), and the Independent Research Project of the State Key Laboratory of Traction Power (Grant No. 2022TPL_T03).

ORCID iD

Shuancheng Wang

References

Salunkhe

Desavale

Kumbhar

. Vibration analysis of deep groove ball bearing using finite element analysis and dimension analysis. J Lubr Technol 2022; 144: 081202.

Toumi

Murer

Bogard

, et al. Numerical simulation and experimental comparison of flaw evolution on a bearing raceway: case of thrust ball bearing. J Comput Des Eng 2018; 5: 427–434.

Torun

Swaminathan

. High-dimensional global optimization method for high-frequency electronic design. IEEE Trans Microw Theory Tech 2019; 99: 1–15.

Liu

Qiu

Chen

. Mechanical properties of coupled bending and contact to hollow cylindrical roller contacts. Tribol Trans 2022; 65: 963–976.

Reddy

Hemasunder

Ramana

, et al. State of art on FEM approach in inverse heat transfer problems for different materials. Mater Today Proc. Epub ahead of print 1 July 2023. DOI: 10.1016/j.matpr.2023.06.323.

Pascual

Ortega

Plaza

, et al. Analysis of the influence of L-PBF porosity on the mechanical behavior of AlSi10Mg by XRCT-based FEM. J Mater Res Technol 2023; 22: 958–981.

Cheng

Thornton

Luding

, et al. Concurrent multi-scale modeling of granular materials: role of coarse-graining in FEM-DEM coupling. Comput Method Appl Mech Eng 2023; 403: 115651.

Zhao

Feng

, et al. Service life prediction method for wheel-hub-bearing under random multi-axial wheel loading. Eng Fail Anal 2021; 122: 105211.

Song

Zeng

, et al. Analysis of the service life of railway hydraulic dampers considering temperature and loading. Proc IMechE, Part F: J Rail and Rapid Transit 2021; 235: 3–11.

10.

Lang

Yang

Cao

, et al. Thermal effect analysis of sliding bearing and reasonable determination of the equivalent temperature of flow field. Lubr Eng 2012; 37: 70–73+98.

11.

Liu

. Structural deterioration knowledge ontology towards physics-informed machine learning for enhanced bridge deterioration prediction. J Comput Civil Eng 2023; 37: 04022051.

12.

Ren

Cui

, et al. Optical imaging spectroscopy coupled with machine learning for detecting heavy metal of plants: a review. Front Plant Sci 2022; 13: 1007991.

13.

Nguyen

Ionescu

Géhin

, et al. Predicting the opening state of a group of windows in an open-plan office by using machine learning models. Build Environ 2022; 225: 109636.

14.

Wei

Chu

Sun

, et al. Machine learning in materials science. Info Mat 2019; 1: 338–358.

15.

Ruijie

Chunlei

, et al. Research progress on rolling contact fatigue damage of bainitic rail steel. Eng Fail Anal 2023; 143: 106875.

16.

Nguyen-Le

Tao

Nguyen

, et al. A data-driven approach based on long short-term memory and hidden Markov model for crack propagation prediction. Eng Fract Mech 2020; 235: 107085.

17.

Liao

Yang

Qin

, et al. Research on evolution law of short fatigue cracks based on genetic wavelet neural network. J Chin Rai Soc 2018; 40: 66–72.

18.

Zheng

Zhang

, et al. Performance prediction of perovskite materials based on different machine learning algorithms. Chin J Nonferrous Met 2019; 29: 803–809.

19.

Schmidt

Lipson

. Distilling free-form natural laws from experimental data. Science 2009; 324: 81–85.

20.

Chen

Liang

Huang

, et al. Symbolic discovery of optimization algorithms. Adv Neural Inf Process Syst 2024; 36: 1–30.

21.

Yang

Luo

. A Survey on deep learning for symbolic music generation: representations, algorithms, evaluations, and challenges. ACM Comput Surv 2023; 56: 1–39.

22.

Tenachi

Ibata

Diakogiannis

. Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws. arXiv preprint 2303:03192, 2023.

23.

Antoniadis

Lambert

Poggi

. Random forests for global sensitivity analysis: a selective review. Reliab Eng Syst Saf 2021; 206: 107312.

24.

Zhou

Yang

Xiao

, et al. Crack growth rate model derived from domain knowledge-guided symbolic regression. Chin J Mech Eng 2023; 36: 1–16.

25.

Harris

Akotzalas

. Advanced concepts of bearing technology: rolling bearing analysis. New York, NY: CRC Press Taylor & Francis Group, 2006.

26.

Udrescu

Tegmark

. AI Feynman: a physics-inspired method for symbolic regression. Sci Adv 2020; 6: 2631.

27.

Gan

Zhong

. Integration of symbolic regression and domain knowledge for interpretable modeling of remaining fatigue life under multistep loading. Int J Fatigue 2022; 161: 106889.

28.

Matchev

Matcheva

Roman

. Analytical modeling of exoplanet transit spectroscopy with dimensional analysis and symbolic regression. Astrophys J 2022; 930: 1–13.

29.

Weng

Song

Zhu

, et al. Simple descriptor derived from symbolic regression accelerating the discovery of new perovskite catalysts. Nat Commun 2020; 11: 1–8.

30.

Cranmer

Alvaro

. Discovering symbolic models from deep learning with inductive biases. Adv Neural Inf Process Syst 2020; 33: 17429–17442.

31.

Jiang

Wang

Chu

, et al. Analysis of internal contact characteristics of high-speed train axle box bearings considering temperature. J Chongqing Univ Technol 2023; 37: 109–118. (in Chinese)

32.

Han

Chen

, et al. Orthogonal test design for optimization of synthesis of Bi₂WO₆ superstructure with high photocatalytic activity by hydrothermal method. J Chem Technol Biotechnol 2021; 97: 943–949.

33.

Zhang

Lin

. Optimization calculation of single stage coil launcher based on prediction model. IEEE Trans Plasma Sci 2021; 49: 2176–2182.

34.

Ren

QCR

Liu

, et al. Frost-heave characteristics of saturated lean clay fillers by orthogonal test design considering the static surcharge. Case Stud Constr Mater 2022; 17: e01526.

35.

Zheng

. Experimental design and data analysis: based on R language application. Guangzhou: South China University of Technology Press, 2016, pp.98–136.

36.

Brutger

Kertzer

Renshon

, et al. Abstraction and detail in experimental design. Am J Pol Sci 2023; 67: 979–995.

37.

Nie

, et al. Numerical simulation analysis of a combined wind-fog dust removal device in return air roadways based on an orthogonal test. Powder Technol 2023; 413: 117890.

38.

Dong

Qiu

, et al. Analysis and experimental research of transient temperature rise characteristics of high-speed cylindrical roller bearing. Sci Rep 2024; 14: 711.

39.

Kim

Ryu

, et al. On heat, temperature, and cavity oil volume fraction of an under-race lubricated angular contact ball bearing. Tribol Int 2023; 187: 108715.

40.

Liang

, et al. Motion stability analysis of cage of rolling bearing under the variable-speed condition. Nonlinear Dyn 2023; 111: 11045–11063.

41.

Xie

Jiao

Yang

, et al. Experimental and numerical exploration on the nonlinear dynamic behaviors of a novel bearing lubricated by low viscosity lubricant. Mech Syst Signal Process 2023; 182: 109349.

42.

Liu

Wen

. Experimental design and data process. Chongqing: Chonqqing University Press, 2021, pp.102–130.

43.

Wang

Sui

. Experimental design and MATLAB data analysis. Beijing: Tsinghua University Press, 2012, pp.153–175.

44.

Hossain

. Spectral simulation and materials design for camouflage textiles coloration against materials of multidimensional combat backgrounds in visible and near infrared spectrums. MRS Commun 2023; 13: 306–319.

45.

Hamada

Khalil

Abou Al-Sood

, et al. An experimental investigation of nanofluid, nanocoating, and energy storage materials on the performance of parabolic trough collector. Appl Therm Eng 2023; 219: 119450.

46.

Wang

Zhang

Wang

, et al. Experimental study on the influence of temperature and humidity on the thermal conductivity of building insulation materials. Energy Built Environ 2023; 4: 386–398.

Number	Factor1	Factor 2	Factor 3	Number	Factor1	Factor 2	Factor 3
1	1	1	1	14	3	4	5
2	1	2	3	15	3	5	2
3	1	3	5	16	4	1	3
4	1	4	2	17	4	2	5
5	1	5	4	18	4	3	2
6	2	1	5	19	4	4	4
7	2	2	2	20	4	5	1
8	2	3	4	21	5	1	2
9	2	4	1	22	5	2	4
10	2	5	3	23	5	3	1
11	3	1	4	24	5	4	3
12	3	2	1	25	5	5	5
13	3	3	3

Number	Factor1	Factor 2	Factor 3	Number	Factor1	Factor 2	Factor 3
1	1	1	1	14	3	4	5
2	1	2	3	15	3	5	2
3	1	3	5	16	4	1	3
4	1	4	2	17	4	2	5
5	1	5	4	18	4	3	2
6	2	1	5	19	4	4	4
7	2	2	2	20	4	5	1
8	2	3	4	21	5	1	2
9	2	4	1	22	5	2	4
10	2	5	3	23	5	3	1
11	3	1	4	24	5	4	3
12	3	2	1	25	5	5	5
13	3	3	3

Number	Factor1	Factor 2	Factor 3	Number	Factor1	Factor 2	Factor 3
1	1	1	1	14	3	4	5
2	1	2	3	15	3	5	2
3	1	3	5	16	4	1	3
4	1	4	2	17	4	2	5
5	1	5	4	18	4	3	2
6	2	1	5	19	4	4	4
7	2	2	2	20	4	5	1
8	2	3	4	21	5	1	2
9	2	4	1	22	5	2	4
10	2	5	3	23	5	3	1
11	3	1	4	24	5	4	3
12	3	2	1	25	5	5	5
13	3	3	3