Generation mechanism and evolution of five-state meshing behavior of a spur gear system considering gear-tooth time-varying contact characteristics

Teeth disengagement or back-side teeth engagement induced by backlash reduces the transmission quality and dynamic performance of gear systems, and the accurate interpretation of multi-state meshing behavior can provide guidance for structural optimization and performance evaluation. Therefore, the multi-state meshing behavior of the gear system is elaborated. A new nonlinear dynamic model of a spur gear system with five-state meshing behavior is established based on the time-varying backlash and contact ratio. The time-varying meshing stiffness and time-varying backlash considering the elastic contact of gear teeth, gear temperature rise and lubrication are included in the model. The five-state meshing behavior is clearly characterized by constructing five Poincaré maps, and its generation mechanism is revealed using dynamic meshing force time history, teeth relative displacement time history and phase portrait. The bifurcation and evolution of five-state meshing behavior are analyzed under the effects of load factor, meshing frequency and error coefficient. The results show that the mutation in the direction of dynamic meshing force leads to teeth disengaging and back-side single- or double-tooth contact, forming multi-state meshing behavior. Bifurcation caused by parameter changes greatly affects the evolution of five-state meshing behavior, particularly grazing bifurcation can decrease the number of teeth disengagement. Chaotic behavior or trajectory expansion inspires multi-state meshing vibration of the system. Previous gear system models could not reveal these phenomena due to ignoring the multi-state meshing behavior.

of the gear system, such as single-tooth drive-side meshing, double-tooth drive-side meshing, singletooth back-side meshing and double-tooth back-side meshing, is possible due to the interaction between single and double teeth alternately meshing as well as backlash. This exacerbates the contact impact vibration, dynamics instability, and gear failure of the gear system [8,9]. Therefore, this paper is dedicated to studying the possible multi-state meshing behavior of gear systems, analyzing its generation mechanism and evolution process, and exploring the inherent law between multi-state meshing behavior and dynamic response.
The study on the characteristics of single-and double-tooth alternating meshing mainly focuses on the time-varying meshing stiffness [10]. The contact areas of single and double teeth along the tooth profile or line of action were accurately identified according to contact ratio [11]. Yang et al. [12] considered the influence of single-and double-tooth meshing characteristics and derived an analytically model of timevarying meshing stiffness of the spur gear system including Hertzian contact stiffness, bending stiffness and axial compressive stiffness based on the energy method. Subsequently, Tian [13] added shearing stiffness on the base of the Yang's model. Wang [14] and Sainsot [15] evaluated the effect of the filletfoundation deflection and added the fillet-foundation stiffness based on Tian's model. The jump of timevarying meshing stiffness is observed in these studies when the single-tooth and double-tooth meshing are switched, which causes periodic meshing impact between the gear teeth. Thus, the studies on the tooth profile modification are carried out by many scholars in order to reduce shock and vibration during the gear mesh [16]. Recently, various models of the timevarying meshing stiffness under faulty gears are developed, such as cracked gears [17][18][19], tooth surface pitting [20], and tooth surface wear [21]. Although diverse time-varying meshing stiffness models have been obtained considering various influencing factors, the effects of tooth surface temperature rise and lubrication are rarely reported in the current literature. Gear temperature rise and lubrication are indispensable in gear transmission. Therefore, it is necessary to establish a comprehensive time-varying meshing stiffness model including elastic deformation, thermal deformation, and oil film thickness deformation in the single-and double-tooth areas for the analyses of time-varying meshing characteristics and dynamics of gear systems. In addition, it is worth noting that single-tooth and double-tooth meshes are rarely considered in the nonlinear dynamic modeling and analysis of gear systems in the previous studies.
At present, the meshing vibration of the gear system is equivalently simplified to the coupling vibration between two rotating bodies connected by a spring and a damper [22][23][24], which is convenient for modeling and analysis. In these studies, the gear contact characteristics are only described by time-varying meshing stiffness with single and double teeth [25][26][27]. However, the real time-varying contact characteristics are ignored, such as the time-varying geometric properties of the meshing point along the tooth profile and the contact behaviors of gear teeth induced by backlash and contact ratio (single-tooth drive-side meshing, double-tooth drive-side meshing, single-tooth back-side meshing and double-tooth back-side meshing).
The constant backlash can clearly describe the drive-side engagement, disengagement and back-side contact of the gear systems [12,28], but it cannot effectively distinguish between single and double teeth during drive-side or back-side tooth meshing. Nevertheless, this can be accurately identified by the time-varying backlash [29]. Due to gear installation and machining errors, tooth surface wear and tooth deformation, the backlash is time-varying along the tooth profile [30][31][32][33][34][35]. Therefore, numerous calculation models of time-varying backlash involving various influencing factors were proposed for different purposes. Li et al. [30] used a first-order sine function to characterize the time-varying backlash caused by the deformation of the support bearing and the transmission shaft. Chen et al. [31] established a random fluctuation model of backlash including the influence of random factors. Chen et al. [32] proposed the dynamic backlash model by deriving the relationship between the gear center distance and backlash. Based on the finite element method, Wang et al. [33] established a time-varying backlash model considering the geometric error of the gear and the deviation of the gear center distance. Wu et al. [34] established a time-varying backlash model that accumulates with the amount of wear according to the static tooth surface wear prediction model. Chen et al. [35] developed a fractal model of the time-varying backlash of a gear system based on the fractal theory.
Although these models describe the time-varying properties of backlash, they cannot distinguish between single and double teeth during drive-side and back-side meshing since the gear-tooth contact characteristics are not considered. As such, the timevarying contact behavior and multi-state meshing behavior of the gear system cannot be precisely elaborated. The authors proposed a new calculation model of the time-varying backlash of the spur gear system considering the influence of the elastic deformation, thermal deformation, oil film thickness deformation at the meshing point [29], which can accurately identify the multi-state meshing behavior.
Based on the time-varying backlash model, this paper proposes a new nonlinear dynamic model of the involute spur gear system including five-state meshing behavior, such as single-tooth drive-side meshing, double-tooth drive-side meshing, teeth disengaging, single-tooth back-side meshing and double-tooth back-side meshing. The time-varying backlash and the time-varying meshing stiffness considering elastic contact, gear temperature rise and lubrication are calculated in detail, and are introduced into the model of the gear system established in this work. Furthermore, the mechanism of multi-state meshing behavior is analyzed, and the bifurcation and evolution of fivestate meshing behavior under the influence of meshing frequency, load factor and transmission error factor are investigated in detail via phase diagram, Poincaré map, bifurcation plot and time history chart of the dynamic meshing force. Some new and interesting results that were impossible to obtain in previous studies were observed. This research can provide a useful reference for dynamic performance optimization, parameter design and further study of gear transmission system, such as gear-tooth meshingimpact dynamics.
The remainder of this paper is structured as follows. In Sect. 2, a nonlinear dynamic modeling of a spur gear pair based on time-varying meshing characteristics is derived strictly. Classification of multi-state meshing behavior is presented in Sect. 2.1, and the    Figure 2 is a schematic diagram of drive-side and back-side tooth meshes for a spur gear pair. R i and R ai are the pitch radius and addendum radius, respectively. e m is the contact ratio of the gears, which is 1:0\e m \2:0. Thus, the alternate meshing between single tooth and double teeth is obtained along Driveside LoA or back-side LoA. Let x ¼ R bp h p À R bg h g À eðsÞ be the relative displacement of meshing gear teeth. eðsÞ is the comprehensive transmission error, which is eðsÞ ¼ e d ðsÞ for drive-side meshing and eðsÞ ¼ e k ðsÞ for back-side meshing. Based on the moving position of the meshing point and the geometric relationship between x and DðsÞ, five meshing states of the spur gear system and their corresponding boundary conditions are summarized as follows.
(i) Double-tooth drive-side meshing state (AB and CD areas in Fig. 2). The boundary condition is x ! DðsÞ and nT m s ðe m À 1Þ nT m ðn ¼ 0; 1; 2; Á Á ÁÞ; (ii) Single-tooth drive-side meshing state (BC area in Fig. 2). The boundary condition is x ! DðsÞ and ðe m À 1ÞnT m \s ðn þ 1ÞT m ; (iii) Double-tooth back-side meshing state (A 0 B 0 and C 0 D 0 areas in Fig. 2). The boundary condition is x À DðsÞ and nT m s ðe m À 1ÞnT m ; (iv) Single-tooth back-side meshing state (B 0 C 0 area in Fig. 2). The boundary condition is x À DðsÞ and ðe m À 1ÞnT m \s ðn þ 1Þ T m ; (v) Teeth disengaged state. The boundary condition is jxj\DðsÞ and nT m s ðn þ 1Þ T m .
Herein, T m ¼2p ðz p x p Þ is a complete meshing cycle including single-tooth and double-tooth meshes (pitch AC or A 0 C 0 ). z p and x p , respectively, represent the tooth number and rotational angular velocity of pinion. According to the above five meshing states, the corresponding dynamic models are, respectively, established below.
where F Np1 , F Np2 , F Ng1 and F Ng2 can be written as: where L d1 ðsÞ and L d2 ðsÞ are the load sharing ratio (see Ref. [36]). F m is the dynamic meshing force, which can be calculated by Eq. (3).
In Eq. (1), F fp1 , F fp2 , F fg1 and F fg2 can be expressed as: where l di ðsÞ is the friction coefficient, and its direction related to the relative slip velocity of two gears [37,38], see Sect. 3.4 for details. S dp1 ðsÞ, S dp2 ðsÞ, S dg1 ðsÞ and S dg2 ðsÞ in Eq. (1) are the friction moments of the tooth pairs, respectively, which can be written as: Therefore, the relative torsion dynamic equation of the spur gear pair under the double-tooth drive-side meshing state can be got by subtracting the two formulas in Eq. (1), as expressed in Eq. (6).
where m e ¼ I p I g . ðR 2 bp I g þ R 2 bg I p Þ is the equivalent mass of the gear pair. F h ðsÞ ¼ Àm e € e d ðsÞ is the error excitation force of the gear pair. g dj ðsÞ ¼ 2Þ is the equivalent friction arm of the tooth pair.
average torque of the gear pair.  and gear under single-tooth drive-side engaging state can be obtained by Eq. (7).

Single-tooth drive-side meshing
Therefore, the relative torsion dynamic equation of the spur gear pair under the single-tooth drive-side meshing state can be written as:

Double-tooth back-side meshing
When the back-side tooth engages, the rotation direction of the gears remains unchanged, and the back-side tooth profile of gear pushes the drive-side tooth profile of pinion to transmit motion and power along back-side LoA. Figure 5 is a schematic diagram of the force analysis of the gear pair at the meshing points under double-tooth back-side meshing state. Based on 2nd Newtonian law, the absolute torsional vibration equation of the pinion and gear under double-tooth back-side meshing state can be obtained by Eq. (10).
where F Np1 , F Np2 , F Ng1 and F Ng2 can be expressed as: where L k1 ðsÞ and L k2 ðsÞ are the load distribution coefficients of the first and second meshing tooth pairs under double-tooth back-side meshing state [36]. F m is the total dynamic meshing force along back-side LoA, which can be calculated by Eq. (12).
Also, F fp1 , F fp2 , F fg1 and F fg2 can be expressed as: where l k1 ðsÞ and l k2 ðsÞ are the friction coefficients along back-side LoA. In Eq. (10), S kp1 ðsÞ, S kp2 ðsÞ, S kg1 ðsÞ and S kg2 ðsÞ are the friction moments of the tooth pairs at the meshing points under back-side tooth meshing, respectively, which can be expressed as: Thus, the relative torsion dynamic equation of the spur gear pair under the double-tooth back-side meshing state can be achieved by Eq. (15).
where g kj ðtÞ ¼ ðR bp I g S kpj ðtÞ þ R bg I p S kgj ðtÞÞ . R 2 bg I p Þ: ðR 2 bp I g þ R 2 bg I p Þ ðj ¼ 1; 2Þ is the equivalent friction arm of the tooth pairs at the meshing point under backside tooth meshing state.

Single-tooth back-side meshing
Therefore, the relative torsion dynamic equation of the spur gear pair under the single-tooth back-side meshing state can be expressed by Eq. (17).

Teeth disengaged state
The pinion and gear are separated from each other under teeth disengaged state, and there is no movement and power transmission between the gear pair. The two meshing gears can be used as two mutually independent rotating bodies. Figure 7 is a schematic diagram of the gear pair under teeth disengaged state. According to 2nd Newtonian law, the absolute torsional vibration equation of the pinion and gear under teeth disengaged state can be written as: Therefore, the relative torsion dynamic equation of the spur gear pair under teeth disengaged state can be obtained by Eq. (19).
Some dimensionless parameters are defined as: is the natural frequency of the gear pair, where k av is the average meshing stiffness.
x ¼ x h =x n is the dimensionless meshing frequency, where x h is the dimensional meshing frequency. t ¼ s=x n is dimensionless time. DðtÞ ¼ DðtÞ D c is the dimensionless time-varying backlash, D c is the characteristic size. The remaindered dimensionless parameters are The dimensionless dynamic meshing force F m is calculated by Eq. (30).
Therefore, the dimensionless normalized expression of nonlinear dynamic model of the spur gear pair including five states (such as double-tooth drive-side meshing, single-tooth drive-side meshing, doubletooth back-side meshing, single-tooth back-side meshing, and teeth disengaging) and time-varying meshing parameters (such as time-varying meshing stiffness, time-varying meshing damping, and time-varying backlash) is obtained, as written in Eq. (25). This model provides a basis for the study of nonlinear dynamics and multi-state meshing behavior of a spur gear system. Due to the significant importance of the time-varying meshing parameters to the nonlinear dynamics and meshing characteristics of the gear system, the time-varying backlash, the time-varying meshing stiffness, and the time-varying meshing damping considering the elastic deformation and thermal deformation at the meshing point for gear teeth as well as the deformation of the oil film thickness are calculated, as seen from Sect. 3'f[.

Time-varying backlash
Time-varying backlash can identify multi-state meshing behavior of spur gear system. The elastic deformation of gear tooth and the deformation of the oil film thickness increase backlash, and the thermal deformation of tooth reduces backlash. It is assumed that the gear pair is installed as standard, thus the timevarying backlash, D(t), is expressed as [29]: where d 0 is the amount of static backlash, which refers to the tooth side clearance reserved in the design and is considered a constant. dðtÞ is a time variable caused by the elastic deformation and thermal deformation of gear teeth as well as the oil film thickness deformation at the meshing point. dðtÞ can be obtained by Eq. (31) [29].
where a is the pressure angle. d e ðtÞ, d t ðtÞ and h oil ðtÞ are the elastic deformation amount, thermal deformation amount and oil film thickness deformation at the meshing point, respectively. The detailed calculation of these deformations can be found in Ref. [29]. The parameter properties of the studied gear pair are shown in Table 1. The static backlash amount is d 0 ¼ 200lm in this work. The distribution characteristics of the time-varying backlash along the action line are shown in Fig. 8. AB and CD represent the double-tooth meshing areas, and BC is the single-tooth meshing area. The time-varying backlash in the singletooth area is greater than the one in the double-tooth areas. The backlash jumps when switching between single and double teeth, such as B and C points. Additionally, the time-varying backlash can be regarded as a constant in the single-tooth and double-tooth areas due to its very small change, as shown in Fig. 8. Therefore, half of the dimensionless time-varying backlash is Ds = 0.89 in the single-tooth area and Dd = 0.79 in the double-tooth areas. This helps to identify and analyze the single-tooth and double-tooth meshing behavior of the gear pair.

Time-varying meshing stiffness
The accuracy of the time-varying meshing stiffness model reflects the precision and validity of the dynamic model of the gear system. The mesh stiffness is mainly characterized by the elastic stiffness in the dynamic modeling and analysis of the gear systems [10][11][12][13][14][15]. Recently, many improved time-varying meshing stiffness models have been proposed to consider different influencing factors. Ma et al. [39] established an improved time-varying mesh stiffness considering the effects of ETC, nonlinear contact stiffness, revised fillet-foundation stiffness, and tooth profile modification. Also, various time-varying stiffness calculation models including the effects of crack propagation path [40], the realistic spalling morphology [41], and roughness [42] are obtained. However, the influence of gear temperature is not considered in these time-varying mesh stiffness models. The thermal deformation caused by the tooth temperature is of great significance for the meshing stiffness of spur gear systems during gear meshing. This thermal deformation is different from the elastic deformation caused by the contact load, such as Hertzian contact deformation, shear deformation, bending deformation and axial compression deformation, etc. In this paper, the stiffness caused by thermal deformation is defined as temperature stiffness. As the reviewer said, the temperature stiffness should be an elastic stiffness considering the effect of temperature. However, there is a certain difference between elastic stiffness and temperature stiffness. That is, the elastic stiffness does not depend on the change of the load, but the temperature stiffness is positively correlated with the load when the thermal deformation is constant. Therefore, in order to highlight this difference, the temperature stiffness is defined additionally. The temperature stiffness and the oil film stiffness caused by the oil film thickness deformation have an important influence on the time-varying meshing stiffness of the gear system. Thus, a comprehensive time-varying meshing stiffness model, k m ðtÞ, including elastic stiffness, temperature stiffness and oil film stiffness is developed in this paper, as expressed as Eq. (33).
where k e ðtÞ, k t ðtÞ and k oil ðtÞ are elastic stiffness, temperature stiffness and oil film stiffness, respectively. The elastic stiffness, k e ðtÞ, includes bending stiffness k eb ðtÞ, axial compression stiffness k ea ðtÞ, shear stiffness k es ðtÞ, fillet-foundation stiffness k ef ðtÞ and Hertzian contact stiffness k eh , which can be written as in Eq. (34) [12][13][14][15].
The temperature stiffness, k t ðtÞ, can be calculated by Eq. (35).
where F m and d t ðtÞ are the dynamic meshing force and thermal deformation amount at the meshing point, respectively. The oil film stiffness, k oil ðtÞ, can be got by Eq. (36) [43].
where DF m is the change in load, and DH oil ðtÞ is the change in oil film thickness caused by DF m . The studied gear parameters are seen from Table 1. The oil film stiffness k oil ðtÞ, temperature stiffness k t ðtÞ, elastic stiffness k e ðtÞ as well as comprehensive timevarying meshing stiffness k m ðtÞ along the action line (LoA) for one tooth pair are shown in Fig. 9a and b. It is thus clear that the introduction of oil film stiffness and temperature stiffness reduces the comprehensive time-varying meshing stiffness. In contrast, the temperature stiffness has a greater effect on the comprehensive time-varying meshing stiffness, while the effect of the oil film stiffness is small due to its larger value.
According to the boundary conditions of singletooth engagement and double-tooth engagement in the time domain: single-tooth engagement ðe m À 1ÞT m t\ðn þ 1ÞT m and double-tooth engagement nT m t\ðe m À 1ÞT m ðn ¼ 0; 1; 2 Á Á ÁÞ, the expression of the comprehensive time-varying meshing stiffness calculation model under single-and double-tooth engaging conditions is expressed as Eq. (37).
A comparison of the time-varying meshing stiffness in this paper with the one in the published Fig. 9 The changes of the time-varying meshing stiffness along the action line (LoA): a The oil film stiffness k oil ðtÞ; b Elastic stiffness k e ðtÞ, temperature stiffness k t ðtÞ and comprehensive timevarying meshing stiffness k m ðtÞ literatures is carried out in order to verify the proposed model including the temperature stiffness and oil film stiffness, as illustrated in Fig. 10. The dark blue curve is the time-varying meshing stiffness calculated by the proposed model including the temperature stiffness and oil film stiffness. The red curve represents the time-varying meshing stiffness calculated from the published literature models [10][11][12][13][14][15], which is only the elastic stiffness. The stiffness proposed in this paper is less than the stiffness in the literature because the temperature stiffness and oil film stiffness are considered. However, the stiffness changes of the two models are consistent and within the same order of magnitude. Notably, the thermal deformation caused by temperature changes and oil film thickness deformation can reduce the meshing stiffness of the gear system. Therefore, the influence of temperature and lubrication cannot be ignored in the analysis of gear meshing vibration.

Time-varying meshing damping
The meshing damping of the spur gear pair is related to the meshing stiffness. Thus, the time-varying meshing damping c m ðtÞ of the spur gear pair can be calculated by Eq. (38).
where 1 is the meshing damping rate and is 1 ¼ 0:07 in this work.

Time-varying friction coefficient under EHL
The oil film stiffness under the elastohydrodynamic lubrication (EHL) condition is calculated in this work, so the time-varying friction coefficient under EHL condition is referred to in the dynamic modeling and analysis of the gear system. Based on the studies of He et al. [44], the time-varying friction formula for the ith meshing tooth pair are expressed as Eq. (39).
where R aavg ¼ ðR a1 þ R a2 Þ=2 is the averaged surface roughness and g M ¼ 0:058 is the dynamic viscosity of the oil entering. The maximum Hertzian pressure P hi ðtÞ of the ith meshing tooth pair is written as Eq. (40).
Herein, f e is the unit normal load, and q ri ðtÞ is the composite relative radius of curvature of the ith meshing tooth pair.
In Eq. (39), SR i ðtÞ ¼ 2v si ðtÞ=v ei ðtÞ is the dimensionless slide-to-roll ratio of the ith meshing tooth pair. v si ðtÞ is the sliding velocity of the ith meshing tooth pair. v ei ðtÞ is the entraining velocity of the ith meshing tooth pair. k i ðtÞ is the direction coefficient of the friction force and can be obtained by the Eq. (41).
Herein, sgn [•] is the sign function, used to determine the direction of friction force.

Characterization of five-state meshing behavior
Based on the analysis in Sect. 2.1, five-state meshing behavior of the spur gear system, such as double-tooth drive-side meshing, single-tooth drive-side meshing, double-tooth back-side meshing, single-tooth backside meshing as well as teeth disengaging, is obtained. Therefore, the existence areas of the five-state meshing behavior in the phase plane can be observed according to the numerical relationship between the teeth relative displacement and time-varying backlash, as shown in Fig. 11. Ds is half of the dimensionless backlash in the single-tooth area and Dd is half of the dimensionless backlash in the double-tooth areas, where Ds [ Dd.
As shown in Fig. 11, single-tooth drive-side meshing only occurs in area I (x 3 C Ds). Single-tooth disengagement occurs in II [ III [ IV area (|x 3 | \ Ds).
Single-tooth back-side meshing only occurs in V area (x 3 B -Ds). Also, double-tooth drive-side meshing occurs in I [ II area (x 3 C Dd). Double-tooth disengagement only occurs in III area (|x 3 | \ Dd). Doubletooth back-side meshing occurs in IV [ V area (x 3 B -Dd). Additionally, area I only appears single-tooth and double-tooth drive-side meshes, where the phase orbit is characterized by the symbol C1. area II appears single-tooth disengagement and double-tooth driveside meshing, where the phase orbit is characterized by C2. area III appears single-tooth disengagement, where the phase orbit is characterized by C3. area IV appears single-tooth disengagement and double-tooth back-side meshing, where the phase orbit is characterized by C4. area V appears single-tooth and doubletooth back-side meshes, where the phase orbit is characterized by C5.
In order to clearly solve and identify the five-state meshing behavior of the gear system, five different Poincaré mapping sections can be constructed based on Fig. 11, such as time-period mapping section, single-tooth disengagement mapping section, singletooth back-side contact mapping section, double-tooth disengagement mapping section, and double-tooth back-side contact mapping section below.
(1) Time-period mapping section: R n ¼ ðx 3 ; _ x 3 ; tÞ f 2 R 2 Â R þ ; modðt; 2p=xÞ ¼ 0g; (2) Single-tooth disengagement mapping section: (3) Single-tooth back-side contact mapping section: Double-tooth back-side contact mapping section: According to the five different Poincaré mapping sections above, the multi-state meshing behavior of the gear system can be characterized by the symbol np-q-r-s. n represents the number of system periodic motion; p represents the number of single-tooth disengagement; q represents the number of singletooth back-side meshing; r represents the number of double-tooth disengagement; s represents the number of double-tooth back-side meshing. If p, q, r or s is zero, it indicates that the there is no single-tooth disengagement, single-tooth back-side engagement, double-tooth disengagement, or double-tooth back- Fig. 11 The existence area of five meshing states in the phase plane side engagement occurring. Thus, the multi-state meshing behavior and nonlinear dynamics of the gear system are obtained by analyzing the characteristics of n-p-q-r-s.

Mechanism of multi-state meshing behavior
The nonlinear dynamic equation, Eq. (25), of the spur gear system including five-state meshing behavior and time-varying meshing parameters is solved numerically by using fourth-order Runge-Kutta method. The mechanism of five-state meshing behavior is analyzed in this section by means of phase diagrams, Poincaré maps and dynamic meshing force time history diagrams.

Drive-side tooth meshing state
The system dimensionless parameters x = 0.3, F = 0.1 and e = 0.15 are fixed. Figure 12a and b shows the phase diagram and the time history diagram of the corresponding dynamic meshing force, respectively. The relative displacement of the system is always greater than half of the single-tooth and double-tooth backlash, x 3 [ Ds = 0.89 and x 3-[ Dd = 0.79. Also, the dynamic meshing force of the system is always greater than zero, F m [ 0.0. It means that the drive-side tooth meshing persists under this parameter condition, and there are no teeth disengagement and back-side tooth meshing. The teeth relative displacement and relative speed of the gear system jump between single tooth and double teeth more obviously because of their small vibration amplitude, as shown in Fig. 12a.

Drive-side tooth meshing and teeth disengaged states
x = 0.3 and F = 0.1 are fixed and other parameters are kept unchanged. Figure 13a is the phase diagram and represents Poincaré mapping point. The system performs a stable period-1 behavior as shown in Fig. 13a, and the phase orbit is distributed in I, II and III regions, as in red marked by C1, blue marked by C2 and purple marked by C3. Figure 13b is the time history diagram of the corresponding dynamic meshing force F m . F m changes periodically between greater than zero (F m-[ 0) and equal to zero (F m = 0). This means that the drive-side tooth mesh and teeth disengagement occur periodically. The disengagement state occurs in the single-tooth area and the double-tooth area. As shown in Fig. 13b, when the drive-side tooth mesh (F m [ 0) transitions to the teeth disengagement (F m = 0), the dynamic meshing force has two sudden changes, namely: In order to reveal the mechanism of teeth disengagement, it is necessary to analyze the mutation process of dynamic meshing force in detail. Figure 14 shows the changes of the dynamic meshing force and teeth relative displacement in the time domain when the drive-side tooth mesh is transferred to the teeth disengagement. The dynamic meshing force of the system jumps periodically, and the force of the double-tooth area (double-tooth F m ) is greater than that of the single-tooth area (single-tooth F m ).
The dynamic meshing forces are greater than zero before point F1, and the system performs drive-side tooth meshing state. At point F1, the single-tooth F m decreases to zero, and then increases in the opposite direction (see area F1F2). In the F1F2 area, the system is still in contact with the single-tooth drive-side because the teeth relative displacement is greater than half of the single-tooth backlash (x 3 [ Ds). However, the meshing gear teeth tend to separate from each other due to the change in the direction of the meshing force. At point F2, the relative displacement is equal to half of the single-tooth backlash (x 3 = Ds) and the single-tooth F m suddenly becomes zero. It indicates that the meshing gear teeth are about to separate. Subsequently, teeth relative displacement is smaller than half of the single-tooth backlash, and the singletooth meshing force is equal to zero, and thus the single-tooth disengagement occurs in the area F2F3. Double-tooth drive-side meshing persists in F2F3 since double-tooth F m is greater than zero and the teeth relative displacement is higher than half of the doubletooth backlash (x 3 [ Dd). The double-tooth meshing force decreases to zero at point F3 and then increases in the opposite direction (see area F3F4). In the F3F4 area, the system is in double-tooth drive-side contact state because the relative displacement is still greater than half of the double-tooth backlash (x 3 [ Dd). But the meshing teeth tend to separate from each other due to the change of the direction of the meshing force. At F4, x 3 = Dd and double-tooth F m = 0 mean that the meshing gear teeth are about to separate. And then, the relative displacement is less than half of the doubletooth backlash (x 3 \ Dd) and the dynamic meshing force is equal to zero (F m = 0). The system is completely disengaged after F4.
Therefore, the transition process from drive-side tooth mesh to teeth disengagement occurs in the F1F4 area. The mechanism of gear teeth disengagement is described as follows. The single-tooth dynamic Fig. 14 Mechanism of the transition from drive-side tooth meshing to teeth disengagement meshing force first decreases to zero and then increases in the opposite direction until the relative teeth displacement is equal to the half value of the single-tooth backlash (point F2). Then the teeth disengagement occurs in the single-tooth area (such as F2F3 area). Subsequently, the double-tooth dynamic meshing force decreases to zero, and then increases in the opposite direction. Until the relative teeth displacement is equal to the half value of the double-tooth backlash (point F4), the double-tooth dynamic meshing force suddenly decreases to zero, and then the system is completely in a state of teeth disengagement. Thus, the dynamic meshing force experienced two sudden changes in the direction, respectively, in the single-tooth and double-tooth areas, resulting in the occurrence of gear teeth disengagement.

Back-side tooth meshing and teeth disengaged states
Take x = 1.4, other parameters are consistent with those in Fig. 13. Figure 15a is the phase diagram and Poincaré map of the gear system where double-tooth drive-side mesh (C1 and C2), single-tooth drive-side mesh (C1), double-tooth back-side mesh (C4 and C5), single-tooth back-side mesh (C5) and teeth disengagement (C3) are observed. The system performs a stable period-3 behavior and its phase orbit are distributed in I, II, III, IV and V regions. The meshing state of the system becomes extremely complicated. Figure 15b shows the time history of the corresponding dynamic meshing force F m , which changes periodically between greater than zero (F m [ 0), equal to zero (F m = 0) and less than zero (F m \ 0).
It indicates that the single-tooth drive-side meshing, double-tooth drive-side meshing, single-tooth backside meshing, double-tooth back-side meshing and teeth disengagement appear periodically as shown in Fig. 15a. Additionally, when the drive-side tooth meshing (F m [ 0) transitions to teeth disengaging (F m = 0) or the back-side tooth meshing (F m \ 0) transitions to teeth disengaging (F m = 0), the dynamic meshing force has undergone two sudden changes, such as local magnifications in the picture. The mutation process from drive-side tooth engagement to teeth disengagement is the same as that in Fig. 14. However, the mutation process from back-side tooth meshing to teeth disengagement is analyzed in detail below. Figure 16 shows the changes of dynamic meshing force F m and teeth relative displacement x 3 in the time domain when the back-side tooth meshing transitions to teeth disengagement. The dynamic meshing forces are less than zero, and the system performs the backside tooth meshing state before point F1. At F1, the absolute value of single-tooth F m is reduced to zero and then increases in the opposite direction (see area F1F2). Between F1 and F2, the single-tooth drive-side contact persists because the teeth relative displacement is less than half of the negative value of singletooth backlash (x 3 \ -Ds). However, the contact teeth tend to separate from each other in the singletooth area because of the change in the direction of the dynamic meshing force (single-tooth F m ). At F2, the relative displacement is equal to half of the negative value of single-tooth backlash (x 3 = -Ds) and the single-tooth F m suddenly becomes zero. It means that the back-side meshing teeth are about to separate. Subsequently, the relative displacement is greater than , and the single-tooth meshing force is equal to zero, so the single-tooth disengagement appears between F2 and F3. However, double-tooth back-side meshing persists in the F2F3 area because double-tooth F m is less than zero and the teeth relative displacement is lower than half of the negative value of double-tooth backlash (x 3 \ -Dd). The doubletooth meshing force decreases to zero at point F3 and then increases in the opposite direction (see area F3F4). The double-tooth drive-side contact state persists because the relative displacement is less than half of the negative value of double-tooth backlash (x 3 \ -Dd) between F3 and F4. But the double-tooth meshing teeth tend to separate from each other due to the change of the meshing force in direction. At F4, x 3 = -Dd and single-tooth F m = 0 indicate that double-tooth meshing teeth are about to separate. Subsequently, the relative displacement is larger than half of the negative value of double-tooth backlash (x 3 [ -Dd) and the dynamic meshing force is equal to zero (F m = 0). The completely teeth disengagement is obtained after F4. Therefore, the transition process from back-side tooth meshing to teeth disengagement appears in the F1F4 area. The single-tooth dynamic meshing force first changes in direction between F1 and F2. And then single-tooth back-side teeth disengagement occurs between F2 and F3. Subsequently, the double-tooth dynamic meshing force changes in direction between F3 and F4. Finally, the system enters a completely teeth disengaged state after point F4. Similarly, the mutation in the direction of the dynamic meshing force leads to the occurrence of teeth disengagement.

Bifurcation and evolution of five-state meshing behavior
Due to time-varying meshing stiffness, time-varying backlash, dynamic transmission error, and switching between different meshing states, the gear system becomes a non-smooth nonlinear strong time-varying system. The five-state meshing behavior and nonlinear dynamics caused by the system parameters have a significant impact on the transmission efficiency and motion quality of the system. The effects of meshing frequency, load factor, and transmission error coefficient on the nonlinear dynamics of five-state meshing behavior of the gear system are discussed in detail as follows.

Effect of meshing frequency
Take F = 0.13 and e = 0.25 and the meshing frequency x is used as the control parameter. Based on five different Poincaré mapping sections established in Sect. 3.1, the bifurcation diagram of the five-state meshing behavior of the system and the corresponding top Lyapunov exponent (TLE) are shown in Fig. 17ad increasing in x. Figure 17a is the bifurcation diagram in time-period mapping section R n in blue. In Fig. 17b, the bifurcation diagrams in single-tooth disengagement mapping section R p and single-tooth back-side contact mapping section R q are in magenta and green, respectively. In Fig. 17c, the bifurcation diagrams in double-tooth disengagement mapping section R r and double-tooth back-side contact mapping section R s are in red and grass green, respectively. Figure 17d shows the corresponding top Lyapunov Fig. 16 Mechanism of the transition from back-side tooth meshing to teeth disengagement exponent (TLE), which is used to analyze the dynamic stability of the system. TLE value less than zero indicates a stable periodic behavior, and TLE value greater than zero means an unstable chaotic behavior. TLE value is approximately equal to zero or jumps indicate the bifurcation or sudden change of the dynamic behavior of the system, such as the black dot ''•'' (A, B, C, D, E and F). When the meshing frequency x is small (left side of point A), the system shows a stable 1-0-0-0-0 behavior, where no teeth disengagement and backside tooth meshing occur. The corresponding TLE values are lower than zero. Figure 18a shows the phase portrait and Poincaré map calculated for x = 0.2 where the symbol ''•'' is Poincaré mapping point. The vibration amplitude between single-tooth and double-tooth is obviously jumping due to the slight vibration. The time history of the corresponding dynamic meshing force is shown in Fig. 18b, which is always greater than zero. With the increase in x, single-tooth disengagement appears at point A. The system performs a stable 1-1-0-0-0 behavior between A and B points. The phase portrait and Poincaré map are displayed in Fig. 18c as x = 0.5. The phase portrait crosses the half of the value of single-tooth backlash (Ds = 0.89), where C1 is the portrait of single-tooth and double-tooth drive-side meshing portrait and C2 is the portrait of single-tooth disengagement and double-tooth drive-side mesh. Figure 18d shows the corresponding meshing force in time domain. It is observed that the dynamic meshing force is equal to zero in the single-tooth area (see partial enlargement), and the dynamic meshing force in the double-tooth area is always greater than zero. This means that the single-tooth disengagement occurs periodically.
At B, 1-1-0-0-0 behavior jumps and transitions to a 1-1-1-1-1 behavior and then the system performs a stable 1-1-1-1-1 behavior between B and C points. The multi-state dynamics of the system become more complicated. The phase trajectory of the system passes through half of the value and half of the negative value of the single-tooth and double-tooth backlash, as marked by C1, C2, C3, C4 and C5, as shown in Fig. 18e. The phase trajectory of the system becomes smoother than that of Fig. 18a and c due to the expansion of its topology, and thus the jump of vibration amplitude between single and double teeth is not obvious. The corresponding dynamic meshing force changes periodically between greater than zero (drive-side meshing), equal to zero (disengagement) and less than zero (back-side meshing), as illustrated in Fig. 18f. Also, the mutation of the meshing force is observed when the drive-side meshing transitions to disengagement or back-side meshing transitions to disengagement.
At C, 1-1-1-1-1 behavior shifts to a stable 1-1-0-1-0 behavior. Single-tooth and double-tooth back-side meshes disappear for 1-1-0-1-0, but single-tooth and double-tooth disengagements persist. Phase trajectories of different meshing states are marked by the symbols C1, C2 and C3 as shown in Fig. 18g. The dynamic meshing force of the system appears periodically between greater than zero and equal to zero, as seen from Fig. 18h. 1-1-0-1-0 behavior enters a short chaotic behavior via period-doubling bifurcation sequence and then degenerates to a stable 2-2-0-2-0 behavior by inverse period-doubling bifurcation sequence. The phase orbit and dynamic meshing force in time domain of 2-2-0-2-0 behavior are shown in Fig. 18i and j. 2-2-0-2-0 behavior degenerates to 2-1-  Fig. 18k and the corresponding dynamic meshing force is shown in Fig. 18l. The meshing force changes periodically between F m [ 0, F m = 0 and F m \ 0. Thus, the single-tooth and double-tooth disengagements as well as single-tooth and double-tooth back-side meshes are observed in chaos. As x increases further, the chaos degenerates into a 2-1-1-1-1 behavior, and then 2-1-1-1-1 behavior enters chaos. Subsequently, the chaos degenerates into a 3-1-1-1-1 behavior. The phase trajectory of the system crosses Ds, -Ds, Dd and -Dd, as shown in Fig. 18m. The five-state meshing behaviors are observed as marked by C1, C2, C3, C4 and C5. The time history of the corresponding dynamic meshing force is shown in Fig. 18n, which appears periodically between F m [ 0, F m = 0 and F m \ 0. 3-1-1-1-1 behavior transfers to a 1-1-0-1-0 behavior at point F. The system performs drive-side tooth mesh, singletooth and double-tooth disengagements after F. Figure 18o is the phase portrait and Poincaré map of 1-1-0-1-0. Therefore, the system performs a stable drive-side meshing state when the meshing frequency x is small. With the increase in x, the system gradually appeared multi-state meshing behavior such as single-tooth and double-tooth disengagement, double-tooth back-side meshing and single-tooth back-side meshing. The multi-state meshing behavior of the system become complicated decreasing x. In addition, the bifurcation, periodic jump, and chaos are found between B and F points (0.5 \ x \ 2.0). Teeth disengagements and back-side meshes caused by bifurcation or chaos are also obtained in this area. Bifurcations such as B, C, D, E and F points not only modify the type of motion, but also change the multi-state meshing behavior of the system. Also, chaotic behavior stimulates tooth disengagement and back-side tooth contact, forming a multi-state meshing behavior. When x is large (to the right of point F), although the system exhibits a stable period-1 behavior, the single-tooth and doubletooth disengagement still occurs, as seen from 1-1-0-1-0 behavior.

Effect of load factor
x = 1.3 and e = 0.12 are fixed, and other parameters remain unchanged. Figure 19 shows the bifurcation diagrams of the gear system in different Poincaré sections and TLE diagram with the decrease in load factor F. When F is large (on the right of point A), the system exhibits a 1-1-0-1-0 behavior, which is a period-1 motion and the corresponding TLE values are less than zero. Only the single-and double-tooth disengagements are observed and no back-side teeth b Fig. 18 Fig. 20e where the curve C2 is tangent to the straight line Dd, and then doubletooth disengaging is reduced once. The phase trajectory crosses x 3 = Ds twice and crosses x 3 = Dd once, as indicated in Fig. 20f. The corresponding dynamic meshing force is illustrated in Fig. 20g.
At C point, 2-2-0-1-0 jumps and transitions to a 4-2-0-2-0 behavior. Figure 20h and i shows the phase trajectory and dynamic meshing force in time domain for 4-2-0-2-0. Subsequently, 4-2-0-2-0 enters into a chaotic behavior by short period-double bifurcations and TLE values of the chaotic behavior are higher than zero. Figure 20j is Poincaré mapping diagram of chaos calculating for F = 0.09, and Fig. 20k is the corresponding dynamic meshing force. The dynamic meshing force changes periodically between F m [ 0 and F m = 0, and no F m \ 0 appearing. Thus, there is not any back-side teeth mesh occurring for the chaos. However, the chaotic attractor enlarges at E point and then single-tooth and double-tooth back-side meshes appear. The enlarged chaotic attractor is depicted in Fig. 20l and the corresponding dynamic meshing force is shown in Fig. 20m where F m \ 0 is observed.
As F is further decreased to point F, the chaotic behavior degenerates to a 3-1-1-1-1 behavior. The system performs a stable period-3 behavior, but the teeth disengagements and back-side tooth mesh persist as shown in Fig. 20n. The corresponding dynamic meshing force changes periodically between F m [ 0, F m = 0 and F m \ 0 as shown in Fig. 20o.
Therefore, when the load factor is large, the system behaves as a stable period-1 behavior. The multi-state meshing behavior is simple, and only single-tooth and double-tooth disengagements occur. As decreasing in load factor F, the system undergoes period-doubling and grazing bifurcations, and five-state meshing behavior is gradually observed. In addition, grazing bifurcation (B point) reduces the number of doubletooth disengagement. Bifurcation changes the multistate meshing behavior of the system. Additionally, the single-tooth and double-tooth back-side meshes are gradually observed. Both the multi-state meshing behavior and nonlinear dynamics become complicated decreasing in F.

Effect of transmission error coefficient
Take F = 0.05 and x = 1.3, Fig. 21 shows the bifurcation diagrams of the five-state meshing behavior of the system with the increase in transmission error coefficient e. When e is small (left of point A), the system behaves as a stable 1-0-0-0-0 behavior, and the corresponding TLE value is less than zero. There is a complete drive-side meshing state is observed, and no disengagement and back-side meshing appearing. The system exhibits a flutter motion in the range of extremely small vibration amplitude, as shown in Fig. 22a calculating for e = 0.012.
With the decreasing in e, single-tooth disengagement is observed at A point and then 1-1-0-0-0 behavior is obtained between A and B points. The phase portrait and Poincaré map are depicted in Fig. 22b, and the corresponding dynamic meshing force is shown in Fig. 22c. The phase portrait crosses x 3 = Ds, as marked by C1 and C2. At B, 1-1-0-0-0 behavior jumps to 2-1-0-1-0 behavior. The system performs a period-2 motion and appears double-tooth disengagement between B and C points. Figure 22d is the phase portrait and Poincaré map of 2-1-0-1-0 b Fig. 20  behavior. The phase portrait crosses x 3 = Ds and x 3 = Dd, as marked with C1, C2 and C3. The corresponding dynamic meshing force changes periodically between F m [ 0 and F m = 0. Drive-side mesh, single-tooth disengagement and double-tooth disengagement occur periodically during gears operating. 2-1-0-1-0 behavior jumps at C and D points, and then enters chaos. The corresponding TLE value is higher than zero. Poincaré map and dynamic meshing force of the chaos are shown in Fig. 22f and g as e = 0.07. The meshing force changes between F m [ 0 and F m = 0. Single-tooth and double-tooth disengagements persist, but no back-side meshing occurs for the chaotic behavior. However, the chaotic behavior degenerates into 6-2-2-2-2 behavior at point E. Single-tooth and double-tooth back-side meshes are observed. The phase portrait passes through x 3 = Ds, -Ds, Dd and -Dd, as shown in Fig. 22h. The phase trajectories under different meshing states are marked with C1, C2, C3, C4 and C5. The dynamic meshing force changes periodically between F m [ 0, F m = 0 and F m \ 0, as shown in Fig. 22i.
Subsequently, 6-2-2-2-2 behavior enters to chaos via a short period-doubling bifurcation sequence. Teeth disengagements and back-side tooth meshes persist for the chaotic behavior. Furthermore, the chaotic behavior degenerates into a 6-1-2-1-2 behavior at point F. 6-1-2-1-2 behavior exists in a larger range of e and then enters chaos after a series of bifurcation and evolution. Figure 22j and k shows the Poincaré map and dynamic meshing force of the chaotic behavior. The dynamic meshing force alternates between F m [ 0, F m = 0 and F m \ 0. Teeth disengagements and back-side meshes always occur.
Therefore, when e is small, the system exhibits a complete drive-side meshing state without disengagement and back-side meshing. With the increase in e, however, single-tooth disengagement occurs at point A, and then double-tooth disengagement appears at point B. Single-tooth and double-tooth back-side meshes are observed at point E. Subsequently, teeth disengagements and back-side tooth meshing persist after point F. In addition, the system behaves as a stable periodic behavior when the transmission error coefficient e is small. Periodic jumps and bifurcation points are obtained increasing in e, such as B, C, D, E and F points. The system performs chaotic and periodic behaviors where teeth disengagements and back-side tooth engagements are observed as e is large. Bifurcation leads to gear tooth disengagement and back-side tooth contact, resulting in multi-state meshing behavior. The type of multi-state meshing behavior is also changed through bifurcation or periodic jump, such as B and E point. The transmission error coefficient e greatly affects the nonlinear dynamics as well as multi-state meshing behavior of the gear system. When the motion trajectory is enlarged or chaotic, the system behaves as a multi-state meshing behavior.
In order to demonstrate the advantages of the proposed model, the nonlinear vibration obtained from the traditional method is compared with the proposed method. A bifurcation diagram with the increase in the transmission error coefficient e is calculated by the traditional model, as shown in Fig. 23, which corresponds to Fig. 21a.
Comparing Fig. 21a and Fig. 23, when e is small (at the left of F), the two bifurcation graphs have a certain similarity, and the evolutionary trend of dynamic  behavior is basically the same. However, the difference is that the period-2 behavior near point B cannot be revealed in Fig. 23, and the position of the bifurcation point E (chaotic behavior degenerates to period-6 behavior) is advanced, and the parameter range of the period-6 behavior is shortened in Fig. 21a compared to Fig. 23. When e is large (at the right of point F), the bifurcation diagrams obtained by the two methods are quite different. The dynamic behavior obtained by the traditional method mainly shows periodic behavior when e is large, while the dynamic behavior obtained by the proposed model mainly shows period-6 behavior and chaotic behavior, and the nonlinear vibration of the gear system is unstable and complex. Therefore, compared with the traditional model, the proposed model in this work can better reflect the dynamical characteristics of the influence of the transmission error coefficient e on the gear system. The proposed model in this work has certain advantages.
Notably, based on the actual engaging process of the gear pair, a dimensionless normalized model of the gear pair system including five-state mesh is established. Then five different Poincaré mapping sections are obtained according to the switching conditions between different meshing states. Based on this, the multi-state meshing behavior is clearly characterized, and the evolutionary mechanisms of single-tooth or double-tooth disengagement and back-side tooth contact are well revealed. However, these behaviors and phenomena cannot be obtained in the traditional model because the multi-state meshing characteristics are not considered.

Conclusions
The time-varying contact characteristics of an involute spur gear system considering backlash and contact ratio are elaborated. A new nonlinear dynamic model of the gear system with five-state meshing behavior, such as single-tooth drive-side mesh, double-tooth drive-side mesh, teeth disengagement, single-tooth back-side mesh, and double-tooth back-side mesh, is established in this work. The time-varying meshing stiffness and time-varying backlash considering the elastic contact of gear teeth, gear temperature rise and lubrication are calculated and introduced into this model. The generation mechanism of multi-state meshing behavior is revealed by using dynamic meshing force time history and phase portrait. The bifurcation and evolution of five-state meshing behavior are analyzed separately under the effects of meshing frequency, load factor and transmission error coefficient. The nonlinear vibration obtained from the traditional method has been compared with the proposed method to demonstrate the advantages of the model. Some beneficial conclusions are obtained as follows: (1) Five different meshing states of a spur gear system and their boundary conditions are clearly described. A dimensionless normalized nonlinear dynamic model of the gear system including five-state meshing behavior together with the effects of time-varying backlash and timevarying meshing stiffness are elaborated. The time-varying backlash and time-varying meshing stiffness are calculated in detail with the considerations of the elastic contact of the gear teeth, and the gear temperature rise and lubrication. (2) The mutation in the direction of the dynamic meshing force results in the disengagement of gear teeth, and then five-state meshing behavior is observed. This mutation occurs successively in the single-and double-tooth meshing regions. The five-state meshing behavior is well characterized by constructing five different Poincaré mapping sections, such as time-period mapping section R n , single-tooth disengagement mapping section R p , single-tooth back-side contact mapping section R q , double-tooth disengagement mapping section R r , and double-tooth back-side contact mapping section R s . (3) The meshing frequency, load factor, and transmission error coefficient greatly affect the nonlinear dynamics and multi-state meshing behavior of the gear system. When meshing frequency or load factor is small, the system behaves as a complete drive-side tooth meshing state and a stable period-1 behavior. The multistate meshing behavior, bifurcation, and chaotic behavior gradually occur as meshing frequency or load factor increase. The multi-state meshing behavior stimulates the complexity of the nonlinear dynamics of the gear system. In addition, grazing bifurcation reduces the number of double-tooth disengagement. These phenomena, however, cannot be obtained in the traditional method due to the neglect of the multistate meshing behavior.
This paper provides a scientific basis for the parameter optimization design, dynamic performance improvement and reliable application of the gear transmission system, and also lays a foundation for further research of the gear system in meshing-contact impact characteristics.