## Abstract

Rotate vector (RV) reducer is an essential mechanical transmission device in industrial machinery, robotics, aerospace, and other fields. The dynamic transmission characteristics and strength of the cycloidal pin gear and turning-arm bearing significantly affect the motion accuracy and reliability of RV reducer. Uncertainties from manufacturing and assembly errors and working loads add complexity to these effects. Developing effective methods for uncertainty propagation and reliability analysis for the RV reducer is crucial. In this work, the mail failure modes of RV reducer are studied, and an effective reliability analysis method for RV reducer considering the correlation between multi-failure modes is proposed by combining polynomial chaos expansion (PCE) and saddlepoint approximation method (SPA). This paper develops an uncertainty propagation strategy for RV reducer based on dynamic simulation and PCE method with high accuracy. On this basis, a surrogated cumulant generating function (CGF) and SPA are combined to analyze the stochastic characteristics of the failure behaviors. Then, based on the probability density function (PDF) and cumulative distribution function (CDF) calculated by SPA, the copula function is employed to quantify the correlations between the multi-failure modes. Further, the system reliability with multi-failure modes is estimated by SPA and optimal copula function. The validity of the proposed approach is illustrated by RV-320E reducer reliability estimation, and the results show that the proposed method can provide an effective reliability assessment technology for complex system under unknown physical model and distribution characteristics.

## 1 Introduction

Rotate vector (RV) reducer is a two-stage closed differential gear mechanism, which is composed of a planetary gear reducer on the high-speed stage, and a cycloidal pin gear drive on the low-speed stage [1]. With the development of intelligent equipment, RV reducer is more widely applicated, and the performance requirements and reliability of RV reducer have become increasingly strict. Complicated structure and multi-failure modes interaction make the uncertainty propagation and reliability analysis for RV reducer a challenge problem [2,3].

Dynamic analysis is a key step to study the performance of RV reducer, and many factors affect the dynamic transmission performance and strength of RV reducer, such as the manufacturing and assembly errors, working loads, and tooth modification [4,5]. The dynamic analysis of RV reducer is a highly nonlinear problem, and the addition of uncertainties further complicates the dynamic transmission analysis.

Many researches have been developed from different aspects to study the dynamic transmission performance. Ren et al. [6] analyzed the dynamic transmission accuracy of RV reducer based on a new tooth modification method. Li et al. [7] proposed a meshing contact analysis method for cycloidal pin gear transmission considering the influence of manufacturing error. Wang et al. [8] analyzed the transmission accuracy considering the influence of multi-tooth contact stiffness and profile modification. The effect of geometrical design parameters on the dynamic transmission performance of RV reducer was also studied [9,10]. Cao et al. [11] studied the dynamic transmission performance of a RV reducer considering the rigid–flexible coupled effect. Zhang et al. [12] studied the transmission error and torsional stiffness of RV reducer under wear. Han [13] proposed a global sensitivity analysis method to analyze the effects of manufacturing errors, assembly errors, bearing clearances, etc., on the transmission accuracy of RV reducer. The dynamic performance of RV reducer has been extensively studied based on virtual simulation method [14,15]; however, the influence of uncertainties on the dynamic transmission performance is not well studied.

Uncertainties from manufacturing and assembly errors and working loads make the dynamic transmission behavior of RV reducer is stochastic, and performance estimation and design optimization of RV reducer considering the effects of uncertainties are important issues. Sun et al. [16] analyzed the reliability of RV reducer based on fault tree and Monte Carlo simulation (MCS) method. Qian et al. [17,18] proposed a time-variant reliability modeling and estimation method for RV reducer based on Kriging method, and the strength of RV reducer was studied. Yang et al. [19] analyzed the fatigue strength of RV reducer under uncertainties and proposed a reliability-based design optimization method.

The available reliability estimation methods mostly consider fatigue strength; however, few researches focus on dynamic transmission accuracy reliability. Both the strength and dynamic transmission accuracy of RV reducer needed to be focused on Refs. [20] and [21]: (1) fatigue pitting between the cycloidal gear teeth and the pin tooth surface; (2) the bending strength and stiffness of the needle tooth pins; (3) the fatigue life of the turning-arm bearing; and (4) the dynamic transmission accuracy. Therefore, in this paper, a reliability estimation method is proposed considering both the strength and dynamic transmission accuracy.

Reliability modeling and complex uncertainties propagation analysis are two important issues for reliability estimation of complex mechanical system. For such system, the limit state function is usually not available, and performance simulation and physical tests should be conducted. To improve the reliability modeling efficiency and accuracy, sampling strategy combining surrogate modeling methods has been employed [22–27]. PCE methods approximate the limit state function as a linear combination of orthogonal polynomial functions [28,29], which can accurately describe the randomness of random variables and propagate uncertainty. In this paper, considering the effects of uncertainties from working loads and manufacturing and assembly errors, a system reliability model of RV reducer is developed based on PCE method. The dynamic transmission error of RV reducer and the maximal contact force are analyzed by adams simulation [30]. On this basis, the dynamic transmission accuracy reliability and strength reliability estimation models are formulated.

Reliability analysis methods including sampling-based methods [31], most probable failure point (MPP)-based methods [32,33], moment-based methods [34,35], probability density evolution-based method [36], and direct probability integral method [37] have been deeply developed to analyze the failure probability of mechanical system. To estimate the reliability problem with highly nonlinear response function or unknown response function, moment-based methods and saddlepoint approximation method (SPA)-based methods have been developed. Du et al. [38–41] developed a first-order and second-order SPA-based reliability analysis method by linearizing the limit state function at the MPP point. Meng et al. [42] proposed a mean-value first-order SPA-based reliability estimation method. To address highly nonlinear problem, high order moments-based saddlepoint approximation (TMSA) methods were proposed, and an approximated cumulant generating function (CGF) was employed to deal with highly nonlinear response function without requiring the existence of the CGF [43,44].

In this work, to improve the efficiency and accuracy of reliability modeling, data-driven PCE method is employed to analyze the multi-failure behavior of RV reducer. On this basis, the first four order moments of the failure behavior are calculated. A reliability analysis method based on an approximated CGF combining SPA is developed to analyze the probability of failure. To deal with the correlation between the multi-failure modes, standard normal distribution transformation method [40] and copula function-based method [45–47] were thoroughly investigated. Copula functions provide different forms to quantify the nonlinear correlation of random variables and link marginal distributions with a joint distribution. In this paper, the proper copula function is chosen to describe the correlation between the multi-failure modes, and the reliability of RV reducer under multi-failure modes is estimated by SPA and copula function.

The remainder of this paper is organized as follows: in Sec. 2, the physical model of RV reducer and the main failure modes are analyzed and discussed based on dynamic simulation. In Sec. 3, an efficient reliability estimation method is developed based on data-driven PCE and SPA, and the reliability estimation process is validated based on a numerical problem. In Sec. 4, the effects of uncertainties from load and manufacturing and assembly errors on the failure behaviors are analyzed based on PCE, and a reliability estimation method is developed by combining an approximated CGF, SPA method, and copula function for RV 320E with multi-failure modes. In Sec. 5, some conclusions are presented.

## 2 Failure Modes Analysis of Rotate Vector Reducer Based on Dynamic Simulation

Rotate vector reducer is a closed differential gear mechanism, as shown in Fig. 1. The first stage is a planetary gear transmission reducer, while the second stage is a cycloid pin-wheel drive. The crankshafts connected with the planetary gear rotate at the same speed. The output mechanism is driven by $nc$ pairs of crankshafts supported by bearings.

The main failure modes of cycloidal pin gear planetary transmission are: (1) fatigue pitting between the cycloidal gear teeth and the pin tooth surface. (2) Insufficient bending strength and stiffness of the needle tooth pins, which can result in their bending deformation or breakage. (3) The fatigue failure of the turning-arm bearing due to the large force on the turning-arm bearing and the relatively high relative speed of the inner and outer rings. (4) The dynamic transmission accuracy does not satisfy the requirement.

In this paper, RV-320E reducer is taken as an object, and the structure parameters and dynamic performance parameters of RV-320E reducer are, respectively, listed as in Tables 1 and 2.

Description | Symbols | Value |
---|---|---|

Number of crankshaft | $nc$ | 2 |

Tooth number of sun gear | $zs$ | 20 |

Tooth number of planetary gear | $zsp$ | 60 |

Module of planetary gear/mm | $msp$ | 0.6 |

Meshing pressure angle of planetary gear/deg | $\alpha sp$ | 20 deg |

Basic circle radius of the planetary gear/mm | $rgp$ | 36 |

Central distance of the sun–planet/mm | $asp$ | 48 |

Tooth number of cycloid gear | $zc$ | 39 |

Tooth number of pin gear | $zp$ | 40 |

Radius of pin-wheel/mm | $rp$ | 145.5 |

Eccentric distance of cycloid gear/mm | $e$ | 2.2 |

Pitch radius of the cycloid gear/mm | $rcp$ | 5 |

Description | Symbols | Value |
---|---|---|

Number of crankshaft | $nc$ | 2 |

Tooth number of sun gear | $zs$ | 20 |

Tooth number of planetary gear | $zsp$ | 60 |

Module of planetary gear/mm | $msp$ | 0.6 |

Meshing pressure angle of planetary gear/deg | $\alpha sp$ | 20 deg |

Basic circle radius of the planetary gear/mm | $rgp$ | 36 |

Central distance of the sun–planet/mm | $asp$ | 48 |

Tooth number of cycloid gear | $zc$ | 39 |

Tooth number of pin gear | $zp$ | 40 |

Radius of pin-wheel/mm | $rp$ | 145.5 |

Eccentric distance of cycloid gear/mm | $e$ | 2.2 |

Pitch radius of the cycloid gear/mm | $rcp$ | 5 |

Performance parameters | Value |
---|---|

Rated output rotating speed | $15\u2009rpm$ |

Rated output torque | $3136\u2009N\u22c5m$ |

Start–stop allowable torque | $7840\u2009N\u22c5m$ |

Instantaneous maximum allowable torque | $15,680\u2009N\u22c5m$ |

Maximum output rotating speed | $50\u2009rpm$ |

Performance parameters | Value |
---|---|

Rated output rotating speed | $15\u2009rpm$ |

Rated output torque | $3136\u2009N\u22c5m$ |

Start–stop allowable torque | $7840\u2009N\u22c5m$ |

Instantaneous maximum allowable torque | $15,680\u2009N\u22c5m$ |

Maximum output rotating speed | $50\u2009rpm$ |

To calculate the performance response, the first step is to calculate the dynamic contact force in the cycloidal pin gear.

### 2.1 Theory Analysis of Transmission Principle.

where $\Delta rrp$ is the displacement modification amount, $\Delta rrp$ is the equidistant modification amount, $\phi $ is the intersection angle between *X-*axis and the line connecting the point on the tooth profile and the origin $Ocg$, $e$ is the eccentric distance of the cycloid wheel, $iH$ is the transmission ratio of cycloidal pin stage transmission, $Sr=1+K12\u22122K1\u2009cos\u2009\phi $, and $K1$ is the short amplitude coefficient of the cycloid gear.

The second step is to determine how many tooth pairs are engaged. Assume that output torque is $Tout$, the torque on each cycloid gear of the RV reducer is $Tout/2$ based on the transmission principle. Further, based on the contact model as displayed in Fig. 2, the number of meshing tooth pairs can be determined by the algorithm in Fig. 3.

The loads acting on the cycloid gear must be in load equilibrium with the turning-arm bearing reaction forces as shown in Fig. 4.

where $Tcg$ is the torque transmitted by the cycloidal gear, and $Fpk$ is the meshing force between the $kth$ pair of teeth. $lk=rcp(1+K12\u22122K1\u2009cos(\phi pk))\u22121/2\u2009sin(\phi pk)$ is the moment arm of force, and the force arm corresponding to the maximum load $Fmax$ is $rcp$.

Given the design parameters of RV reducer and the output torque, the time-variant contact deformation can be calculated as shown in Table 3.

Algorithm 1: The contact deformation calculation |

Step 1: Set the transmission gear parameters and input torque $Tcg$, manufacturing and assembly errors. |

Step 2: Analysis the curvature of cycloid tooth profile with and without gear profile modification based on Eq. (1). |

Step 3: Transform the gear profiles into global coordinate system. |

Step 4: Determine the meshing tooth pairs based on Algorithm 1 in Fig. 3, and calculate the contact force on the $ith$meshing tooth pair: $\u2211Fpk=Tcg2eZc,Fmax=Tcg\u2211k=1nk(lkrc\u2212\Delta \phi pk\delta max)lk$ |

Step 5: Contact deformation calculation $\delta k=lk\beta =\u2009sin\u2009\phi pk1+K12\u22122K1\u2009cos\u2009\phi pk\delta max$ |

Algorithm 1: The contact deformation calculation |

Step 1: Set the transmission gear parameters and input torque $Tcg$, manufacturing and assembly errors. |

Step 2: Analysis the curvature of cycloid tooth profile with and without gear profile modification based on Eq. (1). |

Step 3: Transform the gear profiles into global coordinate system. |

Step 4: Determine the meshing tooth pairs based on Algorithm 1 in Fig. 3, and calculate the contact force on the $ith$meshing tooth pair: $\u2211Fpk=Tcg2eZc,Fmax=Tcg\u2211k=1nk(lkrc\u2212\Delta \phi pk\delta max)lk$ |

Step 5: Contact deformation calculation $\delta k=lk\beta =\u2009sin\u2009\phi pk1+K12\u22122K1\u2009cos\u2009\phi pk\delta max$ |

### 2.2 Dynamics Simulation Analysis.

where $Js$, $Jp$, $Jh$, $Jc$, and $Jo$ are the rotational inertia of sun gear, planetary gear, crankshaft, cycloid gear, and planet carrier, respectively. $kb$ and $kq$ denote the torsional stiffness of input shaft and crank shaft, respectively. $\theta a$, $\theta s$, $\theta hi$, $\theta pi$, $\theta cj$, and $\theta o$ are the torsion angles of input shaft, sun gear, planetary gear, crankshaft, cycloid gear, and planet carrier, respectively. $mp$, $mh$, and $mc$ denote the mass of the planetary gear, crankshaft, and cycloid gear, respectively. $rh$ is the radius of the distribution circle of the crankshaft. $xhi$ and $xcj$ are the displacement of crankshaft and cycloid gear, respectively. $F=[Fspi,Fhicj,Fckj]$ denote the force vector of components in RV reducer system, which is composed of the contact force between the sun gear and the planetary gear, contact force between the planetary carrier and the crankshaft, and contact force between the cycloidal gear and the needle gear.

Dynamic analysis is a highly nonlinear complex problem; to solve this dynamic equation under the effects of uncertainties, dynamics simulation analysis by adams is used. Based on the input uncertainties, sampling points are chosen by Latin hypercube sampling [48], and the corresponding dynamics analysis is carried out.

To establish the dynamics simulation analysis model displayed in Fig. 5, the maximum payload of the RV-320E reducer is set as 3136 $N\u22c5m$, the maximum output speed is 15 $rpn$, transmission ratio is 118.5, and the simulation time is 5 s. Define the step function as the input for the sun gear, the dynamic transmission force on one needle tooth is analyzed as shown in Fig. 6, and the transmission angular velocity and displacement of cycloidal gear are analyzed as shown in Fig. 7.

The meshing force on the needle teeth participating in meshing gradually increases from 0 to maximum and then decreases to 0 with the change of the contact area.

where $\theta in$ is the rotation angle of input gear shaft, $\theta out,p$ is the rotation angle of output shaft by dynamics simulation analysis, $\theta out,th$ denotes the theoretical rotation angle of output shaft, and $iR$ is the transmission ratio of RV reducer. The dynamic transmission error is calculated based on the dynamics simulation results as shown in Fig. 8.

### 2.3 Failure Modes Analysis of Cycloid Pin Gear.

where $\sigma H$ denotes the contact stress on the cycloidal gear teeth, and $\sigma Hp$ is the allowable contact stress. $Ee$ is the equivalent elastic modulus, $Ee=E1E2E1+E2$, $E1$ and $E2$ denote the elastic modulus of the needle teeth and the cycloidal gear, respectively, $bc$ is the width of the cycloidal gear, and $\rho e=|\rho rzp\u2212rz|$ is the equivalent curvature radius.

where $\sigma F$ is the bending stress of needle tooth pin, and $\sigma Fp$ denote the allowable bending stress. $Fmax$ is the maximum pressure applied on the needle tooth, $L$ is the span of the needle tooth pin, and $dz$ is the diameter of the needle tooth pin.

where $Lh$ is the fatigue life, and $n$ is the speed of turning-arm bearing. $C$ denotes the basic rated dynamic load of turning-arm bearing. $Fe=KdFr$ is the equivalent dynamic load, and $Fr$ is the force subjected on the turning-arm bearing by the cycloidal gear.

where $\theta errorlim$ is the allowable dynamic transmission error.

## 3 Reliability Estimation by Combing Polynomial Chaos Expansion and Saddlepoint Approximation Method

In this paper, the PCE and SPA are combined to deal with the failure physics modeling and reliability evaluation of complex system with complex structure and complicated random response.

### 3.1 Reliability Modeling Based on Polynomial Chaos Expansion Method.

where $a=[a0,ai1,ai1i2,ai1i2i3,\cdots ],i=1,2,\cdots ,n$ is the expansion coefficients to be calculated. $\psi =[\psi 1,\psi 2,\psi 3,\cdots ]$ is an *p*-order generalized chaos polynomial of the *n*-dimensional uncertainty variable $\xi =[\xi 1,\xi 2,\cdots ,\xi n]$. The orthogonal polynomials corresponding to different distributions are listed in Table 4.

Distribution types | Orthogonal polynomials | Weight functions | Variables range |
---|---|---|---|

Normal distribution | Hermite | $e\u2212x2/2$ | $(\u2212\u221e,+\u221e)$ |

Uniform | Legendre | 1 | $(\u22121,1)$ |

Exponential | Laguerre | $e\u2212x$ | $(0,+\u221e)$ |

Beta | Jacobi | $(1\u2212x)\alpha (1+x)\beta $ | $(\u22121,1)$ |

Gamma | General Laguerre | $xae\u2212x$ | $(0,+\u221e)$ |

Distribution types | Orthogonal polynomials | Weight functions | Variables range |
---|---|---|---|

Normal distribution | Hermite | $e\u2212x2/2$ | $(\u2212\u221e,+\u221e)$ |

Uniform | Legendre | 1 | $(\u22121,1)$ |

Exponential | Laguerre | $e\u2212x$ | $(0,+\u221e)$ |

Beta | Jacobi | $(1\u2212x)\alpha (1+x)\beta $ | $(\u22121,1)$ |

Gamma | General Laguerre | $xae\u2212x$ | $(0,+\u221e)$ |

*P*, and the prediction model can be approximated as

### 3.2 Reliability Estimation by Saddlepoint Approximation Method and the First Four Moments.

*x*is available, the moment generating function (MGF) of a random variable is defined by the following equation:

*x*, and the CGF is the logarithm of MGF as shown in the following equation:

### 3.3 Results Validation and Discussion.

Step 1: Performance modeling based on PCE method

Based on the random distribution of random inputs, the Laguerre polynomial with two orders is employed to establish the PCE modelThe coefficients are solved based on least square method, the first fourth statistical moments are calculated based on Eq. (28), and the results are displayed in Table 5.$gpre\u2248a0+a11(1\u2212x1)+a12(x12\u22124x1+2)+a21(1\u2212x2)+a22(x22\u22124x2+2)$(27)- Step 2: CGF approximation: The CGF $g(x)$ is approximated, and the coefficients are solved based on the following equation:${Kg(t)=q1t+q2t2\u2212q3ln{(1\u2212q4t)2}Kg(1)(t)=q1+2q2t+2q3q4/(1\u2212q4t)=M1,gKg(2)(t)=2q2+2q3q42/(1\u2212q4t)2=M2,g2Kg(3)(t)=4q3q43/(1\u2212q4t)3=M3,gM2,g3Kg(4)(t)=12q3q44/(1\u2212q4t)4=M4,gM2,g4$(28)
- Step 3: Saddlepoint $tg$ calculation according to the equation as$Kg\u2032(tg)\u2212glim=0$(29)
- Step 4: Reliability estimation: The reliability of $g\u2264glim$ can be evaluated based on saddlepoint $tg$ by the following equation:$R=Pr(g\u2264glim)\u2248\Phi (w)+\phi (w)(1w\u22121v)$(30)

Statistical moments | Mean value | Standard deviation | Skewness | Kurtosis |
---|---|---|---|---|

MCS ($106$ samples) | 1.9977 | 1.4134 | 1.4162 | 6.0079 |

PCE method (20 samples) | 1.9977 | 1.4134 | 1.4162 | 6.0079 |

Statistical moments | Mean value | Standard deviation | Skewness | Kurtosis |
---|---|---|---|---|

MCS ($106$ samples) | 1.9977 | 1.4134 | 1.4162 | 6.0079 |

PCE method (20 samples) | 1.9977 | 1.4134 | 1.4162 | 6.0079 |

$a$ | FORM | SORM | FOSPA | MCS ($106$samples) | PCE + FOSPA |
---|---|---|---|---|---|

0.5 | 0.1795 | 0.2301 | 0.2482 | 0.2474 | 0.2481 |

1 | 0.0990 | 0.1393 | 0.1459 | 0.1452 | 0.1452 |

2.0 | 0.0286 | 0.0466 | 0.0469 | 0.0466 | 0.0466 |

3.0 | 0.0079 | 0.0145 | 0.0142 | 0.0141 | 0.0141 |

4.0 | 0.0021 | 0.0043 | 0.0041 | 0.0041 | 0.0042 |

5.0 | 0.0006 | 0.0012 | 0.0012 | 0.0012 | 0.0012 |

6.0 | 0.0001 | 0.0003 | 0.0003 | 0.0003 | 0.000,335 |

$a$ | FORM | SORM | FOSPA | MCS ($106$samples) | PCE + FOSPA |
---|---|---|---|---|---|

0.5 | 0.1795 | 0.2301 | 0.2482 | 0.2474 | 0.2481 |

1 | 0.0990 | 0.1393 | 0.1459 | 0.1452 | 0.1452 |

2.0 | 0.0286 | 0.0466 | 0.0469 | 0.0466 | 0.0466 |

3.0 | 0.0079 | 0.0145 | 0.0142 | 0.0141 | 0.0141 |

4.0 | 0.0021 | 0.0043 | 0.0041 | 0.0041 | 0.0042 |

5.0 | 0.0006 | 0.0012 | 0.0012 | 0.0012 | 0.0012 |

6.0 | 0.0001 | 0.0003 | 0.0003 | 0.0003 | 0.000,335 |

From Table 6, it can be found that the combination of PCE and SPA can effectively solve the problem of structural uncertainty propagation and reliability evaluation under unknown failure physical models.

## 4 Reliability Analysis of Rotate Vector Reducer Under Multifailure Modes

The uncertainties from profile modification, displacement modification, payload, and assembly error have great influence on the strength and dynamic transmission accuracy of cycloid gear. Therefore, the effects are analyzed based on simulation and theory model. The input stochastic characteristics are listed in Table 7.

Parameters | Distribution | Mean value | Standard deviation |
---|---|---|---|

Profile modification $\Delta rp/mm$ | Normal | 0.0833 | 0.0083 |

Displacement modification $\Delta rrp/mm$ | Normal | −0.0532 | 0.0053 |

Payload/$N\u22c5m$ | Normal | 3136 | 31.36 |

Assembly error/mm | Normal | 0.002 | 0.0002 |

Parameters | Distribution | Mean value | Standard deviation |
---|---|---|---|

Profile modification $\Delta rp/mm$ | Normal | 0.0833 | 0.0083 |

Displacement modification $\Delta rrp/mm$ | Normal | −0.0532 | 0.0053 |

Payload/$N\u22c5m$ | Normal | 3136 | 31.36 |

Assembly error/mm | Normal | 0.002 | 0.0002 |

Copula function | Copula model |
---|---|

Gaussian copula | $C(u,v)=\Phi 2(\Phi \u22121(u),\Phi \u22121(v),\theta k),\u2009\u22121\u2264\theta k\u22641$ |

Frank copula | $C(u,v)=1\theta kln(1+(e\theta ku\u22121)(e\theta kv\u22121)e\theta k\u22121),\u2009\u2212\u221e<\theta k<\u221e$ |

Gumbel copula | $C(u,v,\theta k)=exp(\u2212[((\u2212lnu)1\theta k+(\u2212lnv)1\theta k)]\theta k)$ |

Clayton copula | $C(u,v,\theta k)=(u\u2212\theta k+v\u2212\theta k\u22121)\u22121/\theta k$ |

Copula function | Copula model |
---|---|

Gaussian copula | $C(u,v)=\Phi 2(\Phi \u22121(u),\Phi \u22121(v),\theta k),\u2009\u22121\u2264\theta k\u22641$ |

Frank copula | $C(u,v)=1\theta kln(1+(e\theta ku\u22121)(e\theta kv\u22121)e\theta k\u22121),\u2009\u2212\u221e<\theta k<\u221e$ |

Gumbel copula | $C(u,v,\theta k)=exp(\u2212[((\u2212lnu)1\theta k+(\u2212lnv)1\theta k)]\theta k)$ |

Clayton copula | $C(u,v,\theta k)=(u\u2212\theta k+v\u2212\theta k\u22121)\u22121/\theta k$ |

Algorithm 2: PCE combining SPA method to calculate the system reliability of RV reducer |

Step 1: Quantify the uncertainties of RV reducer considering the random load and manufacturing and assembly errors. Choose representative points in the stochastic variables space. |

Step 2: Establish the dynamic model of RV reducer considering the random inputs, and analyze the four-performance response on each representative point. |

Step 3: Based on the random input, Hermite polynomial is employed to establish PCE models for the four-failure modes. ${\phi 1(\xi )=a0+a11\xi 1+a21\xi 2+a12(\xi 12\u22121)+a22(\xi 22\u22121)+a3\xi 1\xi 2\phi 2(\xi )=b0+b11\xi 1+b21\xi 2+b12(\xi 12\u22121)+b22(\xi 22\u22121)+b3\xi 1\xi 2\phi 3(\xi )=c0+c11\xi 1+c21\xi 2+c12(\xi 12\u22121)+c22(\xi 22\u22121)+c3\xi 1\xi 2\phi 4(\xi )=d0+d11\xi 1+d21\xi 2+d12(\xi 12\u22121)+d22(\xi 22\u22121)+d3\xi 1\xi 2$ |

Step 4: Calculate the first fourth central moments, and the approximated CGF is formulated for each failure mode by Eq. (20). |

Step 5: Saddlepoint calculation and reliability estimation for each failure mode ${Kg1(1)(tg1)=g1lim,Kg2(1)(tg2)=g2limKg3(1)(tg3)=g3lim,Kg4(1)(tg4)=g4limRi=Pr(gi\u2264gilim)\u2248\Phi (wi)+\phi (wi)(1wi\u22121vi)$ |

Step 6: Analyze the correlation between multi-failure modes based on copula function |

Step 7: System reliability estimation for RV reducer considering the correlation of the four-failure mode $Rsys=1\u2212Pf,sys$ |

Algorithm 2: PCE combining SPA method to calculate the system reliability of RV reducer |

Step 1: Quantify the uncertainties of RV reducer considering the random load and manufacturing and assembly errors. Choose representative points in the stochastic variables space. |

Step 2: Establish the dynamic model of RV reducer considering the random inputs, and analyze the four-performance response on each representative point. |

Step 3: Based on the random input, Hermite polynomial is employed to establish PCE models for the four-failure modes. ${\phi 1(\xi )=a0+a11\xi 1+a21\xi 2+a12(\xi 12\u22121)+a22(\xi 22\u22121)+a3\xi 1\xi 2\phi 2(\xi )=b0+b11\xi 1+b21\xi 2+b12(\xi 12\u22121)+b22(\xi 22\u22121)+b3\xi 1\xi 2\phi 3(\xi )=c0+c11\xi 1+c21\xi 2+c12(\xi 12\u22121)+c22(\xi 22\u22121)+c3\xi 1\xi 2\phi 4(\xi )=d0+d11\xi 1+d21\xi 2+d12(\xi 12\u22121)+d22(\xi 22\u22121)+d3\xi 1\xi 2$ |

Step 4: Calculate the first fourth central moments, and the approximated CGF is formulated for each failure mode by Eq. (20). |

Step 5: Saddlepoint calculation and reliability estimation for each failure mode ${Kg1(1)(tg1)=g1lim,Kg2(1)(tg2)=g2limKg3(1)(tg3)=g3lim,Kg4(1)(tg4)=g4limRi=Pr(gi\u2264gilim)\u2248\Phi (wi)+\phi (wi)(1wi\u22121vi)$ |

Step 6: Analyze the correlation between multi-failure modes based on copula function |

Step 7: System reliability estimation for RV reducer considering the correlation of the four-failure mode $Rsys=1\u2212Pf,sys$ |

The Pearson linear correlation [52] between the four-failure modes is also analyzed as shown in Fig. 10, and the results from copula function is displayed in Fig. 11. Based on the correlation coefficients, the joint CDF is calculated as shown in Fig. 12.

Failure modes | $g1,g2$ | $g1,g3$ | $g1,g4$ | $g2,g3$ | $g2,g4$ | $g3,g4$ |
---|---|---|---|---|---|---|

Copula function | Gaussian copula | Gaussian copula | Clayton copula | Gaussian copula | Clayton copula | Frank copula |

Correlation coefficients | 0.99 | −0.463 | 0.081 | −0.4635 | 0.081 | 0.0923 |

Failure modes | $g1,g2$ | $g1,g3$ | $g1,g4$ | $g2,g3$ | $g2,g4$ | $g3,g4$ |
---|---|---|---|---|---|---|

Copula function | Gaussian copula | Gaussian copula | Clayton copula | Gaussian copula | Clayton copula | Frank copula |

Correlation coefficients | 0.99 | −0.463 | 0.081 | −0.4635 | 0.081 | 0.0923 |

From Figs. 10 and 11, it can be observed that the greater the fatigue and bending stress, the more significant their influence on the dynamic transmission error, which demonstrates a strong negative correlation existing between the fatigue life of the turning-arm bearing and the fatigue stress and bending stress of the cycloid pin gear.

The reliability of the four-failure modes is computed as $[R1=0.997,R2=0.937,R3=0.998,R4=0.917]$, which indicates that the bending strength and dynamic transmission accuracy are the weak links of the RV reducer. Based on the joint CDF, the probability of system failure can be calculated.

## 5 Conclusions

In this paper, a reliability assessment method is developed to deal with multi-failure modes of RV reducer by combining PCE and SPA method. The main conclusions are listed as follows:

The dynamic analysis under uncertainties for RV reducer is complex; to deal with this problem, dynamic simulation combining PCE method is developed to analyze the uncertainties propagation problem of RV reducer. Further, SPA method is employed to estimate the reliability of RV reducer with high accuracy.

The bending strength and dynamic transmission accuracy are the wear links of RV reducer. Therefore, further study will focus on the reliability design optimization of RV reducer to improve both bending strength and dynamic transmission accuracy.

The interaction mechanisms between the multi-failure modes are analyzed based on Pearson linear correlation method and copula function, and proper copula functions are accurately chosen based on saddlepoint approximated PDF and CDF, which provide an effective reliability method for complex system with multi-failure modes.

## Funding Data

National Natural Science Foundation of China (Contract No. 51875088; Funder ID: 10.13039/501100001809).

Fundamental Research Funds for the Central Universities (No. ZYGX2020ZB024; Funder ID: 10.13039/501100012226).

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## References

**104**