，，，

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

（Received 18 April 2020；revised 9 July 2020；accepted 14 September 2020）

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.

Key words：meshing beyond pitch point；planetary gear system；nonlinear；time-varying bearing stiffness；timevarying 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 zoneACandEC，as shown in Fig.1.In Fig.1，Cis the gear mesh point，F(xiàn)μthe tooth surface friction force，and arrow the direction of the friction force.In the standard meshing gear transmission system，F(xiàn)μwill change the direction at the pointCthat may cause vibration of the system.

As a non-standard meshing gear transmission system，the gear transmission system only meshes in one side of theCpoint，as shown in Figs.2 and 3.In Fig.2，ACis the actual mesh zone.In Fig.3，CEis the actual mesh zone.

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

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［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［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，timevarying 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.［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 adynamicmodelby 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.［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.［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.［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.［15-16］investigated the load sharing characteristics of herringbone planetary gear train and multi-power face gear split flow system.Jin et al.［17］studied the effect of friction on dynamic response of a power split transmission system.Bao et al.［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

1.1 Calculation of phase angle

λspiis the meshing phase coefficient between theith sun-planet meshing pair and the first sun-planet meshing pair，λrpithe meshing phase coefficient between theith ring-planet meshing pair and the first ring-planet meshing pair，andλsrthe phase difference between the sun and the ring.The corresponding formulas can be written as

wherezsandzrare the numbers of teeth for the sun and ring gears，respectively；φiis the position angle of the planet，ppbthe base pitch of the planet，F(xiàn)′2the point of the opposite tooth surface of the planet base circle relative to the meshing starting point of the sun，Br2the 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 lineBpri_1Bpri_2is located on the side of the pitch point，and the friction arms at any moment according to the geometric relationship are given as follows

wherehspi_1andhspi_2are the friction arms of the sun in the sun-planet meshing pairs，andhpsi_1andhpsi_2the friction arms of the planet in the sun-planet meshing pairs；εspis the contact ratio of the sun-planet meshing pairs，ωscthe relative angular velocity of the sun relative to the carrier，pbthe base pitch；andNspi_1Nspi_2the theoretical meshing line of the sunplanet meshing pairs.

wherehpri_1andhpri_2are the friction arms of the planet in the ring-planet meshing pairs，andhrpi_1andhrpi_2the friction arms of the ring gear in the ringplanet meshing pairs；εpris the contact ratio of the ring-planet meshing pairs，ωpcthe relative angular velocity of the planet relative to the carrier，andNpri_1Npri_2the theoretical meshing line of the sunplanet 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.

wherek0is 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

1.4 Calculation of time-varying friction coefficient of tooth surface

Theelastohydrodynamiclubrication （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

wherePhis the maximum Hertz contact stress（GPa），Rthe comprehensive radius of curvature at the contact point（m），SRthe slip ratio at the contact point，andVethe convolution rate（m/s）.f（SR，Ph，ηm，Savg）can be written as

The maximum Hertz contact stress is defined as

whereW′is the unit normal load（GN/m）.E′is the comprehensive elastic modulus（GPa），which can be calculated as

whereμ1andμ2are Poisson ratios of the driving and driven wheels，respectively，andE1andE2the elastic moduli（GPa）of the driving and driven wheels，respectively.

Various princesses were proposed to him, and the fairy, who was anxious to get the affair over before she left the Court for ever, gave it as her opinion that the Princess Diaphana would make the most suitable wife

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

whereρpandρgare the radii of curvature（m）of the driving and driven wheels，respectively，andsis 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

The root-mean-square value of the roughness isSavg=0.6 μm.The values ofb1，b2，… ，b9are shown in Table 2.

Table 2 Coefficients of the EHL friction model

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

whereDwis the outer diameter of rolling body，Zthe number of rolling body，αthe contact angle，andδrthe 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

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

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，φspi=αspi-φi，whereαspiis the actual meshing angle of the sun-planet meshing pairs.φrpi=αrpi+φi，whereαrpiis the actual meshing angle of the ring-planet meshing pairs.

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

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

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 theith planetary-planet carrier along theX，Y，and tangential directions can be expressed as

To facilitate the analysis，the vibration angular displacements（θs，θpi，θr，θc）of the sun，planet，inner ring gear，and carrier are transformed into line displacements（us，upi，ur，uc）.

whererbs，rbpi，andrbrare the radii of the base circle of the sun，planet，and ring gears，respectively；andrbcis 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

whereIsis the moment of inertia of the sun；Msthe quality of the sun；andTDthe input torque.FspiandFfspiare the dynamic meshing force and dynamic friction force of the sun in theith sun-planet meshing pair，respectively.

（2）Dynamic equations for the planets

whereIpiis the moment of inertia of theith planet andMpithe quality of theith planet.FpsiandFfpsiare the dynamic meshing force and dynamic friction force of the planet in theith sun-planet meshing pair，respectively；andFpriandFfprithe dynamic meshing force and dynamic friction force of the planet in theith planet-ring meshing pair，respectively.

（3）Dynamic equations for the ring gear

whereIris the moment of inertia of the ring andMrthe quality of the ring.FrpiandFfrpiare the dynamic meshing force and dynamic friction force of ring in theith planet-ring meshing pair，respectively.

（4）Dynamic equations for the carrier

whereIcis the moment of inertia of the carrier，Mcthe quality of the carrier，Tcthe output torque，kcuthe torsional stiffness，andccuthe torsional damping.

crepresents bearing damping；sx，sy，piy，pix，ry，rx，cx，cyare 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δspiandδrpiare introduced，and the concrete expressions are Eqs.（13）and（14）.Thus，we can conclude the following

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

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（bc=10-5mm）and the natural frequencyωsnare introduced into the dynamic differential equations.ωsnis defined as

wherekspis the average meshing stiffness of the sunplanet meshing pairs.Then，the dimensionless time and dimensionless displacement are

The dimensionless acceleration and velocity can be written as

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

The specific expressions forare

3 DynamicAnalysisofPlanetary 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 Influenceofinputspeed 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 ofthe 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.

Table 4 Operating conditions

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.

Fig.12 Bifurcation along the external meshing line direction

Fig.13 Bifurcation along the internal meshing line direction

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.

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.

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.

Fig.14 Bifurcation with different damping ratios

Fig.15 Bifurcation with different damping ratios

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.

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 backlashis less than 1.26，the system is in the single-period state，and with the increase of back-lash，system enters the period three motion state from the bifurcation state.At=1.7，the system is in the chaotic state.Whenis 2.4—3.2，the system is in the multi-period state.Whenis 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.

Fig.16 Bifurcation with respect to the backlash

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 clearanceoris 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.