Contact stress reliability analysis based on first order second moment for variable hyperbolic circular arc gear

Abstract

Aiming at the contact strength reliability of variable hyperbolic circular arc gear, a reliability analysis method for contact strength of variable hyperbolic circular arc gear based on Kriging model and advanced first-order and second-moment algorithm is proposed. Kriging model was used to establish the limit state equation of the contact stress reliability analysis of variable hyperbolic circular arc gear, and the advanced first-order second-moment method was used to analyze the contact stress reliability of variable hyperbolic circular arc gear based on the limit state equation of the contact stress. In order to verify the effectiveness of the proposed algorithm, a Markov Chain Monte Carlo reliability analysis method based on Important Sampling was proposed. Markov Chain and Important Sampling were exploited to improve the accuracy of contact reliability analysis based on Monte Carlo method for variable hyperbolic circular arc gear. The comparison between the analysis results of Markov Chain Monte Carlo with Important Sampling method and first order second moment shows that it is feasible to analyze the reliability of variable hyperbolic circular arc gear by first-order second-moment method.

Introduction

Gear is the most commonly used power transmission and motion transmission device. Installation error, manufacturing error, external load and other factors are random, which will lead to gear vibration, noise and pitting failure. These random factors seriously affect the reliability of gear meshing transmission. ¹ When study the reliability analysis of gear transmission, due to the limitations of human, material and financial factors, it is difficult to carry out real reliability test to obtain a large amount of useful data. Therefore, virtual simulation with sampling method has become one of the commonly used tools for reliability analysis instead of real response. ¹ At the same time, the surrogate model method can effectively reduce expensively time-consuming in virtual simulation. ² Kriging, response surface, RBF and so forth have been widely used analytical methods.^3–5 Tong et al. used response surface and Markov Chain Monte Carlo (MCMC) to study the influence of random factors such as installation error, manufacturing error and external load on the reliability of gear transmission. ¹ Yu et al. ⁶ established the reliability analysis model of gear thermal transfer error based on PC-Kriging model and active learning function, through which the reliability of gear transfer error was effectively investigated. Liu et al. ⁷ proposed a corrected-partial least squares regression to study the dynamic reliability analysis of a gear transmission system (GTS) of wind turbine (WT). Hu et al. ⁸ combined Response Surface Methodology (RSM) and first-order second-moment method (FOSM) analysis methods to study the effects of elastohydrodynamic lubrication (EHL) on contact fatigue reliability of spur gear. Bai et al. ⁹ proposed a MDCRSA method based on the response surface methodology and Monte Carlo to study the reliability analysis of the gear transmission errors of the steering mechanism of the gun (GTESMG). Yin et al. ¹⁰ set up the reliability limit state equation of the retraction system by busing a quadratic polynomial that has no cross terms and discussed how the key parameters affect the reliability of the retraction system. Zhang et al. ¹¹ derived a quadratic polynomial function by using the response surface method to describe the quantitative relationship between the tooth modification parameters and the dynamic transmission error fluctuations of helical planetary gear pair. He has established a quantitative relationship between the tooth modification parameters and the transmission error fluctuations, and carried out a reliability sensitivity analysis to demonstrate the effects of the tooth modification parameters on the dynamic transmission error fluctuations of the based on the quantitative relationship. ¹¹ Liang et al. ¹² proposed contact stress reliability limit state equation by using full quadratic RSM to study the effects of size parameters on reliability of the main reduction gears in a car. Zhang et al. ¹³ used Kriging model to simplify the calculated stress-strength reliability mathematics of gear transmission and employed genetic algorithm to globally optimize the volume and reliability of large ball mill gear transmission. Zhou et al. ¹⁴ used the adaptive Kriging model to study the vibration reliability and sensitivity of retraction-extension system under the condition of single failure and multiple failure modes. Yang and Tong ¹⁵ proposed reliability analysis method based on dimensionality reduction visualization and Kriging model to solve the problems of large computation and low precision during gear vibration reliability analysis. Cui et al. ¹⁶ used the Kriging method to simplify the reliability calculation model established on the basis of stress-strength interference theory. Zhang et al. ¹⁷ used the Kriging metamodel to build the approximate the real limit state function of the landing gear shock absorber to analyze the reliability of the multi-mode sensitivity analysis method of the shock absorber.

Variable hyperbolic circular arc gear (VHCAG) is a new type of gear transmission, which, theoretically, has the characteristics of high contact rate, high bearing capacity, no-axial force, high transmission efficiency, low noise, and so forth. ¹⁸

At present, there is no study on the contact strength reliability for the variable hyperbolic circular arc gear. This paper proposes a method for calculating the contact strength reliability of variable hyperbolic circular arc gear based on Kriging model and first-order second-moment. In the proposed method, the limit state equation of contact strength for reliability analysis of variable hyperbolic circular arc gear is established by using the Kriging model, and the first-order second-moment. The proposed reliability analysis method was used to analyze the contact strength reliability of variable hyperbolic circular arc gear based on the established limit state equation of contact strength.

Mathematical model and contact stress finite element model of circular arc gear

Cylindrical gear with arc tooth transmission device, shown in Figure 1, which is a new type of gear transmission has the advantages of good meshing performance, high coincidence, no axial force, stable transmission, and so on.

Figure 1.

Cylindrical gear with arc tooth.

Mathematical model

According to its forming principle, the coordinate system of forming principle of circular arc tooth line cylindrical gear is shown in Figure 2.^18,19

Figure 2.

Coordinate system of forming principle of cylindrical gear with arc tooth, R: dividing circle radius of gear; R_T: cutter head radius of machining tool; R_n: inner edge radius in the tangent direction of dividing circle; R_w: outer edge radius in the tangent direction of dividing circle; R_w = R_T+πm/4; m: gear module; a: pressure angle; ω: cutter head rotation angle speed of machining tool; V_T: cutter head moving speed of machining tool.

In Figure 1, S(O—XYZ) is the static coordinate system, S₁(O₁—X₁Y₁Z₁) is the solidification coordinate of the gear blank, and S_T(O_T—X_TY_TZ_T) is the tool coordinate, which moves relative to S(O—XYZ) coordinate at the speed of V_T = R×ω.

Unit normal vector of tool surface

The tool surface equation is:

$$\{\begin{array}{c}\overrightarrow{r_{\mathrm{dq}}}=\pm \sin \alpha \cos \theta \overrightarrow{i_{\mathrm{T}}}\pm \sin \alpha \sin \theta \overrightarrow{j_{\mathrm{T}}}+\cos \alpha \overrightarrow{k_{\mathrm{T}}}\hfill \\ \multicolumn{1}{c}{\overrightarrow{r_{\mathrm{d}\mathrm{\theta }}}=-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}+\pi m/4)\sin \theta \overrightarrow{i_{\mathrm{T}}}}\\ \multicolumn{1}{c}{+(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}+\pi m/4)\cos \theta \overrightarrow{j_{\mathrm{T}}}}\end{array}$$

(1)

Where, q is the distance from the tool to the X axis of the reference frame along the bus direction

While the unit normal vector of the tool surface is:

$$\begin{array}{c}\multicolumn{1}{c}{\overrightarrow{n}=\frac{\overrightarrow{r_{\mathrm{dq}}}\times \overrightarrow{r_{\mathrm{d}\mathrm{\theta }}}}{|\overrightarrow{r_{\mathrm{dq}}}\times \overrightarrow{r_{\mathrm{d}\mathrm{\theta }}}|}=-\cos \alpha \cos \theta \overrightarrow{i_{\mathrm{T}}}-\cos \alpha \sin \theta \overrightarrow{j_{\mathrm{T}}}}\\ \multicolumn{1}{c}{\pm \sin \alpha \overrightarrow{k_{\mathrm{T}}}}\end{array}$$

(2)

Relative speed of tool and gear at meshing point

The direction vector can be expressed as:

$$\overrightarrow{\lambda }=\overrightarrow{O_1{O}_{\mathrm{T}}}=(\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})\overrightarrow{i_{\mathrm{T}}}+R\overrightarrow{k_{\mathrm{T}}}$$

(3)

Therefore, the relative speed of the tool and the gear at the meshing point is:

$$\begin{array}{c}{\overrightarrow{v}}^{12}={\overrightarrow{\omega }}^{12}\times {\overrightarrow{r}}^1+\frac{d\lambda }{\mathrm{dt}}-{\overrightarrow{\omega }}^2\times \overrightarrow{\lambda }\hfill \\ =-{\omega }_{1}q\cos \alpha \overrightarrow{i_{\mathrm{T}}}-{\omega }_1(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}+\pi m/4)\hfill \\ \multicolumn{1}{c}{\cos \theta \overrightarrow{k_{\mathrm{T}}}+{\omega }_1\mathrm{R}\overrightarrow{i_{\mathrm{T}}}}\\ \multicolumn{1}{c}{-[{\omega }_1\mathrm{R}\overrightarrow{i_{\mathrm{T}}}-{\omega }_1(\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})\overrightarrow{k_T}]}\\ =-q\cos \alpha {\omega }_1\overrightarrow{i_{\mathrm{T}}}-[(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}+\pi m/4)\cos \theta \hfill \\ \multicolumn{1}{c}{+(\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})]{\omega }_1\overrightarrow{k_{\mathrm{T}}}}\end{array}$$

(4)

Meshing function

Based on the meshing principle, the meshing function is expressed as follows:

$$\mathrm{\Gamma }=\overrightarrow{\mathrm{n}}·{\overrightarrow{v}}^{12}=0$$

(5)

As ω1≠0, it can be concluded from the above formula that:

$$q=\mp \sin \alpha \frac{\cos \theta ({\mathrm{R}}_T+\pi m/4)+(\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})}{\cos \theta }$$

(6)

Conjugate surfaces

The contact line equation between the tool and the gear to be machined is as follows:

$$\{\begin{array}{c}{x}_{\mathrm{T}}=-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\cos \theta \hfill \\ \multicolumn{1}{c}{y_{\mathrm{T}}=(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\sin \theta }\\ \multicolumn{1}{c}{z_{\mathrm{T}}=q\cos \alpha }\\ \multicolumn{1}{c}{q=\mp \sin \alpha \frac{\cos \theta ({\mathrm{R}}_{\mathrm{T}}+\pi m/4)+(\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})}{\cos \theta }}\end{array}$$

(7)

By converting the coordinates in S_T (O_T—X_TY_TZ_T) to S₁(O₁—X₁Y₁Z₁), the tooth surface equation of the cut gear can be obtained as follows:

$$\{\begin{array}{c}{x}_1=[-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\cos \theta +\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}}]\hfill \\ \multicolumn{1}{c}{\cos \phi -(q\cos \alpha +\mathrm{R})\sin \phi }\\ \multicolumn{1}{c}{y_1=[-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\cos \theta +\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}}]}\\ \multicolumn{1}{c}{\sin \phi +(q\cos \alpha +\mathrm{R})\cos \phi }\\ \multicolumn{1}{c}{z_1=(\pm q\sin \alpha +{\mathrm{R}}_T\pm \pi m/4)\sin \theta }\\ \multicolumn{1}{c}{q=\mp \sin \alpha \frac{\cos \theta ({\mathrm{R}}_{\mathrm{T}}+\pi m/4)+(\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})}{\cos \theta }}\end{array}$$

(8)

Instantaneous contact line

Meshing function $\Gamma =\Gamma (q,\theta ,\phi )$ is a function of variables $q,\theta ,\phi$ . As ${\phi }_1={\omega }_{1}t$ , ${\phi }_1$ is regarded as a function of entering Г in a certain instant, thus the function expression of the instantaneous contact line between the rack cutter and the gear blank can be obtained. If the value ${\phi }_1$ at any time is substituted into equation (4), the contact line equation of the whole surface can be obtained as follows:

$$\{\begin{array}{c}{x}_T=-[(1+{\sin }^2\alpha )(R_T\pm \pi m/4)+{\sin }^2\alpha \frac{R\phi +R_T}{\cos \theta }]\cos \theta \hfill \\ \multicolumn{1}{c}{y_T=[(1+{\sin }^2\alpha )(R_T\pm \pi m/4)+{\sin }^2\alpha \frac{R\phi +R_T}{\cos \theta }]\sin \theta }\\ \multicolumn{1}{c}{z_T=\mp \sin \alpha \cos \alpha (R_T\pm \pi m/4+\frac{R\phi +R_T}{\cos \theta })}\end{array}$$

(9)

Tooth profile equation

In the axial middle section of the gear, from its extended coordinate system, it can be seen that b = 0, then $\theta$ is also equal to zero. Substituting these two parameters into equation (8), the tooth profile equation of the middle section is obtained:

$$\{\begin{array}{c}{x}_1=[-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)+\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}}]\hfill \\ \multicolumn{1}{c}{\cos \phi -(q\cos \alpha +\mathrm{R})\sin \phi }\\ \multicolumn{1}{c}{y_1=[-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)+\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}}]}\\ \multicolumn{1}{c}{\sin \phi +(q\cos \alpha +\mathrm{R})\cos \phi }\\ \multicolumn{1}{c}{q=\mp \sin \alpha ({\mathrm{R}}_{\mathrm{T}}+\pi m/4+\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})}\end{array}$$

(10)

According to equation (10), it can be seen that the tooth profile of the axial symmetry surface of the gear is in involute.

Similarly, on the axial asymmetric plane, considering that z1 = b, according to the expressions of Z₁ and q in equation (8), we can get:

$$\begin{array}{c}\multicolumn{1}{c}{\phi =\{[-\mathrm{b}/\tan \theta +({\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\cos \theta ]/{\sin }^2\alpha }\\ \multicolumn{1}{c}{-({\mathrm{R}}_{\mathrm{T}}+\pi m/4)\cos \theta -{\mathrm{R}}_{\mathrm{T}}\}/\mathrm{R}}\end{array}$$

(11)

Thus, the expression of the tooth profile of the non intermediate section can be obtained as follows:

$$\{\begin{array}{c}{x}_{\mathrm{T}}=[-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\cos \theta +\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}}]\hfill \\ \multicolumn{1}{c}{\cos \phi -(q\cos \alpha +\mathrm{R})\sin \phi }\\ \multicolumn{1}{c}{y_{\mathrm{T}}=[-(\pm q\sin \alpha +{\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\cos \theta +\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}}]}\\ \multicolumn{1}{c}{\sin \phi +(q\cos \alpha +\mathrm{R})\cos \phi }\\ \multicolumn{1}{c}{q=\mp \sin \alpha \frac{\cos \theta ({\mathrm{R}}_{\mathrm{T}}+\pi m/4)+(\mathrm{R}\phi +{\mathrm{R}}_{\mathrm{T}})}{\cos \theta }}\\ \multicolumn{1}{c}{\phi =\{[-\mathrm{b}/\tan \theta +({\mathrm{R}}_{\mathrm{T}}\pm \pi m/4)\cos \theta ]/{\sin }^2\alpha }\\ \multicolumn{1}{c}{-({\mathrm{R}}_{\mathrm{T}}+\pi m/4)\cos \theta -{\mathrm{R}}_{\mathrm{T}}\}/\mathrm{R}}\end{array}$$

(12)

Contact stress finite element model

The point cloud data of gear 3D model was calculated by numerical solution method based on the mathematical model of gear meshing principle established in sections 2.1, and precise gear 3D model was established by UG. The finite element analysis model of gear is established based on the 3D model. Contact stress finite element model of Variable Hyperbolic Circular Arc Gear is shown in Figure 3.

Figure 3.

Contact stress finite element model of variable hyperbolic circular arc gear: (a) constraint and (b) mesh.

Set contact type of finite element analysis of circular arc tooth line gear pair as “friction free.” The MPC (multi-point constraint) is established by the rotating center of the active and the driven. At the same time, MPC is added to the driving wheel: the direction of its rotation axis is set as free, and other rotation and translation are set as fixed. For the driven wheel, MPC is fixed and omnidirectional fixed. The method of sweeping is used to divide the mesh of gear pair, while the type of mesh is C3D8I. The mesh of gear tooth contact area is subdivided locally to improve the analysis accuracy of contact stress and reduce the calculation time. The calculated stress pattern of driving wheel and driven wheel is shown in Figure 4. The Figure 4(a) shows the dynamic contact stress diagram of driving wheel, and the Figure 4(b) shows the contact stress diagram of driving wheel. As can be seen from Figure 4, the values of contact stress of driven wheel and driving wheel are (a) and (b) respectively, which are basically the same, and the contact form is point contact, which is consistent with the characteristics of gear.

Figure 4.

Stress pattern analysis of driving wheel and driven wheel: (a) driving wheel and (b) driven wheel.

Limit state equation of contact stress of circular arc gear

Stress-strength interference model

Figure 5 shows the stress-strength interference model. The safety margin determined by conventional design methods in parts will attenuates, while the safety margin becomes smaller. The strength and stress distribution function varied from interval to crossover, and the crossover region is the unsafe region (width is $h$ ). When $g(x)>0$ , it indicates that the system is able to achieve the predetermined design function, defined. ²⁰

$$g(x)=g(x_1,x_2,\cdots ,x_n)=\sigma -\delta >0$$

(13)

according to the definition of reliability, the reliability is:

$$R=P(\sigma >\delta )=P(\sigma -\delta >0)$$

(14)

Figure 5.

Stress-strength interference model.

Figure 6 shows the calculation of part strength over stress in the interference region. According to the joint integration method of probability density function, the probability in the interference region is:

$$R=P(\sigma >\delta )={\int }_{-\infty }^{\infty }f(\sigma )[{\int }_{-\infty }^{\infty }f(\delta )d\delta ]d\sigma$$

(15)

Figure 6.

Schematic diagram of reliability calculation in interference region.

Limit state equation

Tooth Width ( $x_1$ ), Modulus ( $x_2$ ), Pressure Angle ( $x_3$ ), Radius Of Tooth Line ( $x_4$ ), Moment ( $x_5$ ) were taken as input variables. The optimal Latin Hyper design method was used to complete the experimental design. Contact stress under different experimental parameters was obtained by using the finite element method. Kriging, SVM and RBF were respectively used to establish gear contact stress prediction models. The accuracy indexes of Kriging, SVR and RBF are shown in Table 1.

Table 1.

Accuracy Indexes of surrogate model.

Model type	R 2	RMSE	RMAE
Kriging	0.9950	3.2887	0.2843
SVR	0.9941	3.5709	0.2506
RBF	0.7584	22.7731	1.4728

After comprehensive analysis of R², RMSE and RMAE, Kriging method has higher accuracy, so Kriging was selected as the method to establish the limit state equation of gear contact strength. In order to further improve its accuracy, an AQPSO-Kriging surrogate model was proposed in Ma et al. ¹⁸ The proposed AQPSO-Kriging surrogate model was used to calculate the gear contact stress. However, the AQPSO-Kriging surrogate model cannot directly provide the explicit expression of the established approximate model function between input and output. ¹⁹ In this paper, Isight and Matlab are used to jointly simulate and construct the binomial expression of the gear contact stress state equation. According to equation (13), it is obtained

$$\begin{array}{c}\multicolumn{1}{c}{\delta =892.858-48.424x_1-129.022x_2+27.193x_3}\\ \multicolumn{1}{c}{-2.389x_4+3.193x_5+0.842x_1^2+7.033x_2^2-0.351x_3^2}\\ +0.002x_4^2-0.002x_5^2+0.5753x_1{x}_2+0.1179x_1{x}_3\hfill \\ \multicolumn{1}{c}{+0.019x_1{x}_4-0.034x_1{x}_5-0.362x_2{x}_3+0.138x_2{x}_4}\\ -0.261x_2{x}_5-0.0005x_3{x}_4+0.0005x_3{x}_5-0.0001x_4{x}_5\hfill \end{array}$$

(16)

Gear selection material is $17\mathrm{CrNiMo}6$ , its allowable strength is $\left[\sigma \right]=520.8M\mathrm{pa}$ . According to equation (13), we have the equation of state of the gear. ²¹

$$\begin{array}{c}g(x)=\left[\sigma \right]-\delta =520.8-(892.858-48.424x_1\hfill \\ \multicolumn{1}{c}{-129.022x_2+27.193x_3-2.389x_4+3.193x_5}\\ \multicolumn{1}{c}{+0.842x_1^2}\\ +7.033x_2^2-0.351x_3^2+0.002x_4^2-0.002x_5^2\hfill \\ \multicolumn{1}{c}{+0.5753x_1{x}_2+0.1179x_1{x}_3+0.019x_1{x}_4}\\ \multicolumn{1}{c}{-0.034x_1{x}_5}\\ -0.362x_2{x}_3+0.138x_2{x}_4-0.261x_2{x}_5\hfill \\ \multicolumn{1}{c}{-0.0005x_3{x}_4+0.0005x_3{x}_5-0.0001x_4{x}_5)}\end{array}$$

(17)

Reliability analysis principle based on improved the first-order second-moment method (AFOSM)

Let the function of the structure be $Z=g(x_1,x_2,\cdots ,x_n)$ and variable $x_i~N(u_{x_i},{\sigma }_{x_i})(i=1,2,\cdots ,n)$ , define $F=\{x:g(x)\le 0\}$ as the failure zone. Under the functional function at that time, the possible failure point under the limit state is $P^{*}(x_1^{*},x_2^{*},\cdots ,x_n^{*})$ , and then $g(x_1,x_2,\cdots ,x_n)=0$ . The functional function is expanded by the Taylor expansion to the nonlinear functional function at the possible failure point, and (18) is taken as the linear part.

$$\begin{array}{c}Z=g(x_1,x_2,\cdots ,x_n)\hfill \\ \approx g(x_1^{*},x_2^{*},\cdots ,x_n^{*})+\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}(x_i-x_i^{*})\hfill \end{array}$$

(18)

According to literature, (18) can be rewritten as follows:

$$\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}(x_i-x_i^{*})=0$$

(19)

$$\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}{x}_i-\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}{x}_i^{*}=0$$

(20)

The exact solution of reliability index $\beta$ and failure probability $P_f$ based on the AFOSM algorithm can be obtained by (21) and (22).

$$\beta =\frac{\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}{u}_{x_i}-\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}{x}_i^{*}}{{\left[\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}^2{\sigma }_{x_i}^2\right]}^{\frac{1}{2}}}=\frac{\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}(u_{x_i}-x_i^{*})}{{\left[\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{p^{*}}^2{\sigma }_{x_i}^2\right]}^{\frac{1}{2}}}$$

(21)

$$P_f=\mathrm{\Phi }(-\beta )$$

(22)

The reliability calculation process of the AFOSM algorithm is shown in the Figure 7.

Figure 7.

Process of solving reliability by AFOSM algorithm.

Reliability analysis principle based on Markov Chain Monte Carlo

Monte Carlo

Suppose $X=(X_1,X_2,\cdots ,X_n)$ is a random variable of the structure, its probability distribution is $f_{\mathrm{x}}\left(x\right)$ , while $g(x)$ is a function, the following multi-dimensional integration can be used to compute the structure failure probability $P_{\mathrm{f}}$ .^22,23

$$\begin{array}{c}{P}_{\mathrm{f}}=\int \cdots {\int }_{g\left(x\right)\le 0}{f}_{\mathrm{x}}(x_1,x_2,\cdots ,x_n)d_{x_1}{d}_{x_2}\cdots d_{x_n}\hfill \\ =\int \cdots {\int }_{R^{\mathrm{n}}}{I}_F\left(x\right)f_{\mathrm{x}}(x_1,x_2,\cdots ,x_n)d_{x_1}{d}_{x_2}\cdots d_{x_n}\hfill \\ =E\left[I_F\left(x\right)\right]\hfill \end{array}$$

(23)

Where, $I_F\left(x\right)=\{\begin{array}{c}0,x\notin F\hfill \\ 1,x\in F\hfill \end{array}$ is the indicating function of the failure domain; $R^{\mathrm{n}}$ is n-dimensional variable space; $E[·]$ is the mathematical expectation operator.

The $j$ -th sample vector of random variable $X$ is expressed as $X_j=(x_{j1},x_{j2},\cdots ,x_{\mathrm{jn}})$ and the expected estimate of structural failure probability $P_{\mathrm{f}}$ is :

$${\hat{P}}_f=\frac{1}{N}\sum _{j=1}^N{I}_F\left(x_j\right)=\frac{N_f}{N}$$

(24)

Flow chart of reliability analysis using Monte Carlo method is shown in Figure 8. The whole process mainly includes the following three steps.

• Step 1. According to the distribution types and parameters of random variables, N groups of samples $x_j=(x_{j1},x_{j2},\cdots ,x_{\mathrm{jm}})(j=1,2,\cdots ,N)$ are randomly generated by a certain sampling method.

• Step 2. Substituting each group of samples $x_j$ into the limit state equation, the indicator function $I_F(x_j)$ value is determined, and the results that meet the relevant conditions are accumulated.

• Step 3. Estimating failure probability ${\hat{P}}_{\mathrm{f}}$ by formula (24).

Figure 8.

Flow chart of reliability analysis using Monte Carlo method.

Markov Chain Monte Carlo with important samples

Important sampling

In order to improve the accuracy of reliability calculation using Monte Carlo method and make the sample points more likely to fall into the failure domain as much as possible, importance sampling is adopted in this paper. The important sampling method makes the sample points more likely to fall into the failure domain by changing the center of random sampling, and reduces the number of simulations appropriately to ensure the accuracy and efficiency of the algorithm. Sample sampling is simulated based on Metropolis-Hastings criterion. The sample points extracted are more likely to fall into the failure domain or the so-called important region, which improves the solution accuracy. From the above, the formula (25) is obtained

$$\begin{array}{c}{P}_{\mathrm{f}}={\int }_{R^{\mathrm{n}}}\frac{I_F[g\left(v\right)]f_{\mathrm{x}}\left(v\right)}{p_V(v)}{p}_V(v)\mathrm{dv}\hfill \\ =E\left[\frac{I_F[g\left(v\right)]f_{\mathrm{x}}\left(v\right)}{p_V(v)}\right]\hfill \end{array}$$

(25)

In equation (25), $p_V(v)$ is the sampling density function in the important sampling, which increases the chance of sample points falling into the failure area through $p_V(v)$ . However, if most samples fall into the failure area, it also increases the possibility of failure invisibly. Therefore, the selection of the center point of sampling is a crucial problem. The general method is to use FOSM or AFOSM to solve the most likely failure point. Sampling was conducted with the most likely failure point as the sampling center.

$$\begin{array}{c}{\hat{P}}_f=E\left[\frac{I_F[g\left(v\right)]f_{\mathrm{x}}\left(v\right)}{p_V(v)}\right]\hfill \\ =\frac{1}{N}\sum _{i=1}^N\frac{I_F[g\left(v_i\right)]f_{\mathrm{x}}\left(v_i\right)}{p_V\left(v_i\right)}\hfill \end{array}$$

(26)

Suppose $h_V(v)=I_F[g\left(v\right)]f_{\mathrm{x}}\left(v\right)/p_V(v)$ , $h_V(v_i)$ are sample values, whatever $h_V(v)$ obeys, according to mathematical statistics theory, there are the following relationships.

$${\mu }_{{\hat{P}}_f}={\mu }_{h_V(v)},{\sigma }_{{\hat{P}}_f}^2={\sigma }_{h_V(v)}^2/N$$

(27)

The variance of ${\hat{P}}_f$ based on (25) and (27) is as follows:

$$\begin{array}{c}{\sigma }_{{\hat{P}}_f}^2=\frac{{\sigma }_{h_V(v)}^2}{N}=\frac{1}{N}E\left[\frac{I_F[g\left(v\right)]{f_{\mathrm{x}}}^2\left(v\right)}{{p_V}^2(v)}\right]-\frac{1}{N}{p}_f^2\hfill \\ =\frac{1}{N}{\int }_{R^{\mathrm{n}}}\frac{I_F[g\left(v\right)]{f_{\mathrm{x}}}^2\left(v\right)}{p_V(v)}\mathrm{dv}-\frac{1}{N}{p}_f^2\hfill \end{array}$$

(28)

According to equation (27), when $p_V(v)=\underset{v\in g\left(v\right)<0}{\int }{f}_{\mathrm{x}}\left(v\right)\mathrm{dv}$ , the variance of ${\hat{P}}_f$ is the smallest. There is:

$$h_V(v)=\frac{I_F[g\left(v\right)]f_{\mathrm{x}}\left(v\right)}{\underset{v\in g\left(v\right)<0}{\int }{f}_{\mathrm{x}}\left(v\right)\mathrm{dv}}$$

(29)

Formula (29) is called the sampling density function of optimal importance sampling.

Process of generating random numbers by rejection technique

The main idea of generating random numbers by rejection technique is whether the selected random number satisfies a criterion. If satisfied, it is an effective random number, otherwise it is generated again. Suppose variable $x\in [a,b]$ , and its probability density function is $f(x)$ . On $x\in [a,b]$ , $f(x)$ has an upper bound. The value is $f_0=\underset{a\le x\le b}{sup}f(x)$ . The specific process of generating random numbers by rejection technique is as follows:

1. Uniform distribution in $[0,1]$ produces two random numbers $R_1$ and $R_2$ ;

2. Let $x=a+(b-a)R_1$ , and calculate the function value $f(x)=f(x=a+(b-a)R_1)$ .

3. Discriminant $R_2{f}_0<f(a+(b-a)R_1)$ is valid. If $x=a+(b-a)R_1$ is set as the newly generated random number, there is a random value $x=a+(b-a)R_1$ at that time. If not, back off at step 1 to repeat the whole process to generate new samples.

Markov Chain simulation sample process

In order to generate simulation samples by Markov Chain, a probability density function $f^{*}(\xi |x)$ , and $\forall \xi ,f^{*}(\xi |x)=f^{*}(x|\xi )$ are selected, and $q(x)=I_F(x)f(x)$ is defined. Then, the Markov Chain simulation sample process can be realized according to the following steps:

• Step 1. The state value $X^{(0)}$ of Markov Chain is initialized (randomly generated in the failure region or determined).

• Step 2. The rejection technique is used to generate the random number $\xi =a+(b-a)R_1$ under state $X^{(i)}$ according to $f^{*}(\xi |x)$ , and $r=q(\xi )/q(X^{(i)})$ is calculated.

• Step 3. The next state value $X^{(i+1)}$ of Markov Chain is determined according to Metropolis-Hastings criterion.

• Step 4. Repeat Step 2, Step 3 until the number of samples reaches the required number.

Markov Chain Monte Carlo reliability analysis based on important sampling

The contact stress reliability analysis process based on important sampling MCMC is shown in Figure 9.

• Step 1. Generated Markov chains containing N samples. According to the state point $X^{(i)}$ generated by the Markov Chain simulation sample process, the mean ${\mu }_i$ and variance ${\sigma }_i$ are calculated.

• Step 2. Kernel density estimation, $f^{*}(\xi |x)$ is used as the sampling density function of the new important sampling, ${\mu }_i$ is used as the sampling center.

• Step 3. According to the principle of importance sampling, the important sampling technology is used for sampling in interval [ ${\mu }_i-3{\sigma }_i$ , ${\mu }_i+3{\sigma }_i$ ] to generated sample points for reliability calculation.

• Step 4. Using equation (23) to calculate reliability by using the sample points generated by the important sampling technology.

Figure 9.

Process of solving reliability by Markov Chain Monte Carlo algorithm.

Contact stress reliability analysis variable hyperbolic circular arc gear

Contact stress reliability analysis variable hyperbolic circular arc gear with AFOSM

According to AFOSM reliability and sensitivity analysis method, the limit state equation between input parameters (pressure Angle, tooth width, modulus, tooth radius, and torque) and output (contact stress) of variable hyperbolic circular arc tooth line cylindrical gear was established in Section 3. The reliability of contact strength of cylindrical gear with variable hyperbolic arc tooth line and its sensitivity to design variables, mean value and variance are calculated. Reliability index, failure probability and corresponding sensitivity and sensitivity coefficient are:

$$\beta =10.7962095888155$$

$$P_f=1.793538326955980\mathrm{e}-27$$

$$\frac{\partial P_f}{\partial \beta }\frac{\partial \beta }{\partial {\mu }^T}={10}^{-27}\times \left[\begin{array}{c}-7.13048846621923\hfill \\ \multicolumn{1}{c}{-61.3644643187168}\\ \multicolumn{1}{c}{0.156137127583956}\\ \multicolumn{1}{c}{-0.693084731665894}\\ \multicolumn{1}{c}{0.343260703833141}\end{array}\right]$$

$$\frac{\partial P_f}{\partial \beta }\frac{\partial \beta }{\partial {\sigma }^T}={10}^{-26}\times \left[\begin{array}{c}1.40555758505874\hfill \\ \multicolumn{1}{c}{20.8196761675733}\\ \multicolumn{1}{c}{0.000673941888692689}\\ \multicolumn{1}{c}{0.663976558991730}\\ \multicolumn{1}{c}{0.105145764493842}\end{array}\right]$$

$$\begin{array}{c}{\lambda }_i=-{\left(\frac{\partial g}{\partial x_i}\right)}_{x_i={\mu }_{\mathrm{xi}}}{\sigma }_{x_i}/\sqrt{\sum _{i=1}^n{\left(\frac{\partial g}{\partial x_i}\right)}_{x_i={\mu }_{\mathrm{xi}}}^2{\sigma }_{x_i}^2}\hfill \\ =\cos {\theta }_i=\left[\begin{array}{c}-0.185673568791041\hfill \\ \multicolumn{1}{c}{-0.392180640139827}\\ \multicolumn{1}{c}{0.00479550691646143}\\ \multicolumn{1}{c}{-0.672604954048888}\\ \multicolumn{1}{c}{0.599415757429034}\end{array}\right]\hfill \end{array}$$

Contact stress reliability analysis variable hyperbolic circular arc gear with Markov Chain Monte Carlo with important samples

In order to verify the correctness of AFOSM method. Markov Chain Monte Carlo with important samples method was used to calculate the reliability of the research object in this paper. Sample space of Markov Chain Monte Carlo with important samples is shown in Figure 10.

Figure 10.

Simulates sample distribution based on Markov Chain sample sampling: (a) x1-x2, (b) x1-x3, (c) x1-x4, (d) x1-x5, (e) x2-x3, (f) x2-x4, (g) x2-x5, (h) x3-x4, (i) x3-x5, and (j) x4-x5.

The sample points of Markov Chain simulation samples based on the Metropolis-Hastings criterion are all near the failure points determined by AFOSM algorithm, which improves the accuracy and reliability of reliability calculation.

Contact stress reliability analysis results of variable hyperbolic circular arc gear with Markov Chain Monte Carlo with important samples was obtained as follows:

$$P_f=1.72838330847567\mathrm{e}-27$$

Compared with AFOSM, the calculation accuracy is close to that of AFOSM. The reliability obtained by MCMC algorithm is close to that obtained by AFOSM algorithm, so it is correct and feasible to analyze the reliability of gears by AFOSM algorithm.

According to the definition of reliability, we have. ²⁴

$$R=1-P_f$$

(30)

The structural failure probability of the gear under the current working condition is close to 0, So R is close to 1, indicating that the gear is reliable under the current working condition.

Conclusions

At present, there is no study on the contact strength reliability for the variable hyperbolic circular arc gear. A method for calculating the contact strength reliability of variable hyperbolic circular arc gear based on Kriging model and AFOSM is proposed in this paper.

A method of construction the explicit expression of limit state equation of contact strength is proposed, and the explicit expression of limit state equation of contact strength for reliability analysis of variable hyperbolic circular arc gear was established.

The AFOSM reliability analysis method was used to analyze the contact strength reliability of variable hyperbolic circular arc gear based on the established limit state equation of contact strength. The structural failure probability calculated by the AFOSM is $1.793538326955980\mathrm{e}-27$ .

The correctness of the reliability algorithm in this paper is studied by Markov Chain Monte Carlo reliability analysis method based on Important Sampling. The structural failure probability calculated by the Markov Chain Monte Carlo reliability analysis method based on Important Sampling is $P_f=1.72838330847567\mathrm{e}-27$ .

The results show that the reliability algorithm in this paper can effectively realize the reliability analysis of variable hyperbolic circular arc gear. The structural failure probability of the gear under the current working condition is close to 0, indicating that the gear is reliable under the current working condition.