Dynamic Modelling Considering Nonlinear Factors of Coupled Spur Gear System and Its Experimental Research

For the nonlinear vibration problem of spur gear transmission system caused by friction and gear backlash, the dynamic model of gearbox with 16-degree-of-freedom (16-DOF) considering time-varying meshing stiffness is established by Lagrange equation. The dynamic equations of the gear system are solved by the Newmark- $\beta $ method. The effects of friction coefficient, error fluctuation and meshing stiffness on the vibration response are investigated. In order to verify the validity of the model, a spur gear transmission rotor test bench was built by simulating the actual working conditions. The root mean square (RMS) is used to verify the accuracy of the model at multiple speeds, and makes up for the shortcomings of many researches without experimental argumentation. The experimental research can provide a reference for the research of gear transmission system.


I. INTRODUCTION
The application range of the gear rotor system is very wide. It plays an important role in the metallurgy, aerospace, chemical industry, electric power system, and machinery manufacturing. With the need of high reliability, high speed and mute, the higher demands are brought to the vibration characteristics of the gear rotor system. The nonlinear vibration characteristics have become one of the hot topics in the research of the machinery industry [1]- [3]. Therefore, it is of great scientific significance and application value to research the nonlinear vibration of gear transmission system by reliable dynamic modelling.
In the recent years, many meaningful researches have been carried out by scholars. With the continuous development of nonlinear dynamics [4], engineering mathematics and computer technology [5], the gear system model has also The associate editor coordinating the review of this manuscript and approving it for publication was Hamid Mohammad-Sedighi . undergone major updates. In the early research, scholars only considered the influence of single factors on the dynamics of the gear system, and then established the single degree of freedom model [6]- [9]. Kahraman and Singh [6] analyzed the nonlinear frequency response characteristics of the spur gear pair with backlash under internal and external excitation. Souza et al. [7] research the nonlinear dynamic behavior of the gear pair model under sinusoidal drive backlash. Rocca and Russo [8] explored the effects of periodic backlash fluctuations on gear vibration. However, the actual gear transmission system is more complicated. This singledegree-of-freedom gear model considering only a single factor does not reflect the rich dynamic behavior of the actual system. Based on this, many scholars research the multi-degree-of-freedom model [10]- [14]. Theodossiades and Natsiavas [10] considered the influence of gear meshing clearance, static transmission error, bearing clearance and contact characteristics on the nonlinearity of the system, established a 2-DOF system model. Oğuz and Fatih [11] proposed a 4DOF cracked gear model based on dynamic transmission error (DTE), and studied the influence of cracks along the tooth thickness direction on the dynamic transmission error of the gear system. Chen [12], [13] et al. developed a 5-DOF of the gear-bearing model, which was used to analyze the effect of friction and dynamic backlash on the system. Xiang et al. [14] studied a 5-DOF of the gear-bearing model with time varying stiffness, gear backlash and surface friction.
In order to better guide engineering practice, researchers need to consider the combined effects of system factors such as friction coefficient [15], meshing stiffness and bending-torsional coupling vibration [16], as well as pitting, spalling and cracking [11], [17], and thermo-structural coupling [18]- [20], on the gear system. This need to create a more reliable gear transmission system model to better explain the complex nonlinear phenomena [21]- [22]. Mohammed [23] established a 12-DOF model for the firststage reducer by introducing the gyro effect and crack fault in the model. Li and Kahraman [24] introduced the mixed elastohydrodynamic lubrication model into the dynamic model considering transverse torsional vibration. Ma [25], [26] introduced time-varying meshing stiffness and crack propagation path into the gear transmission model respectively. Zhou et al. [27] developed the nonlinear dynamics model of gear-rotor-roller bearing system. The effects of speed and bearing clearance on the dynamic characteristics were presented, while the influence of friction was not considered.
In summary, many great achievements on gear system had been made. However, most of the researches are focused on the low degree of freedom model and the process of modeling or analysis is relatively simple. The vibration coupling effect between gear and transmission shaft is not considered. Most importantly, there are no experimental verification in most literatures. For the whole gear transmission system, the coupling relationship has an important influence on the dynamic characteristics of the system, but the relative work is relatively small. Based on this, the paper considering eccentricity, tooth face friction, transmission error and backlash, a nonlinear dynamic model of gearbox with 16 DOF has been established by the theory of gear dynamics and nonlinear dynamics. What's more, the influence law of the tooth surface friction, the error fluctuation and the meshing stiffness on the system is studied through the dynamic model. Last but not least, the model is verified by a gear experiment platform. The analysis results can provide a theoretical basis for the dynamics research of gearbox system. The contribution of this paper is to present an effective 16-degree-of-freedom (16-DOF) gears model. Scholars can use this model, or combine it with engineering data to study the dynamics of a class of gear systems.

II. DYNAMIC MODEL OF GEAR ROTOR BEARING SYSTEM
The dynamic model of gear rotor bearing system takes the effect of time-varying meshing stiffness, the transmission error and the tooth face friction into account as shown in Fig. 1. In the model, there is no axial movement. The rotation centers of gears are defined as O 1 (x 1 , y 2 ), O 2 (x 2 , y 2 ), and the centers of mass are G 1 (x g1 , y g1 ), G 2 (x g2 , y g2 ). The angle displacement of the input/output and gears can be defined as ϕ i (t) (i = d, 1, 2, g), and the corresponding variation displacements are θ i . ω 1 and ω 2 are the angular velocities of the driving shaft and driven shaft respectively. According to the geometrical relationships, the ϕ i can be expressed as: It is shown in Fig. 1(a) that B i (i = 1, 2, 3, 4) are the geometric center of bearing. A 1 , A 2 are the ideal centers of the driving and driven gear, which are located on the ideal centerline of bearing bore. m bi (i = 1, 2, 3, 4) represent the equivalent mass of each bearing. m i (i = 1, 2) are the equivalent mass of the gears, and the moments of inertia corresponding to the center are defined as I i . J j (j = d, g) represent the moment of inertia of the input and output respectively. The distance from the center of gear to the center of bearing is written as l ij , l i are the lengths of the shaft. The meshing gear is simplified as the stiffness k m and damping c m in the direction of the meshing line, as shown in the Fig. 1(b). The eccentricity of the gears is expressed as e 1 , e 2 . r bi (i = 1, 2) are the base radii of gears.

A. GEAR ROTATION CENTER POSITION
The dynamic position relationship between the driving gear and the bearings at both sides can be observed in Fig.2. Where, B 1 A 1 B 2 indicates the initial state of the driving shaft. The driving shaft will be shifted and bent due to the system vibration when the driving gear rotating. The rotating center of support bearings on both sides of the shaft will move from the theoretical centers to the actual rotating centers, and the rotation center of the gear will move from A 1 to O 1 . The B o1 O 1 B o2 is the dynamic location of the driving shaft. The points of B o1 , O 1 , B o2 are projected onto the x-axis and the y-axis. Where, Based on the deformation coordination conditions, the elastic deformation amount of driving gear in x-axis and the y-axis can be expressed as δ 1 x , δ 1 y respectively. Similarly, the elastic deformation amount of driven shaft can be also obtained. Therefore, the elastic deformation of the driving shaft and the driven shaft of the gear system in the x -axis and y -axis can be described as follows: where, τ j = l ij /l i (i = 1, j = 1, 2; i = 2, j = 3, 4). In this paper, the value of l ij are the same, and the lengths of shaft are also same.

B. THE FRICTION FORCE
The friction of the tooth surface is a potential internal excitation, which mainly comes from the sliding and rolling friction of the teeth. The amount and direction of friction will change periodically with the meshing point, which makes it become an important nonlinear factor affecting the dynamic behavior of gear system. Fig. 3 displays the schematic diagram of the meshing process of a pair of gears. The speed of the driving gear is ω 1 . The meshing point maintains a constant motion state along the N 1 N 2 direction at an absolute speed of ω 1 r 1 . Therefore, s 1 and s 2 can be expressed as The direction of the friction force is perpendicular to the meshing line. The sliding speed of gears will be zero when the meshing point is at the pitch line. With the meshing process continuing, the direction of the relative speed of the gears will change. It will also change the direction of the tooth friction force. Therefore, the coefficient of friction direction is timevarying and can be written as Based on the coulomb's law, the friction force [28] can be expressed as where, F m is dynamic meshing force and can be described specifically as where δ is the relative displacement of the meshing point between the teeth, k (t) m is the meshing stiffness, shown as where, k m and k r are mean and amplitude of the meshing stiffness, its specific values are shown in Table 1; ϕ m is initial position angel; ω m is the meshing frequency. In Eq. (6), c m is the meshing damping, shown as where, ξ m is the meshing damping ratio, its specific is shown in Table 1.
In Eq. (6), f (δ) is a nonlinear function representing backlash fluctuation. The backlash is usually a discontinuous and non-differentiable time-varying function. It will cause the system to generate strong nonlinear vibration. It can be expressed as where, b is half of the backlash, its specific is shown in Table 1.

C. GOVERNING EQUATIONS
The generalized coordinate vector of the nonlinear dynamic model can be described as Therefore, the governing equations of the nonlinear dynamic model with 16-DOF can be derived by Lagrange's equation, as shown in Eq. (11).
where, k ti and k si (i = 1, 2) are the torsional stiffness and lateral stiffness of the shaft; c ti and c si (i = 1, 2) are the corresponding damping coefficients. The support bearings are simplified as the radial stiffness and damping. k bxi , k byi , c bxi , and c byi (i = 1, 2, 3, 4) represent the stiffness and damping of the bearings in the x-direction and y-direction [29], [30]. T d / T g is the input / output torque, respectively. The specific parameters of the model are presented in Table 1.

III. NUMERICAL ANALYSIS A. THE EFFECT OF FRICTION COEFFICIENT
Among the parameters affecting gear transmission, the tooth surface friction has become an indispensable research object. It is mainly caused by the sliding and rolling friction of the gear, which has a great influence on the smoothness and vibration of the gear transmission. In this section, the friction coefficient u is as the control parameter, and set ω 1 = 1200 r/min, namely the basic  frequency of driving shaft f 1 = 20 Hz. The meshing frequency f m = z 1 ×f 1 = 55×20 = 1100Hz. In order to analyze the effect of friction coefficient on the gear rotor system specifically, the vibration responses of the system at two friction coefficients of u = 0, u = 0.1 are analyzed, the specific dynamic response of x 1 and y 2 are shown in Fig. 4 and Fig. 5 respectively. It is shown that when considering friction, the vibration amplitude of x 1 increases significantly, but has little effect on y 2 . The vibration response of x 1 and y 2 is obviously different. This is because the gear is also affected by the meshing force in the y direction in addition to the eccentric force.
In order to research the influence of friction u on the dynamic characteristics of the system more accurate, the 3-D frequency spectrums of x 1 and y 2 are obtained, as shown in Fig. 6(a) and Fig. 6(b) respectively.
It can be seen rom Fig. 6(a) that, the main frequency components of x 1 are 8.5f 1 (170Hz), 17f 1 (340Hz), 0.5f m (550Hz), f m (1100Hz), 1.5f m (1650Hz). The f 1 and 0.5f m are in the whole frequency range. The amplitude of 0.5f m is almost linearly increasing with the increase of friction coefficient. There are continuous spectrum on both sides of 8.5f 1 (170Hz). For y 2 , the meshing frequency f m exist in the whole frequency range, and its amplitude does not change significantly in Fig. 6(b). The amplitude of 17f 1 shows a gradual upward trend throughout the interval. However, the amplitude of 19 f 1 is gradually decreasing. When 0.1≤ u ≤0.3, the combined frequency 19f 1 +0.5f m can be observed.

B. THE EFFECT OF ERROR FLUCTUATION
Due to the influence of complex factors such as manufacturing, machining, installation error and tooth surface wear, the uncertainty of the error will have a great influence on  the vibration characteristics and stability of the gear system. In this section, the vibration response with error is analyzed. The rotation speed is set as 1200 r/min and friction coefficient is set as 0.1. Three errors e r conditions are compared as e r = 2.0×10 −5 , e r = 6.0×10 −5 , e r = 2.0×10 −4 . Fig. 7 and Fig. 8 show the vibration displacements in the x 1 and y 2 directions, respectively. As the error increases, the vibrations in the two directions are more intense, and the equilibrium point in the x 1 direction changes greatly. VOLUME 8, 2020  In order to more clearly explain the frequency domain characteristics in the x 1 and y 2 directions. The Fig. 9(a) and Fig. 9(b) display the 3-D frequency spectrums of the x 1 and y 2 . When e r ∈[1.0 × 10 −4 , 3.0 × 10 −4 ], the main frequency components of x 1 are 8.5f 1 (170Hz), 19f 1 (380Hz), 0.5f m (550Hz), f m (1100Hz), 1.5f m (1650Hz). Meanwhile, there is a weak frequency component 19f 1 + f m (1480Hz). Where, 0.5f m occupies the dominant frequency, and its amplitude decreases first and then increases with the error, and the amplitude is the smallest when e r = 3.0 × 10 −5 . The frequency component of y 2 in Fig. 9(b) is different from x 1 , which is due to the different excitations from the two directions. Meshing frequency f m is the main frequency. Other frequency components include 10.5f 1 (210Hz), 0.5f m , 1.5f m (1650Hz) and 8.5f 1 +0.5f m .

C. THE EFFECT OF MESHING STIFFNESS
The meshing stiffness of the gear is one of the main internal excitations, which has a great influence on the transmission precision. In this section, the vibration responses will be analyzed under three meshing stiffness of k m = 2.5 × 10 8 N/m, k m = 4.5 × 10 8 N/m and k m = 9.0 × 10 8 N/m.
The relative dynamic response of x 1 and y 2 is shown in Fig. 10 and Fig. 11 respectively. It can be seen that as the vibration in both directions becomes more intense with meshing stiffness increase.
From the above research, the change of system parameters will have a great influence on the amplitude of the characteristic frequencies. In order to quantitatively study the effect of the meshing stiffness, the amplitudes of rotation frequency of the driving gear f 1 in x 1 direction, set as x 1 (f 1 ) are analyzed when the meshing stiffness change from 2.0 × 10 8 N/m to 9.0 × 10 8 N/m. Similarly, x 1 (f m ), y 2 (f 2 ) and y 2 (f m ) are also analyzed. The corresponding results are normalized, as shown in Fig. 12. It can be seen from Fig. 12(a) that as the meshing stiffness increases, the amplitude of x 1 (f 1 ) increases monotonously, while the amplitude of x 1 (f 1 ) exhibits a trend of increasing first and then decreasing. The amplitude is the largest when meshing stiffness is 6.5 × 10 8 N/m. The amplitudes of x 1 (f m ) and y 2 (f m ) increase monotonically with the increase of meshing stiffness in Fig.12 (b), and the trend of increasing of x 1 (f m ) is more obvious.
Through the above analysis, the meshing stiffness is an important parameter. Therefore it can be concluded that under the condition of meeting the strength requirement, the meshing stiffness k m should be appropriately reduced to decrease the vibration and improve the stability of the system.

IV. EXPERIMENTAL VERIFICATION
The gear equipment and the position of the test point are shown in Fig. 13. The gearbox adopts oil bath lubrication   and the motor provides power to drive the transmission shaft. The drive shaft and input shaft are connected by a coupling, the output shaft of gearbox is connected with the magnetic powder clutch brake. The horizontal and vertical vibration signals of the input and output shafts of the gearbox are respectively measured. The laser sensor and the eddy current sensor are close to the bearing end covers, and connected with the NI9234 acquisition card. Finally, the experimental data can be read in the LABVIEW program. The main parameters of the experimental gear pair are shown in Table 2.
In order to verify the validity of the model, the displacement of the driving gear in x direction and the driven gear in y direction with respect to simulation and experiment   are compared at 1200 r/min. The relative results are shown in Fig. 14 and Fig. 15 respectively. For the simulation, the vibration of the driven gear in the y direction is larger than that of the driving gear in the x direction. That is similar with the experiment. A periodic ''high and low peak staggering'' phenomenon occurs in the simulation and experiment of the driving gear in x direction.
The displacement signals are performed with FFT, and the relative spectrum are obtained, as shown in Fig. 16 and Fig. 19. It can be seen in Fig. 16(b) that the spectrum in the x-direction of the driving gear includes frequency components 0.5f m (552.3 Hz) and f m (1103 Hz). Near the 0.5f m , there is frequency components 0.5f m -2f 1 (514 Hz) and 0.5f m +2f 1 (594.2 Hz). The reason is that the gears are slightly eccentric during the manufacturing process, and the drive shaft has a slight misalignment. In addition, due to interference from other factors, complex spectral bands are generated in the low-order frequency region. It can be seen in Fig. 17(b) that the spectrum of the driven gear in y direction includes frequencies such as 4f 2 (59.38 Hz), 15f 2 (304.7 Hz),  0.5f m (547.3 Hz) and f m (1102 Hz), which is largely consistent with the simulation results in Fig. 17(a).
From the above analysis, the results of simulation and experiment results are basically the same. In fact, due to some human factors, the amplitude of the experiment and the simulation may different, but the trend can also effectively verify the correctness of the theoretical model.
In dynamic analysis, the root mean square (RMS) is often used to reflect the vibration response. In order to further prove the effectiveness of the theoretical model at multiple speeds, the RMS is used to compare the vibration trend with rotation speed.

RMS
wherex is the average, x(i).  The values of RMS of the driving gear in x direction and the driven gear in y direction at 900 r/min, 1050 r/min, 1200 r/min, 1350 r/min and 1500 r/min are calculated respectively. In order to compare the trends of RMS in simulation and experiment conveniently, the results are normalized, as shown in Fig. 18. The overall change trend of the simulation results and the experimental results is relatively consistent. The amplitudes increase with the increase of the rotation speed. It can further prove the validity of the theoretical model.
In order to better verify the correctness of the model, the characteristic frequencies of the driving gear in x direction and the driven gear in y direction at the above five speeds are further analyzed. The experimental meshing frequencies and simulation meshing frequencies of the driving gear in x direction at the five speeds are compared in Table 3. The maximum error of the meshing frequencies are 2.0 Hz, 2.65 Hz, 3.0 Hz, 2.0 Hz, and 3.0 Hz respectively. The relative maximum error of the meshing frequency of the driven gear in y direction are 1.5 Hz, 3.8 Hz, 2.7 Hz, 5.75 Hz and 5.5 Hz respectively, as shown in Table 4. In summary, the meshing frequency of the simulation and the experiment are basically the same, which can further verify the correctness of the simulation results.

V. DISCUSSION
Through the above analysis, the trend of the simulated waveform and the experimental waveform has a good similarity, and the main frequency components are also generally consistent, which verifies the correctness of the theoretical model. However, the experimental signal contains more information, and there is a certain gap between the theoretical value and the experimental value. The main sources of these errors can be classified as: 1) There is a certain deviation between the theoretical calculation point and the actual measurement point. The theoretical calculation point is at the gear node, but the experiment point is to take a position closer to the gear.
2) The gear manufacturing error and the randomness of the installation accuracy make the quality and eccentricity in the experiment different with the theoretical model. Moreover, due to the wear of the teeth and the lubrication of the oil film between the teeth, the gears are subject to some variable excitation.
3) The influence of the components such as the box and the bearing housing is not considered, and the vibration generated by these components during the experiment is bound to be transmitted to the gear.

VI. CONCLUSIONS
In this paper, the nonlinear dynamic model of gearbox with 16 DOF is established through the Lagrange equation, which takes tooth face friction, eccentricity and transmission error into account. The influence law of friction coefficient, error fluctuation and meshing stiffness on the nonlinear dynamics of spur gearbox system is analyzed by numerical simulation. Firstly, with the increase of the friction coefficient, the amplitude of 0.5f m of driving gear in x direction has increased. However, the amplitude of meshing frequency f m of driven gear in y direction does not change significantly. Secondly, with the increase of error, the frequency components are relatively complex and the amplitudes of characteristic frequency also increase. Furthermore, the amplitude of x 1 (f 1 ) increase first and then decrease with the meshing stiffness.
Lastly, the dynamic model of gear system is verified by the experiment. The difference between the experiment and simulation is compared with the RMS method and characteristic frequencies. Results show a good similarity between the theoretical simulation and the experiment, which proves the validity of the theoretical model. The theoretical model, the test bench and the measurement method, have certain guiding significance for the subsequent research of gear transmission system.