节点外啮合行星齿轮传动系统非线性动力学分析

唐鑫; 鲍和云; 陆凤霞; 靳广虎; TANG Xin; BAO Heyun; LU Fengxia; JIN Guanghu

Transactions of Nanjing University of Aeronautics & Astronautics

+Advanced Search

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Nonlinear Dynamic Analysis of Planetary Gear Train System with Meshing Beyond Pitch Point PDF

- ORCID：
TANG Xin ^1,2
- ORCID：
BAO Heyun ¹
- ORCID：
LU Fengxia ¹
- ORCID：
JIN Guanghu ¹

1. National Key Laboratory of Science and Technology on Helicopter Transmission， Nanjing University of Aeronautics and Astronautics， Nanjing 210016， P. R. China； 2. National Key Laboratory of Science and Technology on Helicopter Transmission， AECC Hunan Aviation Powerplant Research Institute， Zhuzhou 412002， P. R. China

CLC： TH132.4

Updated：2021-01-29

DOI：http://dx.doi.org/10.16356/j.1005-1120.2020.06.006

OUTLINE

Abstract

Nonlinear dynamic analysis was performed on a planetary gear transmission system with meshing beyond the pitch point. The parameters of the planetary gear system were optimized， and a two-dimensional nonlinear dynamic model was established using the lumped-mass method. Time-varying meshing stiffness was calculated by the energy method. The model consumes the backlash， bearing clearance， time-varying meshing stiffness， time-varying bearing stiffness， and time-varying friction coefficient. The time-varying bearing stiffness was calculated according to the Hertz contact theory. The load distribution among the gears was computed， and the time-varying friction coefficient was calculated according to elastohydrodynamic lubrication （EHL） theory. The dynamical equations were solved via numerical integration. The global bifurcation characteristics caused by the input speed， backlash， bearing clearance， and damping were analyzed. The system was in a chaotic state at natural frequencies or frequency multiplication. The system transitioned from a single-period state to a chaotic state with the increase of the backlash. The bearing clearance of the sun gear had little influence on the bifurcation characteristics. The amplitude was restrained in the chaotic state as the damping ratio increased.

Keywords

meshing beyond pitch point; planetary gear system; nonlinear; time-varying bearing stiffness; time-varying meshing stiffness; multiple clearances; bifurcation; time-varying friction coefficient

0 Introduction

As a standard meshing gear transmission system， the gear transmission meshes in the actual mesh zone AC and EC， as shown in Fig.1. In Fig.1， C is the gear mesh point， F_μ the tooth surface friction force， and arrow the direction of the friction force. In the standard meshing gear transmission system， F_μ will change the direction at the point C that may cause vibration of the system.

Fig.1 Tooth surface friction of a standard meshing gear transmission system

As a non-standard meshing gear transmission system， the gear transmission system only meshes in one side of the C point， as shown in Figs.2 and 3. In Fig.2， AC is the actual mesh zone. In Fig.3， CE is the actual mesh zone.

Fig.2 Tooth surface friction of a gear transmission system meshed upside of the pitch point

Fig.3 Tooth surface friction of a gear transmission system meshed underside of the pitch point

The non-standard meshing gear transmission system with meshing beyond the pitch point can avoid the change of the tooth surface friction force direction.

As a non-standard meshing gear transmission system， the gear transmission system with meshing beyond the pitch point can avoid the change of the tooth surface friction force and has been investigated by scholars in recent years. In 1997， Gao and Zhou^［

1］ defined the coefficient of the pitch point and proved the feasibility of the theory after optimizing the design of the parameters by using the modification coefficient as the design variable. In 2012， Liu^{［Reference 2

Baidu Scholar}2］ adopted an equal modulus and equal pressure angle to achieve meshing transmission beyond the pitch point by changing the gear-modification coefficient and the coefficient of tooth depth. And the resulting strength and system dynamic response were determined. In 2013， Sun^{［Reference 3

Baidu Scholar}3］ established the dynamic model of the planetary gear transmission system with meshing beyond the pitch point （PGTSMPP） by using non-adopted equal modulus and non-equal pressure three-degree of freedom （DOF） gears and planetary gears. Considering the influence of nonlinear factors such as the time-varying friction， time-varying stiffness， and angle， the system dynamic response and load-sharing coefficient were determined by using the numerical integration method while ignoring the translational vibration of the planets， carrier， and ring gear. In 2016， Bao et al.^{［Reference 4

Baidu Scholar}4］ analyzed the inherent characteristics and dynamic characteristics of the system under the influence of a flexible ring on the PGTSMPP.

In the process of planetary gear transmission， which is inevitably affected by nonlinear factors such as the backlash， bearing clearance， and single and double teeth alternately meshing， the system may be in a multi-period or even a chaotic state， aggravating the vibration and noise and affecting the stability of the system. Domestic and foreign scholars have performed substantial research on the nonlinear dynamics of gears. The scholars^［

5‑9］ studied the nonlinear dynamics and backlash of single-DOF gears and developed a dynamic model by using the lumped-mass method and finite-element method. The nonlinear dynamic characteristics were solved by using the harmonic balance method， inverse Fourier transform method， and Newton-Rapson method.

Taking the spur gear as the research object and considering the influence of the tooth surface friction， backlash， and time-varying meshing stiffness on the dynamic characteristics of the system， Wang et al.^［

10］ established a pure torsional two-DOF dynamic model in 2002. The time-domain graph， frequency-domain graph， phase diagram， Poincaré section， bifurcation diagram， and maximum Lyapunov exponent for the system were obtained， and the effects of the friction on the system periodic response， chaotic response， and bifurcation were analyzed comprehensively. In 2008， considering the relative sliding velocity and the effect of the single and double teeth alternately meshing， Tang et al.^{［Reference 11

Baidu Scholar}11］ adopted a period-expansion method to establish a nonlinear dynamic model in consideration of the friction of the tooth surface and the time-varying stiffness. In 2015， taking a single-stage planetary gear system with a non-equal modulus and non-equal pressure angle as the research object， Ye et al.^{［Reference 12

Baidu Scholar}12］ investigated the effects of the clearance， modulus， and pressure angle on the load-sharing coefficient in consideration of the backlash， comprehensive transmission error， and pressure angle. In the same year， considering the effects of the time-varying stiffness， tooth surface friction， backlash， and bearing clearance on a spur-gear system， Sheng et al.^{［Reference 13

Baidu Scholar}13］ solved the nonlinear dynamic equations by using the four-order Runge‑Kutta method.

The effects of the friction coefficient， damping ratio， and clearance on the bifurcation characteristics were examined by using the Poincaré section. In 2015， taking a higher-contract ratio planetary gear system as the research object， Li^［

14］ established a dynamic model by using the energy method in consideration of the effects of the tooth surface friction， backlash， bearing clearance， time-varying meshing stiffness， and comprehensive transmission error on the nonlinear characteristics of the system. The effects of the main parameters affecting the contact ratio on the dynamic load-sharing coefficient of the high-contact ratio planetary gear transmission system were analyzed， and the theoretical analysis was validated by experiments. Mo et al.^{［Reference 15-16}15-16］ investigated the load sharing characteristics of herringbone planetary gear train and multi-power face gear split flow system. Jin et al.^{［Reference 17

Baidu Scholar}17］ studied the effect of friction on dynamic response of a power split transmission system. Bao et al.^{［Reference 18

Baidu Scholar}18］ analyzed the dynamic of external gear system with meshing beyond pitch considering time-varing friction coefficient.

Thus， the dynamic response， inherent characteristics and load-sharing characteristics of the gear transmission system with meshing beyond the pitch point have been studied. However， domestic and foreign scholars mainly studied the nonlinear dynamics of the standard gear transmission system； nonlinear research on PGTSMPP has rarely been reported and is therefore the focus in this paper.

1 Parameters of System

According to Refs.［

3-4］， the parameters of PGTSMPP are optimized. The results are shown in Table 1.

Table 1 Coefficients of PGTSMPP

Parameter	Sun	Planet	Ring
Tooth number	27	25	81
Modulus/mm	2.63	2.63	2.57
Pressure angle/ (°)	25	25	22.3
Modification coefficient	0.623	0.775	1
Tooth width/mm	60	60	60
Accuracy	6	6	6

1.1 Calculation of phase angle

λ_sp_i is the meshing phase coefficient between the ith sun‑planet meshing pair and the first sun‑planet meshing pair， λ_rp_i the meshing phase coefficient between the ith ring‑planet meshing pair and the first ring‑planet meshing pair， and λ_sr the phase difference between the sun and the ring. The corresponding formulas can be written as

\{\begin{array}{l} λ_{s p i} = d e c (z_{s} φ_{i} / 2 π) \\ λ_{r p i} = d e c (z_{r} φ_{i} / 2 π) \\ λ_{s r} = d e c (F_{2}^{'} B_{r 2} / p_{p b}) \end{array}

(1)

where z_s and z_r are the numbers of teeth for the sun and ring gears， respectively； φ_i is the position angle of the planet， p_pb the base pitch of the planet， $F_{2}^{'}$ the point of the opposite tooth surface of the planet base circle relative to the meshing starting point of the sun， $B_{r 2}$ the meshing starting point of the ring gear， and dec the decimal part.

1.2 Calculation of friction arm

A schematic of the friction arm of the planet is shown in Fig.4.

Fig.4 Friction arm of PGTSMPP

In the model established in this paper， the external meshing involves the standard gear meshing pairs， and the inner meshing involves the meshing pairs after the pitch point. The actual meshing line B_pr_{i_}₁B_pr_{i_}₂ is located on the side of the pitch point， and the friction arms at any moment according to the geometric relationship are given as follows

\{\begin{array}{l} h_{s p i_1} = N_{s p i_1} B_{s p i_2} + ω_{s c} r_{b s} m o d (t_{s p i}, T) & m o d (t_{s p i}, T) \leq (ε_{s p} - 1) * T \\ h_{s p i_2} = N_{s p i_1} B_{s p i_2} + ω_{s c} r_{b s} m o d (t_{s p i}, T) + p_{b} & m o d (t_{s p i}, T) \leq (ε_{s p} - 1) * T \\ h_{p s i_1} = N_{s p i_1} N_{s p i_2} - h_{s p i_1} \\ h_{p s i_2} = N_{s p i_1} N_{s p i_2} - h_{s p i_2} \end{array}

(2)

where h_sp_{i_}₁ and h_sp_{i_}₂ are the friction arms of the sun in the sun‑planet meshing pairs， and h_ps_{i_}₁ and h_ps_{i_}₂ the friction arms of the planet in the sun-planet meshing pairs； ε_sp is the contact ratio of the sun-planet meshing pairs， ω_sc the relative angular velocity of the sun relative to the carrier， p_b the base pitch； and N_sp_{i_}₁N_sp_{i_}₂ the theoretical meshing line of the sun-planet meshing pairs.

\{\begin{array}{l} h_{p r i_1} = N_{p r i_1} B_{p r i_2} + ω_{p c} r_{b p} m o d (t_{r p i}, T) & m o d (t_{r p i}, T) \leq (ε_{p r} - 1) * T \\ h_{p r i_2} = N_{p r i_1} B_{p r i_2} + ω_{p c} r_{b p} m o d (t_{r p i}, T) + p_{b} & m o d (t_{r p i}, T) \leq (ε_{p r} - 1) * T \\ h_{r p i_1} = N_{p r i_1} N_{p r i_2} + h_{s p i_1} \\ h_{r p i_2} = N_{p r i_1} N_{p r i_2} + h_{s p i_2} \end{array}

(3)

where h_pr_{i_}₁ and h_pr_{i_}₂ are the friction arms of the planet in the ring-planet meshing pairs， and h_rp_{i_}₁ and h_rp_{i_}₂ the friction arms of the ring gear in the ring-planet meshing pairs； ε_pr is the contact ratio of the ring-planet meshing pairs， ω_pc the relative angular velocity of the planet relative to the carrier， and N_pr_{i_}₁N_pr_{i_}₂ the theoretical meshing line of the sun-planet meshing pairs.

1.3 Calculation of meshing stiffness and load distribution of planetary gear train

The time-varying meshing stiffness of the internal and external meshing pairs was determined using the energy method^［

19］， as shown in Figs.5 and 6.

Fig.5 Time-varying stiffness of the external meshing gears of planetary gear transmission system

Fig.6 Time-varying stiffness of the internal meshing gears of planetary gear transmission system

Because of the periodicity of the time-varying meshing stiffness， the Fourier series is used. To simplify the calculation of the dynamical equations， the higher-order terms are usually ignored， and the second-order Fourier series is taken. The formula is shown in Eq.（4）， and the results are shown in Figs.2 and 3.

\begin{array}{l} k (t) = k_{0} + A_{1} c o s (ω t + φ) + B_{1} s i n (ω t + φ) + \\ A_{2} c o s (2 ω t + φ) + B_{2} s i n (2 ω t + φ) \end{array}

(4)

where k₀ is the average stiffness of the gear pairs， ω the meshing frequency， and φ the initial phase of the meshing stiffness.

In the process of gear transmission， there will be single and double teeth alternating meshing. In the single-tooth meshing area， the load is borne by a pair of teeth. In the double-teeth meshing area， the load is shared by two pairs of teeth. Because of the different meshing positions in the double-teeth meshing area， the distribution of the load between the two pairs of teeth differs. In the double-teeth meshing zone， the total deformation of each pair of meshing teeth is considered to be equal. The load-distribution ratio is the ratio of the maximum load to the total load between the simultaneous meshing teeth. The results of the load distribution calculated by MATLAB are shown in Fig.7.

Fig.7 Load distribution among the teeth of every meshing gear of the planetary gear transmission system

The correct calculation of the load distribution among the teeth lays the foundation for the calculation of the time-varying friction coefficient， which is described in the next section.

1.4 Calculation of time‑varying friction coefficient of tooth surface

The elastohydrodynamic lubrication （EHL） model comprehensively considers the effects of the load distribution， the relative sliding velocity， the rolling speed， the surface morphology， and the lubrication condition of the gear teeth during the meshing process. Comparing the results obtained by calculation models with different friction coefficients with the experimental values^［

20］ reveals that the values calculated using the EHL friction coefficient are the closest to the experimental values. Therefore， in the paper， the calculation model with the EHL friction coefficient is employed.

This calculation model is expressed as

μ = e^{f (S_{R}, P_{h}, η_{m}, S_{a v g})} P_{h}^{b_{2}} {|S_{R}|}^{b_{3}} V_{e}^{b_{6}} η_{m}^{b_{7}} R^{b_{8}}

(5)

where P_h is the maximum Hertz contact stress （GPa）， R the comprehensive radius of curvature at the contact point （m）， S_R the slip ratio at the contact point， and V_e the convolution rate （m/s）. f（S_R，P_h，η_m，S_avg ） can be written as

\begin{array}{l} f (S_{R}, P_{h}, η_{m}, S_{a v g}) = b_{1} + b_{4} |S_{R}| P_{h} l o g_{10} η_{m} + \\ b_{5} e^{- |S_{R}| P_{h} l o g_{10} η_{m}} + b_{9} e^{S_{a v g}} \end{array}

(6)

The maximum Hertz contact stress is defined as

P_{h} = \sqrt[]{\frac{W' E'}{2 π R}}

(7)

where $W_{}^{'}$ is the unit normal load （GN/m）. $E^{'}$ is the comprehensive elastic modulus （GPa）， which can be calculated as

E' = \frac{2}{\frac{(1 - μ_{1}^{2})}{E_{1}} + \frac{(1 - μ_{2}^{2})}{E_{2}}}

(8)

where μ₁ and μ₂ are Poisson ratios of the driving and driven wheels， respectively， and E₁ and E₂ the elastic moduli （GPa） of the driving and driven wheels， respectively.

The instantaneous velocity of the two gears at any meshing point can be written as

\{\begin{array}{l} V_{p} = ω_{p} ρ_{p} = ω_{p} \times (s + N_{p} B_{1}) \\ V_{g} = ω_{g} ρ_{g} = ω_{g} \times (N_{p} N_{g} - s - N_{p} B_{1}) \end{array}

(9)

where ρ_p and ρ_g are the radii of curvature （m） of the driving and driven wheels， respectively， and s is the distance from the instantaneous meshing point to the actual starting point.

The formulas for the slip ratio， relative sliding velocity， convolution rate， and rolling speed are as follows

\{\begin{array}{l} S_{R} = \frac{V_{s}}{V_{e}} = \frac{2 V_{s}}{V_{r}} = \frac{2 |V_{p} - V_{g}|}{V_{p} + V_{g}} \\ V_{s} = |V_{p} - V_{g}| \\ V_{e} = (V_{p} + V_{g}) / 2 \\ V_{r} = V_{p} + V_{g} \end{array}

(10)

The root-mean-square value of the roughness is S_avg = 0.6 μm. The values of b₁， b₂，…， b₉ are shown in Table 2.

Table 2 Coefficients of the EHL friction model

Parameter	b₁	b₂	b₃
Value	-8.916 4	1.033 0	1.036 0
Parameter	b₄	b₅	b₆
Value	-0.354 0	2.812 0	-0.100 6
Parameter	b₇	b₈	b₉
Value	0.752 7	-0.390 9	0.620 3

1.5 Calculation of time‑varying bearing stiffness

Bearings are important supporting elements in the gear transmission system and have the function of transferring motion and force. According to analysis and application of rolling bearing^［

21］， the stiffness of the bearing is not a constant and changes with the displacement （or load）. Therefore， when referring to the stiffness of the bearing， the corresponding displacement or load state should be specified. According to the Hertz contact theory， the calculation formula for the bearing stiffness of a rolling bearing can be expressed as

k_{b} = 32 900 Z (D_{w} δ_{r} c o s^{5} {α)}^{1 / 2}

(11)

where $D_{w}$ is the outer diameter of rolling body， Z the number of rolling body， $α$ the contact angle， and $δ_{r}$ the elastic displacement of bearing rings.

According to the actual conditions， the sun selects a 6213 radial ball bearing. The carrier selects a 6020 radial ball bearing. The unknown quantities in Eq.（11） are shown in Table 3.

Table 3 Coefficients of the radial ball bearing

Bearing type	6213	6020
D_w/mm	16.67	16
Z	10	14
α/ (°)	0	0

2 Dynamic Model of Gear Transmission System

The dynamic model was established via the lumped-mass method， as shown in Fig.8.

Fig.8 Dynamic model of the planetary gear transmission system

As shown in Fig.8， the bending‑torsional coupled dynamic model of the planetary gear transmission system was established in consideration of the backlash and bearing clearance of the sun. The system has 3N+9 DOFs （N=3）， which can be expressed as

\{x_{s}, y_{s}, θ_{s}, x_{p i}, y_{p i}, θ_{p i}, x_{r}, y_{r}, θ_{r}, x_{c}, y_{c}, θ_{c}\}

(12)

2.1 Force analysis of planetary gear transmission system

The gear meshing force is caused by the relative displacement along the direction of the meshing line. Therefore， it is necessary to analyze the relative displacement and the force and then derive the dynamic equations of the system. The relative position relationship are shown in Figs.9 and 10.

Fig.9 Relative position relationship of the external meshing gears

Fig.10 Relative position relationship of the internal meshing gears

In the figures， $φ_{s p i} = α_{s p i} - φ_{i}$ ， where $α_{s p i}$ is the actual meshing angle of the sun-planet meshing pairs. $φ_{r p i} = α_{r p i} + φ_{i}$ ， where $α_{r p i}$ is the actual meshing angle of the ring-planet meshing pairs.

The projection of the relative displacement of the ith sun-planetary gear along the direction of the meshing line can be written as

\begin{array}{l} δ_{s p i} = (x_{s} - x_{p i}) c o s (90 ° - φ_{s p i}) + \\ (y_{s} - y_{p i}) c o s (- φ_{s p i}) + u_{s} + u_{p i} - e_{s p i} (t) = \\ (x_{s} - x_{p i}) s i n φ_{s p i} + (y_{s} - y_{p i}) c o s φ_{s p i} + \\ u_{s} + u_{p i} - e_{s p i} (t) \end{array}

(13)

The projection of the relative displacement of the ith ring-planetary gear along the direction of the meshing line can be written as

\begin{array}{l} δ_{p r i} = (x_{p i} - x_{r}) c o s (φ_{r p i} - 90 °) + (y_{p i} - y_{r}) \\ c o s (180 ° - φ_{r p i}) + u_{r} - u_{p i} - e_{r p i} (t) = \\ (x_{p i} - x_{r}) s i n φ_{r p i} - (y_{p i} - y_{r}) \cdot \\ c o s φ_{r p i} + u_{r} - u_{p i} - e_{p r i} (t) \end{array}

(14)

The relative position relationship of the planetary‑planet carrier is shown in Fig.11.

Fig.11 Relative position relationship of the planetary‑planet carrier

Therefore， the projection of the relative displacement of the ith planetary-planet carrier along the X， Y， and tangential directions can be expressed as

\{\begin{array}{l} δ_{p i c x} = x_{p i} - x_{c} + u_{c} s i n φ_{i} \\ δ_{p i c y} = y_{p i} - y_{c} - u_{c} c o s φ_{i} \\ δ_{p i c u} = - (x_{p i} - x_{c}) s i n φ_{i} + (y_{p i} - y_{c}) c o s φ_{i} - u_{c} \end{array}

(15)

To facilitate the analysis， the vibration angular displacements （θ_s， θ_p_i， θ_r， θ_c） of the sun， planet， inner ring gear， and carrier are transformed into line displacements （u_s， u_p_i， u_r， u_c）.

\{\begin{array}{l} u_{s} = r_{b s} θ_{s} & u_{p i} = r_{b p i} θ_{p i} \\ u_{r} = r_{b r} θ_{r} & u_{c} = r_{b c} θ_{c} \end{array}

(16)

where r_bs， r_bp_i， and r_br are the radii of the base circle of the sun， planet， and ring gears， respectively； and r_bc is the center distance between the sun and the planet.

2.2 Dynamic equations of planetary gear train

According to the force relationship of Figs.6 and 7， the dynamic differential equations of the system are obtained using Newton’s law of motion.

（1） Dynamic equations for the sun

\{\begin{array}{l} M_{s} ({\ddot{x}}_{s} - 2 ω_{c} {\dot{y}}_{s} - ω_{c}^{2} x_{s}) + k_{s x} f (x_{s}, b_{s x}) + c_{s x} {\dot{x}}_{s} = \\ - \sum_{i = 1}^{N} F_{s p i} s i n φ_{s p i} - \sum_{i = 1}^{N} F_{f s p i} c o s φ_{s p i} \\ M_{s} ({\ddot{y}}_{s} + 2 ω_{c} {\dot{x}}_{s} - ω_{c}^{2} y_{s}) + k_{s y} f (y_{s}, b_{s y}) + c_{s y} {\dot{y}}_{s} = \\ - \sum_{i = 1}^{N} F_{s p i} c o s φ_{s p i} + \sum_{i = 1}^{N} F_{f s p i} s i n φ_{s p i} \\ I_{s} {\ddot{θ}}_{s} = - \sum_{I = 1}^{N} F_{s p i} r_{b s} + \sum_{i = 1}^{N} F_{f s p i} h_{s p i} + T_{D} \end{array}

(17)

where I_s is the moment of inertia of the sun； M_s the quality of the sun； and T_D the input torque. F_sp_i and F_fsp_i are the dynamic meshing force and dynamic friction force of the sun in the ith sun-planet meshing pair， respectively.

（2） Dynamic equations for the planets

\{\begin{array}{l} M_{p i} ({\ddot{x}}_{p i} - 2 ω_{c} {\dot{y}}_{p i} - ω_{c}^{2} x_{p i}) + k_{p i x} δ_{p i c x} + c_{p i x} {\overset{\cdot}{δ}}_{p i c x} = \\ F_{p s i} s i n φ_{s p i} + F_{f p s i} c o s φ_{s p i} - F_{p r i} s i n φ_{r p i} - F_{f p r i} c o s φ_{r p i} \\ M_{p i} ({\ddot{y}}_{p i} + 2 ω_{c} {\dot{x}}_{p i} - ω_{c}^{2} y_{p i}) + k_{p i y} δ_{p i c y} + c_{p i y} {\overset{\cdot}{δ}}_{p i c y} = \\ F_{p s i} c o s φ_{s p i} - F_{f p s i} s i n φ_{s p i} + F_{p r i} c o s φ_{r p i} - F_{f p r i} s i n φ_{r p i} \\ I_{p i} {\ddot{θ}}_{p i} = - F_{p s i} r_{b p i} + F_{f p s i} h_{p s i} + F_{p r i} r_{b p i} - F_{f p r i} h_{p r i} \end{array}

(18)

where I_p_i is the moment of inertia of the ith planet and M_p_i the quality of the ith planet. F_ps_i and F_fps_i are the dynamic meshing force and dynamic friction force of the planet in the ith sun-planet meshing pair， respectively； and F_pr_i and F_fpr_i the dynamic meshing force and dynamic friction force of the planet in the ith planet‑ring meshing pair， respectively.

（3） Dynamic equations for the ring gear

\{\begin{array}{l} M_{r} ({\ddot{x}}_{r} - 2 ω_{c} {\dot{y}}_{r} - ω_{c}^{2} x_{r}) + k_{r x} x_{r} + c_{r x} {\dot{x}}_{r} = \\ \sum_{i = 1}^{N} F_{r p i} s i n φ_{r p i} + \sum_{i = 1}^{N} F_{f r p i} c o s φ_{r p i} \\ M_{r} ({\ddot{y}}_{r} + 2 ω_{c} {\dot{x}}_{r} - ω_{c}^{2} y_{r}) + k_{r y} y_{r} + c_{r y} {\dot{y}}_{r} = \\ - \sum_{i = 1}^{N} F_{r p i} c o s φ_{r p i} + \sum_{i = 1}^{N} F_{f r p i} s i n φ_{r p i} \\ I_{r} {\ddot{θ}}_{r} = - \sum_{i = 1}^{N} F_{r p i} r_{b r} + \sum_{i = 1}^{N} F_{f r p i} h_{r p i} \end{array}

(19)

where I_r is the moment of inertia of the ring and M_r the quality of the ring. F_rp_i and F_frp_i are the dynamic meshing force and dynamic friction force of ring in the ith planet‑ring meshing pair， respectively.

（4） Dynamic equations for the carrier

\{\begin{array}{l} M_{c} ({\ddot{x}}_{c} - 2 ω_{c} {\dot{y}}_{c} - ω_{c}^{2} x_{c}) + k_{c x} x_{c} + c_{c x} {\dot{x}}_{c} = \\ \sum_{i = 1}^{N} (k_{p i x} δ_{p i c x} + c_{p i x} {\dot{δ}}_{p i c x}) \\ M_{c} ({\ddot{y}}_{c} + 2 ω_{c} {\dot{x}}_{c} - ω_{c}^{2} y_{c}) + k_{c y} y_{c} + c_{c y} {\dot{y}}_{c} = \\ \sum_{i = 1}^{N} (k_{p i c y} δ_{p i c y} + c_{p i c y} {\dot{δ}}_{p i c y}) \\ I_{c} {\ddot{θ}}_{c} + k_{c u} θ_{c} + c_{c u} {\ddot{θ}}_{c} = \sum_{i = 1}^{N} (k_{p i c u} δ_{p i c u} + c_{p i c u} {\dot{δ}}_{p i c u}) r_{b c} - T_{c} \end{array}

(20)

where I_c is the moment of inertia of the carrier， M_c the quality of the carrier， T_c the output torque， k_c_u the torsional stiffness， and c_c_u the torsional damping.

c represents bearing damping； sx， sy， piy， pix， ry， rx， cx， cy are the subscript of the sun gear， the planet， the ring gear and the carrier， respectively. c is determined as Ref.［

22］.

2.3 Rigid-body displacement elimination and dimensionlessness of dynamic equations

Because the dynamic Eqs.（17）—（20） are positive semidefinite， there is rigid-body displacement， and the solution is uncertain. Therefore， the relative coordinates δ_sp_i and δ_rp_i are introduced， and the concrete expressions are Eqs.（13） and （14）. Thus， we can conclude the following

\begin{array}{l} δ_{p r i} = (x_{p i} - x_{r}) s i n φ_{r p i} - (y_{p i} - y_{r}) c o s φ_{r p i} + \\ δ_{p r N} - (x_{p N} - x_{r}) s i n φ_{r p N} + \\ (y_{p N} - y_{r}) c o s φ_{r p N} + e_{r p N} (t) + δ_{s p N} - \\ (x_{s} - x_{p N}) s i n φ_{s p N} - (y_{s} - y_{p N}) c o s φ_{s p N} + \\ e_{s p N} (t) - [δ_{s p i} - (x_{s} - x_{p i}) s i n φ_{s p i} - \\ (y_{s} - y_{p i}) c o s φ_{s p i} + e_{s p i} (t)] - e_{p r i} (t) \end{array}

(21)

After the elimination of the rigid-body displacement， the system has 3N + 8 DOFs， which can be written as

\{x_{s}, y_{s}, x_{p i}, y_{p i}, x_{r}, y_{r}, x_{c}, y_{c}, δ_{s p i}, δ_{p r N}, u_{c}\}

(22)

Because the numerical gap of the stiffness value and the vibration micro-displacement value is too large， the calculated results cannot converge under the numerical integration method. To obtain the ideal results， the dimensionless displacement （b_c= 10^-5 mm） and the natural frequency ω_sn are introduced into the dynamic differential equations. ω_sn is defined as

ω_{s n} = \sqrt[]{k_{s p} \frac{I_{p i} r_{b s}^{2} + I_{s} r_{b p i}^{2}}{I_{s} I_{p i}}}

(23)

where k_sp is the average meshing stiffness of the sun-planet meshing pairs. Then， the dimensionless time and dimensionless displacement are

\begin{matrix} \bar{X} = X / b_{c} & τ = ω_{s n} t \end{matrix}

(24)

The dimensionless acceleration and velocity can be written as

\{\begin{array}{l} \dot{X} = b_{c} ω_{s n} \dot{\bar{X}} & \ddot{X} = b_{c} ω_{s n}^{2} \ddot{\bar{X}} \\ \ddot{\bar{e}} (τ) = \ddot{e} (t) / b_{c} ω_{s n}^{2} & \dot{\bar{e}} (τ) = \dot{e} (t) / b_{c} ω_{s n} \end{array}

(25)

After the elimination of the rigid-body displacement and dimensionlessness， the dynamic differential equations can be expressed as

\{\begin{array}{l} {\ddot{\bar{x}}}_{s} = 2 \frac{ω_{c}}{ω_{s n}} {\dot{\bar{y}}}_{s} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} f ({\bar{x}}_{s}, {\bar{b}}_{s x}) - \frac{k_{s x}}{M_{s} ω_{s n}^{2}} f ({\bar{x}}_{s}, {\bar{b}}_{s x}) - \frac{c_{s x}}{M_{s} ω_{s n}} {\dot{\bar{x}}}_{s} - \frac{1}{M_{s} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{s p i} s i n φ_{s p i} - \frac{1}{M_{s} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{f s p i} c o s φ_{s p i} \\ {\ddot{\bar{y}}}_{s} = - 2 \frac{ω_{c}}{ω_{s n}} {\dot{\bar{x}}}_{s} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} f ({\bar{y}}_{s}, {\bar{b}}_{s y}) - \frac{k_{s y}}{M_{s} ω_{s n}^{2}} f ({\bar{y}}_{s}, {\bar{b}}_{s y}) - \frac{c_{s y}}{M_{s} ω_{s n}} {\dot{\bar{y}}}_{s} - \frac{1}{M_{s} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{s p i} c o s φ_{s p i} + \frac{1}{M_{s} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{f s p i} s i n φ_{s p i} \\ {\ddot{\bar{x}}}_{p i} = 2 \frac{ω_{c}}{ω_{s n}} {\dot{\bar{y}}}_{p i} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} {\bar{x}}_{p i} - \frac{1}{M_{p i} ω_{s n}^{2}} k_{p i x} {\bar{δ}}_{p i c x} - \frac{1}{M_{p i} ω_{s n}} c_{p i x} {\dot{\bar{δ}}}_{p i c x} + \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{p s i} s i n φ_{s p i} + \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{f p s i} c o s φ_{s p i} - \\ \begin{matrix}  \end{matrix} \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{p r i} s i n φ_{r p i} - \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{f p r i} c o s φ_{r p i} \\ {\ddot{\bar{y}}}_{p i} = - 2 \frac{ω_{c}}{ω_{s n}} {\dot{\bar{x}}}_{p i} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} {\bar{y}}_{p i} - \frac{1}{M_{p i} ω_{s n}^{2}} k_{p i y} {\bar{δ}}_{p i c y} - \frac{1}{M_{p i} ω_{s n}} c_{p i y} {\dot{\bar{δ}}}_{p i c y} + \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{p s i} c o s φ_{s p i} - \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{f p s i} s i n φ_{s p i} + \\ \begin{matrix}  \end{matrix} \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{p r i} c o s φ_{r p i} - \frac{1}{M_{p i} ω_{s n}^{2}} {\bar{F}}_{f p r i} s i n φ_{r p i} \\ {\ddot{\bar{x}}}_{r} = 2 \frac{ω_{c}}{ω_{s n}} {\dot{\bar{y}}}_{r} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} {\bar{x}}_{r} - \frac{1}{M_{r} ω_{s n}^{2}} k_{r x} {\bar{x}}_{r} - \frac{1}{M_{r} ω_{s n}} c_{r x} {\dot{\bar{x}}}_{r} + \frac{1}{M_{r} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{r p i} s i n φ_{r p i} + \frac{1}{M_{r} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{f r p i} c o s φ_{r p i} \\ {\ddot{\bar{y}}}_{r} = - 2 \frac{ω_{c}}{ω_{s n}} {\dot{x}}_{r} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} y_{r} - \frac{1}{M_{r} ω_{s n}^{2}} k_{r y} {\bar{y}}_{r} - \frac{1}{M_{r} ω_{s n}} c_{r y} {\dot{\bar{y}}}_{r} - \frac{1}{M_{r} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{r p i} c o s φ_{r p i} + \frac{1}{M_{r} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{f r p i} s i n φ_{r p i} \\ {\ddot{\bar{x}}}_{c} = 2 \frac{ω_{c}}{ω_{s n}} {\dot{\bar{y}}}_{c} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} {\bar{x}}_{c} - \frac{1}{M_{c} ω_{s n}^{2}} k_{c x} {\bar{x}}_{c} - \frac{1}{M_{c} ω_{s n}} c_{c x} {\dot{\bar{x}}}_{c} + \frac{1}{M_{c} ω_{s n}^{2}} \sum_{i = 1}^{N} (k_{p i x} {\bar{δ}}_{p i c x} + c_{p i x} {\dot{\bar{δ}}}_{p i c x}) \\ {\ddot{\bar{y}}}_{c} = - 2 \frac{ω_{c}}{ω_{s n}} {\dot{x}}_{c} + \frac{ω_{c}^{2}}{ω_{s n}^{2}} y_{c} - \frac{1}{M_{c} ω_{s n}^{2}} k_{c y} {\bar{y}}_{c} - \frac{1}{M_{c} ω_{s n}} c_{c y} {\dot{\bar{y}}}_{c} + \frac{1}{M_{c} ω_{s n}^{2}} \sum_{i = 1}^{N} (k_{p i y} {\bar{δ}}_{p i c y} + c_{p i y} {\dot{\bar{δ}}}_{p i c y}) \\ {\ddot{\bar{δ}}}_{s p i} = ({\ddot{\bar{x}}}_{s} - {\ddot{\bar{x}}}_{p i}) s i n φ_{s p i} + ({\ddot{\bar{y}}}_{s} - {\ddot{\bar{y}}}_{p i}) c o s φ_{s p i} - {\ddot{\bar{e}}}_{s p i} (t) - \frac{1}{m_{e s} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{s p i} + \frac{1}{r_{b s} m_{e s} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{f s p i} h_{s p i} + \frac{T_{D}}{r_{b s} m_{e s} b_{c} ω_{s n}^{2}} - \\ \begin{matrix}  \end{matrix} \frac{1}{m_{e p i} ω_{s n}^{2}} {\bar{F}}_{p s i} + \frac{1}{r_{b p i} m_{e p i} ω_{s n}^{2}} {\bar{F}}_{f p s i} h_{p s i} + \frac{1}{m_{e p i} ω_{s n}^{2}} {\bar{F}}_{p r i} - \frac{1}{r_{b p i} m_{e p i} ω_{s n}^{2}} {\bar{F}}_{f p r i} h_{p r i} \\ {\ddot{\bar{δ}}}_{p r N} = ({\ddot{\bar{x}}}_{p i} - {\ddot{\bar{x}}}_{r}) s i n φ_{r p N} - ({\ddot{\bar{y}}}_{p i} - {\ddot{\bar{y}}}_{r}) c o s φ_{r p N} - {\ddot{\bar{e}}}_{r p N} (t) - \frac{1}{m_{e r} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{r p i} + \frac{1}{r_{b r} m_{e r} ω_{s n}^{2}} \sum_{i = 1}^{N} {\bar{F}}_{f r p i} h_{r p i} + \\ \begin{matrix}  \end{matrix} \frac{1}{m_{e p N} ω_{s n}^{2}} {\bar{F}}_{p s N} - \frac{1}{r_{b p N} m_{e p N} ω_{s n}^{2}} {\bar{F}}_{f p s N} h_{p s N} - \frac{1}{m_{e p N} ω_{s n}^{2}} {\bar{F}}_{p r N} + \frac{1}{r_{b p N} m_{e p N} ω_{s n}^{2}} {\bar{F}}_{f p r N} h_{p r N} \\ {\bar{u}}_{c} = \frac{1}{m_{e c} ω_{s n}^{2}} \sum_{i = 1}^{N} (k_{p i c u} {\bar{δ}}_{p i c u} + c_{p i c u} ω_{s n} {\dot{\bar{δ}}}_{p i c u}) - \frac{k_{c u}}{m_{e c} r_{b c}^{2}} {\bar{u}}_{c} - \frac{c_{c u}}{m_{e c} r_{b c}^{2}} {\dot{\bar{u}}}_{c} - \frac{T_{c}}{r_{b c} m_{e c} b_{c} ω_{s n}^{2}} \end{array}

(26)

The specific expressions for ${\bar{F}}_{s p i}, {\bar{F}}_{p r i}, {\bar{F}}_{f s p i}, {\bar{F}}_{f p r i}$ are

\{\begin{array}{l} {\bar{F}}_{s p i} = k_{s p i} f ({\bar{δ}}_{s p i}, {\bar{b}}_{s p i}) + c_{s p i} ω_{s n} {\dot{\bar{δ}}}_{s p i} \\ {\bar{F}}_{f s p i} = f_{s p i} {\bar{F}}_{s p i} \\ {\bar{F}}_{p r i} = k_{p r i} f ({\bar{δ}}_{p r i}, {\bar{b}}_{p r i}) + c_{p r i} ω_{s n} {\dot{\bar{δ}}}_{p r i} \\ {\bar{F}}_{f p r i} = f_{p r i} {\bar{F}}_{p r i} \end{array}

（27）

3 Dynamic Analysis of Planetary Gear Trains

According to the above analysis， the effects of the input speed， backlash， bearing clearance， and damping ratio on the global bifurcation characteristics were investigated using the 4—5-order Runge‑Kutta method.

3.1 Influence of input speed on bifurcation characteristics

As shown in the previous section， the vibration and instability of the planetary gear transmission system are mainly related to the change of the meshing force. Therefore， it is necessary to study the bifurcation characteristics of the relative displacement along the meshing line. To determine the effect of the velocity on the global bifurcation characteristics， the gear parameters and operation conditions shown in Tables 1 and 4 were used， respectively. The bifurcation characteristics of PGTSMPP along the meshing direction are shown in Figs.12 and 13.

Fig.12 Bifurcation along the external meshing line direction

Fig.13 Bifurcation along the internal meshing line direction

As shown in Figs.12 and 13， with the increase of the velocity， the bifurcation characteristics along the external and internal meshing lines are similar and rich. At a low speed， the system is mainly in the single-period state. With the increase of the speed， the system enters the states of single-period， violent change， chaos， bifurcation， inverted bifurcation， and multi-period， as well as other motion states. At the frequency of 1/2， 1， 2， and 3 times， the system exhibits a significant resonance phenomenon.

Because of the similarity of the bifurcation characteristics along the external and internal meshing lines， it is only necessary to analyze the system entering the chaotic channel with the change of the speed， as shown in Fig.12. The system enters a chaotic state at the dimensionless meshing frequency of Ω_sp=0.6 and enters the single-period state from the inverted-bifurcation state at Ω_sp=0.75—0.78. At Ω_sp=1.15， the system re-enters the chaotic state， which is accompanied by the multi-period motion state. At Ω_sp=1.45， the system enters single-period state， and when the meshing frequency is close to twice the natural frequency， the system undergoes a sudden change and continues operating for a while. At Ω_sp=2.75， the system transitions from the bifurcation state to the period two motion state. The system transitions from the inverted bifurcation state to the stable single-period state at Ω_sp=3.06—3.11.

Therefore， the appropriate speed can effectively prevent the system from entering the chaotic state， improve the stability， load-sharing properties and lifetime of the system， and reduce the vibration and noise.

3.2 Influence of damping ratio on bifurcation characteristics

To analyze the influence of damping ratio on the global bifurcation characteristics， without losing generality， we consider Ω_sp = 1.21 and Ω_sp = 2.83， with typical nonlinear dynamic characteristics. The operating conditions are shown in Table 4， and the results are presented in Fig.14.

Table 4 Operating conditions

Power/kW	Damping ratio	Backlash b_sp_i /μm	Backlash b_pr_i /μm	Bearing Clearance
Power/kW	Damping ratio	Backlash b_sp_i /μm	Backlash b_pr_i /μm	b_s_x /μm b_s_y /μm
400	0.05	40	40	10	10

Fig.14 Bifurcation with different damping ratios

As shown in Fig.14（a）， with the increase of the damping ratio， the amplitude is significantly suppressed and converges to zero. Fig.14（b） indicates that when the damping ratio is less than 0.35， the system is in the chaotic state， and the system enters the period two motion state with the further increase of the damping ratio. When the damping ratio is greater than 0.075， the system is in the single-period state， and the amplitude decreases significantly.

To confirm that the increase of the damping ratio can suppress the vibration amplitude of the chaos， the damping ratio coefficients ζ of 0.05， 0.07， and 0.09 are selected； the corresponding bifurcation results are shown in Figs.12 and 15. Under large damping ratio， the vibration of the system is suppressed. This is because with the increase of the damping ratio， the amount of dissipated energy increases， which suppresses the movement of the system.

Fig.15 Bifurcation with different damping ratios

Therefore， the damping can be increased by changing the material of the meshing pairs， which can effectively ameliorate the vibration characteristics of the system， reduce the noise， and improve the system stability.

3.3 Influence of backlash on bifurcation characteristics

To analyze the influence of the sun-planet meshing backlash on the global bifurcation characteristics， without losing generality， we consider Ω_sp=0.78 and Ω_sp=1.17， with typical nonlinear dynamic characteristics. The operating conditions are shown in Table 4. Assuming that the change of the backlash of three sun‑planet meshing pairs is the same， the results are shown in Fig.16.

Fig.16 Bifurcation with respect to the backlash

As shown in Fig.16（a）， with the increase of the backlash， the nonlinear characteristics of the system are obvious， and it has the same characteristics no matter which one is analyzed. For a small backlash， the system is in a stable single-period state. For a large backlash， the system is in a chaotic state. When the meshing frequency is changed， only the channel of the system entering the chaotic state is affected.

Fig.16（b） shows that when the dimensionless backlash ${\bar{b}}_{s p i}$ is less than 1.26， the system is in the single-period state， and with the increase of backlash， system enters the period three motion state from the bifurcation state. At ${\bar{b}}_{s p i}$ = 1.7， the system is in the chaotic state. When ${\bar{b}}_{s p i}$ is 2.4—3.2， the system is in the multi-period state. When ${\bar{b}}_{s p i}$ is greater than 3.2， the system is in the chaotic state.

Thus， in the design， manufacture， and installation of the gear system， reasonable backlash can effectively prevent the system from entering the chaotic state， which suppresses the system vibration and noise， improves the system stability， and extends the service life of the product.

3.4 Influence of bearing clearance on bifurcation characteristics

To analyze the influence of the bearing clearance on the global bifurcation characteristics， without losing generality， we consider Ω_sp = 0.85 and Ω_sp = 2.05. The operating conditions are shown in Table 4. Assuming that the bearing clearance in the horizontal direction of the sun is the same as that in the vertical direction. The results are shown in Fig.17， respectively.

Fig.17 Bifurcation in different bearing clearance

As shown in Fig.17， there is no multi-period state in the whole parameter field. In Fig.17（a）， only when the dimensionless bearing clearance ${\bar{b}}_{s x}$ or ${\bar{b}}_{s y}$ is 0—0.8， the system state changes with the change of the bearing clearance of the sun. In other areas， the system is in a stable state. Because the sun floating support structure can be understood as a larger sun bearing clearance， it can more effectively show that the planetary gear transmission system using the sun floating support structure can improve the nonlinear dynamic characteristics and load-sharing properties of the system.

4 Conclusions

（1） With the increase of the damping ratio， the amplitude of the system can be significantly restrained in the chaotic state. The system near the resonant frequency region is often accompanied by an unstable motion state， such as violent change and chaos， and the proper selection of the system speed parameters can effectively prevent the chaotic motion and improve the stability of the system.

（2） The backlash has a great impact on the bifurcation characteristics of the system. With the increase of the backlash， the system transitions from the single-period state to the multi-period and chaotic states.

（3） Compared with the backlash， the bearing clearance of the sun has little effect on the bifurcation characteristics. The influence of the nonlinear characteristics of the system is obvious only in the segment with small parameters of the sun bearing clearance， which indicates the superiority of the planetary gear transmission system.

The results can provide theoretical support for the parameter selection and operating conditions of PGTSMPP.

Contributions Statement

Mr. TANG Xin designed the study. Dr. BAO Heyun wrote the manuscript. Dr. LU Fengxia established the dynamic model. Dr. JIN Guanghu solved the dynamic equations. All authors commented on the manuscript draft and approved the submission.

Acknowledgements

This work was supported by the National Natural Science Foundation of China（No.51975274） and National Key Laboratory of Science and Technology on Helicopter Transmission（Nanjing University of Aeronautics and Astronautics）（No. HTL-A-19K03）.

Conflict of Interest

The authors declare no competing interests.

References

GAO Ming， ZHOU Ying. Study of addendum modification coefficient optimization of spur gear engaging beyond pitch point‍［J］. Machinery， 1997， 24（2）： 2-4. (in Chinese) [百度学术]

LIU Jingjing. Design method and dynamics analysis of node external gear-pair‍［D］. Nanjing： Nanjing University of Aeronautics and Astronautics， 2012.（in Chinese） [百度学术]

SUN Yongzheng. Design methods and dynamic analysis of planetary gear transmission drive with meshing beyond pitch point‍［D］. Nanjing： Nanjing University of Aeronautics and Astronautics， 2013.（in Chinese） [百度学术]

BAO Heyun， ZHOU Xingjun， ZHU Rupeng. Load sharing analysis of planetary gear with meshing beyond pitch point considering flexible deformable ring‍［J］. Journal of Central South University （Science and Technology）， 2016， 47（9）： 3005-3010. (in Chinese) [百度学术]

KAHRAMAN A. Planetary gear train dynamics‍［J］. Journal of Mechanical Design， 1994， 116（3）： 713-720. [百度学术]

AL-SHYYAB A， KAHRAMAN A. A non-linear dynamic model for planetary gear sets‍［J］. Proceedings of the Institution of Mechanical Engineers Part K Journal of Multi-body Dynamics， 2007， 221（4）： 567-576. [百度学术]

AMBARISHA V K， PARKER R G. Nonlinear dynamics of planetary gears using analytical and finite element models‍［J］. Journal of Sound and Vibration， 2007， 302（3）： 577-595. [百度学术]

BLANKENSHIP G W， KAHRAMAN A. Steady state forced response of a mechanical oscillator with combined parametric excitation and clearance type non-linearity‍［J］. Journal of Sound and Vibration， 1995， 185（5）： 743-765. [百度学术]

HE S， GUNDA R， SINGH R. Effect of sliding friction on the dynamics of spur gear pair with realistic time-varying stiffness‍［J］. Journal of Sound and Vibration， 2007， 301（3/4/5）： 927-949. [百度学术]

WANG Sanming， SHEN Yunwen， DONG Haijun. Chaos and bifurcation analysis of spur gear pair with combined friction and clearance‍［J］. Chinese Journal of Mechanical Engineering， 2002， 38（9）： 217-223. (in Chinese) [百度学术]

TANG Jinyuan， CHEN Siyu， ZHONG Jue. A improved nonlinear model for a spur gear pair system‍［J］. Engineering Mechanics， 2008， 25（1）： 217-223. (in Chinese) [百度学术]

YE Fumin， ZHU Rupeng， JIN Guanghu. Load sharing characteristics of a planetary gear train system with non-equivalent modulus and pressure angle considering backlashes‍［J］. Journal of Vibration and Shock， 2015， 34（11）： 206-211. (in Chinese) [百度学术]

SHENG Dongping， ZHU Rupeng， BAO Heyun. Bifurcation characteristics and modeling of transverse-torsional spur gear vibration considering friction and multiple clearances‍［J］. Journal of Aerospace Power， 2015， 30（2）： 498-506. (in Chinese) [百度学术]

LI Fajia. The strength and dynamics analysis of hight contact ratio planetary gear transmission‍［D］. Nanjing： Nanjing University of Aeronautics and Astronautics， 2015.（in Chinese） [百度学术]

MO Shuai， ZHANG Ting. Analytical investigation on load sharing characteristics of herringbone planetary gear train with flexible support and floating sun gear‍［J］. Mechanism and Machine Theory， 2020， 144 （2）： 1-27. (in Chinese) [百度学术]

MO Shuai， YUE Zhongxiang. Analytical investigation on load-sharing characteristics for multi-power face gear split flow system‍［J］. Journal of Mechanical Engineering Science， 2020， 234（2）： 676-692. (in Chinese) [百度学术]

JIN Guanghu，YANG Haoyun，LU Fengxia，et al. Effect of friction on dynamic response of a power split transmission system‍［J］. Transactions of Nanjing University of Aeronautics and Astronautics， 2019， 36（3）： 387-400. [百度学术]

BAO Heyun， ZHANG Yayun， ZHU Rupeng， et al. Dynamic analysis of external gear system with meshing beyond pitch considering time-varing friction coefficient‍［J］. Journal of Nanjing University of Aeronautics & Astronautics， 2016， 48（6）： 815-821.（in Chinese） [百度学术]

ZHU Xiaolu， E Zhongkai. Analysis of load capacity of gears［M］. Beijing： Higher Education Press， 1992. (in Chinese) [百度学术]

XU Hai. Development of a generalized mechanical efficiency prediction methodology for gear pairs‍［D］. Columbus： The Ohio State University， 2005. [百度学术]

LUO Jiwei， LUO Tianyu. Rolling bearing analysis and calculation and application［M］. Beijing： China Machine Press， 2009. (in Chinese) [百度学术]

WU Hao. Research on the dynamic characteristics of rolling element bearings and the dynamic model of bearing rotor system‍［D］.Wuhan： Huazhong University of Science and Technology， 2011.（in Chinese） [百度学术]