Coupled Bending-Torsional Nonlinear Vibration and Bifurcation Characteristics of Spiral Bevel Gear System

A spiral bevel gear system supported on thrust bearings considering the coupled bending-torsional nonlinear vibration is proposed and an eight degrees of freedom (8DOF) lumped parameter dynamic model of the spiral bevel gear system combined with timevarying stiffness, static transmission error, gear backlash, and bearing clearances is investigated. The spiral bevel gear system is analyzed with the equations of motion and the dynamic response is solved using the Runge-Kutta method. The effects of mesh frequency, mesh damping coefficient, load coefficient, and gear backlash are revealed, which describe the true mesh characteristics of the spiral bevel gear system.The bifurcation characteristics as jump discontinuities, periodic windows, and chaos are obtained by studying time histories, phase plane portraits, Poincaré maps, Fourier spectra, and global bifurcation diagrams of the gear system. The results presented in this study provide some useful information for engineers in designing and controlling such gear systems.


Introduction
As a typical gear transmission system, the spiral bevel gear system has been used in main reducer of automobile driving axle; particularly dynamic response is one of the most important factors affecting NVH performance; its characteristic research has attracted the attention of domestic and foreign scholars.Accordingly, there are many related literatures studying the spiral bevel gear system in the past several decades.
In earlier years, many fundamental researches were focused on the linear analysis of the spiral bevel gear system [1].A two degrees of freedom vibration model of a pair of bevel gears was established by Kiyono et al. [2]; the model was applied to conduct a stability analysis, in which the line of action vector was simulated by a sine curve.Kahraman and Singh [3] derived a single degree-of-freedom model considering the constant stiffness; the dynamic equations with backlash and transmission error were presented and solved in the involute gear model.Gosselin et al. [4] analyzed the static transmission error of spiral bevel gears, and the effects of the shape and amplitude of the unloaded transmission error curve on the loaded dynamic behaviors were demonstrated.
In order to get a basic understanding of dynamic behavior, a lot of researches had been focusing on nonlinear vibration analysis.Litvin et al. [5] presented tooth contact analysis and stress analysis in spiral bevel gears by means of finite element method.Tang et al. [6] studied the effect of static transmission error on nonlinear dynamic response of the spiral bevel gear system combining with time-varying stiffness, gear backlash, and observed various nonlinear phenomena including periodic solutions, bifurcations, and chaos.Yang et al. [7] established a single degree-of-freedom hypoid gear pair dynamic equation, which included the time-varying stiffness, transmission error, and backlash, and obtained the FFT responses.Wang et al. [8] described a generalized nonlinear time-varying (NLTV) dynamic model of a hypoid gear pair with backlash nonlinearity and time-dependent mesh point, line of action, mesh stiffness, and kinematic transmission error.Periodic motions were obtained by the incremental harmonic balance method (IHBM).Chang-Jian [9] performed dynamic analysis of bevel-geared rotor system supported on a thrust bearing and journal bearings under nonlinear suspension.Theodossiades and Natsiavas [10] established a simplified dynamic model of motordriven gear pair, considering the gear backlash and bearing clearance; then numerical results are presented in the form of classical frequency response diagrams, revealing the effect of the system parameters on its dynamics.Cheng and Lim [11] pointed out that the gear kinematic transmission error was the primary source of the vibratory energy excitation; a new analytical derivation of the hypoid gear-meshcoupling mechanism based on the simulation of tooth contact assuming idealized gear geometry was proposed.However, only in recent years did the related studies in finding nonlinear behaviors of spiral bevel gear system gain some attention.The dynamic analysis of a spiral bevel-geared rotorbearing system was studied by Li and Hu [12].The modeling of coupled axial-lateral-torsional vibration of the rotor system geared by spiral bevel gears was discussed.Different degrees of freedom (DOF) gear dynamic models were implemented by Mohammed et al. [13]; their limitations were evaluated by simulating different DOF for each gear disc.Wang et al. [14] established a nonlinear dynamic model of the spiral bevel gear with the dynamic relative transmission error, backlash, and time-varying stiffness.The vibration displacement and velocity in the torsional, horizontal, and vertical directions in the spiral bevel gear model under different conditions were depicted, and the dynamical responses of the geared system with harmonic internal excitation and parameter excitation were obtained.From the literatures above, the nonlinear characteristics of gear system such as stability, periodic solutions, bifurcations, and chaos have become the most interesting research areas.Different nonlinear parameters will generate an obvious change on the dynamic response.However, the bifurcation characteristics researches of nonlinear dynamic parameters as gear backlash and bearing clearances seem a little deficient; the dynamics analysis combined with thrust bearings clearances considering the coupled bending-torsional vibration is also rarely seen.In this paper, a nonlinear dynamic model of the spiral bevel gear system is formulated, where the timevarying stiffness, static transmission error, gear backlash, and bearing clearances are included.The dimensionless equations of the system are then solved using the Runge-Kutta numerical method.The influence of the nonlinear parameters on the spiral bevel gear system is studied and the nonlinear mesh characteristics are detected and analyzed by construction of the time histories, phase plane portraits, Poincaré maps, Fourier spectra, and global bifurcation diagrams.

Mathematical Modeling and Equations of Motion
Considering the supported stiffness of a spiral bevel gear system is large, so the twist vibration can be neglected.The complex spiral bevel gear system is simplified by the lumped masses method; Figure 1 shows a generalized dynamic model for eight degrees of freedom (8DOF) considering the coupled bending-torsional vibration.The gear system is modeled with rotational and translational displacements as their coordinates.The generalized coordinates vector of the nonlinear dynamic model can be expressed as where  1 ,  1 , and  1 are the translations of the pinion along the axes , , and ;  2 ,  2 , and  2 are the translations of the gear along the axes , , and ;  1 and  2 are the torsional displacements of the pinion and gear, respectively.
Here, the term "pinion" refers to the smaller gear, which is a driver gear connected to the input shaft, and the term "gear" refers to the larger gear, which is a driven gear connected to the output shaft.
The pinion has mass  1 and moment of inertia  1 ; the gear has mass  2 and moment of inertia  2 .
Static transmission error () can be expressed in the form where  0 is constant amplitude of static transmission error,   is variable amplitude of static transmission error, Ω ℎ is excitation frequency, and Φ  is phase angle.Since the stiffness is periodically time-varying with the mesh frequency, its analytical formulation can be obtained by means of a Fourier expansion [15]: Here,   is mean value of the mesh stiffness,   is stiffness fluctuation amplitude, and Φ  is phase angle.
The backlash function () and the bearing clearances functions (  ), (  ), and (  ) ( = 1, 2) can be written as where 2 represents the total backlash; 2 1 and 2 2 represent the bearing clearances in horizontal plane; 2 3 represents the bearing clearances in vertical plane.
The dynamic mesh force along the line of action   can be expressed as The dynamic mesh force along the coordinate directions   ( = , , ) can be expressed as where  1 = sin   cos  1 + cos   sin  sin  1 ,  2 = sin   sin  1 − cos   sin  cos  1 ,  3 = cos   sin ,   is mean value of the mesh damping,   is pressure angle,  is helix angle, and  1 is cone angle of the pinion.
Here,  is damping ratio.
From the proposed concept, the equations of motion of the coupled bending-torsional vibration model can be derived as where   ( = , , ;  = 1, 2) are damping coefficients of the supported structure for the thrust bearings;   ( = , , ;  = 1, 2) are stiffness coefficients of the supported structure for the thrust bearings;   ( = 1, 2) are mean load torques on the pinion and gear, respectively;   ( = , , ) are dynamic loads along the axes , , and  for the pinion and gear, respectively;  1 is base radius of the pinion and   is distance between the acting point of the normal force and the center of the rotation of the gear.The motion equations of the spiral bevel gear system in matrix form can be expressed as where  is lumped mass matrix;  is damping matrix;  is stiffness matrix;  is external excitation force vector.
And the equivalent mass and static load of gear transmission are Next, introducing the characteristic length   and assuming the following set of dimensionless parameters, , Input torque fluctuation is not considered in the current study, but only static torque load, which means  V = 0.In order to reduce the complexity of the computation, higher orders of the time-varying stiffness and static transmission error are truncated, without losing generality.So the dimensionless stiffness and the dimensionless static transmission error are calculated by the first-order component: where  is stiffness coefficient;  is amplitude of static transmission error.Equations ( 8)-( 15) can be expressed in a matrix form: ) ) . (21)

Numerical Solutions
Due to the complexity of the spiral bevel gear system and also the difficulty and limitation of the analytical methods, the numerical method is commonly used to analyze the gear system.In this paper, the nonlinear dynamic equations are solved using the fourth order Runge-Kutta method which is generally applicable to strong nonlinearity [17].The algorithm is implemented in MATLAB that is a widely used matrix and numerical analysis program [18].Without losing generality, the characteristic length is set to   = 10m, and the same value can be found in Sun's thesis [19].The other parameters of the gear transmission system are given in Table 1.

Effect of Mesh Frequency.
A bifurcation diagram summarizes the essential dynamics of the gear system and is therefore a useful means of observing its nonlinear dynamic response [20].For our subsequent numerical study, set half of the backlash as  = 40 m, damping ratio as  = 0.08, amplitude of static transmission error as  = 0.3, and stiffness coefficient as  = 0.2.In this system, the dimensionless mesh frequency  ℎ is commonly used as an important control parameter.Figure 2(a) shows the bifurcation diagram for the spiral bevel gear system displacement against the dimensionless mesh frequency  ℎ .Figures 2(b) and 2(c) show the enhanced bifurcation diagrams with bifurcation parameters  ℎ = 0.65∼1.2and  ℎ = 1.4∼2.15,respectively.The bifurcation results show that the system behaves as 1T-periodic motion at low mesh frequency and the periodic motion persists until  ℎ > 0.425.A jump phenomenon can be observed and 2T-periodic motion appears at the region  ℎ = 0.425∼0.8.After a long 1T-periodic motion, the system enters into a transient chaotic motion after  ℎ > 1.15; then the dynamic behavior behaves as nonperiodic motion until  ℎ = 1.5.
With the increase of mesh frequency, a number of periodic windows appear.At the region  ℎ = 1∼2, 1T-periodic motion, 2T-periodic motion, and T-periodic motion are shown in the chaotic area.It can be observed that T-periodic bifurcation enters into chaos and 7T-periodic bifurcation enters into T-periodic bifurcation obviously in Figure 2(c).Finally, 2T-periodic motion transits to 1T-periodic motion until  ℎ = 2.05.For a better clarity, the other analytical methods for observing nonlinear dynamic responses are necessary.Figures 3-5 illustrate the time histories, phase plane portraits, Poincaré maps, and Fourier spectra for the spiral bevel gear model at various values of the dimensionless mesh frequency  ℎ = 1,  ℎ = 1.6, and  ℎ = 1.77.At  ℎ = 1, the dynamic behavior behaves as 1T-periodic motion.At  ℎ = 1.6, chaotic behavior is clearly visible.At  ℎ = 1.77, the system enters into 6T-periodic motion from chaotic motion.

Effect of Damping Ratio.
In the following section, when the damping ratio  is changed and all the other parameters keep fixed, the bifurcation diagrams are shown in Figure 6.The chaotic motions are shown through the numerical simulation for three different values of  = 0.01,  = 0.04, and  = 0.12 ( = 0.08 see Figure 2(a)).According to the comparison of the bifurcation diagrams with different damping ratio values, it can be observed that chaotic region becomes narrower as  is increased; chaotic region gets smaller and smaller.Therefore, the increase in  tends to decrease the nonlinearity dynamic responses.

Effect of Load Coefficient.
Different vibration characteristics will be shown under different loads.For heavy loaded condition, the value of the load coefficient parameter is fixed as  =   / = 2 (  = 0.1,  = 0.05).For light loaded condition, the value of the load coefficient parameter is fixed as  =   / = 1/3 (  = 0.1,  = 0.3).The time histories, phase plane portraits, Poincaré maps, and Fourier spectra are shown in Figures 7 and 8.These figures show dynamic responses for heavy loaded condition and light loaded condition at the dimensionless mesh frequency  ℎ = 1.4,respectively.As illustrated in Figure 7, the system is linear system; the steady-state response is harmonic response.There is obviously no tooth impact, and the system behaves as 1T-periodic motion.It is illustrated that the dynamic characteristics of the gear are not changed much under loaded condition when the dynamic behavior behaves as periodic motion.As illustrated in Figure 8, compared with the heavy loaded condition, the magnitude of the spectrum becomes larger.The displacement along the line of action is greater than 0 and is less than 0, so the case of tooth impact can be observed in the time history plot.Figures 9 and 10 show dynamic responses for heavy loaded condition and light loaded condition at the dimensionless mesh frequency  ℎ = 1.8, respectively.Nonlinear vibration characteristics can be seen occurring in these two situations.However, the intensities are not the same.As illustrated in Figure 9, the motion state is 1T-periodic harmonic motion under heavy loaded condition.It can be seen that the spiral bevel gear system is in chaotic motion under light loaded condition and the response of the system is irregular in the chaotic state from Figure 10.Furthermore, chaotic phenomenon is more likely to occur in light loaded condition compared with heavy loaded condition without changing other conditions.

Effect of Gear Backlash.
Gear backlash is another important parameter which affects the dynamic responses substantially.Figure 11 presents the bifurcation diagrams for three  different values of half of the backlash  = 10 m,  = 80 m, and  = 120 m ( = 40 m; see Figure 2(a)).A substantial region of jump phenomenon, subharmonic, quasi-periodic, and chaotic responses are contained in Figure 11.With increasing b, the system becomes more and more unstable, and consequently chaotic behavior expands.
For the further study on the effect of gear backlash, the time histories, phase plane portraits, Poincaré maps, and Fourier spectra of the dynamic responses at  ℎ = 2.23 for numerous gear backlash values are compared in Figures 12-15.According to these figures, a transition from quasiperiodic motion to the chaotic dynamics is seen with the increase of .As observed in Poincaré maps, the various forms of the system from the quasi-periodic motion are seen, and with the increase of the bifurcation parameter, the unstable attractive region is enlarged, eventually leading to chaos.

Conclusion
In the present paper, a spiral bevel gear system supported on thrust bearings considering the coupled bending-torsional nonlinear vibration is proposed and an 8DOF lumped parameter dynamic model of the spiral bevel gear system combining with time-varying stiffness, static transmission error, gear backlash, and bearing clearances is investigated.The dimensionless equations of the system are solved using the Runge-Kutta numerical method.The dynamics of the system are analyzed with reference to its bifurcation diagrams, time histories, phase plane portraits, Poincaré maps, and Fourier spectra, and jump phenomena, periodic motions, and chaotic motions are found in this study.
In this study, mesh frequency is first expressed to study the effects of their variations on dynamic response.With the change of mesh frequency, the repeated bifurcation phenomena and the complex chaotic motions are observed, and the evolution process and rule of the system are revealed.The effect of damping ratio is secondly analyzed.It is shown that chaotic region becomes narrower as damping ratio is increased.Other important parameters like load coefficient and gear backlash are also considered in the current analysis.Load coefficient affects light loaded condition more than heavy loaded condition.The increase of gear backlash will expand the effect of nonlinearity effect and increase the chaotic region of the gear system.
The selection of the dynamic parameters plays an important role in designing and controlling such gear systems; the suitable values should be chosen so that chaotic behavior will 14 Shock and Vibration be decreased.A more precise model to better describe the spiral bevel gear system will be proposed and deep study will be analyzed in the future.

Figure 1 :
Figure 1: Dynamic model of a spiral bevel gear system.

Table 1 :
The parameters of the spiral bevel gear system.