Plastic deformation-based rolling bearing reliability and sensitivity analysis under incomplete probability information

Abstract

It is currently difficult to adequately characterize the probability distribution information of rolling bearings, and research on the reliability of rolling bearings based on plastic deformation is lacking. To address this problem, the Latin hypercube design, BP neural networks, and higher-order moment method are combined, fully considering the randomness of the load, geometry, and material parameters of the rolling bearing. The dynamics model of the rolling bearing was constructed by using the Hertz contact theory and Jones model, combined with a Latin hypercube design to obtain neural network training samples. Based on the rational construction of the neural network structure, the mapping relationship between the equivalent stress and the design variables is obtained. The state function of the rolling bearing is constructed according to the stress–strength interference model, and then, the reliability index and reliability of the rolling bearing are obtained using the higher-order moment method, followed by a reliability sensitivity analysis. Finally, the proposed method is applied to the reliability analysis of a certain type of angular contact ball bearing, and the calculation results are compared with the Monte Carlo method results to demonstrate the correctness and effectiveness of the proposed method.

Introduction

Rolling bearings are one of the most widely used components in industrial applications, and their failure is one of the most common causes of mechanical failure. If an excessive load is applied to a bearing at rest, the elastic deformation between the rolling element and the raceway of the rolling bearing will be converted to plastic deformation. The indentation inside the bearing due to plastic deformation will cause vibration, noise, and a change in the friction torque during rotation, making the bearing unable to work properly and even become the cause of early fatigue failure. In a rolling bearing, the load acting between the rolling element and the raceway can only exist over a small contact area, so the stress generated on the surface of the rolling element and the raceway is usually very high even if the rolling element is subject to a moderate load, easily leading to plastic yielding. Therefore, there is an urgent need to improve the reliability of rolling bearings based on plastic deformation to prevent catastrophic mechanical failure.

At present, the basic reliability-based design conditions of rolling bearings are not described well, due to lack of effective data, complex working environment and many influencing factors. Hence, the theoretical research on the reliability-based design of rolling bearings still mainly relies on experience and is in its infancy. Nevertheless, considerable research achievements have been made in the reliability engineering of rolling bearings in recent years in terms of modeling methods, testing techniques, engineering applications, and many other aspects. Weibull ¹ was the first to speculate that the probability of fatigue strength failure at any point on the S–N curve is equal to the probability of fatigue life failure and described two methods for determining the fatigue strength probability distribution. Yudong et al. ² solved the three parameters of the Weibull distribution using the correlation coefficient optimization method and established a model of bearing life reliability based on the Weibull distribution. Yoon and Choi ³ analyzed the effect of geometric randomness on the reliability of air-bearing structures using the mean-value first-order second-moment method and carried out reliability-based design optimization. Rackwitz and Flessler^4,5 proposed an equivalent normal variable transformation algorithm for cases with nonnormal input parameters and then calculated the reliability using an advanced first-order second-moment method. Liu et al. ⁶ used the modified first-order and second-moment method to calculate the reliability and reliability sensitivity of computerized numerical control lathe spindle deformation model considering thermal effect. Nazir et al. ⁷ used three-dimensional finite element analysis and surrogate models to determine equivalent stress intensity factors for the evaluation of the propagation uncertainty of fatigue cracks and predicted the failure probability of rolling ball bearings. Li et al. ⁸ presented a method for fault diagnosis of motor rolling bearings based on neural networks and time-frequency domain bearing vibration analysis. Zhang and Gu ⁹ constructed a nonlinear dynamic model of angular contact ball bearings installed in pairs to analyze the motion error of the bearings and obtained the reliability and reliability sensitivity using a stochastic perturbation method and Edgeworth series. Li et al. ¹⁰ proposed a reliability assessment method based on the degenerate Markov model to successfully assess the reliability of wind turbine gearbox bearings based on small samples. Guo et al. ¹¹ proposed a method that enables the reliability assessment of axial cracking in individual wind turbine bearings through statistical analysis of historical data and connected reliability prediction with wind turbine design and operation. Yucesan and Viana ¹² developed a hybrid model for spindle bearing fatigue prognosis based on physics and machine learning. Chudzik and Warda ¹³ used the finite element method to determine the distribution of maximum equivalent subsurface stresses and their depths needed in fatigue life calculation and analyzed the impact of the combined load on the measured fatigue life of radial cylindrical roller bearings. Nguyen et al. ¹⁴ proposed the optimization of the crankpin bearing parameters considering effect of the high-speed dynamic load and micro asperity contact to improve the lubrication and tribology performance of the crankpin bearing. Zheng et al. ¹⁵ presented a comprehensive quasi-static model for calculating internal load and contact pressure distributions in tapered double-inner ring bearing. Hou et al. ¹⁶ focused on the variation in contact load on the most loaded position of the outer raceway of a gearbox bearing in high-speed train. Tiwari and Waghole ¹⁷ presented the single objective optimization of spherical roller bearings by considering the maximization of dynamic capacity as objective function and carried out the sensitivity analysis of bearing parameters to see the effect of manufacturing tolerance on the objective function.

However, research on the reliability-based design of rolling bearings based on plastic deformation under incomplete probability information is still very limited. In engineering practice, it is often difficult to accurately determine the probability distribution of the design parameters of rolling bearings. Reliability analysis and design under incomplete probability information inevitably produces calculation errors and leads to results that differ substantially from actual solutions. In recent years, the development of rolling bearings in multiple directions (e.g. high speed, heavy load, and high reliability) has put forward higher requirements for the reliability of rolling bearings. Therefore, it is necessary to break through the reliability bottleneck that restricts the development of rolling bearings and form a theory and method for the reliability-based design of rolling bearings based on plastic deformation.

In this study, the geometrical dimensions and material properties of rolling bearings and the effect of the combined load on rolling bearings are fully taken into account, and the Hertz contact theory, Jones model, backpropagation (BP) neural network, and reliability-based design method are combined. Considering the high nonlinearity of the static equilibrium equation of rolling bearings, the nonlinear function is fitted using the BP neural network optimized by a particle swarm algorithm. Using probability, statistics, and matrix analysis theory, the theoretical formula for the reliability-based design of rolling bearings based on plastic deformation is derived, and the reliability model and the reliability sensitivity analysis model for rolling bearings based on plastic deformation under the knowledge of only the first four moments of the basic random parameters are established. Finally, the reliability information of rolling bearings based on plastic deformation under missing probability information was used for verification. Numerical examples showed that the method proposed in this study is convenient and practical for reliability-based design and reliability sensitivity analysis of rolling bearings based on plastic deformation and reveals the degree of influence of basic random parameters have on the plastic deformation of rolling bearings, which has far-reaching implications for guiding the improvement in parameters in the design of rolling bearings.

Mechanical model of rolling bearings

The rolling bearing load is transmitted through a very small contact area, where a very large pressure is generated, easily causing plastic yielding. Therefore, it is necessary to first analyze the stress state of the contact area when verifying the contact strength. In this section, based on the Hertz contact theory and the Jones model, a static analysis model for the three-degree-of-freedom angular contact ball bearing is established by considering the elastic contact between the ball and the raceway in addition to the force of the ball and the external load on the inner ring of the bearing, and this model is used to solve the maximum contact stress. Furthermore, the equivalent stress in the case of a point contact is obtained according to the fourth strength theory.

The relative displacements of the inner ring relative to the outer ring of an angular contact ball bearing in the directions of the axial force $F_{\mathrm{a}}$ , radial force $F_{\mathrm{r}}$ , and moment $M$ , which form the combined external load, are the relative radial displacement ${\delta }_{\mathrm{r}}$ , relative axial displacement ${\delta }_{\mathrm{a}}$ , and relative angular displacement $\theta$ , respectively, as shown in Figure 1. When the load is applied to the bearing, if the outer ring is assumed to be fixed, the inner ring will be displaced, and the groove curvature center radii of the inner and outer rings will also be displaced, as shown in Figure 2. In the figure, $m_0$ and $n_0$ are the curvature centers of the inner and outer raceway grooves at the most loaded position $\psi =0$ , respectively, $m_0^{\prime }$ is the position of $m_0$ after displacement, $m_{\Psi }$ and $n_{\Psi }$ are the raceway groove curvature centers at any position $\psi$ , respectively, and $m_{\Psi }^{\prime }$ is the position of $m_{\Psi }$ after displacement. According to the Jones theory, ¹⁸ the force equilibrium equation of the rolling bearing is established to solve the maximum load of the ball.

Figure 1.

Displacement of the inner ring under the combined action of radial, axial, and moment loads (with the outer ring fixed).

Figure 2.

Relative positions of the groove curvature center radii before and after displacement.

When the bearing is in a contact state under no load, the distance between the raceway groove curvature centers at any position is always $A$ :

$$A=(f_{\mathrm{i}}+f_{\mathrm{e}}-1)D_{\mathrm{w}}$$

(1)

where $f_{\mathrm{i}}$ and $f_{\mathrm{e}}$ are the inner and outer raceway groove curvature coefficients, respectively, and $D_{\mathrm{w}}$ is the ball diameter.

After loading, any position $\psi$ is the angle between the maximum loaded ball and any other ball. $0\le \psi \le \pi$ because of the symmetry. At position $\psi$ , the total contact deformation ${\delta }_{\Psi }$ between the ball and the inner and outer rings is

$${\delta }_{\Psi }=S_{\Psi }-A$$

(2)

where the distance $S_{\Psi }$ between the raceway groove curvature centers at any position after displacement is

$$S_{\Psi }=[(A\sin {\alpha }^0+{\delta }_{\mathrm{a}}+R_{\mathrm{i}}\theta \cos \psi {)}^2+(A\sin {\alpha }^0+{\delta }_{\mathrm{r}}\cos \psi {)}^2{]}^{\frac{1}{2}}$$

(3)

For the free contact angle ${\alpha }^0$ and radial clearance $p_{\mathrm{d}}$ of the ball bearing, $\cos {\alpha }^0=1-\frac{p_{\mathrm{d}}}{2A}$ and $P_{\mathrm{d}}=d_{\mathrm{e}}+d_{\mathrm{i}}-2D_{\mathrm{w}}$ , where $d_{\mathrm{i}}$ and $d_{\mathrm{e}}$ are the bearing inside and outside diameters, respectively. Substitution of equation (3) into equation (2) gives

$$\begin{array}{c}{\delta }_{\Psi }=A\hfill \\ \{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}-1\}\hfill \end{array}$$

(4)

where the loci of the centers of the inner and outer raceway groove curvature radii $R_{\mathrm{i}}$ and $R_{\mathrm{e}}$ , respectively, are

$$R_{\mathrm{i}}=\frac{d_{\mathrm{m}}}{2}+(r_{\mathrm{i}}-\frac{D_{\mathrm{w}}}{2})\cos {\alpha }^0$$

(5)

$$R_{\mathrm{e}}=\frac{d_{\mathrm{m}}}{2}-(r_{\mathrm{e}}-\frac{D_{\mathrm{w}}}{2})\cos {\alpha }^0$$

(6)

where $d_{\mathrm{m}}$ is the bearing pitch diameter; $r_{\mathrm{i}}$ and $r_{\mathrm{e}}$ are the inner and outer raceway groove curvature radii, respectively. Furthermore, a dimensionless displacement is introduced as follows.

$$\begin{array}{c}\overline{{\delta }_{\mathrm{a}}}=\frac{{\delta }_{\mathrm{a}}}{A}\hfill \\ \overline{{\delta }_{\mathrm{r}}}=\frac{{\delta }_{\mathrm{r}}}{A}\hfill \\ \overline{\theta }=\frac{\theta }{A}\hfill \end{array}$$

(7)

At any position $\psi$ , the load on the contact between the rolling element and the raceway is

$$Q_{\psi }=K_{\mathrm{n}}{\delta }_{\Psi }^n$$

(8)

where n = 3/2 for ball bearings and n = 10/9 for roller bearings. $K_{\mathrm{n}}$ is load-deflection factor, and its specific expression can be found elsewhere. ¹⁹ Substituting equation (4) into equation (8) gives the rolling element load at any position

$$\begin{array}{c}{Q}_{\Psi }=K_{\mathrm{n}}{A}^{\frac{3}{2}}\hfill \\ {\{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}-1\}}^{\frac{3}{2}}\hfill \end{array}$$

(9)

After loading, the contact angle changes, and the operating contact angle $\alpha$ at any position can be expressed as

$$\sin \alpha =\frac{\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi }{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}}$$

(10)

or

$$\cos \alpha =\frac{\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi }{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}}$$

(11)

Under the combined load, the static equilibrium equation of the rolling bearing is

$$F_{\mathrm{a}}=\sum _{\psi =0}^{\psi =\pm \pi }{Q}_{\psi }\sin \alpha$$

(12)

$$F_{\mathrm{r}}=\sum _{\mathrm{\psi }=0}^{\psi =\pm \pi }{Q}_{\psi }\cos \psi \sin \alpha$$

(13)

$$M=\frac{d_{\mathrm{m}}}{2}\sum _{\psi =0}^{\psi =\pm \pi }{Q}_{\psi }\cos \psi \sin \alpha$$

(14)

Combining equations (9), (10), and (11) with equations (12), (13), and (14), the force and moment equilibrium equations of the rolling bearing are obtained as follows:

$$F_{\mathrm{a}}=K_{\mathrm{n}}{A}^{\frac{3}{2}}\sum _{\psi =0}^{\psi =\pm \pi }\frac{{\{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}-1\}}^{\frac{3}{2}}(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}}$$

(15)

$$F_{\mathrm{r}}=K_{\mathrm{n}}{A}^{\frac{3}{2}}\sum _{\psi =0}^{\psi =\pm \pi }\frac{{\{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}-1\}}^{\frac{3}{2}}(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )\cos \psi }{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}}$$

(16)

$$M=\frac{1}{2}{d}_{\mathrm{m}}{K}_{\mathrm{n}}{A}^{\frac{3}{2}}\sum _{\psi =0}^{\psi =\pm \pi }\frac{{\{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}-1\}}^{\frac{3}{2}}(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )\cos \psi }{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta }\cos \psi )}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}}\cos \psi )}^2]}^{\frac{1}{2}}}$$

(17)

Equations (15)–(17) are a system of simultaneous nonlinear equations with unknowns ${\delta }_{\mathrm{a}}$ , ${\delta }_{\mathrm{r}}$ , and $\theta$ and are calculated using numerical methods such as the Newton–Raphson method. Once the values of ${\delta }_{\mathrm{a}}$ , ${\delta }_{\mathrm{r}}$ , and $\theta$ are obtained, the maximum load of the ball can be obtained using equation (9) at $\psi =0$ :

$$Q=K_{\mathrm{n}}{A}^{\frac{3}{2}}{\{{[{(\sin {\alpha }^0+\overline{{\delta }_{\mathrm{a}}}+R_{\mathrm{i}}\overline{\theta })}^2+{(\cos {\alpha }^0+\overline{{\delta }_{\mathrm{r}}})}^2]}^{\frac{1}{2}}-1\}}^{\frac{3}{2}}$$

(18)

According to Hertz contact theory in the analysis of the spatial elastic point contact problem of rolling bearings, the maximum contact stress at the geometric center of the elliptical contact area is ²⁰

$${\sigma }_0=\frac{3Q}{2\pi \mathrm{ab}}$$

(19)

where the major and minor semi-axes of the elliptical contact area, $a$ and $b$ , are

$$a={\left[\frac{3Q\Pi \left(e\right)}{\pi (1-e^2)\tilde{E}\sum \rho }\right]}^{\frac{1}{3}}$$

(20)

$$b={\left[\frac{3Q\sqrt{1-e^2}\Pi \left(e\right)}{\pi \tilde{E}\sum \rho }\right]}^{\frac{1}{3}}$$

(21)

where $\sum \rho$ is the bearing curvature sum, $\tilde{E}$ is the contact modulus, $\tilde{E}={(\frac{1-{{\xi }_1}^2}{E_1}+\frac{1-{{\xi }_2}^2}{E_2})}^{-1}$ , with ${\xi }_1$ and ${\xi }_2$ being the Poisson’s ratios of the raceway and ball, respectively, and $E_1$ and $E_2$ being the elastic moduli of the raceway and ball, respectively.

Brewe and Hamrock ²¹ gave a set of approximate calculation equations by curve fitting, and the errors between the calculated results and the exact values are no more than 3%.

$$e^2=1-1/k^2$$

(22)

$$\Pi \left(e\right)\approx 1.0003+\frac{0.5968}{R_y/R_x}$$

(23)

Subscript $x$ and $y$ represent the directions of major and minor semi-axes of the contact ellipse, respectively. The equivalent radius in each direction is defined as follows: $R_x=\frac{D_{\mathrm{w}}}{2}(1-\gamma )$ and $R_y=\frac{D_{\mathrm{w}}{f}_{\mathrm{i}}}{2f_{\mathrm{i}}-1}$ for the inner raceway contact, and $R_x=\frac{D_{\mathrm{w}}}{2}(1+\gamma )$ and $R_y=\frac{D_{\mathrm{w}}{f}_{\mathrm{e}}}{2f_{\mathrm{e}}-1}$ for the outer raceway contact. In addition, $\gamma =\frac{D_{\mathrm{w}}}{d_{\mathrm{m}}}\cos \alpha$ .

To analyze the yield deformation of rolling bearings and obtain the criterion for the contact failure of the rolling element and the raceway, Harris and Kotzalas ²² gave the estimated stress state curve results. To facilitate the application in practical engineering, the stress state at the center point of the contact area is a triaxial compressive stress state under point contact conditions, which is denoted by engineering principal stress symbols ( ${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3$ ), where the negative sign indicates that the principal stress is a compressive stress. The triaxial static compressive stress can be expressed as ²³

$${\sigma }_1=-{\sigma }_0[0.505+0.255{(b/a)}^{0.6609}]$$

(24)

$${\sigma }_2=-{\sigma }_0[1.01-0.250{(b/a)}^{0.5797}]$$

(25)

$${\sigma }_3=-{\sigma }_0$$

(26)

For a triaxial static compressive stress state, the fourth strength theory needs to be used to determine whether the rolling bearing material undergoes plastic deformation, and the von Mises equivalent stress is

$${\sigma }_{\mathrm{VM}}=\sqrt{\frac{1}{2}[{({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_2-{\sigma }_3)}^2+{({\sigma }_3-{\sigma }_1)}^2]}$$

(27)

Substitution of equations (24)–(26) into equation (27) gives

$${\sigma }_{\mathrm{VM}}={\sigma }_0\sqrt{0.25-0.379{(b/a)}^{0.62}+0.192{(b/a)}^{1.24}}$$

(28)

State function of the rolling bearing based on plastic deformation

Due to the complexity of the static equilibrium equation of rolling bearings, the relationship between the basic random parameters (e.g. the combined load, geometry, and material properties) of the rolling bearing and the equivalent stress is an implicit function. Therefore, the relationship between the two needs to be fitted by a BP neural network and the optimal initial weights and thresholds of the network are obtained using the particle swarm optimization (PSO) algorithm. Then, the stress–strength interference model is used to establish the rolling bearing state function which prepare for the model of the reliability-based design. The flow chart is shown in Figure 3.

• Step 1: Dataset construction. Parameters such as the geometric dimensions, the applied load, and the material strength of the rolling bearing are highly random and regarded as basic random parameters. The Latin hypercube design works by evenly equating the sample intervals, generating random numbers in each small interval, and then disrupting the order of the random numbers to obtain the necessary samples, hence avoiding the problem of sample over-aggregation. Here, we use the Latin hypercube design for sampling according to the probability distribution law of these basic random parameters to obtain the input data of the BP neural network. The sampling results are grouped and brought into the static equilibrium equation of rolling bearing to get the corresponding n maximum contact loads ${\sigma }_0$ . According to Hertz contact theory, for the elliptical contact region formed by the contact between rolling body and raceway the fourth strength theory is used to get the equivalent stress ${\sigma }_{\mathrm{VM}}$ . In this way the basic random parameters and equivalent stress ${\sigma }_{\mathrm{VM}}$ constitute n sets of data sets, and the data sets are partitioned into training set, test set and validation set.

• Step 2: Data normalization. To avoid the influence of the large difference in the order of magnitude between the input and output data, the normalization function is used tonormalize the data by mapping the row minimum and maximum values of the data matrix to [−1, 1] so that the indicators are of the same order of magnitude and, hence, suitable for comprehensive comparison and evaluation. ²⁴ The convergence speed and accuracy of the model can also be further improved by normalization.

Figure 3.

Flow chart for the reliability analysis model of rolling bearings based on plastic deformation.

The random variable data of the training set is normalized by equation (29) as follows:

$$y_k=\frac{(x_k-x_{min})}{(x_{max}-x_{min})}$$

(29)

where $x_k$ is the original data and $y_k$ is the corresponding normalized data. $x_{min}$ and $x_{max}$ are the smallest and largest number in the data series, respectively. Also, the maximum and minimum values of $y_k$ are recorded as $y_{\max }$ and $y_{min}$ , respectively.

The random variable data of the test set are normalized and denormalized by equation (30).

$$y_k=\frac{(y_{max}-y_{min})·(x_k-x_{min})}{(x_{max}-x_{min})}+y_{min}$$

(30)

• Step 3: Optimization by PSO. The main feature of BP neural network which is a multilayer feed-forward neural network is that the signal is transmitted forward and the error is propagated backward, allowing it to construct any form of complex mapping relationship between input and output. ²⁵ In order to solve the defects of local optimum and premature convergence in BP neural networks, PSO is introduced in the penalty function to optimize the BP neural networks. Each particle represents the weight and threshold of the neural network, and the optimal initial weight and threshold of the network are found through particle optimization so that the optimized BP neural network can better predict the output of the function. ²⁶

In this process, the particle positions and velocities are initialized with random values. During each iteration, the particle updates its velocity and position by individual and global extremes, with the following update equations

$$V_{\mathrm{id}}^{k+1}=\omega V_{\mathrm{id}}^k+c_1{r}_1(P_{\mathrm{id}}^k-X_{\mathrm{id}}^k)+c_2{r}_2(P_{\mathrm{gd}}^k-X_{\mathrm{id}}^k)$$

(31)

$$X_{\mathrm{id}}^{k+1}=X_{\mathrm{id}}^k+V_{\mathrm{id}}^k$$

(32)

where $\omega$ is the inertia weight; $d=1,2,\cdots ,D$ ; $i=1,2,\cdots ,n$ ; $k$ is the number of current iterations; $V_{\mathrm{id}}$ is the velocity of the particle; $c_1$ and $c_2$ are non-negative constants that become acceleration factors; $r_1$ and $r_2$ are random numbers distributed between $[0,1]$ . To prevent the blind search of the particle, it is generally recommended to restrict its position and velocity to a certain interval of $[-{\mathit{X}}_{max},{\mathit{X}}_{max}]$ , $[-{\mathit{V}}_{max},{\mathit{V}}_{max}]$ .

The BP neural network is trained with the training data according to the fitness function, and the training data errors are used as individuals’ fitness values F.

$$F=t\left(\sum _{i=1}^n\mathrm{abs}(y_i-o_i)\right)$$

(33)

where each particle constitutes a node, that is, $n$ is the number of network output nodes in the neural network, $y_i$ is the expected output of the i-th node of the BP neural network, $o_i$ is the predicted output of the i-th node, and $t$ is a coefficient.

• Step 4: Construction of the BP neural network. The 3-layer structure BP neural network is selected, which is a one-way propagation forward network consisting of input layer, hidden layer, and output layer, and the network structure is connected as shown in Figure 4. The number of nodes in the input layer is p, the number of nodes in the hidden layer is q, the number of nodes in the output layer is r. The internal layers are passed by neurons to derive their functional relationships as follows.

$$h_k=f_1({\alpha }_k+\sum _{i=1}^p{w_i}_k{x}_i),(i=1,2,\cdots ,p;k=i=1,2,\cdots ,q)$$

(34)

$$y_j=f_2({\beta }_j+\sum _{k=1}^q{v}_{\mathrm{kj}}{h}_k)$$

(35)

where, $h_k$ is the output of the hidden layer, $y_j$ is the output of the output layer, $x_i$ is the input variable, ${w_i}_k$ and $v_{\mathrm{kj}}$ are the weights between the 3 layers respectively, ${\alpha }_k$ and ${\beta }_j$ are the threshold values of the hidden and output layers respectively, $f_1(•)$ is the S-type nonlinear function tansig, and $f_2(•)$ is the linear function purelin.

Figure 4.

Topology of the BP neural network.

The mapping relationship of the neural network shown in equations (34) and (35) is used to fit the input and output of the data within the error allowance, then the fitting function relationship is

$$\widehat{{\sigma }_{\mathrm{VM}}}(X)=f_2[{\beta }_j+\sum _k^q{v}_{\mathrm{kj}}{f}_1({\alpha }_k+\sum _{i=1}^p{w_i}_k{x}_i)]$$

(36)

The predicted results and actual results are compared in terms of the MSE, which is used as the training error of the BP neural network to evaluate the degree of data variation.

$$\mathrm{MSE}=\frac{1}{n}\sum _1^n({{\sigma }_{\mathrm{VM}}}_i-\hat{{{\sigma }_{\mathrm{VM}}}_i})$$

(37)

where ${\sigma }_{\mathrm{VM}}$ is the actual equivalent stress during bearing operation, $\widehat{{\sigma }_{\mathrm{VM}}}$ is the equivalent stress predicted by the BP neural network.

• Step 5: Establishment of a stress–strength interference model. In general, the equivalent stress of a rolling bearing under load should be less than or equal to the allowable yield stress. If the equivalent stress of the rolling bearing is greater than the allowable yield stress, plastic deformation will start to occur in the rolling element and the raceway on the contact surface, thereby shortening the service life of the rolling bearing and, in severe cases, even leading to breakdown and failure. The state function under the stress–strength interference model can be expressed as

$$g(\mathit{X})={\sigma }_{\mathrm{s}}-\widehat{{\sigma }_{\mathrm{VM}}}(\mathit{X})$$

(38)

where ${\sigma }_{\mathrm{s}}$ is the allowable yield stress of the bearing material.

Bringing equation (36) into equation (38), we have

$$g(\mathit{X})={\sigma }_{\mathrm{s}}-f_2[{\beta }_j+\sum _{k=1}^{10}{v}_{\mathrm{kj}}{f}_1({\alpha }_k+\sum _{i=1}^6{w_i}_k{x}_i)]$$

(39)

Reliability design and reliability sensitivity analysis

In this section, the matrix equations of reliability and reliability sensitivity are derived, and the reliability-based design model and reliability sensitivity analysis model of rolling bearings based on plastic deformation are established. The matrix forms of the models can be used to calculate the reliability and reliability sensitivity through computer simulation, determine the variation trend of the reliability with a change in design parameters, and reveal the influence of changes in design parameters on the reliability in a convenient and rapid manner.

Reliability design

The main goal of reliability design of rolling bearings is to determine the bearings’ reliability. By denoting the reliability as $R$ , we have

$$R={\int }_{g\left(\mathit{X}\right)>0}{f}_{\mathit{X}}\left(\mathit{X}\right)\mathrm{d}\mathit{X}$$

(40)

where $f_{\mathit{X}}\left(\mathit{X}\right)$ is the joint probability density of the basic random parameter vector $\mathit{X}={(X_1{X}_2\cdots X_n)}^{\mathrm{T}}$ of the rolling bearing. According to the stress–strength interference model of structural reliability, solving the reliability requires the establishment of a state function $g\left(\mathit{X}\right)$ for reliability analysis, which is used to represent the two states of the reliability of the rolling bearing during operation, that is,

$$\begin{array}{c}g(\mathit{X})\le 0\mathrm{failure}\mathrm{state}\hfill \\ \multicolumn{1}{c}{g(\mathit{X})>0\mathrm{reliable}\mathrm{state}}\end{array}\}$$

(41)

where the state function equation $g(\mathit{X})=0$ is an n-dimensional limit state surface.

By applying the perturbation technique and mechanical reliability theory, the mean ${\mu }_g$ , standard deviation ${\sigma }_g$ , third moment ${\theta }_g$ , and fourth moment ${\eta }_g$ of the state function for reliability design are obtained by mathematical derivation as follows ²⁷ :

$${\mu }_g=\mathrm{E}\left[g\left(\mathit{X}\right)\right]=g\left({\mathbf{\mu }}_{\mathit{X}}\right)=g\left(\overline{\mathit{X}}\right)$$

(42)

$${\sigma }_g^2=\mathrm{Var}\left[g\left(\mathit{X}\right)\right]=\mathrm{E}\left[{(g\left(\mathit{X}\right)-{\mu }_g)}^2\right]=\frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial {\mathit{X}}^{\mathrm{T}}}{\mathbf{C}}_2\left(\mathit{X}\right)\frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial \mathit{X}}$$

(43)

$${\theta }_g=\mathrm{E}\left[{(g\left(\mathit{X}\right)-{\mu }_g)}^3\right]=\frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial {\mathit{X}}^{\mathrm{T}}}{\mathbf{C}}_3\left(\mathit{X}\right)\frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial \mathit{X}}\otimes \frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial \mathit{X}}$$

(44)

$$\begin{array}{c}\multicolumn{1}{c}{{\eta }_g=\mathrm{E}\left[{(g\left(\mathit{X}\right)-{\mu }_g)}^4\right]=\frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial {\mathit{X}}^{\mathrm{T}}}\otimes \frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial {\mathit{X}}^{\mathrm{T}}}{\mathbf{C}}_4\left(\mathit{X}\right)}\\ \frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial \mathit{X}}\otimes \frac{\partial g\left(\overline{\mathit{X}}\right)}{\partial \mathit{X}}\hfill \end{array}$$

(45)

where ${\mathbf{\mu }}_{\mathit{X}}$ , ${\mathbf{C}}_2\left(\mathit{X}\right)$ , ${\mathbf{C}}_3\left(\mathit{X}\right)$ , and ${\mathbf{C}}_4\left(\mathit{X}\right)$ represent the mean matrix, variance and covariance matrix, third-moment matrix, and fourth-moment matrix of the basic random parameters X , respectively. For ease of derivation and calculation, the diagonal terms of the variance and covariance matrix ${\mathbf{C}}_2\left(\mathit{X}\right)$ are taken to form the standard deviation vector ${\mathit{\sigma }}_{\mathit{X}}$ , regardless of whether the basic random parameters X is independent variable.

When the probability distribution of the basic random parameters X cannot be determined, if the first four moments of X (i.e. ${\mathbf{\mu }}_{\mathit{X}}$ , ${\mathbf{C}}_2\left(\mathit{X}\right)$ , ${\mathbf{C}}_3\left(\mathit{X}\right)$ , and ${\mathbf{C}}_4\left(\mathit{X}\right)$ ) are known, the reliability index ${\mathbf{C}}_2^{\mu }\left(\mathit{X}\right)=\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\left[\mathit{\rho }\right]\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)$ can be defined according to the higher-order moment method as follows ²⁸ :

$${\beta }_{\mathrm{FM}}=\frac{3(3{\alpha }_{4g}+1){\beta }_{\mathrm{SM}}+5{\alpha }_{3g}({\beta }_{\mathrm{SM}}^2-1)}{\sqrt{9{(3{\alpha }_{4g}+1)}^2-5{\alpha }_{3g}^2(13{\alpha }_{4g}+11)}}$$

(46)

where ${\beta }_{\mathrm{FM}}=\frac{{\mu }_g}{{\sigma }_g}$ is the reliability index when the basic random parameters X follows a normal distribution; ${\alpha }_{3g}=\frac{{\theta }_g}{{\sigma }_g^3}$ and ${\alpha }_{4g}=\frac{{\eta }_g}{{\sigma }_g^4}$ are the skewness and kurtosis coefficients of the state function $g\left(\mathit{X}\right)$ , respectively. An alternative expression for the reliability index can be derived by the higher-order moment method as follows:

$${\beta }_{\mathrm{FM}}=\frac{3(3{\eta }_g+{\sigma }_g^4){\mu }_g+5{\theta }_g({\mu }_g^2-{\sigma }_g^2)}{\sqrt{9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)}}$$

(47)

After obtaining the reliability index ${\beta }_{\mathrm{FM}}$ of the rolling bearing, using the higher-order moment method, an approximate estimate of the reliability $R_{\mathrm{FM}}$ can be determined as

$$R_{\mathrm{FM}}=\mathrm{\Phi }\left({\beta }_{\mathrm{FM}}\right)$$

(48)

where $\mathrm{\Phi }(·)$ represents the standard normal distribution function. Therefore, the reliability of the rolling bearing can be obtained for the corresponding reliability design.

Reliability sensitivity analysis

Reliability sensitivity reveals the degree of influence of design parameters on rolling bearings, and has far-reaching significance for guiding parameter improvement in the design of rolling bearings.

Using the reliability index equation (47) and the reliability equation (48) obtained by the higher-order moment method for reliability design, the sensitivities of the reliability $R_{\mathrm{FM}}$ to the mean vector ${\mathbf{\mu }}_{\mathit{X}}$ and standard deviation vector ${\mathit{\sigma }}_{\mathit{X}}$ of the basic random parameters X are derived as follows

$$\begin{array}{c}\frac{\partial R_{\mathrm{FM}}}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}=\frac{\partial R_{\mathrm{FM}}}{\partial {\beta }_{\mathrm{FM}}}\hfill \\ [\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\mu }_g}\frac{\partial {\mu }_g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}+\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\sigma }_g}\frac{\partial {\sigma }_g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}+\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\theta }_g}\frac{\partial {\theta }_g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}+\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\eta }_g}\frac{\partial {\eta }_g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}]\hfill \end{array}$$

(49)

$$\begin{array}{c}\frac{\partial R_{\mathrm{FM}}}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}=\frac{\partial R_{\mathrm{FM}}}{\partial {\beta }_{\mathrm{FM}}}\hfill \\ [\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\mu }_g}\frac{\partial {\mu }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}+\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\sigma }_g}\frac{\partial {\sigma }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}+\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\theta }_g}\frac{\partial {\theta }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}+\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\eta }_g}\frac{\partial {\eta }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}]\hfill \end{array}$$

(50)

where

$$\frac{\partial R_{\mathrm{FM}}}{\partial {\beta }_{\mathrm{FM}}}=\varphi \left({\beta }_{\mathrm{FM}}\right)$$

(51)

Here, φ(.) is the standard normal probability density function.

$$\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\mu }_g}=\frac{3(3{\eta }_g+{\sigma }_g^4)+10{\mu }_g{\theta }_g}{\sqrt{9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)}}$$

(52)

$$\begin{array}{c}\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\sigma }_g}=\frac{(12{\mu }_g{\sigma }_g^3-10{\theta }_g{\sigma }_g)[9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)]}{{\left(\sqrt{9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)}\right)}^3}\hfill \\ -\frac{[3{\mu }_g(3{\eta }_g+{\sigma }_g^4)+5{\theta }_g({\mu }_g^2-{\sigma }_g^2)][9{\sigma }_g(3{\eta }_g+{\sigma }_g^4)(3{\eta }_g+5{\sigma }_g^4)-110{\theta }_g^2{\sigma }_g^3]}{{\left(\sqrt{9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)}\right)}^3}\hfill \end{array}$$

(53)

$$\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\theta }_g}=\frac{5({\mu }_g^2-{\sigma }_g^2)[9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)]+5{\theta }_g(13{\eta }_g+11{\sigma }_g^4)[3(3{\eta }_g+{\sigma }_g^4){\mu }_g+5{\theta }_g({\mu }_g^2-{\sigma }_g^2)]}{{\left(\sqrt{9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)}\right)}^3}$$

(54)

$$\frac{\partial {\beta }_{\mathrm{FM}}}{\partial {\eta }_g}=\frac{18{\mu }_g[9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)]-[3(3{\eta }_g+{\sigma }_g^4){\mu }_g+5{\theta }_g({\mu }_g^2-{\sigma }_g^2)][54{\sigma }_g^2(3{\eta }_g+{\sigma }_g^4)-65{\theta }_g^2]}{{\left[2\sqrt{9{(3{\eta }_g+{\sigma }_g^4)}^2{\sigma }_g^2-5{\theta }_g^2(13{\eta }_g+11{\sigma }_g^4)}\right]}^3}$$

(55)

$$\frac{\partial {\mu }_g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}=[\frac{\partial g\left(\mathit{X}\right)}{\partial {\mu }_{X_1}}\frac{\partial g\left(\mathit{X}\right)}{\partial {\mu }_{X_2}}\cdots \frac{\partial g\left(\mathit{X}\right)}{\partial {\mu }_{X_n}}]$$

(56)

$$\frac{\partial {\sigma }_g^2}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}=2{\left({\mathbf{C}}_2\left(\mathit{X}\right)\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)}^{\mathrm{T}}\frac{{\partial }^{2}g}{\partial {\mathbf{\mu }}_{\mathit{X}}\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}=2\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}{\mathbf{C}}_2\left(\mathit{X}\right)\frac{{\partial }^{2}g}{\partial {\mathbf{\mu }}_{\mathit{X}}\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}$$

(57)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial {\theta }_g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}=\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}{\mathbf{C}}_3\left(\mathit{X}\right)\frac{\partial }{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\otimes \frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}})}\\ +\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}\otimes \frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}{\mathbf{C}}_3{\left(\mathit{X}\right)}^{\mathrm{T}}\frac{\partial }{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\hfill \end{array}$$

(58)

$$\frac{\partial {\eta }_g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}=2(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}\otimes \frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}){\mathbf{C}}_4\left(X\right)\frac{\partial }{\partial {\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}}(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\otimes \frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}})$$

(59)

$$\frac{\partial {\mu }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}=0$$

(60)

$$\frac{\partial {\sigma }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}=\left(\frac{1}{2{\sigma }_g}\right)\frac{\partial {\sigma }_g^2}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}=\frac{\partial {\sigma }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}=\left(\frac{1}{{\sigma }_g}\right){\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}{\mathbf{C}}_2^{\mu }\left(\mathit{X}\right)$$

(61)

Here, ${\mathbf{C}}_2^{\mu }\left(\mathit{X}\right)=\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\left[\mathit{\rho }\right]\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)$ , in which the correlation coefficient matrix $\left[\mathit{\rho }\right]$ and the diagonal matrix $\mathrm{diag}(·)$ are respectively expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\left[\mathit{\rho }\right]=\left[\begin{array}{cccc}1\hfill & {\rho }_{X_1{X}_2}\hfill & \cdots \hfill & {\rho }_{X_1{X}_n}\hfill \\ \multicolumn{1}{c}{{\rho }_{X_2{X}_1}}& 1\hfill & \cdots \hfill & {\rho }_{X_2{X}_n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\rho }_{X_n{X}_1}}& {\rho }_{X_n{X}_2}\hfill & \cdots \hfill & 1\hfill \end{array}\right],\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)}\\ =\left[\begin{array}{cccc}\frac{\partial g}{\partial {\mu }_{X_1}}\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{\partial g}{\partial {\mu }_{X_2}}\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & \frac{\partial g}{\partial {\mu }_{X_n}}\hfill \end{array}\right]\hfill \end{array}$$

(62)

$$\frac{\partial {\theta }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}={\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}{\mathbf{C}}_3^{\mu }\left(\mathit{X}\right)\frac{\partial }{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}({\mathit{\sigma }}_{\mathit{X}}\otimes {\mathit{\sigma }}_{\mathit{X}})+{({\mathit{\sigma }}_{\mathit{X}}\otimes {\mathit{\sigma }}_{\mathit{X}})}^{\mathrm{T}}{\mathbf{C}}_3^{\mu }{\left(\mathit{X}\right)}^{\mathrm{T}}$$

(63)

Here, ${\mathbf{C}}_3^{\mu }\left(\mathit{X}\right)=\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\left[{\mathit{\rho }}_3\right]\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\otimes \mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)$ , in which $\left[{\mathit{\rho }}_3\right]$ is the third-order correlation coefficient matrix.

$$\frac{\partial {\eta }_g}{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}=2({\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}\otimes {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}){\mathbf{C}}_4^{\mu }\left(X\right)\frac{\partial }{\partial {\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}}({\mathit{\sigma }}_{\mathit{X}}\otimes {\mathit{\sigma }}_{\mathit{X}})$$

(64)

Here, ${\mathbf{C}}_4^{\mu }\left(\mathit{X}\right)=\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\otimes \mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\left[{\mathit{\rho }}_4\right]\mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)\otimes \mathrm{diag}\left(\frac{\partial g}{\partial {\mathbf{\mu }}_{\mathit{X}}}\right)$ , in which $\left[{\mathit{\rho }}_4\right]$ is the fourth-order correlation coefficient matrix.

The values of mean and standard deviation reliability sensitivities $\mathrm{D}{R}_{\mathrm{FM}}/\mathrm{D}{\mathbf{\mu }}_{\mathit{X}}^{\mathrm{T}}$ and $\mathrm{D}{R}_{\mathrm{FM}}/\mathrm{D}{\mathit{\sigma }}_{\mathit{X}}^{\mathrm{T}}$ for the rolling bearing based on plastic deformation are obtained by substituting the known conditions and relevant data into the reliability sensitivity equations (49) and (50), respectively.

Reliability sensitivity is used to evaluate the magnitude of the influence of design parameters on the reliability of rolling bearings. To unify the description of the degree of influence of each parameter on reliability, the reliability to the mean and standard deviation of the basic random parameters are considered comprehensively and expressed in the form of reliability sensitivity gradient as

$$\mathrm{grad}\left[R_{\mathrm{FM}}({\mu }_{X_i},{\sigma }_{X_i})\right]=\sqrt{{\left(\frac{\partial R_{\mathrm{FM}}}{\partial {\mu }_{X_i}}\right)}^2+{\left(\frac{\partial R_{\mathrm{FM}}}{\partial {\sigma }_{X_i}}\right)}^2}$$

(65)

where $\mathrm{grad}[·]$ denotes the rate of change in reliability.

The nondimensionalization processing of the reliability sensitivity avoids the lack of comparability between reliability sensitivities due to the inconsistency of the units of the random parameters. The reliability sensitivity to the mean and standard deviation in basic random parameters are nondimensionalized as ²⁹

$${\tau }_i=\frac{\partial R_{\mathrm{FM}}}{\partial {\mu }_{X_i}}\frac{{\sigma }_{X_i}}{R_{\mathrm{FM}}}$$

(66)

$${\eta }_i=\frac{\partial R_{\mathrm{FM}}}{\partial {\sigma }_{X_i}}\frac{{\mu }_{X_i}}{R_{\mathrm{FM}}}$$

(67)

The dimensionless reliability sensitivity gradient $s_i$ is expressed in the form of

$$s_i=\sqrt{{{\tau }_i}^2+{{\eta }_i}^2}$$

(68)

$s_i$ is normalized to obtain the sensitivity factor ${\lambda }_i$ as follows

$${\lambda }_i=\frac{s_i}{\sum _{k=1}^n{s}_k}\times 100\%$$

(69)

Numerical examples

Among the many types of bearings, angular contact ball bearings can withstand radial and axial forces simultaneously and have good stability and good lubrication characteristics. Therefore, they are widely used in the rotating system. In this case study, a certain model of an angular contact rolling bearing is selected, with the following known parameters: number of rolling elements $\mathrm{Z}=16$ , free contact angle $\alpha =40°$ , Poisson’s ratio and elastic modulus of the rolling element and raceway ${\xi }_1={\xi }_2=0.3$ and $E_1=E_2=208\mathrm{GPa}$ . In the basic random parameters $\mathit{X}={\left(d_{\mathrm{m}}{D}_{\mathrm{w}}{r}_{\mathrm{i}}{r}_0{F}_{\mathrm{a}}{F}_{\mathrm{r}}{\sigma }_{\mathrm{s}}\right)}^{\mathrm{T}}$ , the allowable yield stress ${\sigma }_{\mathrm{s}}$ is the only random parameter whose probability distribution is difficult to determine, but its first four moments are known: $({\mu }_{{\sigma }_{\mathrm{s}}},{\sigma }_{{\sigma }_{\mathrm{s}}},{\theta }_{{\sigma }_{\mathrm{s}}},{\eta }_{{\sigma }_{\mathrm{s}}})$  = (1100 MPa, 55 MPa, 0.0298 (MPa), ³ 0.3522 (MPa) ⁴ ). Other basic random parameters all follow independent normal distributions, with their means and standard deviations shown in Table 1.

Table 1.

Basic random parameters of rolling bearings and their statistical characteristics.

Basic random parameter	Mean value	Standard deviation
Bearing pitch diameter $d_{\mathrm{m}}/\mathrm{mm}$	125.25	0.6250
Ball diameter $D_{\mathrm{w}}/\mathrm{mm}$	22.225	0.1111
Inner raceway groove curvature radius $r_{\mathrm{i}}/\mathrm{mm}$	11.446	0.0572
Outer raceway groove curvature radius $r_{\mathrm{e}}/\mathrm{mm}$	11.668	0.0583
Axial force $F_{\mathrm{a}}/\mathrm{N}$	17,800	890
Radial force $F_{\mathrm{r}}/\mathrm{N}$	17,800	890
Allowable yield stress ${\sigma }_{\mathrm{s}}/\mathrm{MPa}$	1100	55

Lain hypercube design is used to sample $d_{\mathrm{m}}$ , $D_{\mathrm{w}}$ , $r_{\mathrm{i}}$ , $r_{\mathrm{e}}$ , $F_{\mathrm{a}}$ , and $F_{\mathrm{r}}$ according to the normal distribution, and the sampled data are substituted by group into the static equilibrium equation of the rolling bearing to obtain the equivalent stress ${\sigma }_{\mathrm{VM}}$ . In this way, 500 groups of data are obtained, of which 60% are used as the training set, 20% as the validation set, and the remaining 20% as the test set. A three-layer BP neural network is selected and is a unidirectional forward propagation network composed of an input layer, a hidden layer, and an output layer. The numbers of nodes in the input, hidden, and output layers are 6, 10, and 1, respectively. According to the neural network mapping relationship, the training error is set to ${10}^{-5}$ . The input and output of the data are fitted within the allowable range of error to obtain a fitting relationship as shown in the following function:

$$\widehat{{\sigma }_{\mathrm{VM}}}=f_2[{\beta }_j+\sum _{k=1}^{10}{v}_{\mathrm{kj}}{f}_1({\alpha }_k+\sum _{i=1}^6{w_i}_k{x}_i)]$$

(70)

where $x_i$ is the input variable, ${w_i}_k$ and $v_{\mathrm{kj}}$ are the weights between the three layers, ${\alpha }_k$ and ${\beta }_j$ are the thresholds of the hidden layer and the output layer, respectively, $f_1(•)$ is the S-type nonlinear function sigmoid, and $f_2(•)$ is the linear function purelin.

The regression plot of the BP neural network is shown in Figure 5 and gives the correlation between the equivalent stress output by the neural network and all the data that constitute the training set, the validation set, and the test set, along with the mean of the three plot sets (R = 0.99995), indicating that the correlation between the output and the target is 99.995%, suggesting that the trained BP neural network has a good fitting effect.

Figure 5.

Regression plots of the BP neural network.

According to the stress–strength interference model, equation (39) to express the state function for the reliability analysis of rolling bearings based on plastic deformation under incomplete probability information as follows:

$$g(\mathit{X})={\sigma }_{\mathrm{s}}-f_2[{\beta }_j+\sum _{k=1}^{10}{v}_{\mathrm{kj}}{f}_1({\alpha }_k+\sum _{i=1}^6{w_i}_k{x}_i)]$$

(71)

Based on the theory derived earlier, a computer program for the reliability-based design and reliability sensitivity analysis of rolling bearings under incomplete probability information is developed. The data related to rolling bearings are substituted into the expressions for the first four moments of the state function for reliability analysis (equations (42)–(45)) and then substituted into the expressions for the unit reliability index and unit reliability of rolling bearings (equations (47) and (48)) for reliability-based design calculation to obtain the reliability index ${\beta }_{\mathrm{FM}}$ and reliability $R_{\mathrm{FM}}$ of rolling bearings as follows:

$${\beta }_{\mathrm{FM}}=2.7727$$

$$R_{\mathrm{FM}}=0.99722$$

Monte Carlo method has become the benchmark for reliability calculation because it is independent of calculation error and problem dimensionality and does not require the discretization of continuity problems. The reliability $R_{\mathrm{MC}}$ after ${10}^5$ Monte Carlo numerical simulation is obtained as

$$R_{\mathrm{MC}}=0.99805$$

The results obtained using the method proposed in this study are in good agreement with those obtained by the Monte Carlo method.

The sensitivity results calculated according to equations (66)–(69) are shown in Table 2. Based on Table 2 and Figure 6, the following findings are obtained:

1. The analysis results of the reliability sensitivity to mean values are shown in Figure 6(a). It can be seen from mean reliability sensitivity $\partial R_{\mathrm{FM}}/\partial {\mu }_X$ that the reliability of the rolling bearing based on plastic deformation is positively sensitive to the mean values of bearing pitch diameter $d_{\mathrm{m}}$ , ball diameter $D_{\mathrm{w}}$ , axial force $F_{\mathrm{a}}$ , and allowable yield stress ${\sigma }_{\mathrm{s}}$ , indicating that the reliability of the rolling bearing increases with an increase in their mean values; in other words, their influences are positive. The sensitivities to inner raceway groove curvature radius $r_{\mathrm{i}}$ , outer raceway groove curvature radius $r_{\mathrm{e}}$ , and radial force $F_{\mathrm{r}}$ are negative, indicating that their influences are negative. For example, in the engineering design of rolling bearings, improving the strength properties of the material can effectively reduce the occurrence of plastic deformation and increase the corresponding reliability; with an increase in inner raceway groove curvature radius $r_{\mathrm{i}}$ and outer raceway groove curvature radius $r_{\mathrm{e}}$ , the degree of coincidence of the bearing decreases and the corresponding reliability of the bearing based on plastic deformation also decreases.

2. The analysis results of the reliability sensitivity to standard deviations are shown in Figure 6(b). It can be seen from standard deviation reliability sensitivity $\partial R_{\mathrm{FM}}/\partial {\sigma }_X$ that the reliability of the rolling bearing decreases, that is, the bearing tends to fail, as the standard deviations of the seven parameters increase. The reliability is the most sensitive to the standard deviations of ball diameter $D_{\mathrm{w}}$ and inner raceway groove curvature radius $r_{\mathrm{i}}$ but not very sensitive to those of other parameters. Larger standard deviations of the basic random parameters indicate higher random dispersion of the parameters, and the standard deviations of the basic random parameters have a negative impact on the reliability of rolling bearings based on plastic deformation, which also coincides with actual engineering design provisions.

3. The analysis results of the reliability sensitivity gradients are shown in Figure 6(c). The gradients are sorted in descending order as follows: ① ball diameter $D_{\mathrm{w}}$ ; ② inner raceway groove curvature radius $r_{\mathrm{i}}$ ; ③ outer raceway groove curvature radius $r_{\mathrm{e}}$ ; ④ allowable yield stress, ${\sigma }_{\mathrm{s}}$ ; ⑤ bearing pitch diameter $d_{\mathrm{m}}$ ; ⑥ radial force, $F_{\mathrm{r}}$ ; and ⑦ axial force, $F_{\mathrm{a}}$ . On this basis, the degree of influence of the variation in each basic random parameter on the reliability can be evaluated, which can be used as an effective and practical tool for engineering modifications, reanalysis, and redesign.

4. Figure 6(d) shows the analysis results of the proportion of the dimensionless reliability gradient. It is seen from $s_{\mathrm{i}}$ and ${\lambda }_{\mathrm{i}}$ that allowable yield stress ${\sigma }_{\mathrm{s}}$ , radial force $F_{\mathrm{r}}$ , and ball diameter $D_{\mathrm{w}}$ have a large influence and lead to a high sensitivity, accounting for approximately 95.5% of the total, followed by inner raceway groove curvature radius $r_{\mathrm{i}}$ . Therefore, these four parameters should be strictly controlled in the design, production, and use process of rolling bearings. The results also show that the most important thing to determine whether rolling bearings undergo plastic deformation is the bearing material itself and the radial force acting on the bearing, which is also consistent with the actual engineering situation.

Table 2.

Reliability sensitivity of R_FM to μ _X and σ _X.

Basic random parameter	$\frac{\partial R_{\mathrm{FM}}}{\partial {\mu }_X}$	$\frac{\partial R_{\mathrm{FM}}}{\partial {\sigma }_X}$	$\mathrm{grad}[·]$	$s_i$	${\lambda }_i$ (%)
$d_{\mathrm{m}}$	2.6608 × 10−4	−2.6323 × 10−6	2.6609 × 10−4	3.6970 × 10−4	0.07
$D_{\mathrm{w}}$	8.4074 × 10−3	−2.0045 × 10−3	8.6430 × 10−3	4.4683 × 10−2	8.45
$r_{\mathrm{i}}$	−8.3287e × 10−3	−9.4307 × 10−4	8.3819 × 10−3	1.0835 × 10−2	2.05
$r_{\mathrm{e}}$	−7.8127 × 10−4	−1.0796 × 10−5	7.8135 × 10−4	1.3434 × 10−4	0.02
$F_{\mathrm{a}}$	4.2206 × 10−6	−8.3769 × 10−8	4.2214 × 10−6	4.0527 × 10−3	0.77
$F_{\mathrm{r}}$	−9.1981 × 10−6	−2.9562 × 10−6	9.6615 × 10−6	5.3402 × 10−2	10.10
${\sigma }_{\mathrm{s}}$	1.4516 × 10−4	−3.7626 × 10−4	4.0329 × 10−4	4.1512 × 10−1	78.53

Figure 6.

Sensitivity of the reliability to basic random parameters: (a) Mean sensitivity, (b) Standard deviation sensitivity, (c) Reliability sensitivity gradients and (d) The sensitivity of reliability factor level.

In summary, the analysis results reveal that improving the material properties of the bearing material and reducing the magnitude of the applied force are the most effective ways to reduce the failure probability of plastic deformation of rolling bearings. Regarding the bearing geometry, the ball diameter $D_{\mathrm{w}}$ and inner raceway groove curvature radius $r_{\mathrm{i}}$ should be strictly controlled to prevent plastic deformation. Figures 7 and 8 explain the relationship between the plastic deformation-based rolling bearing reliability and the design parameters $D_{\mathrm{w}}$ and $r_{\mathrm{i}}$ , respectively. The reliability increases and decreases with the increase of $D_{\mathrm{w}}$ and $r_{\mathrm{i}}$ , respectively, which is consistent with the reliability sensitivity analysis results.

Figure 7.

Relationship between reliability and ball diameter D_w.

Figure 8.

Relationship between reliability and inner raceway groove curvature radius r_i.

The theory and practice of reliability engineering show that reliability index is closely related to and highly sensitive to statistical characteristics, such as the mean and standard deviation of basic random parameters. The calculation results obtained from the theoretical model for the reliability of rolling bearings based on plastic deformation under incomplete probability information built in this study are consistent with the actual engineering design provisions. Furthermore, the theory established in this study quantifies the influence of basic random parameters on the reliability of rolling bearings based on plastic deformation, which is precisely the value of the theoretical reliability model and reliability sensitivity analysis established in this study.

Conclusion

Currently, due to the complexity of rolling bearings in actual projects and the relative lack of statistical data, the design parameters follow various types of probability distributions, and sometimes it is even impossible to determine their distribution types. For this reason, reliability analysis based on probabilistic and statistical information inevitably lead to results laden with errors or results far from the true solution. In this study, the problem of differential sensitivity analysis of reliability under limited probability information is investigated by using the higher-order moment method for reliability-based design and the plastic deformation analysis technique of rolling bearings. A simulation method for the differential sensitivity analysis of the reliability of rolling bearings under incomplete probability information is proposed, and the pattern of the variation in the reliability index of rolling bearings with changes in basic random parameters is investigated, providing a quantitative basis for the design of rolling bearings. This study has made an important contribution to the innovation and development of the theory of design of rolling bearings based on plastic deformation in engineering practice. In summary, the major results achieved in this study are emphasized as follows:

1. Due to the complexity of the static equilibrium equation of rolling bearings, it is impossible to directly use the higher-order moment method. This study fits the relationship between the basic random parameters and the equivalent stress through the BP neural network. Additionally, to address the deficiency of local optimality and premature convergence in artificial neural networks, PSO is introduced into the penalty function to optimize the BP neural network so that the optimized BP neural network can better predict the output of the function.

2. The matrix form equation for the reliability-based design of rolling bearings based on plastic deformation under incomplete probability information is derived. Additionally, the matrix form equation for the differential sensitivity analysis of the reliability is derived based on the direct differentiation method of sensitivity analysis. Matrix form equations are also derived for the reliability and reliability sensitivity under incomplete probability information. The matrix forms have the advantages of clear expression and easy programing.

3. Based on the stress–strength interference model, a theoretical method is proposed for the reliability-based design and reliability sensitivity analysis of rolling bearings based on plastic deformation under incomplete probability information and based on allowable yield stress. This method takes into account the geometric dimensions and material properties of the rolling bearing itself in addition to the influence of the combined load on the rolling bearing. Additionally, a certain type of angular contact ball bearing is used as an example to demonstrate that the method proposed in this study is correct and effective compared with the traditional Monte Carlo method.

The results of the reliability and reliability sensitivity analysis under incomplete probability information clarify the variation pattern and degree of influence of each basic random parameter on the reliability and reliability sensitivity of rolling bearings based on plastic deformation. Therefore, the engineering design of rolling bearings based on plastic deformation can be carried out in a reasonable and targeted manner. This study provides a novel theoretical method and implementation technique for the reliability-based design of rolling bearings in engineering practice.