Chaos Analysis of Single-Stage Spur Gear System Considering Backlash Fractal

Considering the time-varying pressure angle and dynamic clearance, the effects of rotational speed, tooth surface friction, and tooth surface morphology on the system’s dynamic response are studied. An improved gear nonlinear model is proposed, which considered nonlinear factors such as time-varying pressure angle, position angle, tooth surface morphology, and tooth surface friction. The time-varying dynamic backlash is deduced, and the nonlinear dynamic equation of the gear is established. The nonlinear dynamic response of the gear system is obtained based on Runge–Kutta method. The influence of the rotational speed and tooth surface friction on single-stage spur gear system response is analyzed through the analysis of the bifurcation diagram, the three-dimensional spectrum diagram, the proportion of the meshing state, and the largest dynamic meshing force (LDMF). Dynamic response differences between the three different model are compared. In addition, by changing the tooth surface roughness and fractal dimension, the influence of tooth surface morphology on the dynamic response of the system is studied. Compared with the traditional model, when the dynamic center distance and dynamic pressure angle are considered, the system response may enter a chaotic state earlier. When the tooth surface friction is further considered, the chaotic state of the system response is suppressed. At the same time, the velocity of the dynamic transmission error is significantly reduced, and the fluctuation amplitude of the dynamic pressure angle is increased. The value of LDMF rose overall. The stability of the system response decreases with the increase of tooth surface roughness and fractal dimension. Compared with the fractal dimension, the tooth surface roughness has a more obvious effect on the dynamic response of the system.


Introduction
The increased application of micro-/nano-structured materials in modern nano-systems has stimulated a great deal of interest in size-dependent mechanics [1][2][3].The dynamic characteristics of small-scale structures are studied by solving the recently developed theory of dimensional-dependent elasticity.Since nano-structures may exhibit both stiffening and softening behaviors, the nonlocal elasticity model and modified strain gradient theory cannot solely represent the entire wide spectrum of the size-dependent phenomena at the nanoscale [4].Faghidian [5] mostly investigated the nonlocal strain gradient model, nano-scale effects of the dilatation gradient, deviatoric stretch gradient, and the symmetric rotation gradient tensors along with the nonlocality were suitably considered.The proposed constitutive model of the nonlocal modified gradient beam is demonstrated to be capable of capturing both the higher-order gradients and the nonlocality effects.Feigao et al. [6] evaluated thermoelastic damping (TED) in circular plates by incorporating nonlocal effects within the constitutive and heat conduction frameworks.By considering symmetric time-harmonic vibrations, the size-dependent thermoelastic frequency equation was derived.According to the definition of TED in the framework of the complex frequency approach, a closedform expression characterizing TED in circular nanoplates was introduced.Faghidian et al. [7] proposed a nonlocal modified gradient elasticity theory, enriched with four intrinsic length scales, which can appropriately characterize the effects of nonlocality, dilatation gradient, deviatoric stretch gradient, and symmetric rotation gradient.Faghidian et al. [8] proposed a mixed variational principle associated with the higher-order unified gradient elasticity in the intrinsic form and studied the mechanics of torsion of nano-bars.The conceived stationary variational formulation can realize the higher-order gradient effects while being exempt from the restrictions associated with the nonlocal elasticity model by limited accessible nonlocal kernels mostly devoted to onedimensional structural analysis.An effective practical methodology to calibrate the characteristic parameters associated with the higher-order unified gradient elasticity theory was presented.Thai et al. [9] presented a nonlocal strain gradient mesh-free model based on the nonlocal strain gradient theory and higher-order shear deformation theory and developed it to examine the bending and free vibration behaviors of functionally graded nanoplates.A detailed study to assess the effects of exponential factor, length-to-thickness ratio, geometry, boundary condition, strain gradient parameter, and nonlocal parameter on the deflections and natural frequencies of the FG nanoplates was revealed by numerical results.For the first time, Faghidian et al. [10] realized the dimensional correlation response of porous nanowire of symmetric fiber grating in the framework of high-order non-local gradient elasticity theory.The integrodifferential constitutive laws of the stress resultant fields were established and reinstated with the equivalent differential relations equipped with non-standard boundary conditions.The closed-form solution of the phase velocity of flexural waves was analytically determined.The ensuing numerical results of the flexural wave dispersion detected new benchmarks for numerical analysis.
Gear is the most widely used mechanism in mechanical transmission, the nonlinear characteristics of the gear system are very important through the precision transmission of the gear system.The chaos and the dynamic characteristics of the gear system have been widely studied [11,12].With the development of size-dependent mechanics, the influence of the morphology of the gear surface on the transmission stability of the gear system is emphasized.Chen et al. [13] studied the influence of gear friction and dynamic clearance on nonlinear dynamic gear systems.The results show that the dynamic backlash increases the vibration amplitude of the system and causes the system to respond more easily to a chaotic state.Chen et al. [14] developed a pure torsional model of gears with dynamic backlash based on fractal theory and proposed a method to describe the dynamic backlash, which compared with the backlash generated by applying the normal distribution.The research shows that fractal theory describes the variation of the tooth side clearance more exactly.Liu et al. [15] developed a 16-degreeof-freedom coupled transverse torsional gear-rotor bearing transmission system with the crack faults considering the fractal backlash and the torsional excitation was found to have a significant effect on the forced response spectrum.Huang et al. [16] developed a pure torque model to study the bifurcation and chaos of the system under two clearance cases, fractal backlash and fixed clearance considering the friction case.They found that the fractal backlash would make the system's dynamic characteristics more complex.In addition, the bearing clearance also affects the transmission performance of the gear system seriously.Because of the bearing clearance in the transmission state, the center distance and meshing pressure angle of the gear pair will change.
So far, much research has been operated on bearing and gear clearance, but few found studies on coupled dynamics were searched.Yi et al. [17] set up a gear system model in which the time-varying pressure angles and central distances were considered and the influence of time-varying pressure angles on the dynamic characteristics of the system was studied.Liu et al. [18] analyzed gear-bearing systems with dynamic central distances and backlash using the Newmark-β method.Zhang et al. [19] established a gearbearing model and analyzed the amplitude of dynamic transfer error (DTE) in the time domain considering multiple clearances.
However, due to the simplification of the previous gear dynamics model, nonlinear factors such as position angle, dynamic center distance, and dynamic pressure angle are neglected, the coupling characteristics of backlash with fractal characteristics are not reflected.In this paper, the change in tooth surface morphology which is modified by tooth surface roughness is described by the Weierstrass-Mandelbrot equation.The nonlinear factors such as dynamic pressure angle, dynamic center distance, and position angle in gear meshing are further studied.The meshing state of the gear is described by meshing force and DTE.The effects of rotational speed, tooth surface friction, tooth surface roughness, and fractal dimension in the meshing state and nonlinear dynamic response are studied concurrently.

Centralized Mass Model
Neglecting the effects of friction on the tooth surface, the lumped mass model of the gear system is shown in Fig. 1.
The involute spur gear bending-torsion coupling dynamic model is built using the centralized parameter method.k h and c h are the support stiffness and support damping of each gear in the cartesian coordinate system, respectively.m i ( i = 1, 2 ) is the gear mass., 2b is the bearing clearance of the gears, and 2b 1 (t) is the backlash between the two gear meshing, r bi (i = 1, 2) is the base radius of gear i .The DTE of the sys- tem is as follows: where x i (i = 1, 2) is the horizontal displacement of the main wheel and the driven wheel, y i (i = 1, 2) is the displacement in the vertical direction of both the main and driven wheels, and i (i = 1, 2) represents the gear rotation angle around the drive shaft, is the dynamic pressure angle, is the position angle.

Dynamic Backlash with Fractal Characteristics
There are objective randomness, multi-scale and self-affinity in the machining process of various types of gear surfaces with different precisions so that they have certain fractal characteristics.It is verified that the fractal theory is effective to describe the machined surface profile [13], so the Weierstrass-Mandelbrot function can be used to simulate the trend of the concave and convex height of the gear surface morphology.The height of surface asperities is as follows: where is the characteristic scale coefficient, D is the fractal dimension, and the value of D is between 1.1 and 1.9.
This formula can reflect the variation of concave and convex heights at a certain fractal dimension.But it does not reflect the actual value of the concave and convex heights.Chen [20] proposed a correction to the above formula which is as follows: (1) where is the correction coefficient, R a is the arithmetic mean difference of the evaluated profile, which is the most widely used parameter to evaluate roughness [20], R a is used to evaluate the corrected roughness.So, the value of the correction factor was defined as follows: R ac (D) is the arithmetic mean difference of the contour height corresponding to the target roughness.R ac (D) is the corresponding profile roughness for different fractal dimension simulation cases.The two-dimensional simulation results of the rough surface profile are shown in Fig. 2.
Therefore, the dynamic clearance of the gear pair is derived as follows:

Multi-clearance Coupling
The dynamic backlash caused by bearing clearance may lead to complex changes in the meshing forces.As shown in Fig. 3.
The initial center distance of the gear is expressed as follows [17]: where M i (i = 1, 2) is the gear modulus, Z i (i = 1, 2) is the number of gear teeth.
In the gear system, the gear backlash and bearing clearance will be coupled with each other, and the overall movement of the gear will cause the change of the center distance due to the bearing clearance, which affects the meshing line .
(5) Fig. 2 The two-dimensional simulation results of the rough surface profile and the actual pressure angle in turn.Then the actual center distance of the gear is expressed by [18]: Therefore, the actual pressure angle is as follows: According to the meshing principle, the actual clearance of the two gears is as follows [21]: where t ′ is the dynamic tooth pitch on the gear pitch cir- cle, s � i (i = 1, 2) is the dynamic tooth thickness on the two gear pitch circles, respectively; 0 is the standard pressure angle for gears.
Therefore, the dynamic backlash considering the fractal clearance is: It can be seen that considering the fractal dimension, roughness, bearing axis displacement and other nonlinear functions, the dynamic backlash is expressed.Compared with the traditional gear-bearing collision model, time-varying pressure angle and tooth-side fractal backlash and other nonlinear factors are taken into account for the subsequent study of gear models as a preparation.

Engagement Model Considering Time-Varying Pressure and Position Angles
The gear dynamics parameters include time-varying meshing stiffness, meshing damping, dynamic backlash, etc.The gear meshing is described by a nonlinear stiffness k 1 (t) and viscous damping c 1 (t) .To meet the demand of friction of the gear drive system, the meshing stiffness of the tooth surface is described by using the method proposed in Ref. [22].The meshing stiffness is expressed as: (8) = arccos((r b1 + r b2 )∕a) where k m is the averaged tooth stiffness, k c (t) is the vari- ant part of the stiffness and it can be expanded in the Fourier series generally.
The static transmission error e 1 (t) is a displacement excita- tion related to the manufacturing and installation accuracy of gears.It can usually be expressed by [16]: where A e is the error fluctuation and m is the frequency of gear meshing.
The cycle of gearing meshing T can be expressed by: The meshing damping of the gears is given by [18]: where k m is the average value of the time-varying mesh- ing stiffness, 2 is the gear meshing damping ratio, usually taken as 0.03-0.17and the average value of 0.1 is taken in this paper.
Position angle can be expressed by [17]: q denotes the dynamic transfer error (DTE): Then the backlash nonlinear function of DTE is as follows: 3 Geometrical relationship between gear backlash and center distance 1 3 The bearing clearance nonlinear function is as follows: where represents x i and y i (i = 1, 2).
During gear meshing, the meshing force F p (t) consists of the elastic meshing force due to the time-varying meshing stiffness and the viscous meshing force due to the meshing damping.Then the meshing force is expressed as: Here is the contact ratio of the gear pairs, ′ rounds the elements of to the nearest integers towards plus infinity.The friction forces can be expressed by: where i is the direction coefficient of the friction force which is determined by the nominal sliding velocity.F m is the gear meshing force.i is the nonlinear friction coefficient as proposed in Ref. [22].
The friction arm of the two gears can be represented as follows [13]: where r ai (i = 1, 2) is the top circle radii of gear.S i (t) (i = 1, 2) is the friction arm of the gear surface.i (i = 1, 2) is the angu- lar speed of gear.
During the meshing process, the pressure angle and the position angle of the time change continuously.Therefore, the direction in which the meshing force and the friction force act on the gears are constantly changing.The force balance equation of the two gears in the X and Y directions can be expressed as follows: where is the angle between the gear engagement line and the Y direction.can be expressed by: The equilibrium equation of the two gears in the torsion direction can be expressed as follows:

Comparison and Discussion
The Runge-Kutta method is applied in this section to solve the problem above, by rounding the first 2000 cycles to eliminate the effect of the transient.A pair of parallel shaft involute spur gears are used as an example for comparison and verification of these models.The bearing clearance is 20 µm and the initial gear backlash is 50 µm.Other parameters are listed in Table 1.

Affection of Rotational Speed
Rotational speed is one of the most important parameters that influence gear dynamics.The global characteristics are studied by the bifurcation diagram, the proportion of meshing states, the three-dimension frequency spectrum map and largest dynamic meshing force (LDMF), in which it is assumed that the surface roughness R a is 0.8, and the fractal dimensions D 1 and D 2 are equal to 1.1.
Three models are compared, namely, the traditional model (i.e., the model that considers tooth surface morphology but not nonlinear factors such as time-varying pressure angle and position angle), the compared model (i.e., the compared model that considers tooth surface morphology along with nonlinear factors such as time-varying pressure angle and position angle), and the improved model (i.e., the improved model that further considers friction based on the compared model).Since the improved model takes more factors into account than the traditional model, the dynamic behavior of the gear system could be reflected more realistically.
To study the meshing force of the three models, the LDMF of the tooth surface is defined as follows: where F s is the static meshing force, defined as the ratio of the average input torque to the gear base radius.
Compared with the traditional model, the DTE of the improved model and the compared model have increased in magnitude.Figure 4c shows that with the amplitude of the DTE increases overall, tooth meshing of the improved model and the compared model are both 100% during the range of 2000-2520 rpm.For the traditional model, when the rotational speed is 2000 rpm, the tooth surface meshing state accounts for 99.747% of the entire meshing process.When the rotational speed is 2520 rpm, the tooth surface meshing state accounts for 98.672% of the entire meshing process.And there is no 100% tooth meshing state during the rotational speed from 2000 to 2520 rpm.Compared with the traditional model, the tooth back meshing appears earlier in the improved model and the compared model.In the improved model, tooth back meshing appears at 6080 rpm for the first time, accounting for 0.0234% of the entire meshing process.In the compared model, tooth back meshing appeared at 6040 rpm for the first time, and the proportion of tooth back meshing in the entire meshing process was 0.0260%.In the traditional model, tooth back meshing occurs at 9280 rpm for the first time with a ratio of 0.3906%.
It can be seen that when the rotational speed reaches Ω 3 = 3360 r/min, the response of the traditional model becomes chaotic.The compared model entered chaos at Ω 2 = 2960 r/min.It shows that the response of the sys- tem will enter chaos earlier when the dynamic angle and dynamic center distance are considered in the model.But the improved model enters chaos for the first time when the rotational speed is Ω 1 = 3340 r/min.It shows that chaos of the system is suppressed when the tooth surface friction is considered in the model.In addition, the dynamic response of the conventional model exhibits a multi-cycle response when the rotational speed reaches 13,800 rpm and enters a chaotic state again at 14,040 rpm.Dynamic response of the compared model exhibits a multi-cycle response up to F s , 14,240 rpm and enters a chaotic state again at 14,480 rpm.Interval lags for multiperiod responses.Dynamic response of the improved model is a multi-cycle response when the rotational speed reaches 14,160 rpm and enters chaos again when the speed is 14,480 rpm.After considering friction, the multi-cycle response interval becomes larger.Figure 4d shows that although the magnitude of the dynamic transfer error of the system in the compared model is greater than in the traditional model, the LDMF in the compared model is not always greater than that in the traditional model, which is due to the dynamic pressure angle and the dynamic center distance changing over time throughout the meshing process.The size of the tooth backlash continues to change.Add tooth surface friction to the model.Figure 4d shows that the LDMF size of the improved model improved overall compared to the comparison model.It is indicated that tooth surface friction would increase LDMF of the tooth surface.
To compare the dynamic response of the three models is defined as a dimensionless time.It can be expressed by: Dynamic response of the three models at the rotational speed of 8000 r/min are shown in Fig. 5.The green dots represent the Poincaré mapping points and the Poincaré mapping points of the compared model's system dynamic response are more dispersed, compared with the response of the traditional model.Meanwhile, the dynamic response of the traditional model shows the quasi-periodic response, while the dynamic response of the improved model has started to become chaotic.Pressure angle of the traditional model is constant during the whole meshing process, while the pressure angle of the compared model is changed with time.Dynamical response of the compared model is more complex than the traditional model.It is necessary to consider time-varying pressure angle and time-varying center distance when studying gear system response.The system response transforms into an approximate multi-cycle response after further consideration of friction in the model shown in Fig. 5a, which means that stability of the system response improved.
When the improved model was compared with the first two models, the velocity of DTE decreased significantly because of tooth friction.And amplitude of the dynamic pressure angle increases due to the coupling of friction with the forces in the X and Y directions.
Through the above analysis, we can know that the improved model can reflect the dynamic response of the gear system more refined than the traditional model because several nonlinear factors have been taken into account such as time-varying pressure angle, dynamic center distance, and friction.
The fractal dimension D is varied from 1.1 to 1.9.During the following analysis, only D= 1.1,1.5, 1.9 will be listed for compactness, and we assume that the tooth surface roughness and the fractal dimension of both gears are equal.
As can be seen in Fig. 6.Although the surface roughness of the gear did not change, the amplitude of DTE gradually increase from 41.3107 m , 42.3190 m to 45.4078 m with the increase of fractal dimension and amplitude of the dynamic pressure angle gradually increase from 0.0097 • , 0.0102 • to 0.0105 • .The blue dots represent the Poincaré mapping points.The orbits in the phase diagram change from thin to thick, and as the fractal dimension increases, the Poincaré mapping points gradually disperse.Therefore, it is known that for the same tooth surface roughness, the stability of the gear system becomes worse as the fractal dimension increases.
Based on the conclusion above, we can improve the dynamic stability of the gear system by changing the fractal dimension of the tooth surface without changing the roughness.In practical engineering, when the roughness no longer be reduced, the stability of the dynamic response of the gear system can be improved by reducing the fractal dimension through techniques such as micromachining.Figure 7 shows that as the roughness of the gear surface increasing, the orbits in the phase diagram gradually change from regular to irregular, and the dynamic response of the system changes from quasi-multi-periodic response to chaotic motion.With the increase of roughness, Lyapunov exponent of the system increases from − 2.2143 to 4.2156 and 5.4352.The blue points in Fig. 7a are the Poincaré mapping points and they gradually become dispersed from concentrated with the increase of tooth surface roughness.The number of frequencies of the system response increases significantly, and the stability of the system decreases.The change of dynamic pressure angle also becomes chaotic from the periodic response, and the amplitude of dynamic pressure angle increases with the increase of tooth surface roughness.From Figs. 6c and  7c, it can be easily seen that when the tooth surface roughness increase, the system response would be more chaotic, comparing with the fractal dimension increase of the tooth surface.

Gear Dynamics Influenced by Surface Roughness
It is evident from this analysis that tooth surface roughness has a great influence on the transmission stability of the gear system.The influence on the transmission stability of the gear system is more obvious than the fractal dimension.Therefore, when considering the improvement of gear transmission stability, priority should be given to reducing tooth surface roughness.

Conclusion
Based on the fractal theory, an improved nonlinear model of gear pair with the characteristics of the tooth surface morphology was built, considering the influence of dynamic pressure angle, position angle, and tooth surface friction during gear meshing.Some beneficial conclusions are obtained as follows: 1. Considering the dynamic pressure angle and position angle during gear meshing, the effect of tooth surface morphology on the dynamic response of the system can be better studied when straight cylindrical gears are performing high-precision transmission.2. Considering the tooth surface morphology, if dynamic pressure angle, dynamic center distance, and position angle to the meshing process are considered, the improved model will enter a chaos state earlier and the DTE amplitude of the system will increase synchronously.Once friction is added to the model, the chaotic motion of the system will be suppressed.Meanwhile, LDMF of the gear system will uplift overall and the amplitude of the dynamic pressure angle will increase.In addition, the velocity of DTE drops significantly.We also find that the improved model has 100% tooth surface meshing in the range of 2000-2520 rpm, which is not seen in traditional models.Compared with the traditional model, the tooth-back meshing state of the improved model appears earlier explicitly.3. The stability of the system response would be increased by the decrease of the tooth surface roughness and fractal dimension, but compared with the influence of the fractal dimension, tooth surface roughness would make a greater impact on the stability of the system.

Fig. 1
Fig. 1 Dynamic model of the gear-bearing system Moreover, the main parameters are set as follows:D = 1.1 , b 0 = 10 m , T M = 30N m , 2 = 0.1 , = 7500rpm .We assume separately that R a = 0.8, 1.6, 3.2 .During the tooth surface roughness and the fractal dimension of the two gears the following relationship exists:D 1 = D 2 = D , R a1 = R a2 = R a .

Fig. 4
Fig. 4 Dynamic feature with varying rotational speed

Fig. 5
Fig. 5 Dynamic feature with fractal backlash around n = 8000 r/min

Table 1
Main parameters of the gear pair