# Modelling and Dynamic Analysis of the Spiral Bevel Gear-Shaft-Bearing-Gearbox Coupling System

Authors: Haimin Zhu; Weifang Chen; Rupeng Zhu; Jie Gao; Meijun Liao
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

To accurately study the dynamic characteristics of the spiral bevel gear transmission system in a helicopter tail transmission system, the finite element model of the gear shaft was established by a Timoshenko beam element, and the mechanical model of the spiral bevel gear was created by the lumped mass method. The substructure method is employed to extract the dynamic parameters from the gearbox’s finite element model, and the dynamic model of the spiral bevel gear-shaft-bearing-gearbox coupling system was built according to the interface coordination conditions. In the model, the influences of time-varying stiffness, a time-varying transmission error, gearbox flexibility, unbalance excitation, and a flexible shaft and bearing support on the system vibration were taken into account simultaneously. On this basis, the dynamic differential equations of the full coupling system of the spiral bevel gear were derived, and the effects of the gearbox flexibility, the shaft angle, and the unbalance on the dynamic properties of the system were analysed. The results show that the gearbox flexibility can reduce the gear meshing force and bearing force, in which there is a more significant impact on the bearing force. The shaft angle affects the position, size, and direction of the system’s axis trajectory. Meanwhile, the meshing force and the bearing force of the system are also varied because of the various pitch angles of the driving and driven gears under different shaft angles. The unbalance of the gear shaft has an effect on the vibration of the spiral bevel gear transmission system in all directions, wherein the influence on the torsional vibration is the most significant, and the influence increases as the unbalance rises. The unbalance of the gear shaft also affects the meshing force and bearing force, which increases as the rotational speed rises. This research provides a theoretical basis to optimize dynamic performance and reduce the vibration and noise of a spiral bevel gear full coupling system.


---

## Body

## 1. Introduction

The spiral bevel gear has the characteristics of the firm bearing capacity, a large coincidence degree, and stable operation. Its noise and vibration are relatively small during high-speed operations, and it is widely used in the core transmission systems of helicopters and ships. The most common form of the spiral bevel gear transmission system is the orthogonal shaft system. Because of the limitations of space and shape, the nonorthogonal spiral bevel gear is widely employed in helicopter tail transmission systems, in which the shaft angle is not 90, resulting in the vibration behaviours of the gear pair and the rotor system varying accordingly. In addition, the stiffness of the gearbox is relatively high in general, so the influence of the gearbox flexibility can be ignored when studying the dynamic characteristics of the gear system. However, the wall thickness of the helicopter gearbox is comparatively thin, and the flexible support brought by the thin-wall structure of the gearbox causes a dynamic coupling effect between the gear transmission system and the gearbox. Therefore, in order to obtain more accurate system response and better predict the vibrational noise of such systems, it is necessary to establish a nonorthogonal spiral bevel gear-shaft-bearing- box coupling mechanical model to analyse the dynamic properties of the system.Currently, the dynamics of the spiral bevel gear transmission system has been studied extensively, both domestically and internationally. Wang et al. [1] established a nonlinear vibration equation for the bending-torsional-axial coupling vibration of a 7-degrees-of-freedom spiral bevel gear. The torsional, lateral, and axial vibration displacements and velocities of the spiral bevel gear system under various working conditions are obtained by the operator algorithm. The results show that with the variety of the meshing frequency, the system enters into chaos through period-doubling bifurcation; with the variety of the support stiffness, the system changes into chaos through quasi-period bifurcation; and in the process of the meshing frequency varying, jump phenomena exist in the system. Feng and Song [2] considered the time-varying meshing force of the gear, the position of the meshing point, the transmission error, the direction of the friction force, and the friction coefficient, and then established a dynamic model of concentrated parameters. On this basis, the effects of the dynamic meshing force and dynamic transmission error on the system response are investigated. Wang et al. [3] analysed the torsional vibration of the spiral bevel gear transmission system with backlash and time-varying meshing stiffness. The results showed that when the gear system reciprocates, the backlash and time-varying meshing stiffness will produce a coupling effect, causing the gear system to show obvious nonlinear characteristics. The internal excitation frequency, the meshing stiffness amplitude, the initial condition, and the damping ratio all have an essential impact on the system’s dynamic features. The variation of each parameter makes the system’s response appear in harmonic period solution and chaotic solution. Tang et al. [4], considering such influencing factors as the time-varying meshing stiffness and gear clearance, analysed the effect of the static transmission error on the nonlinear dynamic response of a spiral bevel gear system. The dynamic response, bifurcation diagram, time domain response, phase diagram, and Poincaré map are all obtained by solving the dynamic equations. Fan et al. [5] examined the linear dynamics of the transmission system of a nonorthogonal spiral bevel gear and analysed the system’s inherent characteristics, steady-state response, and axial trajectory. Cheng and Lim [6] built a three-dimensional multi-degrees-of-freedom model of the hypoid gear and analysed the dynamic coupling and vibrational response sensitivity using the response function of the meshing force. Then, the influences of the system parameters on the dynamic response of the hypoid gear were studied, and it was proven that the design parameters and working conditions have significant impacts on the modal and vibration response. Yang and Lim [7] considered the nonlinear clearance and time-varying meshing stiffness and solved the dynamic meshing force and vibration response by the variable step size Runge–Kutta method. The analysis showed that damping could dramatically reduce the dynamic meshing force. Wang et al. [8] established a hypoid gear dynamics model with a time-varying meshing force considering the effect of the nonlinear clearance and studied the impact of the different meshing waveforms on the dynamic response of gear dynamics. Xu et al. [9, 10] constructed the bending-torsion-axial coupling dynamic model of the multi-freedom spiral bevel gear through the lumped mass method and then studied the chaotic and bifurcation properties of the system. Yassine et al. [11] constructed a dynamic model of a two-stage bevel gear transmission system, mainly considering the effect of the time-varying meshing stiffness. The Newmark numerical algorithm was utilised to solve the equations, and the dynamic features of the system were analysed in the time domain and the frequency domain. Peng [12] created a multibody dynamic model of hypoid gears and studied the effects of misalignment, torque, and eccentricity on the dynamic response of the system. Li et al. [13] investigated the effect of an asymmetric time-varying meshing stiffness on the spiral bevel gear transmission system. It was determined that the meshing stiffness has a higher impact on the driving side than that on the driven side. Under a light load, the rise in the meshing stiffness on the drive side deteriorates the vibration stability of the transmission system. Wu and Yang [14] utilised the extended transfer matrix method to analyse the bending-torsional-coupled vibration of the multi-degrees-of-freedom damped shafting caused by external excitation. It was shown that the natural frequency of bending-torsion coupling is close to the natural frequency of torsion. Choi and Mau [15] analysed the dynamic behaviours of the gear-coupled rotor system by using the transfer matrix method. The Timoshenko beam element was employed to model the rotating shaft, and the natural frequency and vibration mode were acquired by studying the bending-torsion-coupled vibration of the system. Neriya et al. [16] adopted the finite element method to simulate the bending-torsion coupling vibration of the gear system comparatively early. The finite element method has been widely utilized by scholars both domestically and abroad because this method can better simulate the influence of a flexible shaft and the bearing parameters on the dynamic response of the gear pair [17–20]. Li and Hu [21] discussed the linear dynamic properties of the bending-torsional coupling of the rotor-bearing system with spiral bevel gears via the finite element method after considering the generalized meshing relationship between the conical teeth. Hua et al. [22] employed the finite element method and the lumped parameter method to analyse the dynamics of the spiral bevel gear-rotor system. Yavuz et al. [23], aking the backlash, time-varying meshing stiffness, and the transmission error into account, built a dynamic model of the spiral bevel gear based on the finite element method, solved differential equations using the harmonic balance method, and studied the effects of backlash, time-varying meshing stiffness, and a bearing stiffness on the system dynamics. Choy et al. [24, 25] examined the effect of the housing’s vibration on the dynamic features of the parallel shaft gear transmission system via the modal synthesis method, but the gearbox was simplified dramatically. Ren et al. [26] established the gear-shaft-bearing-box coupling dynamic model of the parallel shaft via the impedance synthesis method, taking the dynamical coupling between the transmission system and the housing into consideration, and investigated the influence factor of the natural frequency, meshing force, and bearing force. According to the literature retrieved, the modelling methods of the gear transmission system’s dynamics mainly contain the lumped mass method, the transfer matrix method, and the finite element method. Among them, the finite element method not only considers the influence of the shaft flexibility and bearing support but also conveniently couple the forces in the surrounding environment, which is especially suitable for the complex geared rotor coupling systems. Currently, the application of the finite element method is very mature for gear systems with parallel shafts, such as spur gears and helical gears. However, there have been relatively few research works on spiral bevel gear systems with nonparallel shafts, and the effect of thin-walled gearbox flexibility on the dynamics of the transmission system has not been taken into account.In this article, a nonorthogonal spiral bevel gear transmission system supported by a thin-walled gearbox is taken as the research object. First, analysis of the transient dynamics is carried out using finite element software to obtain the time-varying meshing stiffness of the spiral bevel gear pair because the tooth surface of the spiral bevel gear is a complex curved surface. Then, the dynamic equation of the nonorthogonal spiral bevel gear transmission system is constructed via the finite element method and the lumped mass hybrid modelling method. The dynamic parameters of the gearbox system are extracted by the substructure method and are connected by the nodes of the bearing support, thus establishing the dynamic model of the spiral bevel gear-shaft-bearing-box coupling system. On this basis, the dynamic equations of the whole coupling system are solved using the Newmark numerical algorithm, and the effects of the gearbox flexibility, the shaft angle, and the shaft’s unbalance value on the dynamic properties of the coupling system are analysed.


## 2. Mechanical Model of the Spiral Bevel Gear Full Coupling System

The spiral bevel gear full coupling system is a gear-shaft-bearing-box system, which is divided into two parts: the spiral bevel gear transmission system and the gearbox system. As shown in Figure1, it contains a hollow input shaft, a hollow output shaft, a pair of nonorthogonal spiral bevel gears, four ball bearings, and a gearbox connected with the gear transmission system. The coordinates O‐XYZ in Figure 1 are the global coordinate of the entire coupling system.Figure 1
Full coupling system of the spiral bevel gear. (a) The three-dimensional model of the spiral bevel gear fully coupled transmission system. (b) Gear transmission system. (c) Gearbox system.
(a)
(b)
(c)The three-dimensional model of the system is simplified, and a mechanical model of the full coupling transmission system of the spiral bevel gear is constructed. Because of the distinct angle between the driven gear shaft and the driving gear shaft, the vibration displacements of the driven gear and its rotating shaft need to be converted to the global coordinate system through a coordinate transform. The mechanical model of the full coupling system of the spiral bevel gear includes the flexible rotating shaft unit, the meshing unit, the bearing support unit, and the gearbox unit. First, the nodes of the beam elements are divided along the axis at the abrupt change of the shaft section, the bearing support, the gear meshing point, the power point, etc. The gearbox system, because of its complex structure, is more difficult to model directly using the finite element method. The substructure method is employed to convert the dynamic parameters of the box element to the connection node of the bearing support, and the bearing element is applied to connect the gearbox node and the rotating shaft node at the bearing position. Then the finite element method is utilised to assemble each element matrix, and the overall motion equation of the system with the displacement of each node as a generalized coordinate can be obtained. Through the above analysis, the vibration problem of a full coupling system of spiral bevel gears with a continuous mass distribution can be reduced to a vibration problem of a system with limited degrees-of-freedom. Figure2 shows the mechanical model of a full coupling system of the nonorthogonal spiral bevel gear, in which nodes 1–9 are the nodes corresponding to the driving shaft’s beam element, nodes 10–18 are the nodes corresponding to the beam element of the driven axle, and nodes 19–22 are the nodes of the gearbox unit acquired by the substructure method.Figure 2
Mechanical model of the nonorthogonal spiral bevel gear fully coupled transmission system.


## 3. Dynamic Equation of the Spiral Bevel Gear Full Coupling System

### 3.1. Calculation of the Meshing Stiffness of the Spiral Bevel Gear

The gear’s meshing stiffness refers to the synthetic effect of each pair of gear teeth involved in meshing during the meshing process. It is difficult to calculate the meshing stiffness of the gear teeth using the theoretical analysis method because the tooth surface of the spiral bevel gear is a space meshing gear. Meanwhile, the tooth thickness is variable. Therefore, this study utilises finite element software to perform dynamic analysis and obtain an accurate meshing stiffness.The comprehensive elastic deformationδn and contact force Fn are calculated when multiple pairs of gears mesh at the same time, and then the comprehensive meshing stiffness can be obtained directly. The comprehensive elastic deformation of the gear teeth includes the contact elastic deformation due to the local Hertz contact δH, the displacement of the contact position of the gear teeth caused by bending of the gear teeth δb, and the deformation of the bearings, shafts, and supporting structures on the contact positions δf. The deformation of the bearings, shafts, and supporting structures has a relatively small influence on the meshing stiffness. Therefore, this study only takes δH and δb into account, so the deformation of the gear at a certain moment can be expressed as follows:δn=∑i=12δHi+∑i=12δbi.First, the gear’s meshing period is divided inton segments, and the meshing stiffness at each time point Kii=1,2,3,…,n is calculated [27, 28]. Then, a Fourier transform is utilized to convert a series of discrete values Ki into a meshing stiffness’s curve. In this article, the main parameters of the nonorthogonal spiral bevel gear are shown in Table 1, and then the comprehensive meshing stiffness is analysed. In order to improve the efficiency of calculation, nine pairs of gear teeth are selected from the spiral bevel gear pair for the meshing analysis. Figure 3 shows the finite element analysis of the meshing stiffness of the spiral bevel gear, including the meshing of the gear pair and the contact analysis of gear pair.Table 1
Main design parameters of the pinion and gear.
Pinion
Gear
Young’s modulus (Pa)
E = 2.1e9
Poisson’s ratio
N = 0.3
Pressure angle (°)
22.5
Spiral angle (°)
35
Transverse module (mm)
m = 5.6
Width of teeth (mm)
B = 33
Shaft angle (°)
126
Damping ratio
ξ
=
0.05
Mean of the transmission error (m)
e
m
=
2
×
10
−
5
Amplitude of the transmission error (m)
e
r
=
3
×
10
−
5
Pitch radius (mm)
72.5
87.6
Number of teeth
N
1 = 29
N
2 = 35
Hand
Right
LeftFigure 3
Finite element analysis of the meshing stiffness of the spiral bevel gear pair. (a) The meshing of the gear pairs. (b) Contact analysis of the gear pairs.
(a)
(b)Each periodT is divided into ten segments on average, and two meshing periods are taken. Then, each Ki value at its corresponding time ti is calculated. By performing a Fourier transform on a series of discrete stiffness values obtained from the finite element analysis, a stiffness curve expression can be obtained, and the stiffness is transformed into the following form:(1)kht=km+∑i=1NAicosiΩht+BisiniΩht,where km is the average value of the time-varying meshing stiffness, Ai, Bi is the ith-order harmonic amplitude of the time-varying meshing stiffness, Ωh is the excitation frequency, and N represents the order of the Fourier series, and the series was truncated at N=7 in this article.Figure4 shows a graph depicting the comprehensive meshing stiffness of the spiral bevel gear transmission during the two meshing cycles.Figure 4
Curve of the comprehensive meshing stiffness of the spiral bevel gear transmission for two meshing cycles.


### 3.2. Gear Transmission System

The gear transmission system includes the flexible shaft unit, spiral bevel gear meshing unit, and the bearing unit, which are analysed below.


#### 3.2.1. Flexible Shaft Unit

The gear shaft’s finite element model is established by the Timoshenko beam element, taking the shear deformation, gyro moment, and moment of inertia into account, as shown in Figure5, which gives the local coordinate of the beam element. In this article, each node has six degrees-of-freedom, i.e., the movement along the X, Y, and Z directions and the rotation angle around the X, Y, and Z directions. The element node’s displacement vector includes 12 degrees-of-freedom, which can be expressed as follows:(2)Xe=xiyiziθxiθyiθzixi+1yi+1zi+1θxi+1θyi+1θzi+1T.Figure 5
Local coordinate of the beam element.According to the literature [29, 30], the mass matrix Me, the stiffness matrix Ke, and the gyro matrix Ge of the beam elements can be obtained.It is noted that the directions of the local coordinate systems of the beam elements are various in the nonparallel shaft gear system. The beam element in the local coordinate needs to be transformed by the spatial coordinate to the global coordinate system. The relationship of the displacement of theith beam element node in the local coordinate system and its movement in the global coordinate system can be written as follows:(3)xi′,yi′,zi′,θxi′,θyi′,θzi′,xi+1′,yi+1′,zi+1′,θxi+1′,θyi+1′,θzi+1′T=λxi,yi,zi,θxi,θyi,θzi,xi+1,yi+1,zi+1,θxi+1,θyi+1,θzi+1T,where λ=diagT;T;T;T and T=TγTβTα, in which α, β, and γ are the rotation angles of the counterclockwise coordinate axis surrounding the global coordinate system X,Y,andZ. Tγ,Tβ,andTα are projection vectors along X, Y, and Z, respectively, and can be given by the following equations:(4)Tα=1000cosαsinα0−sinαcosα,Tβ=cosβ0sinβ010−sinβ0cosβ,Tγ=cosγsinγ0−sinγcosγ0001.Substituting the coordinate transformation matrix (4) into formula (5), the mass matrix, the stiffness matrix, and the gyro matrix of the ith beam element node in the global coordinate system are calculated, which can be described as follows:(5)Me′=RTMeR,Ke′=RTKeR,Ge′=RTGeR.The matrices of each element of the gear shaft, namely,Me′, Ke′, and Ge′, in the global coordinate system are assembled to obtain the mass matrix of the overall gear shaft Mr, the stiffness matrix of the overall gear shaft Kr, and the gyro matrix of the overall gear shaft Gr, respectively.The damping of the gear shaft is calculated via Rayleigh damping, i.e.,Cr=αMr+βKr, in which α and β are the mass damping coefficient and stiffness damping coefficient, respectively, which are calculated from the following equations:(6)α=2ω1ω2ξ1ω2−ξ2ω1ω22−ω12,β=2ξ2ω2−ξ1ω1ω22−ω12,where ξ1 and ξ2 denote the first and second mode damping ratios, respectively, generally supposing that ξ1=0.005 and ξ2=0.01. ω1 and ω2 are the first- and second-order natural frequencies of the rotor system, respectively, which can be obtained through modal experiments.When each element matrix of the gear shaft is acquired, the motion differential equation of a single gear shaft can also be obtained:(7)MrX¨+Cr+ωGrX˙+KrX=0.


#### 3.2.2. Meshing Unit

The schematic diagram of the dynamical model of the nonorthogonal spiral bevel gear is shown in Figure6. In this model, the local coordinate systems of the pinion and gear are fixed at the center of the two gears, namely, Op−XpYpZp and Og−XgYgZg, respectively. The direction of the local coordinate system of the driving gear is consistent with that of the global coordinate system. The dynamic model of the meshing unit is built based on the lumped mass method, and the lumped mass and the concentrated moment of inertia are utilised to simulate two spiral bevel gears. The model simulates the coupling effect during the meshing process through a time-varying meshing stiffness Kht, a meshing damping cm, and a transmission error et, without taking the backlash clearance, meshing friction force, and gear meshing position into consideration. The bending-torsion-shaft-pendulum coupling vibration of the two gears is considered so that the 12-degrees-of-freedom spiral bevel gear meshing unit can be obtained, namely, xp,yp,zp,θpx,θpy,θpz,xq,yq,zq,θqx,θqy,θqz.Figure 6
Schematic diagram of the dynamical model of the nonorthogonal spiral bevel gear.The spiral bevel gear is a space meshing gear, and the relative displacement of the two gears along the meshing line direction is difficult to obtain directly. By rotating the local coordinate systems of the pinion and gear and performing multiple coordinate projections, the relative displacement of the meshing direction of the spiral bevel gear can be transformed into that of the virtual spur gear’s meshing direction. The transformation process of the relative displacement of the spiral bevel gear in the direction of the meshing line is shown in Figure7.Figure 7
Transformation process of the relative displacement of a spiral bevel gear in the direction of the meshing line. (a) Projection along the pitch angle direction. (b) Projection along the spiral angle direction. (c) Projection along the pressure angle direction.
(a)
(b)
(c)The directional cosine vector of the pinion and gearni i=p,q can be calculated with a series of coordinate transformation as follows:(8)ni=nix,niy,niz=100TTizTiyTix,i=p,g,where Tix, Tiy, and Tiz indicate the transformation matrix about the x, y, and z axes, respectively, and can be expressed as follows:(9)Tix=1000cosδisinδi0−sinδicosδi,Tiy=cos−βi0sin−βi010−sin−βi0cos−βi,Tiz=cospi−αsinpi−α0−sinpi−αcospi−α0001,where δi = δ1 and pi−δ2 i=p,q represent the angles of the pinion and gear rotating in the direction of the pitch angles, respectively. βi = β and −β i=p,q are the midspiral angles of the gears, and α denotes the pressure angle.The directional rotation radii of the meshing force, which can be written as follows:(10)λix=ni⋅uix×ri=riyniz−rizniy,λiy=ni⋅uiy×ri=riznix−rixniz,λiz=ni⋅uiz×ri=rixniy−riynix,where ri=rixriyrizT represents the equivalent position vector of the meshing point, (this article ignores the variety of the meshing point’s position). Referring to Figure 6, the equivalent position vectors of the corresponding meshing point of the pinion and the gear are rp=0rp0T and rg=0rg0T, respectively, in which rp and rg are the pitch radii of the pinion and gear at the midpoint of the tooth width, respectively. uxi=100T, uyi=010T, and uzi=001T are the unit vectors in the X, Y, and Z directions, respectively. In the above equation, i=p,g represents the pinion and the gear, respectively.Also, the relative displacement of the spiral bevel gear on the meshing line can be expressed as follows:(11)Xn=VXm−ent,where V=npx,npy,npz,λpx,λpy,λpz,−ngx,−ngy,−ngz,λgx,λgy,λgz is the projection vector for the displacement of the spiral bevel gear in all directions to the meshing line; Xm=xp,yp,zp,θpx,θpy,θpz,xg,yg,zg,θgx,θgy,θgzT represents the displacement vector of the spiral bevel gear; ent=em+ersinΩt+ϕ, in which em and er represent the average value and the amplitude of the transmission error, respectively; Ω indicates the gear meshing’s angular frequency; and ϕ represents the initial phase angle of the engagement.Combined with the meshing stiffness and relative displacement of the spiral bevel gear obtained, the meshing force of the spiral bevel gear pair can be acquired by the following equation:(12)Fn=khtXn+cmX˙n,where cm represents the meshing damping, which can be calculated by the following:(13)cm=2ξkm1meq,p+1meq,g,where ξ is the meshing damping ratio, generally taken as 0.03∼0.17; meq,i=Izi/ri2 i=p,g denotes the equivalent mass of the gear.According to Newton’s second law, the dynamic equation of the spiral bevel gear meshing unit can be written as follows:(14)mpx¨p=Fnnpx+mpe1Ω12cosΩ1t,mpy¨p=Fnnpy+mpe1Ω12sinΩ1t,mpz¨p=Fnnpz,Ipθ˙px+JpΩ1θ˙py=Fnλpx,Ipθ˙py−JpΩ1θ˙px=Fnλpy,Jpθ˙pz=Fnλpz,mgx¨g=Fnngx+mge2Ω22cosΩ2t,mgy¨g=Fnngy+mge2Ω22sinΩ2t,mgz¨g=Fnngz,Igθ˙gx+JgΩ2θ˙gy=Fnλgx,Igθ¨gy−JgΩ2θ˙gx=Fnλgy,Jgθ¨gz=Fnλgz,where mp and mg represent the mass of the pinion and gear, respectively; Ip and Ig are the diameter rotational inertia of the pinion and gear, respectively; Jp and Jg represent the polar moment of inertia of the pinion and gear, respectively; e1 and e2 denote the pinion and gear’s eccentricity, respectively; and Ω1 and Ω2 are the angular velocity of the pinion and gear, respectively.Substituting formula (14) into the equations in (15) and arranging the equations into a matrix form, the motion differential equations of the gear meshing unit can be obtained as follows:(15)MmX¨m+Cm+Ω1Gm1Ω2Gm2X˙m+KmXm=kmtent+cme˙ntBT,where e˙nt is the derivative of the comprehensive meshing error; Gm1 and Gm2 represent the gyro matrix of the pinion and the gear, respectively; and Mm is the mass matrix of the meshing unit, which can be expressed as follows:(16)Mm=diagmpx,mpy,mpz,Ipx,Ipy,Ipz,mqx,mqy,mqz,Iqx,Iqy,Iqz.The stiffness matrixKm and the damping matrix Cm of the meshing element in the global coordinate system can be calculated as follows:(17)Km=kmBTB,Cm=cmBTB,where B=npx,npy,npz,λpx,λpy,λpz,−npx,−npy,−npz,−λpx,−λpy,−λpzT.


#### 3.2.3. Bearing Unit

The spiral bevel gear transmission system is supported by four angular contact ball bearings. During actual operation, there is a coupling relationship between the displacements of the bearings, and the stiffness in each direction varies with time. The time-varying stiffness of the bearings is usually not considered, and only the main stiffness term is retained, when calculating the dynamics of the gear system. The stiffness matrix of the ball bearing can be expressed as follows:(18)Kb=diagkxx,kyy,kzz,kθx,kθy,0.The bearing damping takes the same form as the stiffness and can also be written as follows:(19)Cb=diagcxx,cyy,czz,cθx,cθy,0.According to Newton’s second law, the dynamic differential equation of the bearing can be given as follows:(20)MekX¨ek+CbX˙ek−X˙Hk+KbXek−XHk=0,where Mek is the mass matrix of the beam element nodes connected to the kth (k = 1, 2, 3, 4) bearing; Xek, X˙ek, and X¨ek are the displacement, velocity, and acceleration vectors of the beam element nodes connected to the kth (k = 1, 2, 3, 4) bearing, respectively; XHk and X˙Hk denote the displacement and velocity vectors of the equivalent nodes of the gearbox unit connected to the kth (k = 1, 2, 3, 4) bearing, respectively.


### 3.3. Gearbox System

In order to establish the dynamic model of the gearbox, the dynamic parameters of the gearbox need to be obtained, including the mass matrix, the stiffness matrix, and the damping matrix. The gearbox nodes are divided into internal nodes and external nodes when the substructure method is employed to extract the gearbox’s matrix information. The external nodes are acquired by coupling the nodes on the bearing hole’s surface to the hole centre, and the remaining finite element nodes are the internal nodes. The gearbox mass, stiffness, and damping matrix can be extracted by a substructure method that condenses the internal nodes into external nodes. After that, the mass, the stiffness, and the damping matrix of the converted main node of the gearbox are assembled with the mass, stiffness, and damping matrix of the spiral bevel gear transmission system. The substructure method can avoid a large-scale finite element matrix, and it can be combined with the dynamic model of the transmission system conveniently.Figure8 shows a finite element model diagram of the gearbox using the substructure method. The gearbox is made of an aluminium alloy with a density of 2700 kg/m3, an elastic modulus of 7.1 × 1010 Pa, and Poisson’s ratio of 0.33. The model includes 828873 nodes and 585756 elements. Taking the degrees-of-freedom of the translation in each direction and rotation around X and Y axes of the bearing hole centre nodes into account, meanwhile, the screw holes are imposed a fixed constraint at the bottom of the gearbox. Therefore, a 20 × 20 mass matrix MH and a 20 × 20 stiffness matrix ΚH are formed; Rayleigh damping is employed, i.e., CH=αHMH+βHKH, where αH and βH are the mass coefficient and the stiffness coefficient, respectively.Figure 8
Finite element model diagram of the gearbox using the substructure method.
(a)
(b)After polycondensation, the degrees-of-freedom of the gearbox is significantly reduced, and the dynamic equation can be expressed as follows:(21)MHx¨H+CHx˙Hk−x˙ek+KHxHk−xek=0,where MH=MH11MH12⋯MH1nMH21MH22⋯MH2n⋮⋮⋱⋮MHn1MHn2⋯MHnn is the equivalent mass matrix of the gearbox substructure, in which n is the number of the bearings; KH=KH11KH12⋯KH1nKH21KH22⋯KH2n⋮⋮⋱⋮KHn1KHn2⋯KHnn, which is the equivalent stiffness matrix of the gearbox substructure; CH=aHMH+bHKH, which is the equivalent damping matrix of the gearbox substructure; XH, X˙H, and X¨H are the displacement, velocity, and acceleration vectors of the gearbox’s main node, respectively.


### 3.4. Global Dynamics Equation of the System

After obtaining the dynamic equations of each unit, the dynamic equations of the whole system can be established. Combined with the differential equation of the gear shaft unit (7), the differential equation of the meshing element (15), the differential equation of the bearing element (20), and the motion equation of the gearbox element (21), the full coupling differential equations of the entire system can be written as follows:(22)Mx¨t+C+ΩGx˙t−e˙t+Ktxt−et=Qt,where M, C, G, and K are the generalized mass, damping, gyroscopic, and stiffness matrices, respectively. X, X˙, and X¨ are the generalized displacement, velocity, and acceleration vectors, respectively. Q is the excitation force vector of the system.


## 3.1. Calculation of the Meshing Stiffness of the Spiral Bevel Gear

The gear’s meshing stiffness refers to the synthetic effect of each pair of gear teeth involved in meshing during the meshing process. It is difficult to calculate the meshing stiffness of the gear teeth using the theoretical analysis method because the tooth surface of the spiral bevel gear is a space meshing gear. Meanwhile, the tooth thickness is variable. Therefore, this study utilises finite element software to perform dynamic analysis and obtain an accurate meshing stiffness.The comprehensive elastic deformationδn and contact force Fn are calculated when multiple pairs of gears mesh at the same time, and then the comprehensive meshing stiffness can be obtained directly. The comprehensive elastic deformation of the gear teeth includes the contact elastic deformation due to the local Hertz contact δH, the displacement of the contact position of the gear teeth caused by bending of the gear teeth δb, and the deformation of the bearings, shafts, and supporting structures on the contact positions δf. The deformation of the bearings, shafts, and supporting structures has a relatively small influence on the meshing stiffness. Therefore, this study only takes δH and δb into account, so the deformation of the gear at a certain moment can be expressed as follows:δn=∑i=12δHi+∑i=12δbi.First, the gear’s meshing period is divided inton segments, and the meshing stiffness at each time point Kii=1,2,3,…,n is calculated [27, 28]. Then, a Fourier transform is utilized to convert a series of discrete values Ki into a meshing stiffness’s curve. In this article, the main parameters of the nonorthogonal spiral bevel gear are shown in Table 1, and then the comprehensive meshing stiffness is analysed. In order to improve the efficiency of calculation, nine pairs of gear teeth are selected from the spiral bevel gear pair for the meshing analysis. Figure 3 shows the finite element analysis of the meshing stiffness of the spiral bevel gear, including the meshing of the gear pair and the contact analysis of gear pair.Table 1
Main design parameters of the pinion and gear.
Pinion
Gear
Young’s modulus (Pa)
E = 2.1e9
Poisson’s ratio
N = 0.3
Pressure angle (°)
22.5
Spiral angle (°)
35
Transverse module (mm)
m = 5.6
Width of teeth (mm)
B = 33
Shaft angle (°)
126
Damping ratio
ξ
=
0.05
Mean of the transmission error (m)
e
m
=
2
×
10
−
5
Amplitude of the transmission error (m)
e
r
=
3
×
10
−
5
Pitch radius (mm)
72.5
87.6
Number of teeth
N
1 = 29
N
2 = 35
Hand
Right
LeftFigure 3
Finite element analysis of the meshing stiffness of the spiral bevel gear pair. (a) The meshing of the gear pairs. (b) Contact analysis of the gear pairs.
(a)
(b)Each periodT is divided into ten segments on average, and two meshing periods are taken. Then, each Ki value at its corresponding time ti is calculated. By performing a Fourier transform on a series of discrete stiffness values obtained from the finite element analysis, a stiffness curve expression can be obtained, and the stiffness is transformed into the following form:(1)kht=km+∑i=1NAicosiΩht+BisiniΩht,where km is the average value of the time-varying meshing stiffness, Ai, Bi is the ith-order harmonic amplitude of the time-varying meshing stiffness, Ωh is the excitation frequency, and N represents the order of the Fourier series, and the series was truncated at N=7 in this article.Figure4 shows a graph depicting the comprehensive meshing stiffness of the spiral bevel gear transmission during the two meshing cycles.Figure 4
Curve of the comprehensive meshing stiffness of the spiral bevel gear transmission for two meshing cycles.


## 3.2. Gear Transmission System

The gear transmission system includes the flexible shaft unit, spiral bevel gear meshing unit, and the bearing unit, which are analysed below.


### 3.2.1. Flexible Shaft Unit

The gear shaft’s finite element model is established by the Timoshenko beam element, taking the shear deformation, gyro moment, and moment of inertia into account, as shown in Figure5, which gives the local coordinate of the beam element. In this article, each node has six degrees-of-freedom, i.e., the movement along the X, Y, and Z directions and the rotation angle around the X, Y, and Z directions. The element node’s displacement vector includes 12 degrees-of-freedom, which can be expressed as follows:(2)Xe=xiyiziθxiθyiθzixi+1yi+1zi+1θxi+1θyi+1θzi+1T.Figure 5
Local coordinate of the beam element.According to the literature [29, 30], the mass matrix Me, the stiffness matrix Ke, and the gyro matrix Ge of the beam elements can be obtained.It is noted that the directions of the local coordinate systems of the beam elements are various in the nonparallel shaft gear system. The beam element in the local coordinate needs to be transformed by the spatial coordinate to the global coordinate system. The relationship of the displacement of theith beam element node in the local coordinate system and its movement in the global coordinate system can be written as follows:(3)xi′,yi′,zi′,θxi′,θyi′,θzi′,xi+1′,yi+1′,zi+1′,θxi+1′,θyi+1′,θzi+1′T=λxi,yi,zi,θxi,θyi,θzi,xi+1,yi+1,zi+1,θxi+1,θyi+1,θzi+1T,where λ=diagT;T;T;T and T=TγTβTα, in which α, β, and γ are the rotation angles of the counterclockwise coordinate axis surrounding the global coordinate system X,Y,andZ. Tγ,Tβ,andTα are projection vectors along X, Y, and Z, respectively, and can be given by the following equations:(4)Tα=1000cosαsinα0−sinαcosα,Tβ=cosβ0sinβ010−sinβ0cosβ,Tγ=cosγsinγ0−sinγcosγ0001.Substituting the coordinate transformation matrix (4) into formula (5), the mass matrix, the stiffness matrix, and the gyro matrix of the ith beam element node in the global coordinate system are calculated, which can be described as follows:(5)Me′=RTMeR,Ke′=RTKeR,Ge′=RTGeR.The matrices of each element of the gear shaft, namely,Me′, Ke′, and Ge′, in the global coordinate system are assembled to obtain the mass matrix of the overall gear shaft Mr, the stiffness matrix of the overall gear shaft Kr, and the gyro matrix of the overall gear shaft Gr, respectively.The damping of the gear shaft is calculated via Rayleigh damping, i.e.,Cr=αMr+βKr, in which α and β are the mass damping coefficient and stiffness damping coefficient, respectively, which are calculated from the following equations:(6)α=2ω1ω2ξ1ω2−ξ2ω1ω22−ω12,β=2ξ2ω2−ξ1ω1ω22−ω12,where ξ1 and ξ2 denote the first and second mode damping ratios, respectively, generally supposing that ξ1=0.005 and ξ2=0.01. ω1 and ω2 are the first- and second-order natural frequencies of the rotor system, respectively, which can be obtained through modal experiments.When each element matrix of the gear shaft is acquired, the motion differential equation of a single gear shaft can also be obtained:(7)MrX¨+Cr+ωGrX˙+KrX=0.


### 3.2.2. Meshing Unit

The schematic diagram of the dynamical model of the nonorthogonal spiral bevel gear is shown in Figure6. In this model, the local coordinate systems of the pinion and gear are fixed at the center of the two gears, namely, Op−XpYpZp and Og−XgYgZg, respectively. The direction of the local coordinate system of the driving gear is consistent with that of the global coordinate system. The dynamic model of the meshing unit is built based on the lumped mass method, and the lumped mass and the concentrated moment of inertia are utilised to simulate two spiral bevel gears. The model simulates the coupling effect during the meshing process through a time-varying meshing stiffness Kht, a meshing damping cm, and a transmission error et, without taking the backlash clearance, meshing friction force, and gear meshing position into consideration. The bending-torsion-shaft-pendulum coupling vibration of the two gears is considered so that the 12-degrees-of-freedom spiral bevel gear meshing unit can be obtained, namely, xp,yp,zp,θpx,θpy,θpz,xq,yq,zq,θqx,θqy,θqz.Figure 6
Schematic diagram of the dynamical model of the nonorthogonal spiral bevel gear.The spiral bevel gear is a space meshing gear, and the relative displacement of the two gears along the meshing line direction is difficult to obtain directly. By rotating the local coordinate systems of the pinion and gear and performing multiple coordinate projections, the relative displacement of the meshing direction of the spiral bevel gear can be transformed into that of the virtual spur gear’s meshing direction. The transformation process of the relative displacement of the spiral bevel gear in the direction of the meshing line is shown in Figure7.Figure 7
Transformation process of the relative displacement of a spiral bevel gear in the direction of the meshing line. (a) Projection along the pitch angle direction. (b) Projection along the spiral angle direction. (c) Projection along the pressure angle direction.
(a)
(b)
(c)The directional cosine vector of the pinion and gearni i=p,q can be calculated with a series of coordinate transformation as follows:(8)ni=nix,niy,niz=100TTizTiyTix,i=p,g,where Tix, Tiy, and Tiz indicate the transformation matrix about the x, y, and z axes, respectively, and can be expressed as follows:(9)Tix=1000cosδisinδi0−sinδicosδi,Tiy=cos−βi0sin−βi010−sin−βi0cos−βi,Tiz=cospi−αsinpi−α0−sinpi−αcospi−α0001,where δi = δ1 and pi−δ2 i=p,q represent the angles of the pinion and gear rotating in the direction of the pitch angles, respectively. βi = β and −β i=p,q are the midspiral angles of the gears, and α denotes the pressure angle.The directional rotation radii of the meshing force, which can be written as follows:(10)λix=ni⋅uix×ri=riyniz−rizniy,λiy=ni⋅uiy×ri=riznix−rixniz,λiz=ni⋅uiz×ri=rixniy−riynix,where ri=rixriyrizT represents the equivalent position vector of the meshing point, (this article ignores the variety of the meshing point’s position). Referring to Figure 6, the equivalent position vectors of the corresponding meshing point of the pinion and the gear are rp=0rp0T and rg=0rg0T, respectively, in which rp and rg are the pitch radii of the pinion and gear at the midpoint of the tooth width, respectively. uxi=100T, uyi=010T, and uzi=001T are the unit vectors in the X, Y, and Z directions, respectively. In the above equation, i=p,g represents the pinion and the gear, respectively.Also, the relative displacement of the spiral bevel gear on the meshing line can be expressed as follows:(11)Xn=VXm−ent,where V=npx,npy,npz,λpx,λpy,λpz,−ngx,−ngy,−ngz,λgx,λgy,λgz is the projection vector for the displacement of the spiral bevel gear in all directions to the meshing line; Xm=xp,yp,zp,θpx,θpy,θpz,xg,yg,zg,θgx,θgy,θgzT represents the displacement vector of the spiral bevel gear; ent=em+ersinΩt+ϕ, in which em and er represent the average value and the amplitude of the transmission error, respectively; Ω indicates the gear meshing’s angular frequency; and ϕ represents the initial phase angle of the engagement.Combined with the meshing stiffness and relative displacement of the spiral bevel gear obtained, the meshing force of the spiral bevel gear pair can be acquired by the following equation:(12)Fn=khtXn+cmX˙n,where cm represents the meshing damping, which can be calculated by the following:(13)cm=2ξkm1meq,p+1meq,g,where ξ is the meshing damping ratio, generally taken as 0.03∼0.17; meq,i=Izi/ri2 i=p,g denotes the equivalent mass of the gear.According to Newton’s second law, the dynamic equation of the spiral bevel gear meshing unit can be written as follows:(14)mpx¨p=Fnnpx+mpe1Ω12cosΩ1t,mpy¨p=Fnnpy+mpe1Ω12sinΩ1t,mpz¨p=Fnnpz,Ipθ˙px+JpΩ1θ˙py=Fnλpx,Ipθ˙py−JpΩ1θ˙px=Fnλpy,Jpθ˙pz=Fnλpz,mgx¨g=Fnngx+mge2Ω22cosΩ2t,mgy¨g=Fnngy+mge2Ω22sinΩ2t,mgz¨g=Fnngz,Igθ˙gx+JgΩ2θ˙gy=Fnλgx,Igθ¨gy−JgΩ2θ˙gx=Fnλgy,Jgθ¨gz=Fnλgz,where mp and mg represent the mass of the pinion and gear, respectively; Ip and Ig are the diameter rotational inertia of the pinion and gear, respectively; Jp and Jg represent the polar moment of inertia of the pinion and gear, respectively; e1 and e2 denote the pinion and gear’s eccentricity, respectively; and Ω1 and Ω2 are the angular velocity of the pinion and gear, respectively.Substituting formula (14) into the equations in (15) and arranging the equations into a matrix form, the motion differential equations of the gear meshing unit can be obtained as follows:(15)MmX¨m+Cm+Ω1Gm1Ω2Gm2X˙m+KmXm=kmtent+cme˙ntBT,where e˙nt is the derivative of the comprehensive meshing error; Gm1 and Gm2 represent the gyro matrix of the pinion and the gear, respectively; and Mm is the mass matrix of the meshing unit, which can be expressed as follows:(16)Mm=diagmpx,mpy,mpz,Ipx,Ipy,Ipz,mqx,mqy,mqz,Iqx,Iqy,Iqz.The stiffness matrixKm and the damping matrix Cm of the meshing element in the global coordinate system can be calculated as follows:(17)Km=kmBTB,Cm=cmBTB,where B=npx,npy,npz,λpx,λpy,λpz,−npx,−npy,−npz,−λpx,−λpy,−λpzT.


### 3.2.3. Bearing Unit

The spiral bevel gear transmission system is supported by four angular contact ball bearings. During actual operation, there is a coupling relationship between the displacements of the bearings, and the stiffness in each direction varies with time. The time-varying stiffness of the bearings is usually not considered, and only the main stiffness term is retained, when calculating the dynamics of the gear system. The stiffness matrix of the ball bearing can be expressed as follows:(18)Kb=diagkxx,kyy,kzz,kθx,kθy,0.The bearing damping takes the same form as the stiffness and can also be written as follows:(19)Cb=diagcxx,cyy,czz,cθx,cθy,0.According to Newton’s second law, the dynamic differential equation of the bearing can be given as follows:(20)MekX¨ek+CbX˙ek−X˙Hk+KbXek−XHk=0,where Mek is the mass matrix of the beam element nodes connected to the kth (k = 1, 2, 3, 4) bearing; Xek, X˙ek, and X¨ek are the displacement, velocity, and acceleration vectors of the beam element nodes connected to the kth (k = 1, 2, 3, 4) bearing, respectively; XHk and X˙Hk denote the displacement and velocity vectors of the equivalent nodes of the gearbox unit connected to the kth (k = 1, 2, 3, 4) bearing, respectively.


## 3.2.1. Flexible Shaft Unit

The gear shaft’s finite element model is established by the Timoshenko beam element, taking the shear deformation, gyro moment, and moment of inertia into account, as shown in Figure5, which gives the local coordinate of the beam element. In this article, each node has six degrees-of-freedom, i.e., the movement along the X, Y, and Z directions and the rotation angle around the X, Y, and Z directions. The element node’s displacement vector includes 12 degrees-of-freedom, which can be expressed as follows:(2)Xe=xiyiziθxiθyiθzixi+1yi+1zi+1θxi+1θyi+1θzi+1T.Figure 5
Local coordinate of the beam element.According to the literature [29, 30], the mass matrix Me, the stiffness matrix Ke, and the gyro matrix Ge of the beam elements can be obtained.It is noted that the directions of the local coordinate systems of the beam elements are various in the nonparallel shaft gear system. The beam element in the local coordinate needs to be transformed by the spatial coordinate to the global coordinate system. The relationship of the displacement of theith beam element node in the local coordinate system and its movement in the global coordinate system can be written as follows:(3)xi′,yi′,zi′,θxi′,θyi′,θzi′,xi+1′,yi+1′,zi+1′,θxi+1′,θyi+1′,θzi+1′T=λxi,yi,zi,θxi,θyi,θzi,xi+1,yi+1,zi+1,θxi+1,θyi+1,θzi+1T,where λ=diagT;T;T;T and T=TγTβTα, in which α, β, and γ are the rotation angles of the counterclockwise coordinate axis surrounding the global coordinate system X,Y,andZ. Tγ,Tβ,andTα are projection vectors along X, Y, and Z, respectively, and can be given by the following equations:(4)Tα=1000cosαsinα0−sinαcosα,Tβ=cosβ0sinβ010−sinβ0cosβ,Tγ=cosγsinγ0−sinγcosγ0001.Substituting the coordinate transformation matrix (4) into formula (5), the mass matrix, the stiffness matrix, and the gyro matrix of the ith beam element node in the global coordinate system are calculated, which can be described as follows:(5)Me′=RTMeR,Ke′=RTKeR,Ge′=RTGeR.The matrices of each element of the gear shaft, namely,Me′, Ke′, and Ge′, in the global coordinate system are assembled to obtain the mass matrix of the overall gear shaft Mr, the stiffness matrix of the overall gear shaft Kr, and the gyro matrix of the overall gear shaft Gr, respectively.The damping of the gear shaft is calculated via Rayleigh damping, i.e.,Cr=αMr+βKr, in which α and β are the mass damping coefficient and stiffness damping coefficient, respectively, which are calculated from the following equations:(6)α=2ω1ω2ξ1ω2−ξ2ω1ω22−ω12,β=2ξ2ω2−ξ1ω1ω22−ω12,where ξ1 and ξ2 denote the first and second mode damping ratios, respectively, generally supposing that ξ1=0.005 and ξ2=0.01. ω1 and ω2 are the first- and second-order natural frequencies of the rotor system, respectively, which can be obtained through modal experiments.When each element matrix of the gear shaft is acquired, the motion differential equation of a single gear shaft can also be obtained:(7)MrX¨+Cr+ωGrX˙+KrX=0.


## 3.2.2. Meshing Unit

The schematic diagram of the dynamical model of the nonorthogonal spiral bevel gear is shown in Figure6. In this model, the local coordinate systems of the pinion and gear are fixed at the center of the two gears, namely, Op−XpYpZp and Og−XgYgZg, respectively. The direction of the local coordinate system of the driving gear is consistent with that of the global coordinate system. The dynamic model of the meshing unit is built based on the lumped mass method, and the lumped mass and the concentrated moment of inertia are utilised to simulate two spiral bevel gears. The model simulates the coupling effect during the meshing process through a time-varying meshing stiffness Kht, a meshing damping cm, and a transmission error et, without taking the backlash clearance, meshing friction force, and gear meshing position into consideration. The bending-torsion-shaft-pendulum coupling vibration of the two gears is considered so that the 12-degrees-of-freedom spiral bevel gear meshing unit can be obtained, namely, xp,yp,zp,θpx,θpy,θpz,xq,yq,zq,θqx,θqy,θqz.Figure 6
Schematic diagram of the dynamical model of the nonorthogonal spiral bevel gear.The spiral bevel gear is a space meshing gear, and the relative displacement of the two gears along the meshing line direction is difficult to obtain directly. By rotating the local coordinate systems of the pinion and gear and performing multiple coordinate projections, the relative displacement of the meshing direction of the spiral bevel gear can be transformed into that of the virtual spur gear’s meshing direction. The transformation process of the relative displacement of the spiral bevel gear in the direction of the meshing line is shown in Figure7.Figure 7
Transformation process of the relative displacement of a spiral bevel gear in the direction of the meshing line. (a) Projection along the pitch angle direction. (b) Projection along the spiral angle direction. (c) Projection along the pressure angle direction.
(a)
(b)
(c)The directional cosine vector of the pinion and gearni i=p,q can be calculated with a series of coordinate transformation as follows:(8)ni=nix,niy,niz=100TTizTiyTix,i=p,g,where Tix, Tiy, and Tiz indicate the transformation matrix about the x, y, and z axes, respectively, and can be expressed as follows:(9)Tix=1000cosδisinδi0−sinδicosδi,Tiy=cos−βi0sin−βi010−sin−βi0cos−βi,Tiz=cospi−αsinpi−α0−sinpi−αcospi−α0001,where δi = δ1 and pi−δ2 i=p,q represent the angles of the pinion and gear rotating in the direction of the pitch angles, respectively. βi = β and −β i=p,q are the midspiral angles of the gears, and α denotes the pressure angle.The directional rotation radii of the meshing force, which can be written as follows:(10)λix=ni⋅uix×ri=riyniz−rizniy,λiy=ni⋅uiy×ri=riznix−rixniz,λiz=ni⋅uiz×ri=rixniy−riynix,where ri=rixriyrizT represents the equivalent position vector of the meshing point, (this article ignores the variety of the meshing point’s position). Referring to Figure 6, the equivalent position vectors of the corresponding meshing point of the pinion and the gear are rp=0rp0T and rg=0rg0T, respectively, in which rp and rg are the pitch radii of the pinion and gear at the midpoint of the tooth width, respectively. uxi=100T, uyi=010T, and uzi=001T are the unit vectors in the X, Y, and Z directions, respectively. In the above equation, i=p,g represents the pinion and the gear, respectively.Also, the relative displacement of the spiral bevel gear on the meshing line can be expressed as follows:(11)Xn=VXm−ent,where V=npx,npy,npz,λpx,λpy,λpz,−ngx,−ngy,−ngz,λgx,λgy,λgz is the projection vector for the displacement of the spiral bevel gear in all directions to the meshing line; Xm=xp,yp,zp,θpx,θpy,θpz,xg,yg,zg,θgx,θgy,θgzT represents the displacement vector of the spiral bevel gear; ent=em+ersinΩt+ϕ, in which em and er represent the average value and the amplitude of the transmission error, respectively; Ω indicates the gear meshing’s angular frequency; and ϕ represents the initial phase angle of the engagement.Combined with the meshing stiffness and relative displacement of the spiral bevel gear obtained, the meshing force of the spiral bevel gear pair can be acquired by the following equation:(12)Fn=khtXn+cmX˙n,where cm represents the meshing damping, which can be calculated by the following:(13)cm=2ξkm1meq,p+1meq,g,where ξ is the meshing damping ratio, generally taken as 0.03∼0.17; meq,i=Izi/ri2 i=p,g denotes the equivalent mass of the gear.According to Newton’s second law, the dynamic equation of the spiral bevel gear meshing unit can be written as follows:(14)mpx¨p=Fnnpx+mpe1Ω12cosΩ1t,mpy¨p=Fnnpy+mpe1Ω12sinΩ1t,mpz¨p=Fnnpz,Ipθ˙px+JpΩ1θ˙py=Fnλpx,Ipθ˙py−JpΩ1θ˙px=Fnλpy,Jpθ˙pz=Fnλpz,mgx¨g=Fnngx+mge2Ω22cosΩ2t,mgy¨g=Fnngy+mge2Ω22sinΩ2t,mgz¨g=Fnngz,Igθ˙gx+JgΩ2θ˙gy=Fnλgx,Igθ¨gy−JgΩ2θ˙gx=Fnλgy,Jgθ¨gz=Fnλgz,where mp and mg represent the mass of the pinion and gear, respectively; Ip and Ig are the diameter rotational inertia of the pinion and gear, respectively; Jp and Jg represent the polar moment of inertia of the pinion and gear, respectively; e1 and e2 denote the pinion and gear’s eccentricity, respectively; and Ω1 and Ω2 are the angular velocity of the pinion and gear, respectively.Substituting formula (14) into the equations in (15) and arranging the equations into a matrix form, the motion differential equations of the gear meshing unit can be obtained as follows:(15)MmX¨m+Cm+Ω1Gm1Ω2Gm2X˙m+KmXm=kmtent+cme˙ntBT,where e˙nt is the derivative of the comprehensive meshing error; Gm1 and Gm2 represent the gyro matrix of the pinion and the gear, respectively; and Mm is the mass matrix of the meshing unit, which can be expressed as follows:(16)Mm=diagmpx,mpy,mpz,Ipx,Ipy,Ipz,mqx,mqy,mqz,Iqx,Iqy,Iqz.The stiffness matrixKm and the damping matrix Cm of the meshing element in the global coordinate system can be calculated as follows:(17)Km=kmBTB,Cm=cmBTB,where B=npx,npy,npz,λpx,λpy,λpz,−npx,−npy,−npz,−λpx,−λpy,−λpzT.


## 3.2.3. Bearing Unit

The spiral bevel gear transmission system is supported by four angular contact ball bearings. During actual operation, there is a coupling relationship between the displacements of the bearings, and the stiffness in each direction varies with time. The time-varying stiffness of the bearings is usually not considered, and only the main stiffness term is retained, when calculating the dynamics of the gear system. The stiffness matrix of the ball bearing can be expressed as follows:(18)Kb=diagkxx,kyy,kzz,kθx,kθy,0.The bearing damping takes the same form as the stiffness and can also be written as follows:(19)Cb=diagcxx,cyy,czz,cθx,cθy,0.According to Newton’s second law, the dynamic differential equation of the bearing can be given as follows:(20)MekX¨ek+CbX˙ek−X˙Hk+KbXek−XHk=0,where Mek is the mass matrix of the beam element nodes connected to the kth (k = 1, 2, 3, 4) bearing; Xek, X˙ek, and X¨ek are the displacement, velocity, and acceleration vectors of the beam element nodes connected to the kth (k = 1, 2, 3, 4) bearing, respectively; XHk and X˙Hk denote the displacement and velocity vectors of the equivalent nodes of the gearbox unit connected to the kth (k = 1, 2, 3, 4) bearing, respectively.


## 3.3. Gearbox System

In order to establish the dynamic model of the gearbox, the dynamic parameters of the gearbox need to be obtained, including the mass matrix, the stiffness matrix, and the damping matrix. The gearbox nodes are divided into internal nodes and external nodes when the substructure method is employed to extract the gearbox’s matrix information. The external nodes are acquired by coupling the nodes on the bearing hole’s surface to the hole centre, and the remaining finite element nodes are the internal nodes. The gearbox mass, stiffness, and damping matrix can be extracted by a substructure method that condenses the internal nodes into external nodes. After that, the mass, the stiffness, and the damping matrix of the converted main node of the gearbox are assembled with the mass, stiffness, and damping matrix of the spiral bevel gear transmission system. The substructure method can avoid a large-scale finite element matrix, and it can be combined with the dynamic model of the transmission system conveniently.Figure8 shows a finite element model diagram of the gearbox using the substructure method. The gearbox is made of an aluminium alloy with a density of 2700 kg/m3, an elastic modulus of 7.1 × 1010 Pa, and Poisson’s ratio of 0.33. The model includes 828873 nodes and 585756 elements. Taking the degrees-of-freedom of the translation in each direction and rotation around X and Y axes of the bearing hole centre nodes into account, meanwhile, the screw holes are imposed a fixed constraint at the bottom of the gearbox. Therefore, a 20 × 20 mass matrix MH and a 20 × 20 stiffness matrix ΚH are formed; Rayleigh damping is employed, i.e., CH=αHMH+βHKH, where αH and βH are the mass coefficient and the stiffness coefficient, respectively.Figure 8
Finite element model diagram of the gearbox using the substructure method.
(a)
(b)After polycondensation, the degrees-of-freedom of the gearbox is significantly reduced, and the dynamic equation can be expressed as follows:(21)MHx¨H+CHx˙Hk−x˙ek+KHxHk−xek=0,where MH=MH11MH12⋯MH1nMH21MH22⋯MH2n⋮⋮⋱⋮MHn1MHn2⋯MHnn is the equivalent mass matrix of the gearbox substructure, in which n is the number of the bearings; KH=KH11KH12⋯KH1nKH21KH22⋯KH2n⋮⋮⋱⋮KHn1KHn2⋯KHnn, which is the equivalent stiffness matrix of the gearbox substructure; CH=aHMH+bHKH, which is the equivalent damping matrix of the gearbox substructure; XH, X˙H, and X¨H are the displacement, velocity, and acceleration vectors of the gearbox’s main node, respectively.


## 3.4. Global Dynamics Equation of the System

After obtaining the dynamic equations of each unit, the dynamic equations of the whole system can be established. Combined with the differential equation of the gear shaft unit (7), the differential equation of the meshing element (15), the differential equation of the bearing element (20), and the motion equation of the gearbox element (21), the full coupling differential equations of the entire system can be written as follows:(22)Mx¨t+C+ΩGx˙t−e˙t+Ktxt−et=Qt,where M, C, G, and K are the generalized mass, damping, gyroscopic, and stiffness matrices, respectively. X, X˙, and X¨ are the generalized displacement, velocity, and acceleration vectors, respectively. Q is the excitation force vector of the system.


## 4. Solution and Analysis of the Dynamic Equation of the Full Coupling System of the Spiral Bevel Gear

From the derivation process and results, it can be determined that the full coupling system is a complex coupling dynamics question, and it is complicated to obtain satisfactory results using the theoretical analysis method. Meanwhile, the stiffness matrix of the entire dynamic equation of the system is a singular matrix, i.e., the system has rigid body displacement. In this article, the rigid body displacement is eliminated by imposing a fixed constraint on the torsional direction of the output end. Although the absolute displacement of each degrees-of-freedom will deviate from the actual one, the relative microvibration displacement has little effect. Then, the Newmark method is employed to calculate the numerical solution of the full coupling system. In order to obtain an accurate dynamic response, the step size is 1/200 of the gear meshing frequency. The structural parameters of the rotating shaft of the full coupling system are shown in Table2. The relevant parameters of the spiral bevel gear are revealed in Table 1, in which the time-varying meshing stiffness of the gear pair is obtained from Part 3.1; the stiffness and damping of the ball bearings are shown in Table 3, and the structural parameters of the gearbox are acquired from Section 3.3.Table 2
Shaft element parameters.
Driving and driven shaft
Node to node
1-2
2-3
3-4
4–9
10–15
15-16
16-17
17-18
Outer diameter (mm)
48
57
60
75
75
60
57
48
Inner diameter (mm)
44
45
45
52
52
45
45
44
Shaft length (mm)
20
26
28
108
140
28
26
20Table 3
Bearing parameters.
Stiffness
Bearing 1-2
Bearing 3-4
k
xx, kyy (N/m)
2 × 108
2 × 108
k
zz (N/m)
1 × 108
1 × 108
k
θxθx, kθyθy (Nm/rad)
1 × 106
1 × 106
Damping
Bearing 1-2
Bearing 3-4
c (Ns/m)
1 × 103
1 × 103


### 4.1. Effect of the Gearbox Flexibility on the Dynamic Characteristics of the Spiral Bevel Gear Transmission System

#### 4.1.1. Influence of the Gearbox Flexibility on the Meshing Force

In order to investigate the gearbox flexibility’s effect on the meshing force of the gear transmission system, the meshing force of the uncoupling model (only considering the spiral bevel gear transmission system) and the coupling model (the spiral bevel gear-gearbox coupling model considering the gearbox flexibility) are obtained. Figure9 shows the amplitude-frequency response curve of the meshing force of the uncoupling and coupling mode at a rotation speed range of 200∼8000 r/min. As seen in Figure 9, the stiffness of the gearbox is introduced into the spiral bevel gear transmission system, which reduces the overall stiffness of the system, and absorbs and transmits part of the energy, resulting in a reduction in the meshing force at all rotational speeds. However, the magnitude of reduction in the meshing force is different at various rotational speeds. By comparing Figure 9(a) with Figure 9(b), it can be observed that, with the variety of the torque, the effect of the gearbox flexibility on meshing force is also varied. Under torques of 300 Nm and 600 Nm, the meshing forces of the uncoupling model and the coupling model are calculated, respectively, at rotating speeds of 3000 r/min, 4700 r/min, and 7100 r/min. The specific values are shown in Table 4. By analysing the data in Table 4, it is demonstrated that the effect of the gearbox flexibility on the meshing force is also larger when the torque is larger at the same speed.Figure 9
Amplitude-frequency response curve of the meshing forces of the uncoupling and coupling model under different torques. (a)T = 300 Nm. (b) T = 600 Nm.
(a)
(b)Table 4
Change rate of the meshing force under different models.
Torque (Nm)
Speed (r/min)
Meshing force of coupling model (N)
Meshing force of uncoupled model (N)
Rate of change (%)
300
3000
46570
49240
5.4
4700
66190
69000
4.1
7100
42130
44920
6.2
600
3000
71170
76010
6.4
4700
95210
100100
4.8
7100
70070
75880
7.6


#### 4.1.2. Influence of the Gearbox Flexibility on the Bearing Force

Similarly, to study the influence of the gearbox flexibility on the bearing force, the differential equations of the system are solved in order to obtain the amplitude-frequency response curves of the bearing force of the uncoupling model and coupling model in the speed range of 200∼8000 r/min, as shown in Figure10(a). It can be observed that with the variety of the rotating speed, the bearing force decreases in the entire rotating speed range while taking the gearbox stiffness into account, and the most significant decrease of the bearing force occurs at the resonance point. The main reason for this is that the vibration of the gearbox is introduced into the entire dynamic system, and the variety of the system stiffness affects the supporting force significantly, and this also means that the force transmitted by the transmission system to the other structural systems is reduced considerably. The gearbox flexibility in the transmission system helps obtain more accurate dynamic analysis results, which provides a theoretical basis for the optimization design and vibration and noise reduction of the system. Then, the amplitude-frequency response curve of the bearing force of the uncoupling and coupling gearbox in the rotating speed range of 200∼8000 r/min is obtained when the torque increases, which is revealed in Figure 10(b). To facilitate the analysis, under torques of 300 Nm and 600 Nm, the bearing forces of uncoupling and coupling the gearbox at rotating speeds of 3000 r/min, 4700 r/min, and 7100 r/min are calculated, respectively, as shown in Table 5. By observing the data in the table, it can be determined that the change rate of the bearing force is large when the torque is large at the same speed.Figure 10
Amplitude-frequency response curves of the bearing forces of the uncoupling and coupling model under different torques.(a)T = 300 Nm. (b) T = 600 Nm.
(a)
(b)Table 5
Change rate of the bearing force under different models.
Torque (Nm)
Speed (r/min)
Bearing force of coupling model (N)
Bearing force of uncoupling model (N)
Rate of change (%)
300
3000
2213
3307
33
4700
6393
9808
35
7100
1278
1893
32
600
3000
2741
4292
37
4700
7964
11700
40
7100
1545
2431
36According to the above analysis, the flexibility of the gearbox has an impact on the meshing force and the bearing force, of which the influence on the bearing force is prominent. In the analysis of the overall dynamic response of the coupling system supported by the thin-walled gearbox, the effect of the gearbox flexibility needs to be taken into account.


### 4.2. Influence of the Shaft Angle on the Dynamic Characteristics of the Spiral Bevel Gear Full Coupling System

Based on the research in the previous section, the influence of the flexibility of the thin-walled gearbox on the spiral bevel gear transmission system cannot be ignored. Thus, this section directly studies the effect of the shaft angle on the dynamic characteristics of the spiral bevel gear-shaft-bearing-gearbox full coupling system.


#### 4.2.1. Effect of the Shaft Angle on the Axis Trajectory of the Gear

Studying the axis locus of the critical nodes can predict whether the geared rotor transmission system will encounter rub-impact during operation or not, and the investigation of the gear’s axis locus is especially helpful in understanding the meshing status of the gear. It is assumed that the speed of the driving spiral bevel gear is 4200 r/min, the torque is 600 Nm, and the other parameters keep invariant. The variation of the axial trajectory of the driven spiral bevel gear is analysed when the shaft angles of the spiral bevel gear transmission system are 70, 90, 110, and 126, as shown in Figure11. The centre of the trajectory of the driven spiral bevel gear is varied when the shaft angle has diverse values, as shown in Figure 11(a). The centres of the axis locus under each shaft angle are set to the same position to observe the variety of the direction and size of the axis locus. It can be seen in Figure 11(b) that the variation in the shaft angle not only affects the direction of the axis locus but also alters the amplitude of the axis locus. The calculation shows that the axis locus of the other nodes will also change with the variety of the shaft angle. Therefore, the impact of the shaft angle must be taken into account when studying the vibration features of the nonorthogonal spiral bevel gear-shaft-bearing-gearbox full coupling system.Figure 11
Axis trajectory diagram of a driven gear under various shaft angles.
(a)
(b)


#### 4.2.2. Effect of the Shaft Angle on the Meshing Force and Bearing Force

Supposing that the driving gear’s rotation speed is 4200 r/min, the torque is 600 Nm, and other parameters kept invariant, then the meshing force and the bearing force variation in each direction are obtained when the shaft angle is 70, 90, 110, and 126, respectively, as shown in Figure12. It can be demonstrated from the figure that the meshing force and the bearing force of the spiral bevel gear pair are various with the variety of the shaft angles. The meshing force varies because the pitch angles of the pinion and the gear are various under different shaft angles, and the projection results of the vibration displacement of the gear pair in various directions to the meshing line vary. At the same time, the force’s components decomposed into diverse directions are also different, resulting in variations in the bearing force for each direction. In conclusion, the vibration behaviours of the spiral bevel gear pair and the full coupling system will vary with the variety of the shaft angle.Figure 12
Meshing force and bearing force under various shaft angles. (a) Meshing force. (b) Bearing force in theX direction. (c) Bearing force in the Y direction. (d) Bearing force in the Z direction.
(a)
(b)
(c)
(d)


### 4.3. Influence of Unbalance on the Dynamic Characteristics of a Spiral Bevel Gear Full Coupling System

#### 4.3.1. The Influence of the Unbalance on Frequency Response Characteristics of the Coupling System

An unbalance amount, which is equal to5×10−4kg⋅m, is applied to the drive shaft’s node 4, with no imbalance given to the other nodes. The driving shaft speed n1 is in the range of 200∼800 r/min, and the three-dimensional spectrum of the full coupling system at the driving and driven gears are obtained, respectively, as shown in Figure 13. As seen from Figure 13(a), because of the coupling effect between the rotor and the gear, the X-direction displacement vibration spectrum of the driving gear includes the rotation frequency of the driving shaft f1, the gear meshing frequency fm, and its harmonic frequency 2fm. Meanwhile, it can be seen from Figure 13(b) that the torsional vibration spectrum of the driving gear also includes the rotation frequency of the driving shaft f1, the meshing frequency fm, and its harmonic frequency 2fm. As seen from Figure 13(c), the X-direction displacement vibration spectrum of the driven gear includes the meshing frequency fm and its harmonic frequency 2fm. The imbalance of the driving shaft has little effect on the vibration displacement of the driven gear in the X direction, where the rotation frequency amplitude of the driving shaft is small, which is not revealed in the figure. As seen from Figure 13(d), because of the coupling effect of bending and torsion of the spiral bevel gear pair, the driven gear’s torsional vibration not only contains the meshing frequency fm and its harmonic frequency 2fm but also the driving shaft’s rotation frequency f1, where the rotational frequency amplitude increases with the rise of the speed. Thus, the effect of the imbalance of the rotating shaft on the gear’s vibration in the torsional direction is more prominent than that on the lateral vibration.Figure 13
Three-dimensional spectrum of the coupling system as the unbalance is5×10−4kg⋅m. (a) Lateral vibration of the driving gear in the X direction. (b) Torsional vibration of the driving gear. (c) Lateral vibration of the driven gear in the X direction. (d) Torsional vibration of the driven gear.
(a)
(b)
(c)
(d)The unbalance amount applied on the driving shaft’s node 4 is changed to10×10−4kg⋅m, and the other conditions remain unchanged. The speed of the driving shaft n1 varies from 200 r/min to 8000 r/min, and the three-dimensional spectrum of the full coupling system at the pinion and gear are acquired, as shown in Figure 14. Comparing Figures 13(a) and 13(b) and 14(a) and 14(b), it can be demonstrated that the amplitude corresponding to the rotating frequency of the driving gear in the lateral and torsional vibration spectrograms increases significantly when the imbalance on the driving shaft is rising, while the amplitude corresponding to the meshing frequency hardly changes. It can be seen from Figure 14(c) that the driven gear’s vibration spectrum exhibits the driving shaft’s rotating frequency f1 at a high rotational speed, indicating that the coupling effect of the driving shaft’s imbalance on the driven gear’s vibration displacement is also risen with an increase in the unbalance of the driving shaft. Comparing Figures 13(d) and 14(d), it can be revealed that the amplitude of the rotation frequency f1 rises in the driven gear’s torsional vibration spectrum, which also indicates that the increase of the driving shaft’s unbalance increases the effect of the driven gear’s torsional vibration, and the increased amplitude is more obvious than that of the transverse vibration.Figure 14
Three-dimensional spectrum of the coupling system as the unbalance is10×10−4kg⋅m. (a) Lateral vibration of the driving gear in the X direction. (b) Torsional vibration of the driving gear. (c) Lateral vibration of the driven gear in the X direction. (d) Torsional vibration of the driven gear.
(a)
(b)
(c)
(d)


#### 4.3.2. Influence of the Unbalance on the Meshing Force and Bearing Force

In order to analyse the influence of the unbalance of the rotating shaft on the meshing force and the bearing force, the unbalance amount on the driving shaft’s node 4 is taken as 0,5×10−4kg⋅m, and 10×10−4kg⋅m, respectively. The amplitude-frequency response curves of the meshing force and the bearing force within a rotation speed range of 200∼8000 r/min are calculated, as shown in Figure 15. It can be revealed from the figure that the meshing force and the bearing force have little influence at low rotational speeds with the increase of the eccentricity; however, the eccentricity has a great influence on the meshing force and bearing force when speed continuously increase. Therefore, when the speed of the full coupling system is high, the eccentric amount of the rotating shaft should be reduced to diminish its impact on the meshing force and bearing force so as to ensure the stable and reliable operation of the system.Figure 15
Amplitude-frequency response curves of the meshing force and the bearing force. (a) Amplitude-frequency response curve of the bearing force. (b) Amplitude-frequency response curve of the meshing force.
(a)
(b)


## 4.1. Effect of the Gearbox Flexibility on the Dynamic Characteristics of the Spiral Bevel Gear Transmission System

### 4.1.1. Influence of the Gearbox Flexibility on the Meshing Force

In order to investigate the gearbox flexibility’s effect on the meshing force of the gear transmission system, the meshing force of the uncoupling model (only considering the spiral bevel gear transmission system) and the coupling model (the spiral bevel gear-gearbox coupling model considering the gearbox flexibility) are obtained. Figure9 shows the amplitude-frequency response curve of the meshing force of the uncoupling and coupling mode at a rotation speed range of 200∼8000 r/min. As seen in Figure 9, the stiffness of the gearbox is introduced into the spiral bevel gear transmission system, which reduces the overall stiffness of the system, and absorbs and transmits part of the energy, resulting in a reduction in the meshing force at all rotational speeds. However, the magnitude of reduction in the meshing force is different at various rotational speeds. By comparing Figure 9(a) with Figure 9(b), it can be observed that, with the variety of the torque, the effect of the gearbox flexibility on meshing force is also varied. Under torques of 300 Nm and 600 Nm, the meshing forces of the uncoupling model and the coupling model are calculated, respectively, at rotating speeds of 3000 r/min, 4700 r/min, and 7100 r/min. The specific values are shown in Table 4. By analysing the data in Table 4, it is demonstrated that the effect of the gearbox flexibility on the meshing force is also larger when the torque is larger at the same speed.Figure 9
Amplitude-frequency response curve of the meshing forces of the uncoupling and coupling model under different torques. (a)T = 300 Nm. (b) T = 600 Nm.
(a)
(b)Table 4
Change rate of the meshing force under different models.
Torque (Nm)
Speed (r/min)
Meshing force of coupling model (N)
Meshing force of uncoupled model (N)
Rate of change (%)
300
3000
46570
49240
5.4
4700
66190
69000
4.1
7100
42130
44920
6.2
600
3000
71170
76010
6.4
4700
95210
100100
4.8
7100
70070
75880
7.6


### 4.1.2. Influence of the Gearbox Flexibility on the Bearing Force

Similarly, to study the influence of the gearbox flexibility on the bearing force, the differential equations of the system are solved in order to obtain the amplitude-frequency response curves of the bearing force of the uncoupling model and coupling model in the speed range of 200∼8000 r/min, as shown in Figure10(a). It can be observed that with the variety of the rotating speed, the bearing force decreases in the entire rotating speed range while taking the gearbox stiffness into account, and the most significant decrease of the bearing force occurs at the resonance point. The main reason for this is that the vibration of the gearbox is introduced into the entire dynamic system, and the variety of the system stiffness affects the supporting force significantly, and this also means that the force transmitted by the transmission system to the other structural systems is reduced considerably. The gearbox flexibility in the transmission system helps obtain more accurate dynamic analysis results, which provides a theoretical basis for the optimization design and vibration and noise reduction of the system. Then, the amplitude-frequency response curve of the bearing force of the uncoupling and coupling gearbox in the rotating speed range of 200∼8000 r/min is obtained when the torque increases, which is revealed in Figure 10(b). To facilitate the analysis, under torques of 300 Nm and 600 Nm, the bearing forces of uncoupling and coupling the gearbox at rotating speeds of 3000 r/min, 4700 r/min, and 7100 r/min are calculated, respectively, as shown in Table 5. By observing the data in the table, it can be determined that the change rate of the bearing force is large when the torque is large at the same speed.Figure 10
Amplitude-frequency response curves of the bearing forces of the uncoupling and coupling model under different torques.(a)T = 300 Nm. (b) T = 600 Nm.
(a)
(b)Table 5
Change rate of the bearing force under different models.
Torque (Nm)
Speed (r/min)
Bearing force of coupling model (N)
Bearing force of uncoupling model (N)
Rate of change (%)
300
3000
2213
3307
33
4700
6393
9808
35
7100
1278
1893
32
600
3000
2741
4292
37
4700
7964
11700
40
7100
1545
2431
36According to the above analysis, the flexibility of the gearbox has an impact on the meshing force and the bearing force, of which the influence on the bearing force is prominent. In the analysis of the overall dynamic response of the coupling system supported by the thin-walled gearbox, the effect of the gearbox flexibility needs to be taken into account.


## 4.1.1. Influence of the Gearbox Flexibility on the Meshing Force

In order to investigate the gearbox flexibility’s effect on the meshing force of the gear transmission system, the meshing force of the uncoupling model (only considering the spiral bevel gear transmission system) and the coupling model (the spiral bevel gear-gearbox coupling model considering the gearbox flexibility) are obtained. Figure9 shows the amplitude-frequency response curve of the meshing force of the uncoupling and coupling mode at a rotation speed range of 200∼8000 r/min. As seen in Figure 9, the stiffness of the gearbox is introduced into the spiral bevel gear transmission system, which reduces the overall stiffness of the system, and absorbs and transmits part of the energy, resulting in a reduction in the meshing force at all rotational speeds. However, the magnitude of reduction in the meshing force is different at various rotational speeds. By comparing Figure 9(a) with Figure 9(b), it can be observed that, with the variety of the torque, the effect of the gearbox flexibility on meshing force is also varied. Under torques of 300 Nm and 600 Nm, the meshing forces of the uncoupling model and the coupling model are calculated, respectively, at rotating speeds of 3000 r/min, 4700 r/min, and 7100 r/min. The specific values are shown in Table 4. By analysing the data in Table 4, it is demonstrated that the effect of the gearbox flexibility on the meshing force is also larger when the torque is larger at the same speed.Figure 9
Amplitude-frequency response curve of the meshing forces of the uncoupling and coupling model under different torques. (a)T = 300 Nm. (b) T = 600 Nm.
(a)
(b)Table 4
Change rate of the meshing force under different models.
Torque (Nm)
Speed (r/min)
Meshing force of coupling model (N)
Meshing force of uncoupled model (N)
Rate of change (%)
300
3000
46570
49240
5.4
4700
66190
69000
4.1
7100
42130
44920
6.2
600
3000
71170
76010
6.4
4700
95210
100100
4.8
7100
70070
75880
7.6


## 4.1.2. Influence of the Gearbox Flexibility on the Bearing Force

Similarly, to study the influence of the gearbox flexibility on the bearing force, the differential equations of the system are solved in order to obtain the amplitude-frequency response curves of the bearing force of the uncoupling model and coupling model in the speed range of 200∼8000 r/min, as shown in Figure10(a). It can be observed that with the variety of the rotating speed, the bearing force decreases in the entire rotating speed range while taking the gearbox stiffness into account, and the most significant decrease of the bearing force occurs at the resonance point. The main reason for this is that the vibration of the gearbox is introduced into the entire dynamic system, and the variety of the system stiffness affects the supporting force significantly, and this also means that the force transmitted by the transmission system to the other structural systems is reduced considerably. The gearbox flexibility in the transmission system helps obtain more accurate dynamic analysis results, which provides a theoretical basis for the optimization design and vibration and noise reduction of the system. Then, the amplitude-frequency response curve of the bearing force of the uncoupling and coupling gearbox in the rotating speed range of 200∼8000 r/min is obtained when the torque increases, which is revealed in Figure 10(b). To facilitate the analysis, under torques of 300 Nm and 600 Nm, the bearing forces of uncoupling and coupling the gearbox at rotating speeds of 3000 r/min, 4700 r/min, and 7100 r/min are calculated, respectively, as shown in Table 5. By observing the data in the table, it can be determined that the change rate of the bearing force is large when the torque is large at the same speed.Figure 10
Amplitude-frequency response curves of the bearing forces of the uncoupling and coupling model under different torques.(a)T = 300 Nm. (b) T = 600 Nm.
(a)
(b)Table 5
Change rate of the bearing force under different models.
Torque (Nm)
Speed (r/min)
Bearing force of coupling model (N)
Bearing force of uncoupling model (N)
Rate of change (%)
300
3000
2213
3307
33
4700
6393
9808
35
7100
1278
1893
32
600
3000
2741
4292
37
4700
7964
11700
40
7100
1545
2431
36According to the above analysis, the flexibility of the gearbox has an impact on the meshing force and the bearing force, of which the influence on the bearing force is prominent. In the analysis of the overall dynamic response of the coupling system supported by the thin-walled gearbox, the effect of the gearbox flexibility needs to be taken into account.


## 4.2. Influence of the Shaft Angle on the Dynamic Characteristics of the Spiral Bevel Gear Full Coupling System

Based on the research in the previous section, the influence of the flexibility of the thin-walled gearbox on the spiral bevel gear transmission system cannot be ignored. Thus, this section directly studies the effect of the shaft angle on the dynamic characteristics of the spiral bevel gear-shaft-bearing-gearbox full coupling system.


### 4.2.1. Effect of the Shaft Angle on the Axis Trajectory of the Gear

Studying the axis locus of the critical nodes can predict whether the geared rotor transmission system will encounter rub-impact during operation or not, and the investigation of the gear’s axis locus is especially helpful in understanding the meshing status of the gear. It is assumed that the speed of the driving spiral bevel gear is 4200 r/min, the torque is 600 Nm, and the other parameters keep invariant. The variation of the axial trajectory of the driven spiral bevel gear is analysed when the shaft angles of the spiral bevel gear transmission system are 70, 90, 110, and 126, as shown in Figure11. The centre of the trajectory of the driven spiral bevel gear is varied when the shaft angle has diverse values, as shown in Figure 11(a). The centres of the axis locus under each shaft angle are set to the same position to observe the variety of the direction and size of the axis locus. It can be seen in Figure 11(b) that the variation in the shaft angle not only affects the direction of the axis locus but also alters the amplitude of the axis locus. The calculation shows that the axis locus of the other nodes will also change with the variety of the shaft angle. Therefore, the impact of the shaft angle must be taken into account when studying the vibration features of the nonorthogonal spiral bevel gear-shaft-bearing-gearbox full coupling system.Figure 11
Axis trajectory diagram of a driven gear under various shaft angles.
(a)
(b)


### 4.2.2. Effect of the Shaft Angle on the Meshing Force and Bearing Force

Supposing that the driving gear’s rotation speed is 4200 r/min, the torque is 600 Nm, and other parameters kept invariant, then the meshing force and the bearing force variation in each direction are obtained when the shaft angle is 70, 90, 110, and 126, respectively, as shown in Figure12. It can be demonstrated from the figure that the meshing force and the bearing force of the spiral bevel gear pair are various with the variety of the shaft angles. The meshing force varies because the pitch angles of the pinion and the gear are various under different shaft angles, and the projection results of the vibration displacement of the gear pair in various directions to the meshing line vary. At the same time, the force’s components decomposed into diverse directions are also different, resulting in variations in the bearing force for each direction. In conclusion, the vibration behaviours of the spiral bevel gear pair and the full coupling system will vary with the variety of the shaft angle.Figure 12
Meshing force and bearing force under various shaft angles. (a) Meshing force. (b) Bearing force in theX direction. (c) Bearing force in the Y direction. (d) Bearing force in the Z direction.
(a)
(b)
(c)
(d)


## 4.2.1. Effect of the Shaft Angle on the Axis Trajectory of the Gear

Studying the axis locus of the critical nodes can predict whether the geared rotor transmission system will encounter rub-impact during operation or not, and the investigation of the gear’s axis locus is especially helpful in understanding the meshing status of the gear. It is assumed that the speed of the driving spiral bevel gear is 4200 r/min, the torque is 600 Nm, and the other parameters keep invariant. The variation of the axial trajectory of the driven spiral bevel gear is analysed when the shaft angles of the spiral bevel gear transmission system are 70, 90, 110, and 126, as shown in Figure11. The centre of the trajectory of the driven spiral bevel gear is varied when the shaft angle has diverse values, as shown in Figure 11(a). The centres of the axis locus under each shaft angle are set to the same position to observe the variety of the direction and size of the axis locus. It can be seen in Figure 11(b) that the variation in the shaft angle not only affects the direction of the axis locus but also alters the amplitude of the axis locus. The calculation shows that the axis locus of the other nodes will also change with the variety of the shaft angle. Therefore, the impact of the shaft angle must be taken into account when studying the vibration features of the nonorthogonal spiral bevel gear-shaft-bearing-gearbox full coupling system.Figure 11
Axis trajectory diagram of a driven gear under various shaft angles.
(a)
(b)


## 4.2.2. Effect of the Shaft Angle on the Meshing Force and Bearing Force

Supposing that the driving gear’s rotation speed is 4200 r/min, the torque is 600 Nm, and other parameters kept invariant, then the meshing force and the bearing force variation in each direction are obtained when the shaft angle is 70, 90, 110, and 126, respectively, as shown in Figure12. It can be demonstrated from the figure that the meshing force and the bearing force of the spiral bevel gear pair are various with the variety of the shaft angles. The meshing force varies because the pitch angles of the pinion and the gear are various under different shaft angles, and the projection results of the vibration displacement of the gear pair in various directions to the meshing line vary. At the same time, the force’s components decomposed into diverse directions are also different, resulting in variations in the bearing force for each direction. In conclusion, the vibration behaviours of the spiral bevel gear pair and the full coupling system will vary with the variety of the shaft angle.Figure 12
Meshing force and bearing force under various shaft angles. (a) Meshing force. (b) Bearing force in theX direction. (c) Bearing force in the Y direction. (d) Bearing force in the Z direction.
(a)
(b)
(c)
(d)


## 4.3. Influence of Unbalance on the Dynamic Characteristics of a Spiral Bevel Gear Full Coupling System

### 4.3.1. The Influence of the Unbalance on Frequency Response Characteristics of the Coupling System

An unbalance amount, which is equal to5×10−4kg⋅m, is applied to the drive shaft’s node 4, with no imbalance given to the other nodes. The driving shaft speed n1 is in the range of 200∼800 r/min, and the three-dimensional spectrum of the full coupling system at the driving and driven gears are obtained, respectively, as shown in Figure 13. As seen from Figure 13(a), because of the coupling effect between the rotor and the gear, the X-direction displacement vibration spectrum of the driving gear includes the rotation frequency of the driving shaft f1, the gear meshing frequency fm, and its harmonic frequency 2fm. Meanwhile, it can be seen from Figure 13(b) that the torsional vibration spectrum of the driving gear also includes the rotation frequency of the driving shaft f1, the meshing frequency fm, and its harmonic frequency 2fm. As seen from Figure 13(c), the X-direction displacement vibration spectrum of the driven gear includes the meshing frequency fm and its harmonic frequency 2fm. The imbalance of the driving shaft has little effect on the vibration displacement of the driven gear in the X direction, where the rotation frequency amplitude of the driving shaft is small, which is not revealed in the figure. As seen from Figure 13(d), because of the coupling effect of bending and torsion of the spiral bevel gear pair, the driven gear’s torsional vibration not only contains the meshing frequency fm and its harmonic frequency 2fm but also the driving shaft’s rotation frequency f1, where the rotational frequency amplitude increases with the rise of the speed. Thus, the effect of the imbalance of the rotating shaft on the gear’s vibration in the torsional direction is more prominent than that on the lateral vibration.Figure 13
Three-dimensional spectrum of the coupling system as the unbalance is5×10−4kg⋅m. (a) Lateral vibration of the driving gear in the X direction. (b) Torsional vibration of the driving gear. (c) Lateral vibration of the driven gear in the X direction. (d) Torsional vibration of the driven gear.
(a)
(b)
(c)
(d)The unbalance amount applied on the driving shaft’s node 4 is changed to10×10−4kg⋅m, and the other conditions remain unchanged. The speed of the driving shaft n1 varies from 200 r/min to 8000 r/min, and the three-dimensional spectrum of the full coupling system at the pinion and gear are acquired, as shown in Figure 14. Comparing Figures 13(a) and 13(b) and 14(a) and 14(b), it can be demonstrated that the amplitude corresponding to the rotating frequency of the driving gear in the lateral and torsional vibration spectrograms increases significantly when the imbalance on the driving shaft is rising, while the amplitude corresponding to the meshing frequency hardly changes. It can be seen from Figure 14(c) that the driven gear’s vibration spectrum exhibits the driving shaft’s rotating frequency f1 at a high rotational speed, indicating that the coupling effect of the driving shaft’s imbalance on the driven gear’s vibration displacement is also risen with an increase in the unbalance of the driving shaft. Comparing Figures 13(d) and 14(d), it can be revealed that the amplitude of the rotation frequency f1 rises in the driven gear’s torsional vibration spectrum, which also indicates that the increase of the driving shaft’s unbalance increases the effect of the driven gear’s torsional vibration, and the increased amplitude is more obvious than that of the transverse vibration.Figure 14
Three-dimensional spectrum of the coupling system as the unbalance is10×10−4kg⋅m. (a) Lateral vibration of the driving gear in the X direction. (b) Torsional vibration of the driving gear. (c) Lateral vibration of the driven gear in the X direction. (d) Torsional vibration of the driven gear.
(a)
(b)
(c)
(d)


### 4.3.2. Influence of the Unbalance on the Meshing Force and Bearing Force

In order to analyse the influence of the unbalance of the rotating shaft on the meshing force and the bearing force, the unbalance amount on the driving shaft’s node 4 is taken as 0,5×10−4kg⋅m, and 10×10−4kg⋅m, respectively. The amplitude-frequency response curves of the meshing force and the bearing force within a rotation speed range of 200∼8000 r/min are calculated, as shown in Figure 15. It can be revealed from the figure that the meshing force and the bearing force have little influence at low rotational speeds with the increase of the eccentricity; however, the eccentricity has a great influence on the meshing force and bearing force when speed continuously increase. Therefore, when the speed of the full coupling system is high, the eccentric amount of the rotating shaft should be reduced to diminish its impact on the meshing force and bearing force so as to ensure the stable and reliable operation of the system.Figure 15
Amplitude-frequency response curves of the meshing force and the bearing force. (a) Amplitude-frequency response curve of the bearing force. (b) Amplitude-frequency response curve of the meshing force.
(a)
(b)


## 4.3.1. The Influence of the Unbalance on Frequency Response Characteristics of the Coupling System

An unbalance amount, which is equal to5×10−4kg⋅m, is applied to the drive shaft’s node 4, with no imbalance given to the other nodes. The driving shaft speed n1 is in the range of 200∼800 r/min, and the three-dimensional spectrum of the full coupling system at the driving and driven gears are obtained, respectively, as shown in Figure 13. As seen from Figure 13(a), because of the coupling effect between the rotor and the gear, the X-direction displacement vibration spectrum of the driving gear includes the rotation frequency of the driving shaft f1, the gear meshing frequency fm, and its harmonic frequency 2fm. Meanwhile, it can be seen from Figure 13(b) that the torsional vibration spectrum of the driving gear also includes the rotation frequency of the driving shaft f1, the meshing frequency fm, and its harmonic frequency 2fm. As seen from Figure 13(c), the X-direction displacement vibration spectrum of the driven gear includes the meshing frequency fm and its harmonic frequency 2fm. The imbalance of the driving shaft has little effect on the vibration displacement of the driven gear in the X direction, where the rotation frequency amplitude of the driving shaft is small, which is not revealed in the figure. As seen from Figure 13(d), because of the coupling effect of bending and torsion of the spiral bevel gear pair, the driven gear’s torsional vibration not only contains the meshing frequency fm and its harmonic frequency 2fm but also the driving shaft’s rotation frequency f1, where the rotational frequency amplitude increases with the rise of the speed. Thus, the effect of the imbalance of the rotating shaft on the gear’s vibration in the torsional direction is more prominent than that on the lateral vibration.Figure 13
Three-dimensional spectrum of the coupling system as the unbalance is5×10−4kg⋅m. (a) Lateral vibration of the driving gear in the X direction. (b) Torsional vibration of the driving gear. (c) Lateral vibration of the driven gear in the X direction. (d) Torsional vibration of the driven gear.
(a)
(b)
(c)
(d)The unbalance amount applied on the driving shaft’s node 4 is changed to10×10−4kg⋅m, and the other conditions remain unchanged. The speed of the driving shaft n1 varies from 200 r/min to 8000 r/min, and the three-dimensional spectrum of the full coupling system at the pinion and gear are acquired, as shown in Figure 14. Comparing Figures 13(a) and 13(b) and 14(a) and 14(b), it can be demonstrated that the amplitude corresponding to the rotating frequency of the driving gear in the lateral and torsional vibration spectrograms increases significantly when the imbalance on the driving shaft is rising, while the amplitude corresponding to the meshing frequency hardly changes. It can be seen from Figure 14(c) that the driven gear’s vibration spectrum exhibits the driving shaft’s rotating frequency f1 at a high rotational speed, indicating that the coupling effect of the driving shaft’s imbalance on the driven gear’s vibration displacement is also risen with an increase in the unbalance of the driving shaft. Comparing Figures 13(d) and 14(d), it can be revealed that the amplitude of the rotation frequency f1 rises in the driven gear’s torsional vibration spectrum, which also indicates that the increase of the driving shaft’s unbalance increases the effect of the driven gear’s torsional vibration, and the increased amplitude is more obvious than that of the transverse vibration.Figure 14
Three-dimensional spectrum of the coupling system as the unbalance is10×10−4kg⋅m. (a) Lateral vibration of the driving gear in the X direction. (b) Torsional vibration of the driving gear. (c) Lateral vibration of the driven gear in the X direction. (d) Torsional vibration of the driven gear.
(a)
(b)
(c)
(d)


## 4.3.2. Influence of the Unbalance on the Meshing Force and Bearing Force

In order to analyse the influence of the unbalance of the rotating shaft on the meshing force and the bearing force, the unbalance amount on the driving shaft’s node 4 is taken as 0,5×10−4kg⋅m, and 10×10−4kg⋅m, respectively. The amplitude-frequency response curves of the meshing force and the bearing force within a rotation speed range of 200∼8000 r/min are calculated, as shown in Figure 15. It can be revealed from the figure that the meshing force and the bearing force have little influence at low rotational speeds with the increase of the eccentricity; however, the eccentricity has a great influence on the meshing force and bearing force when speed continuously increase. Therefore, when the speed of the full coupling system is high, the eccentric amount of the rotating shaft should be reduced to diminish its impact on the meshing force and bearing force so as to ensure the stable and reliable operation of the system.Figure 15
Amplitude-frequency response curves of the meshing force and the bearing force. (a) Amplitude-frequency response curve of the bearing force. (b) Amplitude-frequency response curve of the meshing force.
(a)
(b)


## 5. Conclusion

A full coupling dynamic model of the spiral bevel gear was built, taking the flexibility of the shaft, the bearing supporting, the unbalance amount, a time-varying meshing stiffness, a comprehensive transmission error, and the gearbox flexibility into account. The Timoshenko beam element was employed to simulate the flexible shaft in the model, and the gears and bearings were treated according to the concentrated mass. Because of the complex structure of the gearbox system, the substructure method was applied to extract its mass, stiffness, and damping matrix. The full coupling dynamic equations of the spiral bevel gear-shaft-bearing-gearbox were constructed, and the coupling dynamic equations were solved via the Newmark numerical algorithm. The following conclusions were obtained.The stiffness of the gearbox affects the meshing force and bearing force of the spiral bevel gear transmission system. The meshing force and bearing forces are smaller than those of the rigid support when considering the gearbox flexibility, wherein the gearbox flexibility has a more significant effect on the bearing force.The shaft angle of the full coupling system of the nonorthogonal spiral bevel gear affects the axis locus, meshing force, and bearing force simultaneously. For various shaft angles, the position, size, and direction of the gear’s axis locus vary. Because the pitch angles of the pinion and gear are different under various shaft angles, the meshing force and bearing force of the system also vary. As a result, the vibrational characteristics of the gear pair and the entire coupling system also change.Because of the coupling of bending and torsion of the system, the unbalance of the rotating shaft exerts an effect on the vibration of the spiral bevel gear transmission system in all directions, where the torsional vibration has the greatest influence, and the effect increase with the variety of unbalance. The unbalance of the rotating shaft also affects the meshing force and the bearing force. The influence is little at a low rotational speed. However, the unbalance exerts a significant effect on the meshing force and bearing force when speed continuously increases.



---

*Source: 9065215-2019-09-08.xml*