Nonlinear dynamics analysis of a gear system considering tooth contact temperature and dynamic wear

Increased temperature and surface wear of high-speed and heavy-load gears are inevitable. Thermal deformation and surface wear modify the position of the action line of the tooth surface and thus influence the dynamic characteristics of the gear mesh. In this study, the elastic modulus and tooth profile thermal deformation were calculated when the tooth contact temperature (TCT) increased. A dynamic wear calculation method was used to combine the dynamic mesh force and dynamic wear coefficient caused by the dynamic mesh force obtained in the nonlinear dynamics model with the quasi-static wear model to obtain the cumulative wear depth. The changed elastic modulus, the tooth profile thermal deformation and the wear depth are considered when calculating the mesh stiffness using the energy method. A nonlinear dynamics model was established by considering the effects of TCT and dynamic wear on the internal dynamic excitation of the gear transmission system. The effects of internal excitation variations such as mesh stiffness, STE and backlash on gear dynamics are analyzed, and the study showed their complex effects on gear dynamics. Comparing the bifurcation diagrams with or without considering TCT and dynamic wear reveals that the system enters chaos earlier after considering TCT and dynamic wear.


Introduction
Gear transmission plays an essential role in aerospace, national defense and military industry, automotive transportation, and other fields. With the high-speed and heavy-load gear system, tooth surface wear often occurs, which will not only alter stress and load distribution, but will also worsen system vibration and noise. 1 Moreover, due to the tooth surface friction, the TCT increases, which induces thermal expansion and modifies the gear meshing characteristics and dynamic characteristics. 2 Consequently, it is essential to research gear dynamics considering TCT and dynamic wear.
Many specialists have studied the influence of tooth surface wear on the dynamics of the gear system. Gao et al. 3 established a dynamics model of spur gear system with 6-degree-of-freedom on a single-stage and proposed a wear fault analysis method on dynamic backlash. Osman and Velex 4 investigated the influence of wide-face solid spur and helical gears on the dynamic performance of gears in quasi-static and dynamic wear. Huangfu et al. 5 suggested a three-dimensional wear prediction model considering the thin-rimmed structure and the modification of the teeth, and analyzed the dynamic performance of the thin-rimmed gear throughout the wear process. Wu et al. 6 combined the static wear calculation model with the dynamic iteration method to study the influence of the clearance from accumulated surface wear on the dynamic characteristics of the composite planetary gear system. Shen et al. 7 considered gear wear as a tooth deflection and investigated the effect of gear wear on dynamic parameters and vibration response of planetary gears. In addition, in the analysis of spur gear system dynamics, the dynamics model was improved by Chen et al. 8 after considering the coupling effect between gear neighboring teeth. The model is able to consider both additional deformation of neighboring teeth caused by loading the gear and the nonlinear teeth contact immediate effects in the dynamics analysis, and the dynamics model is experimentally proven. Wei et al. 9 studied the effect of different uncertain parameters on the dynamic response of a spur gear system. Ambaye and Lemu 10 assumed a uniform distribution of surface wear, modeled the wear as the backlash, and analyzed the gear dynamics using ADAMS. Some scholars have researched the effect of gear dynamics on tooth surface wear. Feng et al. 11 proposed a scheme based on vibration to update wear prediction model. The worn fault tooth profile is regarded as the transmission error and updated into the dynamics model, from which the vibration response and contact force are obtained, and the process is reduplicated to create curves for different wear depths. Ding and Kahraman 12 combined the gear pair dynamics model with the wear model to research the interaction between wear and dynamics under linear and nonlinear conditions. Later, they extended this method to planetary gear system to research the interaction between wear and dynamics of planetary gear system. 13 Liu et al. 1 suggested a dynamic wear calculation method to study the coupling effect between tooth wear and spur gear system dynamics. In this study, the worn surface is represented by TVMS excitation and is applied to the dynamics model to research the influence of tooth wear on the dynamic performance of the system. The Miner rule is also used to uncover the influence of dynamics upon wear behavior. Asymmetric tooth spur gears have the advantages of high bending strength and low vibration. Karpat and Ekwaro-Osire 14 parameterized the tooth tip relief modification and pressure angle and investigated its effect on the wear of asymmetric tooth spur gears. Karpat et al. 15 has developed another new method for calculating the mesh stiffness of asymmetric gears, providing important results for gear researchers to understand involute asymmetric tooth spur gears.
During gear meshing, heat generated by friction causes the TCT to rise and affects the stability of the system. Similarly, many efforts have been made to investigate the dynamic response of the gear system to the TCT. Luo and Li 16 suggested the gear thermal mesh stiffness and brought it into the system dynamics equations in order to research the influence of temperature on the dynamics of gear system. Gou et al. 17 developed a nonlinear dynamics model for a gear-rotor-bearing system to research the effect of tooth flash temperature on the dynamics of a spur gear system. Pan et al. 18 developed a 10 DOF coupled nonlinear dynamics model for a gear-shaft-bearing transmission system and investigated the effects of TCT, and random loading on the dynamic performance. Liu et al. 19 focused on the effect of thermal conditions on the nonlinear vibration of planetary systems under thermal equilibrium. The thermal meshing stiffness is expressed by mathematical expression, the variation of backlash caused by heat is calculated, and the variation trend of nonlinear dynamic characteristics with temperature is obtained. Sun et al. 2 expressed the influence of temperature on the Poisson's ratio and elastic modulus of gear by mathematical expression, and combined it with the improved calculation formula of meshing stiffness to obtain the meshing thermal stiffness. Afterward, the dynamic behaviors of gear system is studied.
At high-speed under heavy-load, rising of tooth surface temperature and increased wear will inevitably occur simultaneously. However, there are few studies on gear dynamics considering both temperature and wear in existing literatures. Zhang et al. 20 used Hertzian theory to obtain the meshing thermal stiffness and calculated the quasi-static wear depth by Archard's equation, and then studied the influence of the changes of meshing stiffness and transmission error on system dynamics. Xu et al. 21 established meshing stiffness model based on flash temperature theory and dynamic wear model by Archard's theory calculation. The effects of initial wear, friction coefficient and damping ratio on system response are also studied. In the above study, gear wear during all operating periods is considered to be the same, without considering the dynamics process of tooth surface wear. Furthermore, the influence of the gear bulk temperature on the modification of the tooth profile thermal deformation is dominant and should be considered. Therefore, in this study, the proposed dynamic wear calculation method is adopted for calculation of the surface wear depth. This method combines the dynamic wear coefficient caused by the dynamic mesh force and the dynamic mesh force obtained from the nonlinear dynamics model with the quasi-static wear model, and finally acquires the accumulated wear depth thanks to a continuous geometry update of the tooth surface. The increase of TCT is accompanied by a modification of the elastic modulus of the tooth and the generation of the tooth profile thermal deformation. This will directly affect internal excitations of the gear system such as TVMS excitation, STE excitation and backlash excitation of gears. The TVMS is calculated using the energy method, and the cumulative wear depth caused by the dynamic surface wear, the tooth profile thermal deformation and the modified elastic modulus caused by the TCT are included. The STE is created by the elastic deformation, wear depth and thermal deformation of gear teeth. Similarly, the wear depth and thermal deformation are taken into account in the backlash. In this study, a nonlinear dynamics model has been established based on the idea that TCT and dynamic wear affect internal dynamic excitation of gear system and then affect nonlinear dynamic characteristics, and the influence of internal excitation changes such as TVMS, STE and backlash on gear dynamics has been studied. Specifically, in Section 2, a method for calculating dynamic wear for a spur gear system is proposed by combining the nonlinear dynamics model with the quasi-static wear model of gear. Section 3 the thermal deformation of the tooth profile along the gear action line in accordance with Blok's theory of flash temperature and the thermal deformation formula was calculated. In Section 4, the research idea is that the TCT and dynamic wear impact the internal excitation of the gear system, and then affect the entire gear system. The impacts of TCT and dynamic wear on TVMS, STE and backlash are analyzed. In the calculation of TVMS, the energy method is used and the impact of temperature on the elastic module of the gear teeth is considered. Then, a nonlinear dynamics model of a spur gear system considering TCT and dynamic wear was established. In Section 5, the influence of changes in internal excitation such as TVMS, STE and backlash on gear dynamics were studied after considering TCT and dynamic wear. Finally, this paper is summarized.

Dynamic wear calculation of a spur gear system
In this section, a method for calculating the dynamic wear of a spur gear system is proposed by combining the nonlinear dynamics model with the quasi-static wear model. On the basis of this calculation method, the depth of surface wear at any point can be obtained, which is more convenient for the study of surface wear of gears.

Nonlinear dynamics model
The dynamics model is established by simplifying the continuous mass-stiffness-damping systems of actual engineering structures into discrete lumped parameter systems.
The nonlinear dynamics model of the spur gear system prepared by Kahraman and Singh 22 was used to develop the motion equation of a spur gear system as follows: where m e , c m , k m , and e(t) denote the equivalent mass, viscous damping, TVMS, and STE. f(x) is a nonlinear description function used to calculate the teeth mesh force with backlash of 2b. F m illustrates the average mesh force, F aT (t) represents the relative fluctuation force of the component of the input external excitation variable. x(t) is defined as: and f(x) can be expressed as: Quasi-static wear model Flodin 23 proposed a wear prediction method for spur gears based on Archard's wear equation. The method can predict the wear depth of any point P on the tooth surface with the following formula: where the wear coefficient k w and the contact pressure p are constants, h P,n denotes the wear depth of the tooth surface at point P with n operation cycles, h P,(n-1) is the wear depth with (n-1) operation cycles, p P,(n-1) indicates the contact pressure, s P illustrates the sliding distance. Based on the above formula, Flodin provided the wear formula for the points on the pinion: and the gear: where a H denotes the half contact width, u p and u g represent peripheral velocities at the two meshing points of the pinion and gear, respectively. According to Hertz contact principle, 24 the contact pressure p P of point P along the action line can be given as follows: where F p is the tangential force of the tooth, B represents the tooth width, y i denotes the distance between the contact point and the Hertz contact center. F p can be calculated as follows: where LSF refers load sharing factor, F denotes the total mesh force tangent to the gear base circle, T illustrates the torque, r b,p represents the base circle radius. a H in equation (7) can be expressed as: where R e and E e indicate the equivalent curvature radius of the contact surface and the equivalent elastic modulus of the material.

Dynamic wear prediction method
In actual gear meshing, k w and p P changed with time due to the internal excitation of the gear. Therefore, in order to calculate more accurately the tooth wear depth of a spur gear system, a dynamic wear calculation method is suggested. The flow chart of dynamic wear calculation method of a spur gear system is demonstrated in Figure 1.
As illustrated in Figure 1, when calculating the dynamic wear depth of a spur gear system, the TVMS model is used to obtain the TVMS k q m(n) (t) after the q-th tooth profile update under different wear depths, firstly. The wear of gear tooth surface will not only affect the TVMS, but also the STE e q (n) (t) and backlash b q (n) (t) will be changed with the aggravation of wear depth. Secondly, the three excitations c hq m(n) (t), e el (n) (t) and b q (n) (t) are brought to the dynamics equation to analyze the dynamic response at this moment. Finally, the dynamics equation is solved to obtain the steadystate response x q (t) after the q-th tooth profile update, and x q (t) refers the relative displacement on the action line during gear meshing, which is used to calculate the dynamic mesh force F q d(n) (t) = k q m(n) (t) f(x q (t)). The load sharing factor is applied to characterize the variation of the load distribution on the tooth surface. LSF represents the ratio of the mesh force of a pair of teeth to the total mesh force, 25 which can be calculated by the following: where F and k demonstrate the mesh force and mesh stiffness of a pair of teeth respectively, superscripts 1 and 2 represent tooth pairs 1 and 2 respectively,Ẽ h refers the error function with the influence of wear, which can be calculated as: where d is the elastic deformation of a gear pair under load. Assumed that the wear depth h w is negative (subscripts p and g represent the pinion and the gear respectively). If two gears contact with an ideal involute profile, then d 1 -d 2 =Ẽ h = 0. Accurate wear calculation depends on reliable load sharing method. Dynamic mesh force should be considered to obtain dynamic load sharing factor LSF d , which can be calculated as: and dynamic contact pressure p Pd can be expressed as: where The wear coefficient 12 is essential for tooth surface wear. The dynamic wear coefficient k wd can be written as: where k wd0 is the wear coefficient under relatively low speed boundary lubrication and it can be given by 26 : where the dimensionless load L wd = FÁLSF d /(bE e R e ), the dimensionless lubricant pressure-viscosity coefficient G w = a pv E e (a pv denote the pressure-viscosity coefficient of the lubrication), the dimensionless composite surface roughness amplitude S w = R rms /E e (R rms represent the composite roughness amplitude), the lambda ratio l = h min /R rms , and the minimum film thickness h min can be given as follows: h min = 3:63U 0:68 G 0:49 w L À0:073 where the dimensionless speed U = h 0 u r /(E e R e ) (h 0 denote the ambient lubricant viscosity and the rolling velocity u r = (u 1 + u 2 )/2), and for spur gears, k is infinity. 12 After LSF d , p Pd and k wd are obtained, the wear depth h Pd of P at any point on the meshing surface can be calculated by the equation for a wear cycle. As illustrated in Figure 1, after n running cycles of the pinion, the wear depth h Pd,(n) at any point on the meshing surface can be calculated. If n is less than the threshold value N 1 , the wear depth is repeated until n reached the set threshold value, under which not only the TVMS will be changed correspondingly, but also the STE and so on will be affected. At this time, it is necessary to reconstruct the gear profile to obtain the dynamic mesh force F q d(n) , dynamic load sharing factor LSF q d(n) , dynamic contact pressure p q Pd(n) and dynamic wear coefficient k q wd(n) after the q-th geometrical update to predict the dynamic wear depth h q Pd , (n) of the tooth surface. Updated the dynamic mesh force and the dynamic load sharing factor can be used to obtained the new dynamic contact pressure and the dynamic wear coefficient to perform new dynamic wear calculations. The process is repeated. After q-th tooth profile update, the pinion had N = qÁN 1 running cycles. Until N reached the threshold value N 2 after Q-th updates. The final accumulated wear depth at all points on the tooth surface can be calculated by summarizing the wear depth at all points updated at different pressures.

Thermal deformation analysis of tooth profile for a spur gear system
In this section, the TCT of the gear under stable running conditions will be calculated and the tooth profile thermal deformation at this temperature will be analyzed. Friction heat occurred during the transmission of the gear system. The friction heat raised the TCT to rise and at the same time accompanies the tooth profile thermal deformation, which will destroy the original involute characteristics of gear and affected the gear transmission.

Calculation of TCT
During the gear meshing process, the TCT T c (t) mainly consists of two parts 27 : the bulk temperature T b and the flash temperature T f (t), as demonstrated below: The bulk temperature of the pinion and the gear is the same and does not change with time under stable working state. However, because of the existence of friction on the surface, and the temperature rises, thereby forming the flash temperature on the surface. The tooth flash temperature may be calculated from the Blok's flash temperature formula 28 : where u = 0.83 indicates the temperature rise coefficient, v i (t)(i = p, g) is the surface tangential velocity, g i shows the heat transfer coefficient, r i means the density, c i illustrates the specific heat, B H (t) demonstrates half of the tooth contact width and the normal load of per tooth width f n =F p /B (F p can be calculated from equation (8) and B refers the tooth width), the friction coefficient f m can be given by 29,30 : where, f mD (t), f mE (t), f mM (t), and f mB (t) demonstrate the friction coefficients under dry friction, elastohydrodynamic lubrication, mixed lubrication and boundary lubrication respectively, and z(t) denotes the oil film thickness.

Calculation of tooth profile thermal deformation
When the gear is in stable operation, the temperature of the gear will reach thermal equilibrium, the temperature field is basically stable, and the bulk temperature will not change any more, but the temperature is non-uniform everywhere. The temperature field inside the base cylinder can be considered as a stable non-uniform steady-state temperature field without internal heat source. Therefore, the temperature spread is constant and only a function of radius, so the temperature changes in a gradient along the radius of the base cylinder. Thus the thermal deformation of the base circle due to the bulk temperature can be expressed as follows 17 : where l denotes the linear expansion coefficient, T(r s,i ) and r s,i (i = p, g) illustrate the temperature of transmission shafts and the shaft radius respectively, T(r b,i ) indicates the temperature of the base circle and takes its value equal to 100°C. 18 As the standard involute spur gear system enters the steady-state running state from the initial state, the base circle radius will change due to the change of TCT. The involute polar coordinate parameter equation of theoretical thermal deformation in the new base circle is as follows: where subscript 1 represents the theoretical value after thermal deformation, r k1,i , u k1,i , and a k1,i are the polar diameter, polar angle and meshing point pressure angle of the theoretical involute after thermal deformation, respectively.
In actual gear transmission, the tooth profile thermal deformation at any meshing point includes two parts: radial direction and tooth thickness direction. Therefore, the polar coordinate parameter equation of actual gear thermal deformation can be calculated as: where subscript 2 represents the actual value after thermal deformation, r k2,i and u k2,i indicate the polar diameter and polar angle after thermal deformation, r k,i and u k,i denote the polar diameter and polar angle before thermal deformation. As demonstrated in Figure 2, Dr k,i and Du k,i can be obtained, which can be calculated as: where Dt = T c (t)2 T 0 (T 0 represents the initial tooth temperature), a k,i indicates the pressure angle before thermal deformation, r i and r k,i denote the pitch circle radius and the radius of any point on the tooth profile before thermal deformation respectively. a and s are pitch circle pressure angle and tooth thickness before thermal deformation respectively.
During gear meshing, the tooth profile thermal deformation due to TCT, resulting in the loss of the involute characteristics of the tooth profile. There is an error between the actual tooth profile and the theoretical tooth profile, which can be expressed by the normal distance d i (t)(i = p, g). As demonstrated in Figure 3 and it can be given by 31 : When d i (t) ł 0, it indicates that the actual tooth profile after thermal deformation is outside the theoretical tooth profile, that is the actual tooth profile envelope theoretical tooth profile, and when d i (t) . 0, the opposite is true.
In the above analysis, d i (t) can also be called noninvolute error after thermal deformation of involute gear. If the standard involute tooth profile is taken as the tooth profile at the TCT at 20°C, the tooth thickness at this time is the tooth thickness s k,i before thermal deformation, while there is a difference Ds kt,i between the tooth thickness s k,i before thermal deformation and the tooth thickness s k1,i after theoretical thermal deformation, and it may be formulated as: where s k1, i = sr k1, i r i À 2r k1, i (inva k, i À inva) Therefore, the tooth profile thermal deformation along the action line d t,i (t) may be calculated as follows: Dynamics modeling of a spur gear system considering TCT and dynamic wear In this section, a dynamics model of spur gear system considering TCT and dynamic wear is proposed, which combines the dynamic wear prediction model of Section 2 with the tooth profile thermal deformation analysis model of Section 3. In the nonlinear dynamics model for a spur gear system, TVMS excitation, STE excitation and backlash excitation are the mainly internal dynamic excitation forms in gear transmission, which have essential influence on the dynamic behavior for spur gear system. In addition, the surface wear occurs only on one side of the teeth, while the tooth profile thermal deformation is bilateral deformation of the teeth.

Calculation of TVMS
In gear transmission system, the TVMS of gears can be calculated according to the energy method. Moreover, the energy saved in meshing gear mainly includes Hertz contact energy, bending potential energy, axial compressive deformation energy, and shear deformation energy. Each energy can calculate its corresponding stiffness separately. The mesh stiffness of the single-tooth-pair may be described as:  and the double-tooth-pair: where k h , k b , k s , and k a represent Hertz contact stiffness, bending stiffness, axial compressive stiffness and shear stiffness respectively. k f is the stiffness of the flexible base, which can be given by 32 : and the parameters of the equation can be found in Chaari et al. 32 Hertz contact stiffness k h is a constant expressed as: where B, E, and n represent tooth width, Elastic modulus and Poisson ratio, respectively. As demonstrated in Figure 4, the energy method for calculating the TVMS of gears is to regard the teeth as cantilever beam on the base circle. Based on energy method, bending stiffness k b , axial compressive stiffness k a and shear stiffness k s are given in the following equation by 33 : where I x and A x denote the moment of inertia and the area of cross-section at a distance of x from the base circle. The geometric parameters h, h x , d, I x , A x are calculated as: Therefore, k b , axial k a and k s are expressed as:

Àa1
(a 2 À a) cos asin 2 a 1 2EB½sin a + (a 2 À a) cos a da 1 k s = ð a2 Àa1 1:2(1 + n)(a 2 À a) cos acos 2 a 1 EB½sin a + (a 2 À a) cos a da where In Section 2, a dynamic wear calculation method of spur gear system was proposed, by which the dynamic wear depth at all points on the tooth profile can be calculated. As illustrated in Figure 5, in actual meshing of spur gear systems, the accumulated wear depth of the tooth surfaces of the pinion and the gear is non-uniformly distributed along the involute tooth profile, the smallest at the pitch point and the largest at addendum. Tooth surface wear changes the geometric parameters h and h x of a cantilever beam, which in turn changes the TVMS.
In Section 3, the tooth profile thermal deformation for a spur gear system was analyzed. Changes in the TCT results in changes in the thickness and height of the teeth, which ultimately results in the profile no longer having involute characteristics after thermal deformation. If the TCT changes Dt, the tooth profile thermal deformation along the action line is d t,i (t). In spur gear system with considering TCT, temperature rise not only causes tooth profile thermal deformation, but also changes elastic modulus of gear material. These two factors together affect TVMS. Li et al. 34 gave a method to predict the elastic modulus at different temperatures: where l(t) denotes the linear expansion coefficient, C v (t) represents the heat capacity at constant volume, T m and T B are the melting point and boiling point of bulk material, respectively. DfusH and DvapH are the molar enthalpy of fusion at the melting point and the molar enthalpy of vaporization at the boiling point, respectively. K is a ratio coefficient. Therefore, the TVMS k 0 considering the TCT and dynamic wear can be calculated as: where Àa 1 3 1 + cos a 1 ½(a 2 À a) sin a À cos a À D f g 2 (a 2 À a) cos a 2E 0 B½sin a + (a 2 À a) cos a + H 3 da Àa 1 (a 2 À a) cos asin 2 a 1 2E 0 B½sin a + (a 2 À a) cos a + H da Àa 1 1:2(1 + n)(a 2 À a) cos acos 2 a 1 E 0 B½sin a + (a 2 À a) cos a + H da where D = (h wd + 2d t cosa 1 )sina 1 /r b , H = (h wd + 2d t cosa)cosa/(2r b ). The wear depth h wd refers the material removed from the ideal involute profile and set to negative, then the tooth profile thermal deformation d t is positive. It is worth noting that d t is the tooth profile thermal deformation along the action line and the wear depth along the action line is h t = h wd /cosa 1 .

Analysis of STE
In dynamics analysis of spur gear system, STE is usually expressed by periodic displacement variation along the action line, which is mainly caused by elastic deformation of loaded teeth and manufacturing error of gears. The manufacturing transmission error of gear pair is composed by the manufacturing transmission error of the pinion and the gear, while the manufacturing transmission error of the pinion and the gear is obtained by the tooth profile manufacturing transmission error of each tooth. Chen and Shao 25 put forward an analytical equation of error along the direction of action line in the mesh of a pair of gear pairs by considering the tooth profile error and spacing error between the teeth. It reads: where d, E, and E s respectively represent the elastic deformation of a gear pair under load, tooth profile error and spacing error between the teeth. Superscripts 1 and 2 represent the first tooth pair and the second tooth pair, and subscripts p and g represent the pinion and the gear respectively. If some material is removed relative to the ideal tooth profile, the tooth profile error E is defined as positive, otherwise it is negative.
In the established dynamics model, the pinion and the gear before meshing are considered as ideal gears without error. In the spur gear system, the cumulative wear depth of surface presents a non-uniform distribution along the involute tooth profile, and the error caused is considered to be the tooth profile error. Although the tooth profile thermal deformation caused by TCT changes the tooth thickness at the gear pitch circle, it will not change the spacing error between the teeth, so the error caused is also considered as tooth profile error.
In this analysis model, the analytical equation of STE could be represented as: where H t = h t, p + h t, g and D t = 2d t, p + 2d t, g.

Calculation of backlash
Backlash in spur gear system can not only prevent thermal expansion caused by rising TCT and stuck, but also store lubricant to reduce surface wear. Therefore, the backlash is very important for the operation of the gear system.
The backlash is caused by the elastic deformation, thermal deformation and wear of the gear teeth. A time-varying backlash model for the gear pair was developed by Shi et al. 35 and the time-varying backlash b 0 (t) can be calculated as: where b 0 denotes the constant part (affected by gear center distance, spacing error between the teeth and installation error etc.), b(t) represents the time-varying part (affected by elastic deformation, thermal deformation and wear of the gear teeth etc.) can be calculated by the following equation: where 2a is the angle between the Back-side action line and the Drive-side action line. 35 In summary, the TVMS k 0 (t), STE e 0 (t) and backlash b 0 (t) considering TCT and dynamic wear are brought into the motion equation (1) to obtain the dynamics model considering TCT and dynamic wear.

Numerical simulation and discussions
In Table 1, the gear pair parameters are summarized, and the prediction of the dynamic wear depth and the calculation of the tooth profile thermal deformation are demonstrated by using the dynamic wear calculation method in Section 2.3 and the tooth profile thermal deformation calculation method in Section 3.2. Then the gear dynamic behaviors considering TCT and dynamic wear is studied., in which TVMS excitation, STE excitation and backlash excitation caused by TCT and dynamic wear are introduced into the gear dynamics model to perform dynamic response analysis.

Calculation of surface wear depth and tooth profile thermal deformation
As described in Section 2, the surface wear depth is primarily influenced from the contact pressure and wear coefficient of the tooth surface. In the dynamic wear model, the tooth contact pressure is related to the dynamic mesh force. The dynamic mesh force  originates from the contribution of the nonlinear dynamic response of the system and varies with the operating speed of the system. Therefore, surface wear is affected by operating speed, and Liu et al. 1 believes that there is not a simple linear relationship between maximum wear depth and operating speed. Surface wear is influenced by the nonlinear dynamic responses of the system, which in turn, as internal excitation, affects the dynamic behaviors of the system. In most literatures, the contact pressure and wear coefficient are considered to be constants. However, in this study, dynamic contact pressure and dynamic wear coefficient are considered. Figure 6 demonstrates the accumulated wear depth of the gear system at all meshing positions on the tooth surfaces of the pinion and the gear with various number of runs N. As demonstrated in Figure 6(a), heavy wear happens in the dedendum area of the pinion due to the high wear coefficient and the high sliding ratio. Conversely, heavy wear happens in the addendum area of the gear (as shown in Figure 6(b)). As there is no relative slip on the pitch of the gear pair, there is no wear. It is worth noting that since the gear pair coincidence is not an integer, there is a sudden change in the tooth surface wear depth during the alternating engagement of single and double teeth. The simulated wear depth demonstrated in Figure 6 is very similar to the experimental results. 36 As described in Section 3, the TCT includes the bulk temperature and the flash temperature. The flash temperature increases monotonically with the increase of friction coefficient, speed and load, which eventually affects the TCT. 17 Tooth profile thermal deformation occurs due to changes in the TCT. When the gear is in stable operation, the tooth profile thermal deformation resulting from the bulk temperature is much larger than that resulting from the flash temperature of the surface, so the tooth profile thermal deformation resulting from the bulk temperature is mainly studied. Figure 7 demonstrates the tooth profile thermal deformation at  all meshing positions on the surfaces of the pinion and the gear with different T b . As demonstrated in Figure  7(a), the tooth profile thermal deformation in the tooth thickness direction decreases as the thickness of the pinion decreases from the dedendum to the addendum of the tooth. Correspondingly, the tooth profile thermal deformation reaches its maximum at the dedendum of the gear (as shown in Figure 7(b)).

Analysis of gear TVMS, STE, and backlash
The influence of TCT and dynamic wear on the TVMS is shown in Figure 8. As the TCT rises, not only the tooth profile thermal deformation but also the elastic modulus changes, which finally results in the monotonous reduction of the TVMS. At the same time, the dynamic wear of tooth surface will also reduce the TVMS of gear. The influence of dynamic wear on stiffness of double-tooth-mesh area is more significant than that of single-tooth-mesh area, which is consistent with other scholars' research results. 37 In fact, the mesh stiffness decreases in the single-tooth-mesh area even if it is not significant. TCT and dynamic wear increase the excitation of TVMS, which in turn leads to the degradation of gear system dynamics.
The influence of TCT and dynamic wear on the STE is shown in Figure 9. During gear meshing cycles, the STE increases monotonously with increasing wear depth and contact temperature. Similar to the change of TVMS, the influence in the double-tooth-mesh area is significantly greater than that in the single-toothmesh area. As one of the main internal excitations, changes in STE due to TCT and dynamic wear also impact the dynamic performance of the gear system.
The variation of the elastic deformation d, wear depth H t and thermal deformation D t of the gear along the action line is illustrated in Figure 10. As observed, the mean value of thermal deformation is relatively large. Therefore, elastic deformation and wear depth have little influence on gear system when wearing slightly, but thermal deformation has great influence on gear system. Consequently, the effect of TCT on spur gear system cannot be neglected. The wear depth jumps and catastrophes in the single-tooth-mesh area and the double-tooth-mesh area, which will lead to backlash jumps. Figure 11 shows the backlash considering the total elastic deformation, total wear depth and total thermal deformation of single and double   tooth meshing. As the TCT and surface wear depth increase, the backlash shows a trend of decreasing first and then increasing along the action line.
TVMS excitation, STE excitation and backlash excitation are the main factors causing vibration and noise for spur gear system. Therefore, in the research of dynamic performance of spur gear systems, the influence of TCT and dynamic wear on TVMS, STE, and backlash must be considered.

Research on nonlinear dynamics of spur gear system
Considering the contributions of TCT and dynamic wear to the internal excitation of the spur gear system, the dynamic behavior of the example system can forecast using a gear nonlinear dynamics model. Based on the established nonlinear dynamics model of spur gear system with TCT and dynamic wear, the impacts of mesh frequency, STE and backlash on system dynamics were analyzed respectively. For convenience of calculation, the periodic TVMS and STE of gears only take the first order component, and then the 4-5 variable-step Runge-Kutta method was applied to solve the nonlinear differential equation of motion of spur gear system.
Effect of mesh frequency on system dynamics. Since the input speed of the pinion directly determines the mesh frequency of the system, it is crucial to research the impact of the mesh frequency on the dynamic behavior of the spur gear system. In dimensionless simulation, the dimensionless meshing frequency v h = v e /v n (v e is the mesh frequency; v n is the natural frequency).
When dimensionless mesh frequency of gear pair is v h = [0.1-1.7], the bifurcation diagrams of mesh frequency without TCT and dynamic wear and with TCT and dynamic wear are demonstrated in Figure 12. In Figure 12  Although it is easy to see that the amplitude of vibration fluctuates with the increase of gear pair mesh frequency in Figure 12, in order to accurately study the system motion mode at a certain mesh frequency, phase portrait is used to find the mesh frequency at which different motion behaviors (single-periodic motion, quasiperiodic motion, period-doubling motion, chaotic motion) are about to occur and cut out the Poincareś ection. Figure 13 demonstrates the phase portraits and Poincare´sections without the TCT and dynamic wear. Figure 14 demonstrates the phase portraits and Poincare´sections with TCT and dynamic wear. By cutting out the Poincare´sections at v h = 0.4 in Figure  12 Based on the above analysis, the bifurcation diagrams (as shown in Figure 12) show the same trend of variation at low and high mesh frequencies in both cases (without TCT and dynamic wear and with TCT and dynamic wear). Compared with the former, the latter system enters into chaos ahead of time. In order to reduce the influence of TCT and dynamic wear on chaotic motion of the system, the mesh frequency should be controlled in the range v h \ 1.23.
Effect of STE on system dynamics. Gear STE is an essential index to measure the dynamic performance of gear pair. The TCT and dynamic wear will cause tooth profile error, so that the STE increases gradually with the transmission of gear pair. Therefore, it is essential to research the effect of STE on the dynamic characteristics of gear system. In dimensionless simulation, dimensionless mesh frequency v h = 1.0 is taken. The mean value of the STE e mm without TCT and dynamic wear is defined as the nominal size, so the mean value of the dimensionless static transfer error e h = e m /e mm (e m is the mean value of STE).
As demonstrated in Figure 15, a bifurcation diagram of dimensionless STE e h = [0.1-3.0] for a spur gear system with TCT and dynamic wear. The system is in single-periodic motion as e h = [0.1-0.5] and quasiperiodic motion as e h = [0.51-1.11]. Subsequently, the system jumps and enters the period-doubling motion. By doubling the period, the system is in period-2 motion as e h = [1.  28-2.29]. The system has experienced alternating transformation areas of complex multi-periodic motion and chaotic motion, one time into period-6 motion, two times into period-2 motion, period-4 motion and three times into chaos  Figure  15, where Figure 16(b) is a point and has a closed ring phase portrait in Figure 16(a), it is in single-periodic motion. With the increase of STE, the phase portrait at e h = 0.68 is a closed ring band and the Poincare´section is a point group, at which time the system performs quasi-periodic motion. In section e h = [1.12-2.29], Poincare´sections are cut at e h = 1.25 and e h = 1.65, as two groups of points, respectively, and the phase portraits shows a closed circle, including many circles of similar shape but different size, which are in period-2 motion. The system begins to double its motion near e h = 1.25 and e h = 1.65, and demonstrates four-point groups on the Poincare´sections at e h = 1.39 and e h = 2.28. The phase portraits are similar to that at e h = 1.25, where the system is in period-4 motion. Then, the Poincare´sections of the system at e h = 1.51, e h = 1.75 and e h = 2.25 is a lot of disordered points, respectively. Phase portraits show many rings which are not compact and the system goes into chaotic motion. The Poincareś ection is doubled to six-point groups at e h = 2.05, which illustrates that the system is moving period-6   In summary, with the increase of STE, the stability of the system gradually decreases and undergoes a complex dynamics process. Single-periodic motion, quasiperiodic motion, period-doubling motion and chaotic motion occur one after another in the system, and complex dynamic behaviors of mutual conversion between period-doubling motion and chaotic motion occurs. Consequently, when considering the TCT and dynamic wear, the spur gear system tends to be unstable due to e h . 1.12 entering the period-doubling motion area when it runs N = 2100 million at stable operating state T b = 100°C.
Effect of backlash on system dynamics. Backlash has strong nonlinear characteristics and its influence on gear dynamics is very complex. Proper backlash is beneficial to the normal working of the gear. Too small backlash will cause the gear to jam due to thermal expansion, while too large backlash will cause collision between high-speed gear teeth and aggravate the vibration and noise of the system. Therefore, it is important to obtain the proper backlash for the gear drive system. In dimensionless simulation, dimensionless mesh frequency v h = 1.0 is taken. The mean value of the initial backlash b mm = 65 mm without TCT and dynamic wear is defined as the nominal size, so the mean value of the dimensionless backlash b h = b m /b mm (b m is the mean value of backlash).
As demonstrated in Figure 17   portrait with a closed ring band in Figure 18(a), it is quasi-periodic motion. As the increase of backlash, it is obvious that the system performs chaotic motion as b h = [0.23-0.46], the phase portrait at b h = 0.35 is not a tight ring and the Poincare´section is a lot of disordered points. In section b h = [0.47-0.91], the Poincare´sections are cut at b h = 0.48, b h = 0.60 and b h = 0.80, respectively, as two groups of points, and the phase portrait illustrates a closed ring band, which includes many rings of similar shape but different size, now in period-2 motion. The Poincare section of the system at b h = 0.50 and b h = 0.69 shows four-point groups and the phase portrait is similar to that at b h = 0.48, where the system is in period-4 motion. Then, the system is transformed into a point group at the Poincare´section of b h = 0.95. The phase portrait shows a closed ring band and the system enters a quasi-periodic motion state. Subsequently, the Poincare´section evolves to a point at b h = 1.25, which indicates a stable singleperiodic motion of the system near it. Finally, the phase portrait at b h = 1.25 and the Poincare´section are similar to those at b h = 0.95, where the system is in quasiperiodic motion.
Thus, the variation of backlash affects the dynamic performance of spur gear systems considering TCT and dynamic wear in a more complex way. With the increase of backlash, there are alternating changes of quasi-periodic motion, chaotic motion, period-doubling motion, quasi-periodic motion, single-periodic motion, quasi-periodic motion. Moreover, the period-doubling motion part changes more frequently. Through perioddoubling and inverse-doubling, the system has the phenomenon of alternating transformation of period-2 motion and period-4 motion. Therefore, when the system has a large backlash b h ø 0.92, it is in a stable quasi-periodic motion or single-periodic motion state.

Conclusions
A nonlinear dynamics model of spur gear system considering TCT and dynamic wear was established, and the effects of mesh frequency, STE and backlash on dynamic behaviors were studied. These conclusions are as follows: (1) By comparing the dynamic characteristics considering TCT and dynamic wear or not, the study shows that with the increase of TCT and dynamic wear of the spur gear system, the TVMS demonstrates a monotonous decrease, STE monotonous increase, and the change trend of backlash first decreases and then increases, which causes the system to enter chaotic motion at the mesh frequency v h = 1.230. And the system enters chaotic motion at the mesh frequency v h = 1.2691 when TCT and dynamic wear are not considered. Therefore, this leads the system to enter chaos in advance. To reduce the influence of TCT and dynamic wear on the chaotic motion of spur gear system, the mesh frequency should be controlled in the range of v h \ 1.230. (2) Errors due to TCT and dynamic wear are considered to be tooth profile errors in STE. With the operation of gear system, the STE increases monotonically, resulting in a complex dynamics process. In the nonlinear dynamic analysis of a spur gear system considering TCT and dynamic wear, the system tends to become unstable due to the STE e h . 1.12 after operating approximately N = 2100 million times at the steady temperature of tooth surface T b = 100°C. (3) In the research of nonlinear dynamic characteristics of a spur gear system, the influence of backlash excitation of the system dynamics is complicated. Due to the change in backlash caused by TCT and dynamic wear, the system will change between stable state and unstable state. Therefore, to reduce the impact of backlash on the system stability, the backlash should be controlled in the range of b h ø 0.92.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.