Nonlinear characteristic analysis of gear rattle based on refined dynamic model

In order to explore the vibration response and the nonlinear behavior of gear rattle, a nonlinear dynamic model is developed in this paper considering the nonlinear oil film force, time-varying meshing stiffness, dragging torque and friction and gear backlash with time-varying characteristics. The time-varying meshing stiffness is calculated using the potential energy method. Considering the convolution and extrusion of the lubricant between the teeth, it is modeled as a nonlinear spring-damped element and a nonlinear oil film force is introduced. The drag torque is calculated numerically using fluid dynamics theory, friction principle, momentum conservation principle, and Bernoulli's principle. The model is verified by comparing the bench test results with the simulation results, and the parametric study is performed to investigate the nonlinear system characteristics by means of fast Fourier transform, Poincare map, and bifurcation diagram. The results show that as the excitation frequency increases, the system enters chaotic state under multiple frequency jumps. Excessive torque fluctuation may lead to severe rattling phenomenon of the gear pair and making the system enter a chaotic state; as the lubricant viscosity increases, the system shows periodic, bifurcation and chaotic behavior alternately, and when the viscosity is large enough, the high-frequency shock intensity decays more rapidly due to large drag torque, which is beneficial to reduce the percussion strength and prevent the system into chaos.

has been carried out. For example, Karagiannis et al. [1] built a single pair of gear rattle vibration test rig and measured the displacement and velocity of gear rattle vibration to study the periodic and chaotic response of gear rattle vibration. Wang et al. [2] developed a linear dynamic model as well as a nonlinear subsystem model considering gear backlash and clutch multistage stiffness to study the system torsional vibration and gear rattle, respectively. Bozca et al. [3] studied the influence of module, number of teeth, axial clearance, and backlash on gear rattle, in order to minimize rattle noise. Farshidianfar [4] developed a generalized nonlinear time-varying dynamic model of a spur gear pair and investigated the transition process from global homogeneous bifurcation to chaos in nonlinear gear systems using Melnikov analysis. Ling et al. [5] developed a sixdegree-of-freedom nonlinear dynamic model using the periodic expansion method to obtain the periodic and chaotic response of the system by observing its dynamic trajectory. Yi et al. [6] proposed a nonlinear dynamic model for spur gears considering the effects of pressure angle, time-varying characteristics of gear backlash, unbalanced mass, and internal and external excitation, and studied the influence of these parameters on the nonlinear behavior and the impact state.
Most of the above models consider the factors such as meshing stiffness, gear backlash, and resistance torque. However, the models developed assume that the gears are in a dry friction state and the lubricant is always present in the meshing gears, even if the two solid surfaces are not in contact, the reaction force between the teeth will be generated due to the squeezed oil film effect, thus ignoring the effect of the nonlinear oil film force generated by the inter-tooth lubrication effect on the gear system response. In the study of gear lubrication effects, Brancati et al. [7,8] proposed an analytical model that accounts for the oil squeeze effect between the impacting teeth of mating gears. The influence of the oil-damping effects on the automotive rattling phenomenon. Theodossiades et al. [9] considered the gear impact surface as a lubricated connection and studied the effect of lubrication on light-load idling gear vibration. It was shown that lubricant viscosity in the right range allows smoother transmission of motion between gears. It was shown that lubricant viscosity in the right range allows smoother transmission of motion between gears. Cruz et al. [10] considered the effect of coiling and squeezing of lubricant between teeth and showed that the oil film force is related to the viscosity of the lubricant, the entrainment speed, and the speed of the squeezed oil film between the teeth of the wheel. Liu et al. [11] described the damping and stiffness of intertooth lubricant during gear rattle by numerical calculations.
It is well known that gear meshing excitation is very important to the dynamic performance of gears, on the base of consideration of oil film forces, it has been argued that the magnitude and direction of tooth friction force will change periodically during gear mesh, thus becomes a non-harmonic type of internal excitation and couples with other nonlinear factors such as gear backlash and time-varying stiffness, which makes the gear train exhibit complex nonlinear vibration characteristics. For example, Vaishya et al. [12,13] developed a direct contact model for Coulomb friction between gear teeth and investigated the dynamic effects caused by friction through linear and nonlinear time-varying systems and showed that friction can reduce large oscillations under certain resonance conditions. Chone [14] constructed a direct contact elastic deformation model considering friction and analyzed the effect of friction on the gear vibration response. The results show that the effect of friction on the nonlinear behavior of gears is small and has a positive effect, but the scholar modeled the friction coefficient as a constant and did not consider the timevarying characteristics of the friction coefficient. In the case of gear backlash studies, most of the existing gear backlash models only consider compensation errors and only consider as a constant backlash function. In fact, the backlash of gear teeth changes during the motion due to wear of gear teeth and deformation of support bearings and shafts. The study of the nonlinear dynamic characteristics of dynamic backlash gear systems has a wide range of applications. For example, Li et al. [15] considered that the presence of lubricant changes the relative displacement of gears and investigated the effect of gear backlash on the dynamic behavior of gears by developing a nonlinear mathematical model of gear backlash. In the study of resistive torques, numerous scholars have proposed different calculation models from different perspectives. For example, Lauster [16] fitted an expression for the churning torque to a churning loss test on a commercial vehicle transmission. Changenet [17] used a single pair of spur gears as the test object and carried out extensive tests considering gear parameters, lubricant oil parameters, and operating parameters to fit a new formula for calculating the churning torque. Gorla et al. [18] used a CFD method to investigate the churning power loss for different parameters at high-speed conditions and verified it experimentally. However, most of the existing calculation methods are based on semiempirical formulas obtained through experimental fitting and have a limited scope of application. The CFD-based simulation usually used simplified mesh, which will lead to reduced efficiency and accuracy. There is still a lack of a general and effective theoretical model for calculating the churning resistance torque.
Based on the summary of the above references, although some dynamic models exist to investigate the nonlinear characteristics of gears, most do not take into account the comprehensive effect of influencing parameters, such as drag torque, nonlinear oil film force, lubricant viscosity, time-varying meshing stiffness, friction force, and gear backlash. Therefore, this paper aims at establishing a nonlinear dynamic model of gear rattling that considers the combined effect of above-mentioned parameters, and investigating the nonlinear behavior and severity of gear rattling. The rest of the paper is structured as follows. Section 2 introduces the gear dynamic model which includes the calculation of time-varying mesh stiffness and dragging torque. The influence of time-varying friction coefficient on the dynamical behavior is also discussed and, Section 2 presents the experimental validation. Section 3 presents the parametric analysis for various parameters. The final conclusions are given in Sect. 4.
2 Dynamic modeling of gear system

Gear dynamic model
In order to explore the vibration response and the nonlinear behavior of gear rattle, this paper takes a certain gear of a dual-clutch transmission as the research object. Equations (1)-(6) describe a strongly nonlinear dynamics model with meshing stiffness, meshing damping, nonlinear backlash, friction, and nonlinear oil film force, as shown in Fig. 1. Here C m and K m represent gear meshing damping and timevarying meshing stiffness, respectively, and h p , h g represent the angular displacement of the two gears, the r bp and r bg represent the base circle radius of the main and driven gears, respectively.
According to Newton's second law, the force analysis of the gear system is carried out, and the mathematical model is shown in Eqs. (1)- (3): where x n is the normal relative displacement of helical gear; r b i ð Þ represent the radius of the base circle of the gear, angle markers p and g, respectively, represent the active and driven gear; T mesh i ð Þ is the torque produced by meshing force F mesh and friction force F f during tooth meshing; F ehl represents the meshing force when the gear is in direct contact, F hdl represents tooth surface fluid lubrication force; T ef i ð Þ and T hf i ð Þ represent the friction torques under two contact states, respectively; b b is spiral Angle of gear base circle, I p and I g represent the rotational inertia of the main and driven gears, respectively; € h p and € h g represent the angular acceleration of the two gears, respectively; T e and T d represent driving torque and towing torque, respectively; the driving torque T e can be expressed as: where T m represents the mean torque of driving torque, T A ¼ I 0 Á a 0 , it represents the amplitude of torque fluctuation, and where I 0 is the rotational inertia of the input and a 0 is the angular acceleration of the input, f m is the excitation frequency, f m ¼ n p =60, n p is the driving gear speed.

Gear contact model
The model established in this paper is a dynamic model with lubrication effect considered. Therefore, the gear meshing force F mesh in lubrication state includes two aspects: That is, the dynamic meshing force F ehl generated by direct contact of gear teeth and the fluid lubrication force generated by extrusion and suction of lubricating oil when the tooth surfaces do not contact. At this point, the lubricating oil film between teeth is equivalent to a spring damping element and presents nonlinear oil film force F hdl . The two contact forms are shown in Fig. 2. However, in the process of gear processing, in order to avoid thermal expansion jam, and to produce ideal lubrication, it is necessary to maintain adequate clearance. According to reference [15], backlash 2c b is symmetrically distributed and lubricants exist between the two tooth surfaces and change the relative displacement of gears. Therefore, considering the backlash, the nonlinear deformation function t ð Þ is used to modify x n , and the piecewise function expression is shown in Eq. (6): where h max is the upper limit of the central oil film thickness at which the meshing force begins, c b À h max is the nonlinear backlash affected by h max . The direct contact condition is regarded as a special case of lubrication condition, there is no oil between the two tooth surfaces, i.e.,h max ¼ 0 mm. Therefore, the meshing force F mesh of gears in lubrication state can be expressed by piecewise function as shown in Eq. (7): where x n is the normal relative displacement of helical gear, c b is half of the backlash, r bp and r bg represent the base circle radius of the main and driven gears, respectively, K m is time-varying meshing stiffness of gear, C m is gear meshing damping [19,20] which can be expressed as: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K m r 2 bp r 2 bg I p I g = r 2 bp I p þ r 2 bg I g r ð8Þ where n represents the damping ratio coefficient, which usually ranges from 0.03 to 0.17, and the value in this paper is 0.1.
For the meshing stiffness, there are different calculation methods such as Weber's method [21], ISO standard [22], and potential energy principle [23][24][25]. Among them, the potential energy method has better accuracy and is widely used in many references as an analytical calculation method for gear mesh stiffness, and the gear teeth are usually considered as non-uniform cantilever beams in the calculation, as shown in Fig. 3. Four kinds of stiffnesses are generated during the meshing of spur gears: Hertzian contact stiffness k h and bending stiffness k b , shear stiffness k s and axial compression stiffness k a . The potential energy method solves them jointly to obtain the single-tooth time-varying meshing stiffness k of the gear, as shown in Eq. (9).
Then, slice the helical gear and simplify it into several spur pieces of gears (as shown in Fig. 3). Finally, the time-varying meshing stiffness K m of the helical gear was obtained considering the changes of contact degree and meshing period, which could be expressed by Fourier series as follows: The time-varying meshing stiffness curve of helical gear is shown in Fig. 4, where the dotted line represents the time-varying meshing stiffness of single tooth in meshing cycle, and the solid line represents the comprehensive time-varying meshing stiffness of multiple tooth meshing.
When the tooth surfaces do not contact, the lubricating oil between the two meshing teeth can be regarded as a nonlinear elastic-damping element with time variation. Under the action of tooth surface extrusion and coiling, the fluid lubrication force is produced, which will affect the dynamic response of the gear meshing. Therefore, According to different conditions of contact between two teeth as shown in Fig. 5, the mathematical model of oil film force is summarized as Eqs. (11) and (12), and the calculation process is explained in detail in the literature [26], and only briefly described here, the reader interested in this procedure is referred to [26].
where the C Wedge and C Squeeze are equivalent damping coefficients of coiling and extrusion, respectively; L represents the meshing curve of helical gear, r eq represents the equivalent radius of curvature of gear, h i represents oil film thickness, dh i =dt is the oil film extrusion speed, g 0 represents dynamic viscosity of lubricating oil, d i represents the distance between tooth profiles, h min is to avoid that the calculated value of Eqs. (11) and (12) tends to be infinite when h i becomes null, so a saturation value was adopted, let h i = 10 R a , where R a represents the mean roughness, u is the suction speed of lubricating oil, which can be expressed as: where v p and v g represent the rolling speed of the tooth surface of the main and driven gears. The oil film thickness can be expressed as:  The drag torque T D of the empty gear mainly includes three aspects: the gear oil stirring resistance torque T g , the roller bearing torque T rb , and the oil film shear resistance torque T sy . Since the torque of the gear power transmission is much larger than the drag torque, the bearing gear in this paper ignores the influence of the drag torque, while ignoring the main wind resistance.
In the calculation of the churning torque, many scholars have proposed different calculation models from different perspectives, but most of the existing calculation models are semi-empirical formulas obtained through experimental fitting. Comparing the calculation results of several models as shown in Fig. 6 and Table 1, and the difference between different churning torque calculation models under the same conditions is more obvious, indicating that there is still a lack of a more unified and effective churning torque theoretical calculation model. Guo et al. [27,28] established a theoretical calculation model of gear churning resistance torque from the mechanism of gear dragging torque, using the theory of fluid attachment layer and fluid momentum theorem, taking into account the dynamic change characteristics of gear end velocity, the dynamic change law of fluid attachment layer and the change rate of lubricant volume in the meshing area, and verified the feasibility of the model through tests. The calculation process is explained in detail in the literature [27], and only briefly described here, the reader interested in this procedure is referred to [27].
The churning torque T g is divided into three parts according to the part of the gear that acts with the lubricant in the oil-immersed state: the friction resistance torque T el generated by the friction between the end face of the gear and the lubricating oil, the friction resistance torque T w generated by the friction between the gear peripheral face and the lubricating oil, and the resistance torque T s generated by the extrusion power loss P s of the lubricating oil in the meshing area. In the formula, T el and T w are calculated as follows: where t and q represent the kinematic viscosity and density of lubricating oil, respectively; h represents the distance between the lubricating oil surface and the gear center line; r represents the radius of gear pitch circle, w g represents the angular velocity of gear, B represents tooth width, R represents the circumference of oil immersion. The gear churning model is shown in Fig. 7.
In the process of gear meshing, the teeth of the active and passive gears will periodically invade the gap area at the root of each other's teeth, causing the lubricant to be squeezed in the gear gap, resulting in a rapid decrease in the volume of the gap area and forcing the lubricant between the gear teeth to be squeezed out from the sides and the gap, thereby generating a resistive torque. Through the Boolean method of calculation to obtain each squeeze side area, and then find the backlash in the lubricant volume change and lubricant speed, the use of Bernoulli's principle to calculate the squeeze lubricant power loss, and then calculate the squeeze resistance torque. The calculation process is as follows: where v p ij is the lubricant speed, DV is the volume change, Dt is the change time, B is the gear width, x i is the angular speed of the driven gear, and Dh is the angular change, S p ij representing the area of the jth squeeze side of gear i at meshing position P (as shown in Fig. 8).
The lube oil extrusion power loss is calculated as follows: The squeeze resistance torque at S p ij is calculated by combining in Eqs. (19) and (20) as in Eq. (21) where T p i is the squeeze resistance torque generated by the jth squeeze side area of gear i at engagement position P, and n i is the speed of gear i.
Similarly, calculate the lubricant squeeze resistance torque of other meshing positions, and finally find the average anti-squeeze resistance torque to obtain the average squeeze oil resistance torque of the gear, as shown in Eq. (22), where N indicates the number of squeeze side areas.
Therefore, the churning torque can be expressed as: In the method of reference [30], the bearing resistive torque is shown in Eq. (24) where f 0 is the lubrication factor, which is related to the type of lubricant and lubrication method, and its value can be referred to the standard ISO TR 14,179-2-2001, DN is the speed difference between the driven gear and the shaft;d m is the average diameter of the bearing.
The oil film shear torque is mainly the resistance torque caused by the relative motion between the inner cone surface of the synchronizer lock ring and the outer cone surface of the gear friction when the synchronizer is not combined, and the mathematical expression is: where l is the lubricant dynamic viscosity, L is the length of the tapered surface, R is the average radius of the tapered surface, j c represents the clearance between the tapered surfaces of the synchronizer ring and the synchronizer gear ring.
In summary, the total no-load gear dragging torque can be expressed as:

Friction and friction torque
Due to the influence of lubricating oil and backlash, there are two kinds of friction in the process of gear meshing, namely direct contact Coulomb friction F ef and hydrodynamic sliding friction F hf obtained through semi-Sommerfeld [31] conditions, both of which affect the dynamic characteristics of gear. According to Coulomb's law of friction, the direct contact with friction F ef and its resulting torque T ef can be expressed as: where i ¼ p; g, and r i ð Þ represents the radius of gear pitch circle, a t represents the pressure angle of end face, l t ð Þ is the friction coefficient between teeth in direct contact, and the friction coefficient model proposed in Ref. [32] is shown in Eq. (29): where w ¼ mod n p ; r bg ; t; p b À Á À L Â Ã , and l avg is the time mean value, / is the regularization factor, and / ¼ 50 is set in this paper, p b is the pitch of base circle, e is gear tooth meshing coincide degree, L is the length of meshing contact line. The hydrodynamic sliding friction F hf obtained by the semi-Sommerfeld condition, which is defined as follows: (where u t is the sliding speed of tooth surface: Therefore, the frictional torque acting on the main and driven gears is: In order to determine the effect of the friction coefficient on the dynamic behavior of the gear, this paper simulates three different analytical models based on frictionless contact, constant friction coefficient, and time-varying friction coefficient. The comparison results of dynamic transmission error (DTE) are shown in Fig. 9. The results show that there are some frequency-dependent discrepancies between these three models, but the influence is limited in general, which is basically consistent with the conclusions of the literature [14,27,32]. The timevarying friction coefficient will still be used for the further analysis in this paper.

Model validation
In order to verify the accuracy of the proposed model, the dual-clutch transmission was tested on the threemotor test stand (Fig. 10) to obtain test data. Considering that the coupling between the power flow gear and the no-load or light-load gear is not strong, the whole dual-clutch transmission is selected for the test. The main motor of the stand is directly connected to the transmission using a rigid shaft, so that the secondorder engine excitation simulated by the N/T control mode control motor to be transferred directly to the transmission input shaft, outputting a continuously fluctuating speed, and the two secondary motors are connected to the two half shafts of the transmission to simulate the load torque. Most transmission bench tests collect vibration and noise signals, which makes it difficult to compare with the simulation results. In this study, a magnetoelectric speed sensor is used to measure the gear speed and angular acceleration can be obtained.
The motor parameters and test conditions are shown in Table 2, and under the condition that the simulation conditions are consistent with the test, the influence of torque fluctuation on the system at 1500 rpm is numerically simulated, and the gear angular acceleration is obtained by post-processing the speed information of the no-load gear obtained from the test and simulation, respectively, and comparing the simulation with the experimental analysis results. It should be noted that the test and simulation results are filtered to ensure that 50 Hz is the main characteristic frequency, and the comparison results are shown in Fig. 11.
From Fig. 11, it can be found that the trend of the simulation results and the test results are in good agreement, and it is worth noting that the simulation results are slightly smaller in value than the test results, and there is a certain error in the amplitude. There are two possible reasons for the error: (1) it is related to the vibration of the table and the installation of the gearbox, and as the fluctuation amplitude increases, it causes the vibration of the test bench and the transmission input shaft to increase, resulting in large test results; (2) the dynamics model established is somewhat simplified compared with the actual one. However, the errors on the results are within the acceptable range. Therefore, the simulation results are basically consistent with the results of the bench test, which verifies the validity of the built single-tooth rattle dynamics model, which can be used for further parametric analysis.

Numerical simulations and analysis
The main factors affecting gear rattle include input excitation and torque fluctuation, drag torque,  Table 3, and in this section, we will use the fourth-order Runge-Kutta method to solve the strongly nonlinear dynamical system with the aforementioned factors, and obtain the bifurcation diagram and relative displacement time domain, frequency domain, phase trajectory and Poincare mapping of the system under different influencing factors. Referring to the literature [33], selecting the rattle evaluation index RI as shown in Eq. (32) is selected as an objective criterion to explore the effect   of different parameters on the rattle strength of gears, and the forward difference method is used to analyze the sensitivity of the rattle index, as shown in Eq. (33).
where _ F mesh is the derivative of gear meshing force to time, oRI=oX i represents the rattle sensitivity, and the integral step is 5e-6. The larger the value of RI, the more obvious the change of tooth meshing force, the more intense the tooth impact, and the more serious the rattling phenomenon.

Excitation frequency
Under different excitation frequencies, the dynamic response of the system exhibits rich vibration characteristics, which have a non-negligible impact on the rattle phenomenon of transmission gears. Therefore, in order to analyze the effect of input excitation frequency on the rattle nonlinear characteristics of the gear vice, the response results at f m of 10, 25, 50, and 60 Hz are compared without loss of generality, and the rattling strength at different input excitation frequencies is visually compared by the rattle index RI, and the simulation analysis results have typical nonlinear dynamic characteristics. In order to highlight the comparability, the numerical calculation is set with a gear backlash of 0.12 mm, a drag torque of 0.1 Nm, an input angular acceleration of 300 rad=s 2 .
As shown in Figs. 12 and 13, when f m is 10 Hz and 25 Hz, the wheel teeth are impacted on both the driven and non-driven sides, the relative displacement x n oscillates between the gear backlash, and the gears produce double-sided rattle, but the contact between the wheel teeth and the non-driven side is relatively less when f m is 25 Hz. The Fourier spectrum diagram shows that the input excitation frequency component has the largest contribution to x n , and the spectrum at 25 Hz is a continuous spectrum with a certain width; the system appears chaotic response in both cases, the time course of the response is a non-periodic response, the phase trajectory is a closed curve with a certain width, the Poincare cross section is distributed in multiple discrete points, and the motion is chaotic motion, which can be It can be seen that the system vibration response increases sharply after the chaos increases.
When f m is 50 Hz, as shown in Fig. 14, the system is a periodic response, and the time courses of the responses are all periodic. It can be seen by the spectrogram that the Fourier spectrum is distributed over discrete points of kf m =2 (k is a positive integer and f m is the excitation frequency), and the periodic solution of the system is quite rich in super-harmonic and subharmonic components, and the phase trajectory is a non-elliptical closed curve. Poincare cross section presents 2 points distributions, respectively, and the system motion is periodic. And when f m is 60 Hz, as shown in Fig. 15, the system appears chaotic response, the time course is non-periodic motion, the Poincare cross section presents multiple discrete point distributions.
Comparing the system response results generated by different excitation frequencies, it shows that the impact frequency of tooth face changes with the increase in excitation frequency; in each cycle under low excitation frequency, multiple impacts occur on the driving face as well as on the non-driving face, and in each cycle under higher excitation frequency, only one impact occurs on the driving face as well as on the non-driving face. It indicates that the excitation frequency is an important factor that causes the tooth surfaces to detach from each other and produce rattle; looking at the whole process of the system response, its motion state experiences chaotic to periodic and chaotic motion, and the nonlinear performance of the rattle dynamics of the gear pair coupled with multiple nonlinear strong factors is obvious. The rattle index RI and rattle sensitivity of the system when the input excitation frequency increases from 0 to 90 Hz are also calculated as shown in Fig. 16. The rattle index shows a monotonically increasing trend, indicating that the higher the input excitation frequency, the higher the high-frequency impact strength increases and the gear rattle intensifies, and the rattle sensitivity is greater than 0, indicating that the RI is positively correlated with the excitation frequency. Therefore, choosing a  suitable input excitation frequency can prevent the system from entering a chaotic state, which helps to reduce the gear rattle strength, and improve the stability and NVH performance of the system.

Input torque fluctuation
When there are periodic fluctuations of input torque, it may lead to different degrees of rattle of the gear pair.  Therefore, in order to analyze the effect of the amplitude of input torque fluctuation on the rattle nonlinear characteristics of the gear pair, the response results at a 0 of 0, 100, 250, and 320 rad/s 2 are compared without loss of generality, and the rattle index RI is used to visually compare the rattle strength at different input torque fluctuation The results of the simulation analysis have typical nonlinear dynamic characteristics. To highlight the comparability, the backlash of 0.12 mm, drag torque of 0.1Nm and dynamic viscosity of 0.01 pa s are set for the numerical calculation.
As shown in Fig. 17, during the change of input angular acceleration amplitude, the whole system motion is chaotic motion. In the region of chaotic motion, there exists periodic motion, i.e., chaotic motion alternates with periodic motion, and there exists a window of periodic motion in chaotic motion. For example, the system exhibits periodic motion in a small region of angular acceleration amplitude of 218-300 rad/s 2 , and these periodic motions are a window of periodic motion in chaotic motion. In terms of the number of periodic solutions, as the angular acceleration increases, the system starts to bifurcate from a single cycle after excitation to a multiplicative cycle, with single-cycle, 3-cycle motion, etc., and the system amplitude decreases significantly. After a 0 is 300 rad/ s 2 the system starts to decay until the generation of chaotic motion.
In Fig. 18, when a 0 is 0 rad/s 2 , the system is a quasi-periodic response, the time course of the response is non-periodic response, due to the time- varying characteristics of the gear mesh stiffness, nonlinear oil film force leads to the existence of oscillations in the relative displacement, but x n is always greater than 0, the amount of oscillation is very small, indicating that the gear is in a normal meshing state in this case, the teeth are only in contact on the drive side, there is no rattle produced. The Fourier spectrum is distributed at discrete points at 425 Hz (It refers to the gear meshing frequency, f ¼ n p =60 Á z 1 ¼ 425), indicating that the system response contains a super-harmonic response. When a 0 is 100 rad/s 2 , as shown in Fig. 19, the time course of the system response is non-periodic response, and the gear oscillates more obviously with respect to the displacement x n , and there is a more obvious phenomenon of tooth disengagement, at this time, the impact-disengagement cycle occurs repeatedly on the driving side of the gear teeth, and there is no contact with the non-driving side, which is shown as singlesided rattle, and the system motion is chaotic.
When a 0 is 250 rad/s 2 , the system response is shown in Fig. 20, and the system is a times-periodic response, and there are high-frequency shocks in the gear wheel teeth on the driving side and the nondriving side, and in one excitation cycle, the gear first produces a series of shocks on the driving side, and under the action of inter-tooth lubricant, the shocks gradually decay to Then transition to the non-driven side, and then produce a series of shocks with gradually decreasing strength, at this time, the relative displacement x n oscillates between the gear backlash, and the gear produces double-sided rattle. The system phase trajectory is a closed curve with a certain width, and the Poincare section is distributed in three discrete points, and the system motion is periodic. Continuing to increase the angular acceleration amplitude, the system response is shown in Fig. 21, and the system again appears chaotic response. Throughout the whole process of the system response, its motion state goes through the process of period, bifurcation, and chaos, and comparing the dynamic response in these previous cases shows that the input fluctuation is an important factor that causes the tooth surfaces to detach from each other and produce rattle.
The rattle index RI and rattle sensitivity of the system when a 0 is increased from 0 to 1000 rad/s 2 are also calculated as shown in Fig. 22. The rattle index shows a monotonically increasing trend, with the increase in a 0 , the vibration amplitude of the system gradually increases and the system torque fluctuation intensifies, making rattle more likely to occur, and the gear rattle sensitivity is greater than 0, indicating that RI is positively correlated with a 0 . Therefore, appropriately reducing the input torque fluctuation amplitude can prevent the system from entering the chaotic state, and help reduce or suppress the gear rattle strength to improve the stability and NVH performance of the system.

Drag torque
The drag torque T D has the effect of impeding the rotation of the hollow-sleeve gear and can weaken the  tendency of tooth surface separation, which is one of the important factors affecting the occurrence of rattle, and the effect of different drag torque on transmission gear rattle will not be negligible. Therefore, without loss of generality, the bifurcation of the system with the variation of the drag torque is calculated, and the response results are compared for T D of 0.1, 0.3, and 0.44 Nm, and the rattle strength under different input torque fluctuation amplitudes is visually compared by the rattle indicator RI. To highlight the comparability, the numerical calculation was set with a backlash of 0.12 mm, an angular frequency of torque fluctuation of 50 Hz (i.e., w = 100 p rad/s), an input angular acceleration of 300 rad/s, and a kinetic viscosity of 0.01 pa s. Figure 23 shows that in the process of towing torque change, the system motion changes from periodic motion to chaotic motion, and the system exhibits periodic motion in the small region of the dragging torque of 0.08-0.16 Nm, these periodic motions are a window of periods in chaotic motion. After T D is 0.16 Nm the system starts to decay until the generation of chaotic motion. When T D is greater than 0.44 Nm, the relative displacement is always greater than the gap value, and the teeth are only in contact on the drive side, and the system tends to stabilize. The overall system with the change of the dragging torque, its periodic motion is not obvious, will be faster transition to the quasi-periodic and chaotic state, in the large dragging torque, the tendency of the empty set gear tooth surface to occur separation is further weakened, and finally the system into a stable single-cycle state. As shown in Fig. 24, when T D is 0.1Nm, the response of the system is a non-simple harmonic periodic second harmonic response, the relative displacement x n oscillates between the backlash, and the wheel teeth are impacted on both the driven and nondriven sides, which behaves as two-sided rattle. The Fourier spectrum is distributed on the discrete points of kf m =2, and the frequency component of 50 Hz has the largest contribution to x n . Poincare cross section is distributed in two points, and the system motion is periodic.
When T D is 0.3 Nm, the system response is shown in Fig. 25, and the chaotic response starts to appear, the time course of the response is non-periodic motion, the gear teeth are impacted on both the driving and non-driving sides, and there is the phenomenon of dislocation, and the meshing process is mixed with one-sided rattle (driving side), because the increase in the dragging torque makes the gear teeth not easy to separate, and the gear teeth are rattled more on the driving side, which reduces the degree of double-sided rattle of the gear vice, but the overall performance is double-sided rattle. When T D is 0.44 Nm, as shown in Fig. 26, the system responds to the quasi-periodic motion, the relative displacement of gears is greater than the clearance value, the gears are in the normal meshing state, and no rattle occurs on the gears. It is found that with the increase in the dragging torque, the gears gradually transition from double-sided rattle to single-sided rattle, and finally stabilize and no more rattle occurs.
The rattle index RI and rattle sensitivity of the system when T D increases from 0 to 0.8 Nm are also calculated as shown in Fig. 27. The rattle index shows a monotonically decreasing trend, with the increase in T D , the vibration amplitude of the system gradually increases, and the teeth are gradually difficult to separate, which makes rattle less likely to occur, and the rattle sensitivity is less than 0, which means that RI is negatively correlated with T D . Therefore, the appropriate dragging torque is helpful to prevent the system from entering the chaotic state, and at the same time, it can effectively reduce or inhibit the occurrence of gear rattle to improve the stability and life of the system.

Dynamic viscosity
When the tooth surfaces are not in contact, the lubricant in the gear backlash can be equated to an elastic-damping element, and the change of lubricant viscosity can be regarded as the change of the stiffness   the rattle strength at different lubricant viscosities is visually compared by rattle index RI, and the simulation analysis results have typical nonlinear dynamic characteristics. It should be noted that the dynamic viscosity also affects the magnitude of the dragging torque, and in order to compare the effect of different dynamic viscosities on the system, therefore, the effect of viscosity on the dragging torque is not considered in this section.
As shown in Fig. 28, during the change in the dynamic viscosity of the lubricant, the system motion evolves from the initial chaotic state to periodic motion and then to chaotic motion, and in the region of chaotic motion, there exists periodic motion, that is, chaotic motion alternates with periodic motion, and there exists a window of periodic motion in chaotic motion. The system exhibits periodic motion in the region to the left of g 0 of 0.0219 and in the small segment of 0.031-0.0318. These periodic motions are a periodic window in the chaotic motion, and there are many such periodic windows throughout the process of dynamic viscosity change. From the number of periodic solutions, as the kinetic viscosity increases, the system starts to bifurcate from a single cycle to a multiple cycle after the excitation, and there are 3cycle and 6-cycle motions, which further show the rich and diverse bifurcation characteristics of the system, and when g 0 is 0.0318, the system reenters the chaotic state, and after g 0 is 0.0336, the system starts to reverse the bifurcation to a stable single-cycle state.
As shown in Figs. 29 and 30, when the lubricant viscosity g 0 is 0.01 pa s, the system has a chaotic response, the time course of the response is nonperiodic motion, the wheel teeth are impacted on both the driving and non-driving sides, the relative displacement x n oscillates between the backlash, the system mainly occurs double-sided rattle, the Poincare cross section is distributed by several discrete points, when g 0 is 0.024 pa s, the system response is times-periodic motion, the phase trajectory is a closed curve with a certain width, and the Poincare section is distributed by six discrete points. In Fig. 31, when g 0 is 0.03 pa s, the system is chaotic response, the time course of the response is    time step used for the numerical calculation, which is 5 9 10 6 in this study.
Numerical calculation of the three viscosities of R dte and O dte comparison results are shown in Fig. 32, while combined with the rattle index RI and rattle sensitivity as shown in Fig. 33 can be found, the greater the viscosity of the lubricant, the stronger the damping effect of the squeeze film, the curve overlap range and peak with the increase in viscosity and decrease, the more rapid the decay of high-frequency impact strength, and RI is monotonically decreasing trend, while rattle sensitivity is greater than 0, RI and g 0 is positively correlated. It shows that lubricant viscosity has a great influence on the dynamic characteristics of gears. Larger dynamic viscosity accelerates the decay of impact force, but excessive dynamic viscosity affects mechanical efficiency, therefore, appropriately increasing the dynamic viscosity can effectively prevent the system from entering a chaotic state, which helps to reduce or inhibit gear rattle strength, and improve the stability and NVH performance of the system.

Summary
In this paper, a nonlinear gear rattle model is proposed, which considers time-varying friction forces and other strong nonlinear factors such as time-varying meshing stiffness, meshing damping, gear backlash, and oil film force. The proposed model is verified by threemotor drive test. The nonlinear dynamic characteristics are studied by analyzing the time history of response, phase trajectory diagram, Poincare cross section, and Fourier spectrum, and based on the gear rattle index, the influence of each parameter is analyzed, and the following conclusions are given: acceleration is in good consistency with the test results, and the numerical error is within the acceptable range. The model can accurately reflect the rattle phenomenon of the gear with good accuracy. 2. As the lubricant viscosity increases, the system shows periodic, bifurcation and chaotic behavior alternately, At the same time, the larger the viscosity, the greater the equivalent damping of the system, the high-frequency shock intensity decays more rapidly due to large drag torque, which can effectively suppress or reduce the occurrence of gear rattle. Therefore, appropriately increasing the dynamic viscosity can effectively prevent the system from entering the chaotic state, reduce the gear rattle strength. 3. With the increase in the torque fluctuation amplitude a 0 , the system has experienced the evolution process of quasi-periodic, bifurcation, period doubling, and chaos, and the system exhibits rich and diverse bifurcation characteristics. Gear rattle index increases with the amplitude of input angular acceleration fluctuation, and reducing the input torque fluctuation of gear pair, which has a positive effect on controlling the occurrence of gear rattle, and can prevent the system from entering a chaotic state. 4. The input excitation frequency and the drag torque likewise have a great influence on the nonlinear behavior of the system, and the system undergoes both periodic and chaotic motions as the parameters change. However, with large drag torque, the tendency of tooth separation is further attenuated and the system eventually enters a stable singlecycle state. Therefore, a larger drag torque is beneficial to prevent the system from entering the chaotic state, which can effectively reduce or inhibit the occurrence of rattle.