Nonlinear Dynamic Analysis and Optimization of Closed-Form Planetary Gear System

Abstract

A nonlinear purely rotational dynamic model of a multistage closed-form planetary gear set formed by two simple planetary stages is proposed in this study. The model includes time-varying mesh stiffness, excitation fluctuation and gear backlash nonlinearities. The nonlinear differential equations of motion are solved numerically using variable step-size Runge-Kutta. In order to obtain function expression of optimization objective, the nonlinear differential equations of motion are solved analytically using harmonic balance method (HBM). Based on the analytical solution of dynamic equations, the optimization mathematical model which aims at minimizing the vibration displacement of the low-speed carrier and the total mass of the gear transmission system is established. The optimization toolbox in MATLAB program is adopted to obtain the optimal solution. A case is studied to demonstrate the effectiveness of the dynamic model and the optimization method. The results show that the dynamic properties of the closed-form planetary gear transmission system have been improved and the total mass of the gear set has been decreased significantly.

1. Introduction

Planetary gear sets have been widely used in engineering including automotive transmissions, aviation transmissions, and crane gearboxes as well as other marine and industrial power transmission systems. Planetary gear trains have many advantages over fixed-center counter-shaft gear systems. The flow of power via multiple-gear meshes increases the power density and helps to reduce the overall size of the transmission train. The ability of multistage planetary sets in providing multiple speed reduction ratios has been the main reason for their extensive use in automatic transmission applications. Closed-form planetary trains are obtained from a number of single-stage differential planetary gear sets and one quasi-planetary stage whose central members are connected according to a given power flow configuration. Input, output, and fixed member assignments are made to certain central members to achieve a given gear ratio.

Because of complexity of structure, most of the earlier published studies on the planetary gear systems were confined to single-stage planetary. In addition, these early models were of linear time-invariant type, so that the eigen solutions and model summation techniques were used to predict the natural modes and the forced response [1–3]. Kahraman [4] employed a purely rotational dynamic model for all possible power flow configurations of complex compound planetary gear sets. He classified the natural model in two categories: asymmetric planet modes and axi-symmetric overall modes. Sun and Hu [5] investigated the frequency response of nonlinear planetary transmission system with multiple clearances using single-term harmonic balance method and focusing only on a single power flow configuration, in which the ring gear was fixed. Al-shyyab and Kahraman [6] developed a rotational single-stage nonlinear dynamic model of a simple planetary gear set and provided a semi-analytical forced response solution using multiterm HBM and showed that these HBM solutions in well agreement with direct numerical integration solution. Also, a recent study by Al-shyyab [7] investigated a compound planetary gear set formed by any number of simple planetary stages, and each planetary stage has a distinct fundamental mesh frequency and any number of planets spaced in any angular positions using multiterm HBM.

Those studies cited previously are mainly aimed at the modeling of planetary gear trains, the analysis of dynamic response as well as the analysis of parameters stability, and so on. Whereas, the studies on nonlinear dynamic optimization design of planetary transmission system with multiple clearances are still very limited. Zeng Bao [8] and Guan Wei [9] investigated multistage helical gear trains and single-stage planetary gear train taking the dynamic properties as objective functions, respectively. But, in their studies the design variables of optimization were limited to these parameters: the number of gears, the pitch-cycle helical-angle, and the modification coefficients. This study aims at providing numerical solutions and analytical solutions for the dynamic response of a closed-form planetary gear train having three planets spaced equally position angle using Runge-Kutta numerical integration and HBM, respectively. Based on the analytical solution, the optimization mathematical model that focused on minimizing the vibration acceleration of structure and the total mass of the gear transmission system is established. Some key design parameters such as the number of each gear, the module of each stage planetary, the transmission ratio of each stage planetary and the pressure angle are chosen as design varies. This study is available for the designing of closed-form planetary gear sets of both minimum weight and best dynamic characteristic for reference.

2. Dynamic Model of System

2.1. Model and Assumptions

The closed-form planetary gear train consists of a single-stage differential planetary (low-speed stage) and a single-stage quasi-planetary (high-speed stage, carrier is fixed) in this study, as shown in Figure 1. Each stage is comprised of three central elements: the sun gear (𝑠₁,𝑠₂), the ring gear (𝑟₁,𝑟₂), and carrier (𝑐₁,𝑐₂). Each stage planetary has 𝑛 planet gears. The parameter 𝑛 representing the number of planet gears is taken as 3 throughout this paper. The planets of each stage are free to rotate with respect to their common carrier. All the gears are mounted on their rigid shafts supported by rolling element bearings. The two rings are connected by torsional linear springs of stiffness 𝑘_𝑟1𝑟2, as shown in Figure 2. Likewise, the other central elements 𝑐₁, 𝑐₂ and 𝑠₁, 𝑠₂ are constrained by torsional linear springs of stiffness 𝑘_𝑐1𝑠2, 𝑘_𝑐2 and 𝑘_𝑠1, 𝑘_𝑐1𝑠2, respectively. In order to establish the mathematical model of the transmission system, a number of simplified assumptions are introduced in the case of speed reduction in the closed-form planetary gear set, as shown in Figure 1.

(a)

(b)

Figure 1 

(a) A closed-form planetary gear system; (b) A discrete model of closed-form planetary gear system.

Figure 2 

Lumped parameter dynamic model of single-stage planetary gear set.

(1)All of the gears in the set are assumed to be rigid and the flexibilities of each gear teeth at the gear mesh interface are modeled by an equivalent spring having time-varying stiffness acting along the mesh directions. These periodically time-varying mesh stiffnesses are subject to piece-linear functions representing gear backlashes.(2)Because the bending stiffness of shafts in the set is very large, the deflection of these shafts can be neglected. Thus, the transverse displacements of gears are not considered.(3)Because the damping mechanisms at the gear meshes and bearings of a planetary gear set are not easy to give a description of mathematical model, viscous gear mesh damping elements are introduced to represent energy dissipation of the transmission system.

2.2. Equivalent Displacements

In order to establish the equations of motion easily, all of torsional angular displacements are unified on the pressure line in terms of equivalent displacements. The equivalent transverse displacements in the mesh line direction caused by rotational displacements are written as follows: 𝑢_𝑐𝑗=𝑟_𝑐𝑗𝜃_𝑐𝑗,𝑢_𝑟𝑗=𝑟_𝑟𝑗𝜃_𝑟𝑗,𝑢_𝑠𝑗=𝑟_𝑠𝑗𝜃_𝑠𝑗,𝑢_𝑝𝑗𝑛=𝑟_𝑝𝑗𝜃_𝑗𝑛,𝑗=1,2;𝑛=1,2,3, (1) where 𝜃 and 𝑟 (subscripts 𝑐𝑗, 𝑟𝑗, 𝑠𝑗, 𝑗𝑛; 𝑗=ℎ,𝑙; 𝑛=1,2,3) are angular displacements and base circle radius of the parts, respectively. While 𝑟 (subscripts 𝑐1, 𝑐2) is the equivalent base circle radius of the carriers defined as follows: 𝑟_𝑠1+𝑟_𝑝ℎ=𝑟_𝑟1−𝑟_𝑝ℎ=𝑟_𝑐1𝑟cos𝛼,_𝑠2+𝑟_𝑝𝑙=𝑟_𝑟2−𝑟_𝑝𝑙=𝑟_𝑐2cos𝛼. (2)

With the symbol 𝑈 is used to represent the relative displacements in the direction of pressure line, the relative displacements are obtained according to the meshing relation and the equivalent displacements as follows. The positive direction of the relative displacements is assumed to be the same direction of the compressive deformation. The relative displacements of the parts in the pressure line direction are written as follows:𝑈_{𝑠1𝑝ℎ𝑖}=𝑟_𝑠1𝜃_𝑠1−𝜃_𝑐1−𝑟_𝑝ℎ𝜃_ℎ𝑖+𝜃_𝑐1=𝑟_𝑠1𝜃_𝑠1−𝑟_𝑐1𝜃_𝑐1cos𝛼−𝑟_𝑝ℎ𝜃_ℎ𝑖,𝑈_{𝑟1𝑝ℎ𝑖}=𝑟_𝑝ℎ𝜃_ℎ𝑖+𝜃_𝑐1−𝑟_𝑟1𝜃_𝑟1+𝜃_𝑐1=𝑟_𝑝ℎ𝜃_ℎ𝑖−𝑟_𝑐1𝜃_𝑐1cos𝛼−𝑟_𝑟1𝜃_𝑟1,𝑈_𝑐1𝑠2=𝑟_𝑐1𝜃_𝑐1−𝜃_𝑠2𝑈cos𝛼,_𝑐2=𝑟_𝑐2𝜃_𝑐2𝑈cos𝛼,_{𝑠2𝑝𝑙𝑖}=𝑟_𝑠2𝜃_𝑠2−𝜃_𝑐2−𝑟_𝑝𝑙𝜃_𝑙𝑖+𝜃_𝑐2=𝑟_𝑠2𝜃_𝑠2−𝑟_𝑐2𝜃_𝑐2cos𝛼−𝑟_𝑝𝑙𝜃_𝑙𝑖,𝑈_{𝑟2𝑝𝑙𝑖}=𝑟_𝑝𝑙𝜃_𝑙𝑖+𝜃_𝑐2−𝑟_𝑟2𝜃_𝑟2+𝜃_𝑐2=𝑟_𝑝𝑙𝜃_𝑙𝑖−𝑟_𝑐2𝜃_𝑐2cos𝛼−𝑟_𝑟2𝜃_𝑟2,𝑈_𝑟1𝑟2=𝑟_𝑟1𝜃_𝑟1−𝜃_𝑟2. (3a)(3b)(3c)(3d)(3e)(3f)(3g)

2.3. Equations of Motion

The closed-form planetary gear system consists of four different kinds of gear pairs, the external gear pair, that is, the sun gear/planet gear-𝑖 pair (subscripts 𝑠1, 𝑝ℎ𝑖 and 𝑠2, 𝑝𝑙𝑖), and the internal gear pair, that is, the ring gear/planet gear-𝑖 pair (subscripts 𝑟1, 𝑝ℎ𝑖 and 𝑟2, 𝑝𝑙𝑖). The mesh of gear 𝑗 (𝑠𝑖 or 𝑟𝑖, 𝑖=1,2) with a planet 𝑝𝑖 (𝑖=ℎ,𝑙) is represented by a periodically time-varying stiffness element 𝑘_𝑗𝑝𝑖(𝑡) subjected to a piecewise linear backlash function 𝑔 that includes a clearance of gap width 2𝑏_𝑗𝑝𝑖. Accordingly, the dynamic model of a closed-form gear set with 𝑛 planets includes 2𝑛 clearances. In this closed-form planetary gear system, damper is described by constant viscous damper coefficient 𝑐_𝑗𝑝𝑖. This is a rather simplified mesh contact model; in reality these contacts are subjected to the hydro-elastic-dynamic regime of lubrication [10]. In this paper, it is supposed that all planets 𝑝𝑖 (𝑖=ℎ,𝑙) and their respective meshes with gear 𝑗 (𝑠1, 𝑠2 or 𝑟1, 𝑟2) are identical so that 𝑘_𝑗𝑝𝑖(𝑡), 𝑏_𝑗𝑝𝑖, and 𝑐_𝑗𝑝𝑖 are the same for each 𝑗𝑝𝑖 mesh, except the phase angles of 𝑘_𝑗𝑝𝑖(𝑡) which differ according to planet phasing conditions.

Thus, the nonlinear dynamic differential equations of motion of closed-form planetary gear set can be established using the Lagrange principle as follows:𝑀_𝑐1𝑢̈_𝑐1+𝐶_𝑐1𝑠2𝑈̇_𝑐1𝑠2−³_𝑖=1𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑖}+𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑖}+𝑘_𝑐1𝑠2𝑈_𝑐1𝑠2−³_𝑖=1𝑘_𝑠1𝑛𝑈_{𝑠1𝑝ℎ𝑖}+𝑘_𝑟1𝑛𝑈_{𝑟1𝑝ℎ𝑖}𝑀=0,_𝑟1𝑢̈_𝑟1+𝐶_𝑟1𝑟2𝑈̇_𝑟1𝑟2−³_𝑖=1𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑖}+𝑘_𝑟1𝑟2𝑈_𝑟1𝑟2−³_𝑖=1𝑘_𝑟1𝑛𝑈_{𝑟1𝑝ℎ𝑖}𝑀=0,_𝑠1𝑢̈_𝑠1+³_𝑖=1𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑖}+³_𝑖=1𝑘_𝑠1𝑛𝑈_{𝑠1𝑝ℎ𝑖}=𝑇_in𝑟_𝑠1,𝑀_𝑝ℎ𝑢̈_ℎ𝑛−𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑛}+𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑛}−𝑘_𝑠1𝑛𝑈_{𝑠1𝑝ℎ𝑛}+𝑘_𝑟1𝑛𝑈_{𝑟1𝑝ℎ𝑛}(),𝑀=0𝑛=1,2,3_𝑐2𝑢̈_𝑐2+𝐶_𝑐2𝑈̇_𝑐2−³_𝑖=1𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑖}+𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑖}+𝑘_𝑐2𝑈_𝑐2−³_𝑖=1𝑘_𝑠2𝑛𝑈_{𝑠2𝑝𝑙𝑖}+𝑘_𝑟2𝑛𝑈_{𝑟2𝑝𝑙𝑖}𝑀=0,_𝑟2𝑢̈_𝑟2−𝐶_𝑟1𝑟2𝑈̇_𝑟1𝑟2𝑟_𝑟1𝑟_𝑟2−³_𝑖=1𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑖}−𝑘_𝑟1𝑟2𝑈_𝑟1𝑟2𝑟_𝑟1𝑟_𝑟2−³_𝑖=1𝑘_𝑟2𝑛𝑈_{𝑟2𝑝𝑙𝑖}=𝑇_out𝑟_𝑟2,𝑀_𝑠2𝑢̈_𝑠2−𝐶_𝑐1𝑠2𝑈̇_𝑐1𝑠2𝑟_𝑐1cos𝛼𝑟_𝑠2+³_𝑖=1𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑖}−𝑘_𝑐1𝑠2𝑈_𝑐1𝑠2𝑟_𝑐1cos𝛼𝑟_𝑠2+³_𝑖=1𝑘_𝑠2𝑛𝑈_{𝑠2𝑝𝑙𝑖}𝑀=0,_𝑝𝑙𝑢̈_𝑙𝑛−𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑛}+𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑛}−𝑘_𝑠2𝑛𝑈_{𝑠2𝑝𝑙𝑛}+𝑘_𝑟2𝑛𝑈_{𝑟2𝑝𝑙𝑛}=0(𝑛=1,2,3), (4a)(4b)(4c)(4d)(4e)(4f)(4g)(4h)where, 𝐼_𝑐1=𝐼^′_𝑐1+3𝑚_𝑝ℎ𝑟²_𝑐1 is the equivalent mass moment of inertia of the carrier including planet gears in high-speed stage; 𝐼^′_𝑐1 is the inertia of the carrier in high-speed stage; 𝑚_𝑝ℎ is the actual mass of planet-gear in high-speed stage; and 𝑟_𝑐1 is the distribution circle radius of the planet gears. 𝐼_𝑐1=𝑀_𝑐1𝑟²_𝑏𝑐1, 𝑀_𝑐1 is the equivalent mass of the carrier in high-speed stage with respect to the its equivalent base circle radius; 𝑟_𝑏𝑐1 is equivalent base circle radius of the carrier in high-speed stage; 𝑟_𝑏𝑐1=𝑟_𝑐1cos𝛼; 𝐼_𝑟1=𝑀_𝑟1𝑟²_𝑟1,  and 𝑀_𝑟1 is the equivalent mass of the ring in low-speed stage with respect to its base circle radius; 𝐼_𝑠1=𝑀_𝑠1𝑟²_𝑠1, and 𝑀_𝑠1 is the equivalent mass of the sun gear in low-speed stage with respect to its base circle radius; 𝐼_𝑝ℎ𝑛=𝑀_𝑝ℎ𝑟²_𝑝ℎ, and 𝑀_𝑝ℎ is the equivalent mass of the planet-gear in high-speed stage with respect to its equivalent base circle radius; 𝐼_𝑐2=𝐼^′_𝑐2+3𝑚_𝑝𝑙𝑟²_𝑐2, and 𝐼^′_𝑐2 is the moment of inertia of the carrier in low-speed stage.

Dynamic model of motion of transmission has considerable difficulties in its solution procedure as follows. (1) As a semidefinite system, its first-order natural frequency is zero corresponding to a rigid body motion. (2) The piecewise-linear function 𝑔 represents the gear backlash, and the number of variables is even different according to the external and internal gear pairs. (3) As both linear and nonlinear restoring forces exist in the equations, it is not possible to write out the governing equation in matrix form, while a general solution technique applicable to the systems of multiple degrees of freedom must be based on matrix form.

Therefore, (4a)–(4h) are simplified further by introducing a set of new variables:𝑈_{𝑠1𝑝ℎ𝑛}=𝑢_𝑠1−𝑢_𝑐1−𝑢_ℎ𝑛,𝑈_{𝑟1𝑝ℎ𝑛}=𝑢_ℎ𝑛−𝑢_𝑐1−𝑢_𝑟1,𝑈_𝑐1𝑠2=𝑢_𝑐1−𝑟_𝑐1cos𝛼𝑟_𝑠2𝑢_𝑠2,𝑈_𝑐2=𝑢_𝑐2,𝑈_{𝑠2𝑝ℎ𝑛}=𝑢_𝑠2−𝑢_𝑐2−𝑢_𝑙𝑛,𝑈_{𝑟2𝑝ℎ𝑛}=𝑢_𝑙𝑛−𝑢_𝑐2−𝑢_𝑟2,𝑈_𝑟1𝑟2=𝑢_𝑟1−𝑟_𝑟1𝑟_𝑟2𝑢_𝑟2. (5a)(5b)(5c)(5d)(5e)(5f)(5g)

The new coordinate variables defined previously not only have intuitional physical meaning, but also eliminate the rigid body motion. Furthermore, the piecewise-linear backlash function 𝑔 can be written as a set of functions with a single variable according to (5a)–(5g).

In addition, nondimensional parameters of (4a)–(4h) can be obtained by a characteristic length 𝑏 and frequency 𝜔_𝑛=𝑘_𝑠1𝑛/𝑀_𝑠1, such that 𝑡=𝜔_𝑛𝑡. (6)

Hence, a set of simplified equations of motions of the transmission system is obtained by substituting (5a)–(5g) and (6) into (4a)–(4h):𝑈̈_{𝑠1𝑝ℎ𝑛}+1𝜔²_𝑛³_𝑖=1𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑖}𝑀_𝑠1𝑐1+𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑖}𝑀_𝑐1+1𝜔²_𝑛𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑛}−𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑛}𝑀_𝑝ℎ−𝐶_𝑐1𝑠2𝑈̇_𝑐1𝑠2𝑀_𝑐1+1𝜔²_𝑛³_𝑖=1𝑘_𝑠1𝑛𝑔𝑈_{𝑠1𝑝ℎ𝑖}𝑀_𝑠1𝑐1+𝑘_𝑟1𝑛𝑔𝑈_{𝑟1𝑝ℎ𝑖}𝑀_𝑐1+1𝜔²_𝑛𝑘_𝑠1𝑛𝑔𝑈_{𝑠1𝑝ℎ𝑛}−𝑘_𝑟1𝑛𝑔𝑈_{𝑟1𝑝ℎ𝑛}𝑀_𝑝ℎ−𝑘_𝑐1𝑠2𝑔𝑈_𝑐1𝑠2𝑀_𝑐1=𝑇_in𝑟_𝑠1𝑀_𝑠1𝑏𝜔²_𝑛(),𝑈̈𝑛=1,2,3_{𝑟1𝑝ℎ𝑛}+1𝜔²_𝑛𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑛}−𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑛}𝑀_𝑝ℎ−𝐶_𝑐1𝑠2𝑈̇_𝑐1𝑠2𝑀_𝑐1−𝐶_𝑟1𝑟2𝑈̇_𝑟1𝑟2𝑀_𝑟1+1𝜔²_𝑛³_𝑖=1𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑖}𝑀_𝑐1+𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑖}𝑀_𝑐1𝑟1+1𝜔²_𝑛³_𝑖=1𝑘_𝑠1𝑛𝑔𝑈_{𝑠1𝑝ℎ𝑖}𝑀_𝑐1+𝑘_𝑟1𝑛𝑔𝑈_{𝑟1𝑝ℎ𝑖}𝑀_𝑐1𝑟1+1𝜔²_𝑛𝑘_𝑟1𝑛𝑔𝑈_{𝑟1𝑝ℎ𝑛}−𝑘_𝑠1𝑛𝑔𝑈_{𝑠1𝑝ℎ𝑛}𝑀_𝑝ℎ−𝑘_𝑐1𝑠2𝑔𝑈_𝑐1𝑠2𝑀_𝑐1−𝑘_𝑟1𝑟2𝑔𝑈_𝑟1𝑟2𝑀_𝑟1(),𝑈̈=0𝑛=1,2,3_𝑐1𝑠2+𝐶_𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛𝑈̇_𝑐1𝑠2−1𝜔²_𝑛³_𝑖=1𝐶_𝑠1𝑛𝑈̇_{𝑠1𝑝ℎ𝑖}+𝐶_𝑟1𝑛𝑈̇_{𝑟1𝑝ℎ𝑖}𝑀_𝑐1+𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑖}𝑟_𝑐1cos𝛼𝑀_𝑠2𝑟_𝑠2+𝐶_𝑐1𝑠2𝑀_𝑠2𝜔²_𝑛𝑟_𝑐1cos𝛼𝑟_𝑠2²𝑈̇_𝑐1𝑠2+𝑘_𝑐1𝑠2𝑔𝑈_𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛+𝑘_𝑐1𝑠2𝑔𝑈_𝑐1𝑠2𝑀_𝑠2𝜔²_𝑛𝑟_𝑐1cos𝛼𝑟_𝑠2²−1𝜔²_𝑛³_𝑖=1𝑘_𝑠1𝑛𝑔𝑈_{𝑠1𝑝ℎ𝑖}+𝑘_𝑟1𝑛𝑔𝑈_{𝑟1𝑝ℎ𝑖}𝑀_𝑐1+𝑘_𝑠2𝑛𝑔𝑈_{𝑠2𝑝𝑙𝑖}𝑟_𝑐1cos𝛼𝑀_𝑠2𝑟_𝑠2𝑈̈=0,_𝑐2+𝐶_𝑐2𝑈̇_𝑐2𝑀_𝑐2𝜔²_𝑛−³_𝑖=1𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑖}+𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑖}𝑀_𝑐2𝜔²_𝑛+𝑘_𝑐2𝑔𝑈_𝑐2𝑀_𝑐2𝜔²_𝑛−³_𝑖=1𝑘_𝑠2𝑛𝑔𝑈_{𝑠2𝑝𝑙𝑖}+𝑘_𝑟2𝑛𝑔𝑈_{𝑟2𝑝𝑙𝑖}𝑀_𝑐2𝜔²_𝑛𝑈̈=0,_{𝑠2𝑝𝑙𝑛}+1𝜔²_𝑛³_𝑖=1𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑖}𝑀_𝑠2+𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑖}+𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑖}𝑀_𝑐2+1𝜔²_𝑛𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑛}−𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑛}𝑀_𝑝𝑙−𝐶_𝑐2𝑈̇_𝑐2𝑀_𝑐2−𝐶_𝑐1𝑠2𝑀_𝑠2𝑟_𝑐1cos𝛼𝑈̇_𝑐1𝑠2𝑟_𝑠2+1𝜔²_𝑛³_𝑖=1𝑘_𝑠2𝑛𝑔𝑈_{𝑠2𝑝𝑙𝑖}𝑀_𝑠2+𝑘_𝑠2𝑛𝑔𝑈_{𝑠2𝑝𝑙𝑖}+𝑘_𝑟2𝑛𝑔𝑈_{𝑟2𝑝𝑙𝑖}𝑀_𝑐2+1𝜔²_𝑛𝑘_𝑠2𝑛𝑔𝑈_{𝑠2𝑝𝑙𝑛}−𝑘_𝑟2𝑛𝑔𝑈_{𝑟2𝑝𝑙𝑛}𝑀_𝑝𝑙−𝑘_𝑐2𝑔𝑈_𝑐2𝑀_𝑐2−𝑘_𝑐1𝑠2𝑀_𝑠2×𝑔𝑈_𝑐1𝑠2𝑟_𝑐1cos𝛼𝑟_𝑠2(),𝑈̈=0𝑛=1,2,3_{𝑟2𝑝𝑙𝑛}+1𝜔²_𝑛𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑛}−𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑛}𝑀_𝑝𝑙−𝐶_𝑐2𝑈̇_𝑐2𝑀_𝑐2+𝐶_𝑟1𝑟2𝑟_𝑟1𝑈̇_𝑟1𝑟2𝑀_𝑟2𝑟_𝑟2+1𝜔²_𝑛³_𝑖=1𝐶_𝑠2𝑛𝑈̇_{𝑠2𝑝𝑙𝑖}+𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑖}𝑀_𝑐2+𝐶_𝑟2𝑛𝑈̇_{𝑟2𝑝𝑙𝑛}𝑀_𝑟2+1𝜔²_𝑛𝑘_𝑟2𝑛𝑔𝑈_{𝑟2𝑝𝑙𝑛}−𝑘_𝑠2𝑛𝑔𝑈_{𝑠2𝑝𝑙𝑛}𝑀_𝑝𝑙−𝑘_𝑐2𝑔𝑈_𝑐2𝑀_𝑐2+𝑘_𝑟1𝑟2𝑟_𝑟1𝑔𝑈_𝑟1𝑟2𝑀_𝑟2𝑟_𝑟2+1𝜔²_𝑛³_𝑖=1𝑘_𝑠2𝑛𝑔𝑈_{𝑠2𝑝𝑙𝑖}+𝑘_𝑟2𝑛𝑔𝑈_{𝑟2𝑝𝑙𝑖}𝑀_𝑐2+𝑘_𝑟2𝑛𝑔𝑈_{𝑟2𝑝𝑙𝑖}𝑀_𝑟2=−𝑇_out𝑀_𝑟2𝑟_𝑟2𝑏𝜔²_𝑛(),𝑈̈𝑛=1,2,3_𝑟1𝑟2+𝐶_𝑟1𝑟2𝑈̇_𝑟1𝑟2𝑀_𝑟1𝑟2𝜔²_𝑛+1𝜔²_𝑛³_𝑖=1𝐶_𝑟2𝑛𝑟_𝑟1𝑀_𝑟2𝑟_𝑟2𝑈̇_{𝑟2𝑝𝑙𝑖}−𝐶_𝑟1𝑛𝑀_𝑟1𝑈̇_{𝑟1𝑝ℎ𝑖}+𝑘_𝑟1𝑟2𝑔𝑈_𝑟1𝑟2𝑀_𝑟1𝑟2𝜔²_𝑛+1𝜔²_𝑛³_𝑖=1𝑘_𝑟2𝑛𝑟_𝑟1𝑔𝑈_{𝑟2𝑝𝑙𝑖}𝑀_𝑟2𝑟_𝑟2−𝑘_𝑟1𝑛𝑔𝑈_{𝑟1𝑝ℎ𝑖}𝑀_𝑟1𝑇=−^out𝑟_𝑟1𝑀_𝑟2𝑟²_𝑟2𝑏𝜔²_𝑛, (7a)(7b)(7c)(7d)(7e)(7f)(7g)where, 1𝑀_𝑠1𝑐1=1𝑀_𝑠1+1𝑀_𝑐1,1𝑀_𝑐1𝑟1=1𝑀_𝑐1+1𝑀_𝑟1,1𝑀_𝑟1𝑟2=1𝑀_𝑟1+𝑟_𝑟1𝑀_𝑟2𝑟_𝑟2. (8) Equations (7a)–(7g) can be written in matrix form as 𝐔̈+𝐂𝐔̇+𝐊𝑔(𝐔)=𝐓, (9) where, the piecewise-linear backlash function is defined as 𝑔𝑈_𝑖𝑗(𝑡=𝑈)^𝑖𝑗(𝑡)−𝑏_𝑖𝑗,𝑈_𝑖𝑗>𝑏_𝑖𝑗,0,−𝑏_𝑖𝑗≤𝑈_𝑖𝑗≤𝑏_𝑖𝑗,𝑈_𝑖𝑗(𝑡)+𝑏_𝑖𝑗,𝑈_𝑖𝑗<−𝑏_𝑖𝑗. (10)

3. Solution of Dynamic Equations

3.1. Dynamic Response Using the Numerical Integration

In this section, the mathematical model will be solved numerically by the variable step-size Runge-Kutta integration method firstly. The parameters of the closed-form planetary gear set shown in Figure 1 are given in Tables 1 and 2. In this work, the values of the gear mesh damping coefficients 𝐶_𝑖𝑗 are assumed to be constants as 0.01 [11].

Table 1 

Geometric parameters of the closed-form planetary gear set.

Table 2 

Physical parameters of the closed-form planetary gear set.

By resolving the dynamic differential equations of motions, the dynamic responses (displacement and speed) of low-stage carrier are gained as shown in Figure 3.

(a) Displacement of low-speed stage carrier

(b) Speed of low-speed stage carrier

Figure 3 

Response of low-speed stage carrier without optimization.

3.2. The Acquisition of Analytical Solutions

In order to establish the optimization mathematical model of the dynamic behavior of the closed-form planetary gear set, it is necessary to obtain the analytical solutions of model using the harmonic balance method. For limiting the number of algebraic balance equations, only the fundamental frequency harmonic of the mesh stiffness functions are considered in this case. Similarly, external torque functions is also considered to be in the form of mesh stiffness functions. And then, attention is paid to the periodic vibrations of system under the harmonic excitation. The procedure of solving (7a)–(7g) using HBM [5] includes some aspects as follows.

(1) According to the assumption previously mentioned, the external excitations that represent fundamental frequency pulsations can be written in the form 𝐓=𝐓_𝑚𝑖+𝐓_𝑎𝑖cosΩ𝑡+𝜙_𝑖, (11) where 𝐓_𝑚𝑖 is the mean component of torque and 𝐓_𝑎𝑖 is the amplitude of the alternating component of the fundamental frequency mesh force. 𝜙_𝑖 is the phase angle.

(2) According to the harmonic excitations given in (11), the harmonic balance method solution to (7a)–(7g) is assumed in the same form 𝐔=𝐔_𝑚𝑖+𝐔_𝑎𝑖cosΩ𝑡+𝜙_𝑖, (12) where 𝐔_𝑚𝑖 and 𝐔_𝑎𝑖 are the mean and alternating components of the steady state response, respectively, and 𝜙_𝑖 is the phase angle.

(3) For relative mesh displacements, the piecewise-linear function in (4a)–(4h) can be written in a unified form: 𝑔𝑈_𝑖=𝑁_𝑚𝑖𝑈_𝑚𝑖+𝑁_𝑎𝑖𝑈_𝑎𝑖cosΩ𝑡+𝜙_𝑖, (13) where𝑁_𝑚𝑖𝑞=1+^𝑎𝑖2𝑞_𝑚𝑖𝐺𝜇₊−𝐺𝜇₋,𝑁_𝑎𝑖1=1−2𝐻𝜇₊−𝐻𝜇₋. (14) The expressions of 𝐺 and 𝐻 are given in the Appendix.

(4) Considering the mean value and the fundamental harmonic value of periodically time-varying mesh stiffness in Figure 2, the elements of stiffness matrix 𝐊 in (9) can be written as 𝑘_𝑖𝑗=𝑘_𝑚𝑖𝑗+𝑘_𝑎𝑖𝑗cosΩ𝜏+𝜑_𝑖𝑗. (15) So, the stiffness matrix 𝐊 is written in terms of two separate matrices for mean stiffness and alternating stiffness as 𝐊=𝐊_𝑚+Δ𝐊, (16) where 𝐊_𝑚=𝑘_{𝑚𝑖𝑗𝑛×𝑛},Δ𝐊=𝑘_𝑎𝑖𝑗cosΩ𝜏+𝜑_{𝑖𝑗𝑛×𝑛}. (17) By substituting (9)–(14) into (7a)–(7g) and balancing the like harmonic terms, the algebraic equations of system can be gained 𝐊_𝑚𝐲_𝑚+𝐊₁𝐲₃2+𝐊₂𝐲₄2−𝐓_𝑚={0}_𝑛×1,𝐊_𝑚𝐲₃+𝐊₁𝐲_𝑚−Ω²𝐲₁−Ω𝐂𝐲₂−𝐓₁={0}_𝑛×1,𝐊_𝑚𝐲₄+𝐊₂𝐲_𝑚−Ω²𝐲₂+𝐂Ω𝐲₁−𝐓₂={0}_𝑛×1, (18) where  𝐲₁={𝑞_𝑎𝑖cos𝜑_𝑖}_𝑛×1, 𝐲₂={𝑞_𝑎𝑖sin𝜑_𝑖}_𝑛×1, 𝐲₃={𝑁_𝑎𝑖𝑞_𝑎𝑖cos𝜑_𝑖}_𝑛×1, 𝐲₄={𝑁_𝑎𝑖𝑞_𝑎𝑖sin𝜑_𝑖}_𝑛×1, 𝐲_𝑚={𝑁_𝑚𝑖𝑞_𝑚𝑖}_𝑛×1, 𝐓₁={𝑇_𝑎𝑖cos𝜙_𝑖}_𝑛×1, 𝐓₂={𝑇_𝑎𝑖sin𝜙_𝑖}_𝑛×1, 𝐊₁={𝑘_𝑎𝑖𝑗cos𝜙_𝑖𝑗}_𝑛×1, and 𝐊₂={𝑘_𝑎𝑖𝑗sin𝜙_𝑖𝑗}_𝑛×1.

The matrices 𝐤_𝑚 and 𝐂 are given in the Appendix.

4. Optimization of the Transmission System

4.1. Variables and Objective Function of Optimization

A number of key design parameters have great influence on the dynamic characteristics of the gear transmissions, so they can be chosen as design variables, such as number of each gear, module, gear width, and pressure angle. To simplify optimization process, the module of each stage 𝑚₁, 𝑚₂, number of sun gears of each stage 𝑧_𝑠1, 𝑧_𝑠2, and pressure angle of each stage 𝛼₁,𝛼₂ are chosen as optimization variables in this study. So, the recurrence variables vector for this optimization procedure can be written as 𝑥=𝑚₁,𝑚₂,𝑧_𝑠1,𝑧_𝑠2,𝛼₁,𝛼₂^𝑇. (19)

The dynamic characteristics standards of the gear transmission include maximum dynamic loads, dynamic load factor, stiffness, and displacement/velocity/acceleration of vibration, each of these standards can be chosen as optimization objective. Considering that the rotation center of the planets in low-stage is unfixed, and reducing the total weight of the transmission system simultaneously, the displacement of rotational vibration of the low-stage carrier and total mass of the gear set are chosen as optimization objective. According to the basic idea of the multiobjective optimization, two optimization objective functions unified using the normalized weighting method [11] can be written in form 𝑓=𝜆₁𝑓₁+𝜆₂𝑓₂, (20) where 𝜆₁ and 𝜆₂ are weighting coefficient, here, 𝜆₁=0.3 and 𝜆₂=0.7 [12]. 𝑓₁ and 𝑓₂ are the ampler of the rotational vibration displacement of the low-stage carrier and total mass of the gear set, respectively. Due to limited space, the expressions of 𝑓₁ and 𝑓₂ are not given in detail.

4.2. Constraints

(1)Distributing the number of teeth of each gear: concentric conditions, adjacent conditions, and assembly conditions;(2)contact fatigue strength and bending fatigue strength constraints 𝑆_𝐻(𝑖)≥𝑆_𝐻min,𝑆_𝐹(𝑖)≥𝑆_𝐹min,𝑖=𝑠₁,𝑠₂,𝑟₁,𝑟₂,𝑝ℎ,𝑝𝑙, (21)(3)contact ratio constraints 𝜀_𝑖𝑗≥1.5,𝑖=𝑠₁,𝑠₂,𝑟₁,𝑟₂,𝑗=𝑝ℎ,𝑝𝑙, (22)(4)transmission ratio without loop-power conditions −32≤𝑖₁+𝑖₂−𝑖₁𝑖₂𝑖≤−30,0<²−𝑖₁𝑖₂𝑖₁+𝑖₂−𝑖₁𝑖₂<1, (23) where 𝑖_𝑛𝑍=−^𝑟𝑛𝑍_𝑠𝑛,𝑛=1,2, (24)(5)minimum tooth thickness constraints 𝑆_𝑎𝑖≥0.4𝑚_𝑗,𝑖=𝑠₁,𝑠₂,𝑟₁,𝑟₂,𝑝ℎ,𝑝𝑙;𝑗=1,2. (25)

5. Results and Discussions

Some more logical parameters of the closed-form planetary gear set are obtained according to optimization constraints conditions previously mentioned: 𝑚₁=2, 𝑚₂=3, 𝑍_𝑠1=17, 𝑍_𝑠2=18, 𝛼₁=𝛼₂=20^∘. Further, other parameters of transmission can be computed as shown in Table 3. Total mass of gear set is 63.5 kg. Comparing the results between before and after the optimization, it is clear that transmission system reduces the total mass of 8.6% and the dynamic characteristic significantly improves as well. It is also evident from the comparison of Figure 3 with Figure 4 that the displacement of each structure significantly reduces at the same time. It is obvious that the scope of the amplitude of displacement of low-stage carrier is −3 to 3 in Figure 4(a), while that in Figure 3(a) is −5 to 3.

Table 3 

Parameters of structure after optimization.

(a)  Displacement of low-speed stage carrier

(b) Speed of low-speed stage carrier

Figure 4 

Response of low-speed stage carrier after optimization.

6. Conclusions

Considering the gear backlash, time-varying mesh stiffness and excitation fluctuation, and so forth, a discrete nonlinear dynamic model of a two stages closed-form planetary set was proposed in this study. In order to facilitate the analysis and comparison of system response between before and after optimization, nonlinear differential equations of motion of the dynamic model were solved using a Runge-Kutta numerical integration method. For optimization of transmission system, the analytical solutions were obtained. The total optimization objective function is obtained using the normalized weights method.

A case shows that the total mass of transmission system has significantly reduced and dynamic characteristic has distinctly improved. Effectiveness of the dynamic model and the dynamic optimization mathematic model is demonstrated.

Appendix

Consider the following: 𝐺𝜇=2𝜋𝜇sin⁻¹𝜇+1−𝜇²,𝜇𝜇,𝜇𝐻𝜇=≤1,>1,−1,𝜇<−1,2𝜇𝜋sin⁻¹𝜇+1−𝜇²,𝜇≤1,1,𝜇>1. (A.1) The elements of the stiffness matrix 𝐤_𝑚 are given as follows 𝑘()=𝑘1,1^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1+1𝑀_𝑝ℎ,𝑘()=𝑘1,𝑖^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1𝑘()=𝑘,𝑖=2,3,1,4^𝑟1𝑛𝜔²_𝑛1𝑀_𝑐1−1𝑀_𝑝ℎ,𝑘()=𝑘1,𝑖^𝑟1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()𝑘,𝑖=5,6,1,7=−^𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛,𝑘()=𝑘2,𝑖^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1𝑘()𝑘,𝑖=1,3,2,7=−^𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛,𝑘()=𝑘2,2^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1+1𝑀_𝑝ℎ,𝑘()=𝑘2,𝑖^𝑟1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()=𝑘,𝑖=4,6,2,5^𝑟1𝑛𝜔²_𝑛1𝑀_𝑐1−1𝑀_𝑝ℎ,𝑘()=𝑘2,𝑖^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1𝑘()=𝑘,𝑖=1,3,2,2^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1+1𝑀_𝑝ℎ,𝑘()=𝑘2,𝑖^𝑟1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()=𝑘,𝑖=4,6,2,5^𝑟1𝑛𝜔²_𝑛1𝑀_𝑐1−1𝑀_𝑝ℎ,𝑘()𝑘2,7=−^𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛,𝑘()=𝑘3,𝑖^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1𝑘()=𝑘,𝑖=1,2,3,3^𝑠1𝑛𝜔²_𝑛1𝑀_𝑠1+1𝑀_𝑐1+1𝑀_𝑝ℎ,𝑘()=𝑘3,𝑖^𝑟1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()=𝑘,𝑖=4,5,3,6^𝑟1𝑛𝜔²_𝑛1𝑀_𝑐1−1𝑀_𝑝ℎ,𝑘()𝑘3,7=−^𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛,𝑘()=𝑘4,1^𝑠1𝑛𝜔²_𝑛1𝑀_𝑐1−1𝑀_𝑝ℎ,𝑘()=𝑘4,𝑖^𝑠1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()𝑘,𝑖=2,3,4,7=−^𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛,𝑘()=𝑘4,4^𝑟1𝑛𝜔²_𝑛1𝑀_𝑝ℎ+1𝑀_𝑐1+1𝑀_𝑟1,𝑘()=𝑘4,𝑖^𝑟1𝑛𝜔²_𝑛1𝑀_𝑐1+1𝑀_𝑟1𝑘()𝑘,𝑖=5,6,4,15=−^𝑟1𝑟2𝑀_𝑟1𝜔²_𝑛,𝑘()=𝑘5,𝑖^𝑠1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()=𝑘,𝑖=1,3,5,2^𝑠1𝑛𝜔²_𝑛1𝑀_𝑐1−1𝑀_𝑝ℎ,𝑘()=𝑘5,𝑖^𝑟1𝑛𝜔²_𝑛1𝑀_𝑐1+1𝑀_𝑟1𝑘()=𝑘,𝑖=5,6,5,5^𝑟1𝑛𝜔²_𝑛1𝑀_𝑝ℎ+1𝑀_𝑐1+1𝑀_𝑟1,𝑘()𝑘5,7=−^𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛,𝑘()𝑘5,15=−^𝑟1𝑟2𝑀_𝑟1𝜔²_𝑛,𝑘()=𝑘6,𝑖^𝑠1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()=𝑘,𝑖=1,2,6,3^𝑠1𝑛𝜔²_𝑛1𝑀_𝑐1−1𝑀_𝑝ℎ,𝑘()=𝑘6,𝑖^𝑟1𝑛𝜔²_𝑛1𝑀_𝑐1+1𝑀_𝑟1𝑘()=𝑘,𝑖=4,5,6,6^𝑟1𝑛𝜔²_𝑛1𝑀_𝑝ℎ+1𝑀_𝑐1+1𝑀_𝑟1,𝑘()𝑘6,7=−^𝑐1𝑠2𝑀_𝑐1𝜔²_𝑛,𝑘()𝑘6,15=−^𝑟1𝑟2𝑀_𝑟1𝜔²_𝑛,𝑘()𝑘7,𝑖=−^𝑠1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()𝑘,𝑖=1,2,3,7,𝑖=−^𝑟1𝑛𝑀_𝑐1𝜔²_𝑛𝑘()=𝑘,𝑖=4,5,6,7,7^𝑐1𝑠2𝜔²_𝑛1𝑀_𝑐1+1𝑀_𝑠2𝑟_𝑐1cos𝛼𝑟_𝑠2²,𝑘()𝑘7,𝑖=−^𝑠2𝑛𝑀_𝑠2𝜔²_𝑛𝑟_𝑐1cos𝛼𝑟_𝑠2𝑘()=𝑘,𝑖=9,10,11,8,8^𝑐2𝑀_𝑐2𝜔²_𝑛,𝑘()𝑘8,𝑖=−^𝑠2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()𝑘,𝑖=9,10,11,8,𝑖=−^𝑟2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()𝑘,𝑖=12,13,14,9,7=−^𝑐1𝑠2𝑀_𝑠2𝜔²_𝑛𝑟_𝑐1cos𝛼𝑟_𝑠2,𝑘()𝑘9,8=−^𝑐2𝑀_𝑐2𝜔²_𝑛,𝑘()=𝑘9,9^𝑠2𝑛𝜔²_𝑛1𝑀_𝑠2+1𝑀_𝑐2+1𝑀_𝑝𝑙,𝑘()=𝑘9,𝑖^𝑠2𝑛𝜔²_𝑛1𝑀_𝑠2+1𝑀_𝑐2𝑘()=𝑘,𝑖=10,11,9,12^𝑟2𝑛𝜔²_𝑛1𝑀_𝑐2−1𝑀_𝑝𝑙,𝑘()=𝑘9,𝑖^𝑟2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()𝑘,𝑖=13,14,10,7=−^𝑐1𝑠2𝑀_𝑠2𝜔²_𝑛𝑟_𝑐1cos𝛼𝑟_𝑠2,𝑘()𝑘10,8=−^𝑐2𝑀_𝑐2𝜔²_𝑛,𝑘()=𝑘10,𝑖^𝑠2𝑛𝜔²_𝑛1𝑀_𝑠2+1𝑀_𝑐2𝑘()=𝑘,𝑖=9,10,10,10^𝑠2𝑛𝜔²_𝑛1𝑀_𝑠2+1𝑀_𝑐2+1𝑀_𝑝𝑙,𝑘()=𝑘10,𝑖^𝑟2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()=𝑘,𝑖=12,14,10,13^𝑟2𝑛𝜔²_𝑛1𝑀_𝑐2−1𝑀_𝑝𝑙,𝑘()𝑘11,7=−^𝑐1𝑠2𝑀_𝑠2𝜔²_𝑛𝑟_𝑐1cos𝛼𝑟_𝑠2,𝑘()𝑘11,8=−^𝑐2𝑀_𝑐2𝜔²_𝑛,𝑘()=𝑘11,𝑖^𝑠2𝑛𝜔²_𝑛1𝑀_𝑠2+1𝑀_𝑐2𝑘()=𝑘,𝑖=9,10,11,11^𝑠2𝑛𝜔²_𝑛1𝑀_𝑠2+1𝑀_𝑐2+1𝑀_𝑝𝑙,𝑘()=𝑘11,𝑖^𝑟2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()=𝑘,𝑖=12,13,11,14^𝑟2𝑛𝜔²_𝑛1𝑀_𝑐2−1𝑀_𝑝𝑙,𝑘()𝑘12,8=−^𝑐2𝑀_𝑐2𝜔²_𝑛,𝑘()=𝑘12,15^𝑟1𝑟2𝑟_𝑟1𝑀_𝑟2𝜔²_𝑛𝑟_𝑟2,𝑘()=𝑘12,9^𝑠2𝑛𝜔²_𝑛1𝑀_𝑐2−1𝑀_𝑝𝑙,𝑘()=𝑘12,𝑖^𝑠2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()=𝑘,𝑖=10,11,12,12^𝑟2𝑛𝜔²_𝑛1𝑀_𝑝𝑙+1𝑀_𝑐2,𝑘()=𝑘12,𝑖^𝑟2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()𝑘,𝑖=13,14,13,8=−^𝑐2𝑀_𝑐2𝜔²_𝑛,𝑘()=𝑘13,10^𝑠2𝑛𝜔²_𝑛1𝑀_𝑐2−1𝑀_𝑝𝑙,𝑘()=𝑘13,𝑖^𝑠2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()=𝑘,𝑖=9,11,13,13^𝑟2𝑛𝜔²_𝑛1𝑀_𝑝𝑙+1𝑀_𝑐2,𝑘()=𝑘13,𝑖^𝑟2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()=𝑘,𝑖=12,14,13,15^𝑟1𝑟2𝑟_𝑟1𝑀_𝑟2𝜔²_𝑛𝑟_𝑟2,𝑘()𝑘14,8=−^𝑐2𝑀_𝑐2𝜔²_𝑛,𝑘()=𝑘14,11^𝑠2𝑛𝜔²_𝑛1𝑀_𝑐2−1𝑀_𝑝𝑙,𝑘()=𝑘14,𝑖^𝑠2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()=𝑘,𝑖=9,10,14,14^𝑟2𝑛𝜔²_𝑛1𝑀_𝑝𝑙+1𝑀_𝑐2,𝑘()=𝑘14,𝑖^𝑟2𝑛𝑀_𝑐2𝜔²_𝑛𝑘()=𝑘,𝑖=12,13,14,15^𝑟1𝑟2𝑟_𝑟1𝑀_𝑟2𝜔²_𝑛𝑟_𝑟2,𝑘()𝑘15,𝑖=−^𝑟1𝑛𝑀_𝑟1𝜔²_𝑛𝑘()=𝑘,𝑖=4,5,6,15,𝑖^𝑟2𝑛𝑟_𝑟1𝑀_𝑟2𝑟_𝑟2𝜔²_𝑛𝑘()=𝑘,𝑖=12,13,14,15,15^𝑟1𝑟2𝜔²_𝑛1𝑀_𝑟1+𝑟_𝑟1𝑀_𝑟2𝑟_𝑟2. (A.2) The elements of 𝐤_𝑚 not listed previously are taken as zeros. The matrix 𝐂 is in the similar form of 𝐤_𝑚 and no longer is listed in detail here.

Notations

HBM:	Harmonic balance method
𝑏:	Half of clearance (backlash)
𝑘:	Gear mesh stiffness
𝑟:	Ring gear
𝑔:	Discontinuous displacement function
𝐼:	Polar mass moment of inertia
𝑖:	Transmission ratio
𝑐:	Carrier
𝑚:	Actual mass
𝑡:	Dimensional time
𝑢:	Actual equivalent transverse displacement in the direction of pressure line
𝑈:	Relative equivalent transverse displacement in the direction of pressure line
𝑠:	Sun gear
𝑛𝑝:	Number of planet
𝐓:	Torque
𝜃∶	Angular displacement
𝐶:	Damping coefficient
𝑍:	Number of gear
𝑟:	Base cycle radius
𝑀:	Equivalent mass
𝑡:	Nondimensional time
𝛼:	Pressure angle.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (no. 51175299), Shandong Province Natural Science Foundation of China (ZR2010EM012) and Shandong University Independent Innovation Foundation of China (no. 2012TS048).