Modelling of time-varying meshing stiffness under tooth wear condition for nonlinear dynamics in compound planetary gear set

Under tooth profile wear condition, with running time increasing, the tooth chord length of the gear tooth is thinner and the tooth profile is not involute profile any more, which results in the changes in amplitude and fluctuation of time-varying meshing stiffness. In this work, in order to reveal the influence of tooth profile wear on time-varying meshing stiffness, an analytical model for calculating the time-varying meshing stiffness under tooth wear condition is built and incorporated with a “translational-rotational” dynamic model of compound planetary gear set. Poincare section, the phase trajectory and bifurcation diagram are employed for analysis of periodic and chaotic motion. The results indicate that, tooth wear accumulation result in the reduction in the amplitude of meshing stiffness at the meshing region of double-pair teeth, but makes no difference in meshing stiffness for the meshing region of single-pair teeth. In order to avoid the system being in Chaos motions, under large tooth wear condition, the input speed should be avoided away from 3700 r/min-6725 r/min, so as to keep the system in a stable periodic motion state.


Introduction
As a strongly nonlinear excitation in gear transmission, time-varying meshing stiffness is a periodic excitation and contributes greatly in the periodic response in gear transmissions. A little change in either amplitude or fluctuation in time-varying meshing stiffness could result in obvious difference in dynamic response, especially in complex gear transmission, such as compound planetary gear set.

Modeling of meshing stiffness
Considering tooth profile modification and tooth root crack, lots of literature described the modeling of meshing stiffness for spur gears, helical gears, bevel gears and so on. Considering the tooth deflections, but neglecting the hertzian deflections, Sanchez [1] evaluated the meshing stiffness of spur gear pairs at any point of the path of contact and approximated by an analytical model. Sun [2] developed a revised time-varying mesh stiffness model for of spur gear pairs, and investigated the impact of the tooth width and torque on mesh stiffness. Assuming that the extended tooth contact is ignored, Ma [3] and Parey [4] established an improved analytical method suitable for gear pairs with tip relief, to determine time-varying mesh stiffness for spur gears. Cui [5] incorporated a finite element model of the spur gear undergoing damaged failure with a simulation calculation program to calculate the time-varying meshing stiffness. By defining the effective thickness of the gear teeth with the occurrence of cracks, a cracked gear-tooth meshing stiffness calculation method is proposed by Cui [6]. Chaari [7] derived an analytical formulation of the time varying mesh stiffness for spur gear tooth crack. Chen [8,9] proposed an analytical mesh stiffness calculation model for non-uniformly distributed tooth root crack along tooth width. By taking the TVMS obtained from the finite element (FE) method as a benchmark, Huangfu [10] obtained the correction coefficient of the gear foundation stiffness of cracked helical gear pairs by the optimization method. By using a combination of finite element method (FEM) and local contact analysis of elastic bodies, Chang and Chen [11][12][13] proposed a model for determining mesh stiffness of cylindrical gears. By comparing results with those obtained by using two-dimensional (2D) finite element (FE) models and specific benchmark software codes [14], Gu [15] gave the time-varying mesh stiffness function for ideal solid spur and helical gears. In order to introduced meshing stiffness into lumped-parameter model for gear transmission, Chen and Wang [16][17][18][19] introduced spring damping pair to describe the meshing stiffness.
According to the literature above, although lots of work is proceeding on modeling of time-varying meshing stiffness under the condition of tooth modifications and tooth crack for spur gear, helical gear, cylindrical gears, the modeling of the meshing stiffness under tooth profile wear condition is absent. In this work, an analytical model for calculating the time-varying meshing stiffness under tooth wear condition is proposed.

Impact of parameter on meshing stiffness
Certain literature focused on the effects of parameter excitation on meshing stiffness. By using three-dimensional boundary element analysis, Saryazdi [20] found that meshing stiffness reduces as crack propagated, but the amount of reduction depended on the position of the contact line and crack propagation path. Sanchez [21] investigated the influence of profile modifications on meshing stiffness, considering both the symmetric and asymmetric reliefs. Huangfu [22] discussed the effects of the tip relief and surface wear on the meshing and the dynamic characteristics. In this work, the impact of tooth profile wear on time-varying meshing stiffness is revealed by using analytical model for meshing stiffness and numerical method for calculation.

Impact of meshing stiffness on dynamics
Works on the impact of meshing stiffness on dynamic behavior have been proceeding on for gear transmission. Lin and Tordion [23,24] identified analytically effects of mesh stiffness parameters, including stiffness variation amplitudes, mesh frequencies, contact ratios, and mesh phasing, on these instabilities in two-stage gear systems. Cui [25,26] introduced a universal gear profile equation reference to the actual manufacturing process to calculate the meshing stiffness and revealed that, the double pulse was generated due to the residual stiffness. Hbaieb [27] listed results of plagiarism investigation on dynamic stability of a planetary gear train under the influence of variable meshing stiffness.
In this work, by incorporating the modeling of the time-varying meshing stiffness under tooth profile wear with a "lumped-parameter" dynamic model of compound planetary gear set, the influence of the meshing stiffness reduced by tooth wear on dynamic behavior of compound planetary gear is revealed, by incorporating analytical model of meshing stiffness and kinematic partial differential equations of compound planetary gear set.

Impact of tooth wear on dynamics
The influence of tooth profile wear on dynamics of gear transmission is also researched before. Zhang [28,29] investigated the influence of backlash generated by tooth profile wear on dynamics of compound planetary gear set, and revealed the coupling relation between tooth wear accumulation and load sharing behavior in compound planetary gear set. By incorporating a torsional dynamic model with a surface wear model, Kahraman [30] also revealed the two-way relationship between nonlinear planetary gear dynamics and tooth surface wear. Shen [31] employed the potential energy method to calculate the mesh stiffness, and adopted the Archard's wear equation to calculate the tooth profile wear, then demonstrated that the tooth wear can result in the reduction of the mesh stiffness [32]. In this work, the influence of time-varying stiffness reduced by tooth profile wear is investigated.
According to literature listed above, lots of work focused on modeling of meshing stiffness under the condition of tooth crack or modification, influence of meshing stiffness on dynamics and influence of wear on backlash. However, the modeling of meshing stiffness under the tooth wear condition is ignored, as well as the impact of meshing stiffness reduced by tooth wear on gear transmission.

Objectives and scope
In this work, "single point observation method" proposed by Flodin [33] is introduced to Weber-Vanaschck [34] model of meshing stiffness to predict the time-varying meshing stiffness along the tooth profile under tooth wear condition. By changing the tooth wear depth, the effects of tooth wear on time-varying meshing stiffness is revealed. Then, the analytical model of the time-varying meshing stiffness under tooth wear condition is introduced into a "translationalrotational" dynamic model [29,35] of compound planetary gear set, so that the impact of meshing stiffness reduced by tooth wear accumulation is investigated by calculating the steady response.

Modeling of time-varying meshing stiffness under tooth wear condition
In order to reflect the influence of tooth profile wear on time-varying meshing stiffness, the Weber-Vanaschck [34] model is employed and incorporated with a static gear wear model proposed by Flodin and Andersson [33], to calculate the changing of meshing stiffness with time increasing.
Firstly, the wear depth distribution along tooth profile is given by Zhang [28] by using "single point observation method", shown as below: where subscript represents the th observation point along tooth profile. In this work, the amount of observation point is 500, which results in a good precision. Subscript represents the driving gear running for rounds. The literature [28,29] gives the definition for other parameters in Eq. (1). ℎ , determines the wear depth distribution along tooth profile.

An analytical model for meshing stiffness under tooth wear condition
In this chapter, two essential questions are discussed and revealed: one is the influence of tooth chord length reduced by tooth profile wear on time-varying meshing stiffness, the other is method of introducing the time-varying meshing stiffness under tooth wear condition into the lumped-parameter dynamic model of compound planetary gear set. The first question is revealed in the modeling of the time-varying meshing stiffness under tooth wear condition, shown in Fig. 1, where , ,

and
, , represent the tooth chord length at th observation point under initial condition of non-wear and the condition after the driving gear running for rounds respectively, and subscript "1" represents the driving gear, which are calculated as: Tooth chord length after first iteration for wear calculation: Tooth chord length after first iteration for wear calculation: Based on assumption that the wear condition of different teeth is same, the gear wear in Fig. 1 is bilaterally symmetrical, and the gear wear in Fig. 1 is similar to the actual wear in planetary gears by Kahraman [30]. Where in Eq. (2), ℎ , , ⋅ sin , is the projection length of tooth wear depth from the th observation point after one wear iteration, on the driving gear. ℎ , , is given in Eq. (1), and , and , is derived in Fig. 2. The beta parameters in Fig. 2 is defined in literature [28], which is described above the Eq. (1).

Fig. 2. Diagram of tooth profile
Where in Fig. 2, is the dedendum radius of driving gear, is the radius of base circle, , is the distance between meshing point and center of rotation. When a pair of mating gears begins to mesh, the driving gear comes into engagement on the point at mesh beginning, shown in Fig. 2. The radius of the point at mesh beginning is . When the driving gear begins to separate from the driven gear on the point at mesh ending, the radius of the point at mesh ending is . One meshing section is the section between the point at mesh beginning and meshing ending on the tooth profile.
is the projection length of the meshing section from tooth profile on meshing line, and is divided into 500 equal pieces to set the observation point for tooth profile wear. The calculation for , and is given as below: Due to nonuniform distribution of the tooth wear depth along tooth profile, in the mesh section, the wear depth ℎ , , at th observation point changes with the position of observation point and tooth wear iteration. In order to reveal the influence of tooth profile wear on meshing stiffness, the tooth chord length under tooth wear condition is given as below: , , where in Eqs. (8)(9)(10), , is the pressure angle of meshing point, or the th observation point on the driving gear, , is the angle between contact force and horizontal direction for meshing point. By introducing Eqs. (10-11) into Eqs. (2-3), the tooth chord length at th observation point after the driving gear running for rounds is calculated. In order to calculate the elastic deformation at meshing point, the tooth chord length of driving gear at dedendum and the distance between the meshing point and horizontal line , is given as below: After the driving gear running for rounds and tooth profile wear accumulation, the meshing deflection consists of tooth deflection , , and matrix deformation , from driving gear, tooth deflection , , and matrix deformation , from driven gear, as well as local contact deformation , : , , = * cos , 10.92 where in Eqs. (15) and (17) and , , respectively, which results in a larger deformation for the tooth deflection of mating gears. Then, the local contact deformation is calculated as below: where in Eq. (18), yℎ , is the distance between meshing point and the pitch point on driving gear, and is calculated in Eq. (19). In consideration of difficulty in convergence of numerical method, ℎ , not changes with tooth wear accumulation, as well as the equivalent radius of curvature * for the th observation point. The time varying meshing stiffness , is calculated as below, and represents the meshing stiffness for the mating gears at the th observation point after the driving gear running for rounds:

Influence of tooth profile wear on time-varying meshing stiffness
In this work, in order to investigate the impact of tooth wear on meshing stiffness, a group of parameters from Eqs.  is assigned in Table 1-2. Influence of tooth profile wear on time-varying meshing stiffness is given as below in Fig. 3. The tooth chord length under tooth wear condition in Fig. 3 is calculated by Eqs. (1-3). By Eqs. (4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17), the time-varying meshing stiffness under the worn tooth chord length after the driving gear running for rounds is calculated by numerical method and displayed in Fig. 3.
As shown in Fig. 3, for the meshing pair formed by a pair of parallel shaft gears, the variation of the variable meshing stiffness in a meshing period varies with tooth profile wear. With the increase of number of revolutions of driving gear, the increasing in gear profile wear reduces the meshing stiffness amplitude of the meshing pair in the meshing region of double teeth meshing region where in Fig. 3 is labeled as "A", and has little influence on the stiffness fluctuation of single tooth meshing region, where in Fig. 3 is labeled as "B". It is indicated that, with the increasing in tooth profile wear, the meshing stiffness of double teeth decreases obviously, and both the stiffness amplitude and fluctuation in single tooth meshing region basically remains unchanged, leading to a decrease in the mean value of time-varying meshing stiffness.  The main cause of this phenomenon is that, the tooth wear depth is higher near the addendum and dedendum than that near the pitch circle [33], resulting in smaller tooth chord length when the mating gears starts and exits meshing than when the mating gears runs in single tooth meshing region. In view of the difference on the tooth accumulation of double teeth and single teeth meshing region, the meshing stiffness of double teeth meshing region decreases with the accumulation of tooth profile wear, and the impact of tooth profile wear on meshing stiffness at single tooth meshing range is not obvious, so the meshing stiffness basically remain unchanged at region "B". Generally speaking, the mean value of time-varying meshing stiffness decreases gradually with the accumulation of tooth profile wear, as well as the fluctuation peak of meshing stiffness.

Model for compound planetary gear set under tooth profile wear
In view that time meshing stiffness is a strongly nonlinear excitation in kinematical equation for gear transmission, especially for complex gear transmission, in this work, a lumped-parameter compound planetary dynamic model is employed and incorporated with the analytical model of meshing stiffness under tooth wear condition. The dynamic model of a compound planetary gear set under tooth wear condition is shown in Fig. 4. Known from the lumped-parameter dynamic model proposed by Zhang [36], in this work, the meshing stiffness of meshing pair , ( = 1 , 2 , , ) changes with the driving gear running for different rounds. The definition of other parameters in Fig. 4 is listed in Table 3.
Due to tooth profile wear increasing, the initial mean value of mesh stiffness ( , , , , , , , ) and coefficient of meshing stiffness fluctuation , decrease with the driving gear running for rounds. By extracting the first harmonic of Fourier transform for the time-varying meshing stiffness after the driving gear running for rounds in Fig. 3, mean value of mesh stiffness , ( = 1 , 2 , , ) and coefficient of meshing stiffness fluctuation , is given as below. By incorporating , and , with the dynamic model in Fig. 4, influence of meshing stiffness under tooth wear condition on dynamics of compound planetary gear set is investigated: where is meshing frequency, ( ) is the relative meshing angle between mating gears.
By increasing the number of revolution of driving gear 1, the tooth profile wear is calculated by Eq. (1) and the time-varying meshing stiffness is calculated by Eqs. . Then the dynamic response in the steady state time domain is solved by fourth-order Runge Kutta method for kinematical Eqs. (24)(25)(26)(27)(28)(29) of the compound planetary gear set.
Planet carrier : Ring : where ( = 1 , 2 , , ) represents the actual penetration depth after the theoretical penetration depth removing the initial backlash .
is the meshing force of meshing pair . Different from the kinematical equation equations in literature [29,30], in this work, the impact of tooth wear on backlash is ignored. The impact of meshing stiffness under tooth wear condition is emphasised on in this dynamic model, and the description of each parameter in Eqs. (24)(25)(26)(27)(28)(29)(30) is listed below.

Chaos and bifurcation
The fluctuation coefficient of meshing stiffness makes key roles in periodic motion of vibration, In order to investigate the influence of meshing stiffness under tooth wear condition on periodic vibration of compound planetary gear set, the diagram of Chaos and Bifurcation under different tooth wear condition is given as below. The max tooth wear depth usually appears on dedendum and addendum and is denoted as  Fig. 5(a-c). In this frequency section, the vibration double period and quasi-period motion occupies a leading position, especially at the section = [0.8-0.1.2], resulting from the influence of running speed increasing on transmission error. With tooth wear increasing, in Fig. 5(a-b), the number of double periods increase slightly at = 0.55, and vibration amplitude increases at = 0.68, which indicates that the tooth wear increasing and time-varying meshing stiffness decreasing cause a larger vibration amplitude. Although a larger amplitude occurs with tooth wear increasing, more regular quasi-periodic and double periodic motion presents, comparing non-wear Fig. 5(a) with tooth wear increasing in Fig. 5(b-c). It is indicated that, with the tooth profile wear increasing, the fluctuation coefficient of meshing stiffness decreasing cause a lower fluctuation of vibration and more regular vibration in periodic motion.
With the running speed and meshing frequency increasing to high frequency, in the section = [1.2-2], Chaos motion occurs in Fig. 5(a). The introduction of tooth wear narrows the Chaos motions section from = [1-1.3] to [1.1-1.2], as well as the scope of vibration amplitude.
The analysis results for Fig. 5 are able to be explained, by illustrating the influence of excitation frequency and meshing stiffness on shock in mesh and vibration. When the gear transmission runs in a low speed, the frequency of internal excitation is also low. Because the periodic motion of vibration is largely depended on time meshing stiffness, the introduction of tooth wear accumulation makes key roles on periodic motion at = [0.4-1.2]. With the meshing frequency increasing, the influence of meshing stiffness decreasing due to tooth profile wear on periodic motion is less and less, compared with the impact of running speed and transmission error on vibration. So the decreasing meshing stiffness due to tooth wear accumulation has little affection on high frequency motion characteristics.

The phase trajectory and Poincare interface
At the low frequency section < 0.4, in Fig. 6(a-c), for one observation point on meshing pair , the track of the vibration at different meshing period is presented as a closed curve. In order to confirm the cycle-index of the closed curve, at the Poincare interface Fig. 6(d-e), one-cycle periodic motion is seen clearly. With tooth wear accumulation, the vibration of one-cycle periodic motion continues, which indicates that the influence of meshing stiffness decreased by tooth wear accumulation is little. That may be explained that, one-cycle periodic motion leads at low running speed section. , the track of the vibration at different meshing period overlaps as a closed curve. With tooth wear increasing and meshing stiffness decreasing, in Fig. 7(b-c), vibration track area increases, resulting a larger vibration displacement and velocity. In order to confirm the cycle-index of the closed curve Fig. 7(a-c), at the Poincare interface Fig. 7(d), many points gather together at non tooth-condition. With tooth wear increasing in Fig. 7(e-f), Two discrete points occurs, which indicates that the meshing stiffness reduction caused by tooth wear accumulation improves the periodic motion from quasi-period to doubling-period, at the medium frequency section. , the different tracks of the vibration at different meshing period overlaps closed. With tooth wear accumulation and meshing stiffness decreasing in Fig. 8(b-c), the number of closed curves decreases sharply, although chaotic motion occurs. At the Poincare interface Fig. 8(d), the points distributes chaotic. With tooth wear increasing in Fig. 8(e-f), the number of points decreases sharply and distributes more regularly, which indicates that the meshing stiffness reduction caused by tooth wear accumulation improves the periodic motion from chaotic motion to quasi-period, at the high frequency section.

Load sharing behavior
Due to tooth wear increasing on tooth profile, the fluctuation and amplitude of meshing stiffness decreasing causes a smaller fluctuation and a larger amplitude of the vibration on gears in compound planetary gear set, resulting in a different load distribution for the planet gears surround the sun gear [37]. In this chapter, the influence of the meshing stiffness reduced by tooth wear accumulation on the load sharing behavior of compound planetary gear set is investigated by calculating the load sharing coefficient cited in [29], under the condition of different tooth wear depth on tooth profile. In order to reflecting the partial load phenomenon in compound planetary gear set, an eccentricity error is applied on sun gear 1, which is the input gear for the gear transmission.
In order to reveal the influence of the time-varying meshing stiffness decreased by tooth wear accumulation on load distribution behavior at certain time, an transient load sharing coefficient is introduced as the function of time ( ) in Eq. (31), which represents the load distribution behavior at certain discrete times: To specify the influence of the time-varying meshing stiffness decreased by tooth wear accumulation on load distribution behavior, a set of different maximum amplitudes (ℎ , , ) of the tooth profile wear on driving gear 1, are introduced into the equations of meshing stiffness in Eqs.  for the calculation. According to a comparison of Fig. 9(a) and (b), a stronger fluctuation in load sharing coefficient appears at planet-planet meshing pair , compared with the load sharing coefficient on meshing pair . Similarly, with a comparison of Fig. 9(b) and 9(c), a larger load sharing coefficient appears on the meshing pair of output stage. That is indicated that the unbalance loading distribution is magnified with the transfer of power from the input ( 1) to the output ( 2). With tooth wear increasing, the decreasing amplitude and fluctuation of time-varying stiffness enlarge the load sharing coefficient and improve the load sharing behavior for compound planetary gear set.

Comparison between the results of this work and other literature
This paper does focus on theoretical analysis and lacks a large number of experimental data to verify the correctness of theoretical results. However, part of the results in chapter 1 and chapter 3 of this paper can be verified by means of literature comparison.
In Fig. 1, the worn tooth profile shape calculated based on 500 discrete points in this work conforms to the worn tooth surface shape obtained by Achard wear theory and single-point observation method. In other words, the wear of tooth tip and tooth root is far greater than that of pitch node, which is verified by [33].
In addition, the influence of tooth wear on meshing stiffness calculated in Fig. 2 shows that, the accumulation of wear on tooth surface reduces the amplitude of meshing stiffness. This result is consistent with the result of reference [31,32], which can prove the correctness of the calculation method of meshing stiffness under tooth wear condition theoretically derived in this JOURNAL OF VIBROENGINEERING work.
Eqs. (23)(24)(25)(26)(27)(28) describe the compound planetary gear kinematics equation derived from Lagrange equation. The equation deriving has been proven in [29,30,36]. The innovation of this paper is that the tooth surface wear is introduced into the kinematics equation in the form of numerical excitation, and the change of meshing stiffness and system dynamic characteristics caused by tooth surface wear is calculated and analyzed. In Fig. 9, due to the transmission error of the input gear 1 ( Table 3 describes the value of the transmission error on input gear 1), partial load distribution occurs in the system, that is, the uniform load coefficient is greater than 1. This is in line with the record in the literature [29,36] that the transfer error is the fundamental factor that affects the uniform loading.

Conclusions
In this work, a numeric model to calculate the time-varying meshing stiffness decreased by tooth wear accumulation on compound planetary gear set is built and incorporated with a "translational-rotational" lumped parameter dynamic model to investigate the influence of the meshing stiffness with tooth profile wear on periodic motion and load sharing behavior of compound planetary gear set. By the analysis of Poincare cross section and phase trajectory, dynamic response in time domain and load sharing behavior for compound planetary gear set is given, it is concluded that: 1) With the accumulation of tooth profile wear, the meshing stiffness of double teeth decreases obviously, and both the stiffness amplitude and fluctuation of single tooth meshing region basically remains unchanged, leading to a decrease in the mean value of time-varying meshing stiffness. The mean value of time-varying meshing stiffness decreases gradually with the accumulation of tooth contact fatigue wear, as well as the fluctuation.
2) The low-frequency motion = 3) The meshing stiffness reduction caused by tooth wear accumulation improves the periodic motion from quasi-period to doubling -period, at the medium frequency section, and the meshing stiffness reduction caused by tooth wear accumulation improves the periodic motion from chaotic motion to quasi-period, at the high frequency section. 4) In order to avoid the system being in Chaos motions, under the condition of large tooth wear, the dimensionless meshing frequency should be avoided [1.3-2], and the input speed of sun gear 1 is converted into the range of 3700 r/min-6725 r/min, so as to keep the system in a stable periodic motion state as far as possible and reduce the noise. 5) Unbalance loading distribution is magnified with the transfer of power from the input to the output. With tooth wear accumulation, the decreasing amplitude and fluctuation of time-varying stiffness enlarge the load sharing coefficient and improve the load sharing behavior for compound planetary gear set.
Purpose of this work is to make a supplement for the research on influence of tooth wear accumulation on dynamic behavior of gear system, based on the work [28][29][30]36]. The calculation model which is used to calculate the decreasing meshing stiffness with tooth profile wear is also appropriate for other gear system, such as parallel-axes gears, simple planetary gear set.
given to this study through the project "Research on Vibration and Noise Reduction Mechanism of High Speed Gear".