Research on the Influence of Backlash on Mesh Stiffness and the Nonlinear Dynamics of Spur Gears

In light of ignoring the effect of backlash on mesh stiffness in existing gear dynamic theory, a precise profile equation was established based on the generating processing principle. An improved potential energy method was proposed to calculate the mesh stiffness. The calculation result showed that when compared with the case of ignoring backlash, the mesh stiffness with backlash had an obvious decrease in a mesh cycle and the rate of decline had a trend of decreasing first and then increasing, so a stiffness coefficient was introduced to observe the effect of backlash. The Fourier series expansion was employed to fit the mesh stiffness rather than time-varying mesh stiffness, and the stiffness coefficient was fitted with the same method. The time-varying mesh stiffness was presented in terms of the piecewise function. The single degree of freedom model was employed, and the fourth order Runge–Kutta method was utilized to investigate the effect of backlash on the nonlinear dynamic characteristics with reference to the time history chart, phase diagram, Poincare map, and Fast Fourier Transformation (FFT) spectrogram. The numerical results revealed that the gear system primarily performs a non-harmonic-single-periodic motion. The partially enlarged views indicate that the system also exhibits small-amplitude and low-frequency motion. For different cases of backlash, the low-frequency motion sometimes shows excellent periodicity and stability and sometimes shows chaos. It is of practical guiding significance to know the mechanisms of some unusual noises as well as the design and manufacture of gear backlash.


Introduction
Gear dynamics is an important method to predict dynamic performance as well as to monitor the status of a gear system.The time-varying mesh stiffness caused by alternately changing the gear pairs in mesh plays a crucial role in dynamic responses.An effective time-varying mesh stiffness model is a basic condition to conduct a dynamic analysis.Thus, it is significant to recognize the mechanism of a gear system's shock and vibration to investigate the influence of time-varying mesh stiffness on vibration responses, especially with regard to the mechanism of the noise, which is hard to identify.
One consistently popular topic in gear dynamics is how to describe the change rule of time-varying mesh stiffness correctly and effectively.The ISO Standard 6336-1 [1] recommends a stiffness formula to calculate the maximum single tooth mesh stiffness for spur and helical gears according to experimental tests.In References [2][3][4], time-varying mesh stiffness was simplified as a rectangular equation which had significant differences from the actual value.The authors in Reference [5] used a polynomial fitting method to smooth the result of the time-varying contact stiffness that was obtained by the Hertz contact algorithm.The Fourier series was used to fit time-varying mesh stiffness function approximatively in Reference [6].Kuang and Yang proposed an equation to calculate the single tooth bending stiffness for the addendum-modified gear and the curve-fitted coefficients took the number of teeth as a variable [7].In Reference [8], the mesh stiffness between an engaged gear pair was regarded as consisting of local Hertzian stiffness and tooth bending stiffness.Cai and Tang revised Kuang's calculation formula and proposed an analytical method to calculate the mesh stiffness by considering the load condition and the tooth profile modification in References [9,10].Chen investigated the differences between the rectangular stiffness and its approximate form on the gear nonlinear dynamic characteristics in Reference [11].
According to the available literature, the main methods to calculate mesh stiffness are Finite Element Method (FEM) and the potential energy method.These two methods have advantages and disadvantages.The potential energy method was proposed by Yang [12] and has been improved several times [13][14][15][16].It is widely used in the failure analysis of a gear system [17][18][19][20][21].The FEM, which can reflect the real contact status, is often applied to figure out engineering problems [22][23][24].Moreover, the FEM is a common method to verify the validity of the potential energy method while both the FEM and potential energy method always ignore the influence of backlash on mesh stiffness.
The backlash is very useful to avoid the jam phenomenon in conditions of elastic deformation and thermal expansion, i.e., adaptive backlash is one of the necessary conditions to ensure the normal operation of gear.One way to achieve backlash is to reduce tooth thickness.The reduction in tooth thickness will cause a decrease in the stiffness, further influencing the dynamic performance.In current gear dynamics, the effect of backlash is limited to the definition of the nonlinear backlash function [25][26][27][28][29][30].
Thus, this work will focus on the influence of backlash on time-varying mesh stiffness, carry out nonlinear dynamic analysis, and investigate the dynamic responses for various values of backlash.This is expected to establish the indirect relationship between backlash and dynamic performance as well as direct the design and control of the gear's backlash.

Mesh Stiffness Model with Backlash
Available mesh stiffness models do not consider the influence of backlash.Precision tooth profile equations play an important role in the analysis method to calculate the mesh stiffness.According to the generating method of gear cutting, the tooth profile equations consist of two parts: the involute and transition curve, as shown in Figure 1.
experimental tests.In References [2][3][4], time-varying mesh stiffness was simplified as a rectangular equation which had significant differences from the actual value.The authors in Reference [5] used a polynomial fitting method to smooth the result of the time-varying contact stiffness that was obtained by the Hertz contact algorithm.The Fourier series was used to fit time-varying mesh stiffness function approximatively in Reference [6].Kuang and Yang proposed an equation to calculate the single tooth bending stiffness for the addendum-modified gear and the curve-fitted coefficients took the number of teeth as a variable [7].In Reference [8], the mesh stiffness between an engaged gear pair was regarded as consisting of local Hertzian stiffness and tooth bending stiffness.Cai and Tang revised Kuang's calculation formula and proposed an analytical method to calculate the mesh stiffness by considering the load condition and the tooth profile modification in References [9,10].Chen investigated the differences between the rectangular stiffness and its approximate form on the gear nonlinear dynamic characteristics in Reference [11].
According to the available literature, the main methods to calculate mesh stiffness are Finite Element Method (FEM) and the potential energy method.These two methods have advantages and disadvantages.The potential energy method was proposed by Yang [12] and has been improved several times [13][14][15][16].It is widely used in the failure analysis of a gear system [17][18][19][20][21].The FEM, which can reflect the real contact status, is often applied to figure out engineering problems [22][23][24].Moreover, the FEM is a common method to verify the validity of the potential energy method while both the FEM and potential energy method always ignore the influence of backlash on mesh stiffness.
The backlash is very useful to avoid the jam phenomenon in conditions of elastic deformation and thermal expansion, i.e., adaptive backlash is one of the necessary conditions to ensure the normal operation of gear.One way to achieve backlash is to reduce tooth thickness.The reduction in tooth thickness will cause a decrease in the stiffness, further influencing the dynamic performance.In current gear dynamics, the effect of backlash is limited to the definition of the nonlinear backlash function [25][26][27][28][29][30].
Thus, this work will focus on the influence of backlash on time-varying mesh stiffness, carry out nonlinear dynamic analysis, and investigate the dynamic responses for various values of backlash.This is expected to establish the indirect relationship between backlash and dynamic performance as well as direct the design and control of the gear's backlash.

Mesh Stiffness Model with Backlash
Available mesh stiffness models do not consider the influence of backlash.Precision tooth profile equations play an important role in the analysis method to calculate the mesh stiffness.According to the generating method of gear cutting, the tooth profile equations consist of two parts: the involute and transition curve, as shown in Figure 1.The subscripts  and  in Figure 1 denote the involute curve and transition curve.Equations of the involute curve are expressed as The subscripts i and t in Figure 1 denote the involute curve and transition curve.Equations of the involute curve are expressed as where x i , y (1) i is defined as the contact point of the rack cutter and gear blank, and x i + y i / tan α.The relationship between the pressure angle of the arbitrary point at the involute curve α ϕi and ϕ i can be deduced by which can be simplified as The value of ϕ for the involute start point is ϕ Equations of the transition curve are expressed as where x Then, the tooth profile, as shown in Figure 2, can be obtained by the new tooth profile equations as follows: Here, δ = π(1 + 1/Z)/2 − b/(4r), where b represents the backlash, Z represents the tooth number, and it is assumed that the gear and pinion have the same value as the reduction in the thickness of the tooth.
Appl.Sci.2018, 8, 2404 3 of 13 where  ( ) ,  ( ) is defined as the contact point of the rack cutter and gear blank, and The relationship between the pressure angle of the arbitrary point at the involute curve  and  can be deduced by which can be simplified as The value of  for the involute start point is  = * , Thus, the range of  is Equations of the transition curve are expressed as Then, the tooth profile, as shown in Figure 2, can be obtained by the new tooth profile equations as follows: Here, , where  represents the backlash, Z represents the tooth number, and it is assumed that the gear and pinion have the same value as the reduction in the thickness of the tooth.According to the potential energy method, the potential energy deposited in the gear set can be divided into Hertzian energy  , bending energy  , shear energy  , axial compressive energy  , and fillet foundation energy  .The relationship between the potential energy and corresponding stiffness can be expressed as follows:  According to the potential energy method, the potential energy deposited in the gear set can be divided into Hertzian energy U h , bending energy U b , shear energy U s , axial compressive energy U a , and fillet foundation energy U f .The relationship between the potential energy and corresponding stiffness can be expressed as follows: where K h , K b , K s , K a , and K f denote the Hertzian contact stiffness, bending stiffness, shear stiffness, axial compressive stiffness, and fillet foundation stiffness, respectively.The stiffness contribution of the gear body is omitted [31].
Then, the total potential energy U for a pair of teeth can be achieved by summing the energy components, i.e., Then, the single-tooth-pair mesh stiffness for a pair of spur gears with a contact ratio between 1 and 2 can be expressed as follows: where the subscripts i = 1, 2 represent the gear and the pinion.According to Reference [32], the Hertzian contact stiffness K h can be expressed as where E is Young's modulus; L is the tooth width; and υ is the Poisson's ratio.The fillet foundation stiffness K f can be calculated according to [13] K f = 1 where δ ϕ = atan(s(ϕ)) is the angle between normal force and axis x in Figure 2. The parameters L * , M * , P * , Q * can be calculated according to the following equation: where h f i = r f /r int , r int and θ f are described in Figure 2, and A i , B i , C i , D i , E i , F i are listed in Table 1.
Table 1.The values of the coefficients of Equation (12).According to cantilever beam theory, axial compressive energy, bending energy, and shear energy can be calculated by where 3 y 3 L .Further, the bending stiffness K b , shear stiffness K s , and axial compressive stiffness K a can be expressed as

Influence of Backlash on Stiffness
The backlash can be obtained by decreasing the thickness of the tooth or by increasing the center distance.Here, only the former situation will be discussed.
According to ISO/TR 10064-2-1996, the value of minimum backlash can be expressed as where a is the center distance.In this paper, three cases will be discussed, i.e., Case 1: b = 0; Case 2: b = b min ; Case 3: b = 2b min .The parameters of an example gear pair are listed in Table 2.According to the above-mentioned mesh stiffness model, the single-tooth-pair mesh stiffness and time-varying mesh stiffness can be calculated, as shown in Figures 3 and 4.
According to the above-mentioned mesh stiffness model, the single-tooth-pair mesh stiffness and time-varying mesh stiffness can be calculated, as shown in Figures 3 and 4.
where  () is the single-tooth-pair mesh stiffness for different values of backlash, a better view can be obtained to present the effect of backlash on stiffness, as shown in Figure 5. Compared to Case 1, Case 2 and Case 3 shows an obvious reduction in single stiffness and mesh stiffness according to Figure 5.Moreover, the reduction is not proportionate.Compared with Case 1, the rate of decline for Case 2 and Case 3 shows a trend of decreasing first and then increasing in a mesh cycle from coming into the mesh to going out of the mesh.

Mesh Stiffness Fitting Method
In many studies, time-varying mesh stiffness is fitted in the form of the Fourier series expansion By introducing a stiffness coefficient η, where K b (t) is the single-tooth-pair mesh stiffness for different values of backlash, a better view can be obtained to present the effect of backlash on stiffness, as shown in Figure 5.
where  () is the single-tooth-pair mesh stiffness for different values of backlash, a better view can be obtained to present the effect of backlash on stiffness, as shown in Figure 5. Compared to Case 1, Case 2 and Case 3 shows an obvious reduction in single stiffness and mesh stiffness according to Figure 5.Moreover, the reduction is not proportionate.Compared with Case 1, the rate of decline for Case 2 and Case 3 shows a trend of decreasing first and then increasing in a mesh cycle from coming into the mesh to going out of the mesh.Compared to Case 1, Case 2 and Case 3 shows an obvious reduction in single stiffness and mesh stiffness according to Figure 5.Moreover, the reduction is not proportionate.Compared with Case 1, the rate of decline for Case 2 and Case 3 shows a trend of decreasing first and then increasing in a mesh cycle from coming into the mesh to going out of the mesh.

Mesh Stiffness Fitting Method
In many studies, time-varying mesh stiffness is fitted in the form of the Fourier series expansion [33,34].Despite being fitted by a high-level Fourier series, the difference between time-varying mesh stiffness and the fitting one is still very obvious in the alternate mesh area for single and double teeth.Here, the single-tooth-pair mesh stiffness, rather than the time-varying mesh stiffness, will be expressed in the form of the Fourier series expansion: where k 0 is the mean value of the single-tooth-pair mesh stiffness.k aj , k bj are the amplitudes of the jth order harmonic.w f is the stiffness fitting frequency when the mesh period is T = 1/ε.Similarly, the stiffness coefficient can also be expressed in the form of Fourier series expansion: where η 0 is the mean value of the stiffness coefficientη aj , η bj are the amplitudes of the jth order harmonic.w ρ is the coefficient fitting frequency when the mesh period is T = 1/ε.For the example gear pair, it can be found that the fitting effect is very well when γ 1 = 1, γ 2 = 3.The corresponding parameters are listed in Table 3.Moreover, η 0 , η aj , and η bj can be further expressed as functions of backlash: The fitting results are shown in Figure 6a: According to Figure 6b, the relative error is less than 0.1%, which means the Fourier fitting has high precision.
Then, the time-varying mesh stiffness can be obtained by defining two pairs of teeth in mesh.In the meshing process of a gear set with a contact ratio between 1 and 2, as shown in Figure 7, the meshing teeth pair with the contact point located between point A and point C is defined as pair #1, and the meshing teeth pair with the contact point located between point C and point D is defined as pair #2.According to Figure 6b, the relative error is less than 0.1%, which means the Fourier fitting has high precision.
Then, the time-varying mesh stiffness can be obtained by defining two pairs of teeth in mesh.In the meshing process of a gear set with a contact ratio between 1 and 2, as shown in Figure 7, the meshing teeth pair with the contact point located between point A and point C is defined as pair #1, and the meshing teeth pair with the contact point located between point C and point D is defined as pair #2.The piecewise stiffness functions can be expressed as Appl.Sci.2019, 9, 1029 9 of 13

Single-Degree-of-Fredom(SDOF) Model and Dynamic Differential Equation
The SDOF model is shown in Figure 8. Two gears are combined with supporting shafts, which are regarded as rigid bodies with large support stiffness.The influences of friction and radial support are ignored.The pinion and gear are represented by their base circle radii of r bp and r bp , respectively.Then, the equations of motion can be expressed as follows: where I p and I g represent the polar mass moments of inertia, respectively.m p and m g represent the masses of pinion and gear.K(t) represents the time-varying mesh stiffness.
c m = 2ζ k m / 1/m p + 1/m g represents viscous damping, and ζ represents the damping ratio.
e(t) represents the static transmission error.T P and T g represent the input and output torques acting on the pinion and gear.To simplify the calculation, the fluctuations of input and output torque are neglected.The backlash function 'h' can be expressed as By employing dynamic transmission error x(t) = r bp θ p − r bg θ g − e(t), Equation ( 27) can be re-formed as m e ..
with m e =

Numerical Results and Discussion
The dynamic Equation ( 29) can be numerically integrated by using the fourth order Runge-Kutta method.The simulation runs for 20,000 periods and only the data of the last 50 periods are plotted to guarantee the data relates to the steady state conditions.Here, a system with  = 500 N • m,  = 2500 rpm,  = 0.11, () = 0 is established.
Figure 9 shows the dynamic response of the system as  =  , the response period of the time history chart equals the excitation period, the phase diagram presents a closed non-circle and nonelliptic curve, and the discrete points are located at the frequencies of  ( = 1,2, ⋯ ) in the FFT spectrogram, which means that the system exhibits a non-harmonic-single-periodic motion, but the periodicity is not strict.Partially enlarged views of the graphics showed that the phase diagram had 50 curves, the resulting trace in the Poincare map was concentrated in 50 points, and the enlarged

Numerical Results and Discussion
The dynamic Equation ( 29) can be numerically integrated by using the fourth order Runge-Kutta method.The simulation runs for 20,000 periods and only the data of the last 50 periods are plotted to guarantee the data relates to the steady state conditions.Here, a system with T p = 500 N•m, n p = 2500 rpm, ζ = 0.11, e(t) = 0 is established.
Figure 9 shows the dynamic response of the system as b = b min , the response period of the time history chart equals the excitation period, the phase diagram presents a closed non-circle and non-elliptic curve, and the discrete points are located at the frequencies of i f m (i = 1, 2, • • •) in the FFT spectrogram, which means that the system exhibits a non-harmonic-single-periodic motion, but the periodicity is not strict.Partially enlarged views of the graphics showed that the phase diagram had 50 curves, the resulting trace in the Poincare map was concentrated in 50 points, and the enlarged spectrum had obvious perturbation peaks.That is to say, the plotted 50 periods were not strictly repetitive, and the situation was the same as the simulation runs for more periods.Thus, we were certain that the system exhibited a non-strictly non-harmonic-single-periodic vibration.In practical engineering applications, the value range of the backlash is  ≤  ≤ 2 .
Here, the research range of the backlash was expanded to 3 to conduct a better study on the influence of backlash on the dynamics.Figure 10 shows the system response in terms of the dynamic transmission error for various backlash values.In the case of  = 1.25 , the dynamic motion is non-harmonic-single-periodic as a whole.Partially enlarged views indicated that the system presented a low-frequency vibration, as shown in Figure 10a.The phase diagram had 16 curves, the resulting trace in the Poincare map was concentrated in 16 points, and the enlarged spectrum had peaks at the points of  /16.As a result, the system mainly exhibited non-harmonic-single-periodic vibration and also exhibited 16T-periodic low-frequency vibration with a small amplitude.
As the backlash increased to 1.5 , the motion was similar to that of  =  where no significant period of low-frequency vibration was presented.As the backlash further increased, the low-frequency vibrations were 7T-periodic in the case of  = 1.75 and 30T-periodic in the case of  = 2 , respectively.According to Figure 10e-h, it was clear that the low-frequency vibrations were 8T-periodic, 16T-periodic, 11T-periodic, and 7T-periodic, respectively.
The dynamic transmission error is a key parameter to evaluate the gear system's vibration and noise.Available studies have often been concerned with the vibration whose minimum response frequency equals that of the excitation frequency, and those small-amplitude and low-frequency motions have received less attention.According to Figures 9 and 10, it is quite clear that the backlash had a great influence on the small-amplitude and low-frequency vibration, which can be regarded as secondary vibration and can guide us to understanding the mechanism of some special noise.In practical engineering applications, the value range of the backlash is b min ≤ b ≤ 2b min .
Here, the research range of the backlash was expanded to 3b min to conduct a better study on the influence of backlash on the dynamics.Figure 10 shows the system response in terms of the dynamic transmission error for various backlash values.In the case of b = 1.25b min , the dynamic motion is non-harmonic-single-periodic as a whole.Partially enlarged views indicated that the system presented a low-frequency vibration, as shown in Figure 10a.The phase diagram had 16 curves, the resulting trace in the Poincare map was concentrated in 16 points, and the enlarged spectrum had peaks at the points of i f m /16.As a result, the system mainly exhibited non-harmonic-single-periodic vibration and also exhibited 16T-periodic low-frequency vibration with a small amplitude.
As the backlash increased to 1.5b min , the motion was similar to that of b = b min where no significant period of low-frequency vibration was presented.As the backlash further increased, the low-frequency vibrations were 7T-periodic in the case of b = 1.75b min and 30T-periodic in the case of b = 2b min , respectively.According to Figure 10e-h, it was clear that the low-frequency vibrations were 8T-periodic, 16T-periodic, 11T-periodic, and 7T-periodic, respectively.
The dynamic transmission error is a key parameter to evaluate the gear system's vibration and noise.Available studies have often been concerned with the vibration whose minimum response frequency equals that of the excitation frequency, and those small-amplitude and low-frequency motions have received less attention.According to Figures 9 and 10, it is quite clear that the backlash had a great influence on the small-amplitude and low-frequency vibration, which can be regarded as secondary vibration and can guide us to understanding the mechanism of some special noise.

Conclusions
An improved mesh stiffness was proposed by taking the effect of backlash into account.The precise tooth profile equations with backlash were established to generate the gears.The potential energy method was employed to calculate the mesh stiffness.The calculation results indicated that when compared with the case of ignoring backlash, the mesh stiffness had an obvious decrease in the cycle of coming into the mesh to going out of the mesh.Moreover, the rate of decline showed a trend of decreasing first before increasing.
To obtain an accurate stiffness function for numerical analysis, the mesh stiffness rather than the time-varying mesh stiffness was fitted in the form of Fourier series expansion, and the time-varying mesh stiffness was presented in terms of a piecewise function.The proposed fitting mesh used a low

Conclusions
An improved mesh stiffness was proposed by taking the effect of backlash into account.The precise tooth profile equations with backlash were established to generate the gears.The potential energy method was employed to calculate the mesh stiffness.The calculation results indicated that when compared with the case of ignoring backlash, the mesh stiffness had an obvious decrease in the cycle of coming into the mesh to going out of the mesh.Moreover, the rate of decline showed a trend of decreasing first before increasing.
To obtain an accurate stiffness function for numerical analysis, the mesh stiffness rather than the time-varying mesh stiffness was fitted in the form of Fourier series expansion, and the time-varying mesh stiffness was presented in terms of a piecewise function.The proposed fitting mesh used a low order Fourier series to obtain very high accuracy.The stiffness coefficient was also fitted by the Fourier series.
The effects of backlash on dynamic transmission error were investigated, and the numerical results revealed that the gear system primarily performed a non-harmonic-single-periodic motion.The partially enlarged views indicated that the system also exhibited small-amplitude and low-frequency motion.
For different cases of backlash, the low-frequency motion sometimes showed excellent periodicity and stability and sometimes showed a non-repetitive and aperiodic situation.
The theory here focused on the impact of backlash on mesh stiffness and dynamic performance and established an indirect relationship between backlash and dynamic transmission error, which is expected in order to understand the mechanism of some unusual noise and to guide the design and manufacture of backlash in the future.

Figure 2 .
Figure 2. The geometric model of the gear profile.

Figure 2 .
Figure 2. The geometric model of the gear profile.

Figure 5 .
Figure 5.  for the different values of backlash.

Figure 5 .
Figure 5.  for the different values of backlash.

Figure 5 .
Figure 5. η for the different values of backlash.

Figure 6 .
Figure 6.The single-tooth-pair mesh stiffness: (a) Comparison of calculated stiffness and fitted stiffness.(b) Relative error curve.

Figure 6 .
Figure 6.The single-tooth-pair mesh stiffness: (a) Comparison of calculated stiffness and fitted stiffness.(b) Relative error curve.

Figure 7 . 26 ) 3 .
Figure 7.A snapshot of the contact pattern (at t = 0) of a spur gear pair.

Figure 7 .
Figure 7.A snapshot of the contact pattern (at t = 0) of a spur gear pair.

13 Figure 9 .
Figure 9.The system response in terms of the dynamic transmission error for  =  .

Figure 9 .
Figure 9.The system response in terms of the dynamic transmission error for b = b min .

Table 2 .
The parameters of the example gear pair.

Table 3 .
The fitting parameters for the example gear pair.