Vibration Response Analysis of a Gear-rotor-bearing System Considering Steady-state Temperature

: As an important factor leading to the failure of gear system, the study of thermal effect is insufficiently deep. Based on the finite element nodal method, a more comprehensive dynamic model of gear-rotor-bearing system is established, which considers the thermal related material properties, time-varying meshing stiffness (TVMS), backlash and friction, gyroscopic effect. The constitutive relation of beam element considering steady-state temperature is reconstructed, and thermal node load is formulated. Considering the influence of temperature on the material properties of flexible shaft and gear, the thermal related TVMS and thermal backlash are obtained. The dynamic response of the system under different steady-state temperature fields is compared, and the influence of hot backlash is studied, then the thermal related vibration characteristics are obtained. Besides, the influence of bearing type on bearing force and axial trajectory is studied. The results show that the system motion changes from period to chaos with the temperature increase in part of the speed range. The appropriate backlash is helpful to restrain the chaotic motion caused by temperature rise. Moreover, the temperature can significantly increase the axial bearing force, and the appropriate bearing can reduce the axial displacement.


Introduction
In the aerospace and other fields, the gear transmission system under high-speed and heavy-load conditions will generate great heat due to friction, which makes the temperature inside the gearbox much higher than the room temperature [1][2][3][4]. The extremely high operating temperature of the gear system changes the material properties of the internal components of the gearbox, and causes even plastic deformation of the gear under load, affects the transmission performance and makes the dynamic response unpredictable [5]. Studying the influence of temperature on the gear transmission system is one of the hot research issues in gear system design and application.
Many researchers have studied the nonlinear dynamics of the gear system. Kahraman and Singh [6,7] studied the nonlinear characteristics of the single-degree-of-freedom spur gear system containing internal and external excitations, and the gear bearing system including clearance through the harmonic balance method. Wang [8] studied the simulation model of the influence of random error on the vibration response of gear system. Vaishya [9] proposed a new gear system model including time-varying meshing stiffness, viscous damping, and sliding friction, then studied the effect of sliding friction on system response and stability. The model excitation includes unloaded transmission error, external torque and periodic sliding friction. Li [10] used on the incremental harmonic balance method to establish a single-degree-of-freedom gear model and analyzed the dynamic characteristics of the system considering dynamic backlash, viscous damping, and periodic external excitation. Cao et al. [11] studied the dynamic characteristics of planetary gear system transmission under non-constant load. Previous papers on gear dynamics researched gear systems with fewer degrees of freedom and generally used the lumped-parameter method, rarely considering the effects of shaft and bearing. However, in the actual gear transmission, the flexibility of the shaft and the stiffness and damping of the bearing have a significant influence on the transmission characteristics of the gear pair and cannot be ignored. Velex [12] studied the influence of mounting error and shape deviation on the gear system dynamics. Kim [13] considered the flexibility of the bearing and proposed a new six-degree-of-freedom gear system model to study the effect of the translational motion due to bearing deformation on the gear system. Ma et al. [14][15][16] studied the influence of tooth tip modification, tooth profile deviation and gear crack on gear vibration response. First, the time-varying meshing stiffness was obtained based on the analytical model, which was then verified by the finite element model, and finally brought to the gear-rotor system model for dynamic response analysis. Inalpolat [17] investigated the effect of indexing error on the dynamic response of the gear system. Hu et al. [18]studied the effects of different tooth profile repair parameters on the dynamic response of the high-speed gear-rotorbearing system. They [19] obtained the time-meshing stiffness of gears with crack faults based on the improved potential energy method, introduced short-term and long-term transmission error excitations, and studied the effects of crack length and transmission error excitation on the dynamic characteristics of the gear system. Chen [20] established the dynamic model of spiral bevel gear-rotor-bearing system to study the influence of elastic ring squeeze film damper on system nonlinearity. Based on the finite element/contact mechanics approach, Liu et al. [21] established the roller bearing model and studied the influence of its stiffness on gear-bearing system vibration.
In high-speed gear systems, the effect of the gyroscopic effect on vibration is pronounced. Peng and Lim [22] studied the influence of the gyroscopic effect on the hyperbolic gear dynamics and found that the gyroscopic effect had little effect at low speeds, but had a significant effect at high speeds. Theoretically, it may be that the higher speed of the rotor due to mass imbalance, eccentricity and other factors causes the rotor to rotate more violently and has a more significant impact on the gear system dynamics. Chen et al. [23] analyzed the dynamic response of the herringbone rotor-bearing system. Hu [24] compared the difference of the dynamic response between the torsional gear dynamics model and the gear-rotor-bearing model coupled with lateral torsion and analyzed the effect of the gyroscopic effect on the lateral-torsional vibration of the gear-rotor-bearing system.
The temperature has become one of the main factors affecting the dynamic performance of the gear system. Since Blok [25] proposed to calculate the flash temperature of the tooth surface during the meshing process, many scholars began to study the temperature field of the gear body and tooth. Li [26] summarized the different temperature fields in the gearbox, and considered that when the lubrication condition of the small gearbox is good, the gear system can be regarded as a steady temperature field. Hohn et al. [27] found that the immersion depth had a great influence on the overall temperature of the gear through experimental research, and the maximum temperature is 300 ℃. Diogo [28] studied the influence of gear speed, load and lubrication conditions on the overall temperature of the gear. In the case of dry friction or lack of lubrication, the maximum temperature of gear teeth can reach 500 ℃. Zhou et al. [29] regarded the gear as a threedimensional model connected by the thermal resistances of heat conduction and heat convection, and established a novel thermal network model to predict the bulk temperature and flash temperature of spur gears.
In order to study the effect of tooth surface temperature on the dynamic characteristics of the gear system, it is necessary to quantify the gear temperature. Many scholars introduced temperature into the meshing stiffness. Gou [30] calculated the flash temperature rise during gear meshing through the Blok formula, and analyzed the influence of the flash temperature rise on the deformation of the tooth profile, and then deduced the expression of the meshing stiffness that changes with temperature according to the Hertz contact theory. Based on the lumped-parameter method, he established a gear-rotor-bearing system that considered factors such as flash temperature, meshing stiffness, friction, backlash, and comprehensive transmission error. He analyzed the nonlinear phenomenon of the system through the bifurcation diagram, phase diagram, and Poincaré diagram. The results showed that the higher tooth surface contact temperature caused the periodic deterioration of the gear system. Shi [31] considered factors such as contact temperature, time-varying meshing stiffness, meshing friction, and transmission error, and established a three-degreeof-freedom nonlinear model of the gear system to analyze the effect of tooth surface flash temperature rise. Pan [32] analyzed the influence of friction coefficient and normal load on the flash temperature, then established a ten degree-of-freedom coupled nonlinear dynamic model and analyzed the influence of flash temperature, random load, and fractal recoil on its dynamic characteristics. The simulation results show that the flash temperature of the tooth surface could significantly change the periodicity of the system movement. Luo and Li [33,34] proposed the concept of thermal stiffness to study the influence of thermal gear stiffness. The thermal stiffness curves were obtained by the analytical method and the finite element method under the thermos-elastic coupling condition, and they were consistent. The mechanism of the influence of heat on the dynamic response was studied experimentally. The results shown that temperature affects the dynamic response of the gear system by changing the involute meshing characteristics and meshing rigidity of the gear. Sun et al. [35] established a single-degree-freedom lumped parameter gear system model, and based on the change of material properties, studied the influence of temperature on the dynamic response of the system. The results shown that the influence of temperature was concentrated in the middle and high speed.
Hence, in this paper, a more comprehensive dynamic model of the gear-rotor-bearing system is established to study the influence of thermal effect on the dynamic response. In this paper, we use the finite element node method to establish a dynamic model of the gear-rotor-bearing system considering nonlinear factors such as steady-state temperature, thermal node load, TVMS, backlash, friction and gyroscopic effect, etc. and analyze its dynamic characteristics. Section 2, the expressions of thermal node load, TVMS and thermal expansion error are proposed under the steady temperature field, and verified by finite element method. Section 3 considers the supporting bearing, flexible shaft and meshing gear, introduces timevarying meshing stiffness, backlash, transmission error excitation, rotor gyroscopic effect and other nonlinear factors, then establishes the gear-rotor-bearing system with the finite element method. Section 4 introduces the parameters of the gear-rotor-bearing system, analyzes the vibration response of the gear system at different temperatures and backlash based on the bifurcation diagram, phase portrait and Poincaré map, and studies the influence of the bearing type on the system. Finally, a conclusion is drawn in Section 5.

Thermal node load
The gear systems under high-speed and heavy-load working conditions generate huge heat due to the friction between the tooth surface and the tooth surface and the bearing. Under loss-of-lubrication condition, the temperature in the gearbox will rise rapidly [36], and temperature changes will cause thermal stress inside the structure, and finally result in the soft surface and deterioration of gear transmission performance. The influence of steady-state temperature field on the flexible shaft element is shown in this section. According to Duamel-Neumann linear thermal stress theory: the stress (or strain) caused by temperature change and the stress (or strain) caused by external force can be superimposed, which can be expressed as: where,  is the total strain,  E and  T are the strain caused by external force and the strain caused by temperature, respectively. When considering temperature, the temperature causes thermal strain due to thermal expansion. The finite element method treats this change as the initial strain, and for isotropic materials, the thermal strain produces only direct strain and no shear strain, as shown in Fig 1 above. For the three-dimensional solid model element, the thermal strain is as follows: here,  i0 and  ii (i=x,y,z)are the initial strain of i direction and the shear strain in the ii plane, respectively. a is the coefficient of linear expansion and ΔT is the amount of temperature change. When considering the temperature factor, the constitutive relation can be expressed as: where  and  are the strain matrix and stress matrix. D is the transformation matrix between  and .
When thermal strain is generated, if the three-dimensional solid is constrained and cannot expand freely, thermal stress will occur. For linear elastic solids, the total potential energy can be written as: where v and s 1 are the element volume and boundary surface, u and r are the displacements vector of this element. P k and u k are the load vector and displacement vector at point k.
By the principle of minimum potential energy, the following expression be obtained as: (5) where N is the total number of displacement degrees, q donates the displacement vector.
Taking Eq. (4) into Eq. (5), we get: where B denotes the transformation matrix between  and q, and N denotes beam element shape function matrix.
Then we can write Eq. (6) as: (9) where K is the element stiffness matrix, and F T is defined as the thermal node load.
The thermal node load vector of the two-node model Timoshenko beam element can be expressed as: where l, I ps ,  and A represent the length of the beam element, polar second area moment of inertia, shear correction factor and sectional area. E * and G * denote elastic modulus and shear modulus obtaining by fitting the values at different temperatures. Table.1 shows the material properties of carbon steel. In this paper, the flexible shaft and gears are made of carbon steel.
The finite element (FE) method is used to verify the correctness of the thermal node load. In the commercial finite element software, firstly, a three-dimensional line model of B31 element with a length of 20 mm is established. The temperature dependent material properties as shown in Table.1 are defined, and the linear expansion coefficient is set to be 1.1710 -5 1/℃. Secondly, create a circular beam section with a radius of 15mm and assign the section. Thirdly, constraints are added to the whole analysis model, and the degrees of freedom θ z is released, while the degrees of freedom in other directions are locked to simulate the actual rotational pairs of the flexible axis. Fourth, add the temperature field. Set the initial steady-state temperature field value as 25 ℃, and set the simulation temperature from 75 ℃ to 425 ℃. Finally, the thermal node load (NFORC3) is shown in Table.2. Because the shear effect is not considered in the analysis of B31 element, only the axial stress and strain can be output. Only axial force (Z direction) is compared here. The results show that the calculated results of thermal node load are very close to those of FM method.

Time-varying meshing stiffness
As one of the most important excitations of gear system, how to accurately calculate the time-varying mesh stiffness (TVMS) is the premise of dynamic analysis. The TVMS model used in this paper is described in Ref. [35], and its solution expressions are as follows: where  1 and  2 denote the correction coefficients of the fillet-foundation stiffness, k h , k a , k b , k s and k f are the non-linear Hertzian contact stiffness, axial compressive stiffness, bending stiffness, shear stiffness and fillet-foundation stiffness, respectively, which can be found in Ref. [35]. n is the number of teeth pairs in engagement. By taking the elastic modulus E * , shear modulus G * into the above formula, the thermal TVMS without considering thermal expansion can be obtained.

Thermal profile error
For the gear system, enough backlash should be left to lubricate and prevent thermal expansion from sticking. In this paper, the thermal expansion of gear is regarded as the change of backlash. Assuming that the temperature rise after thermal deformation is ΔT, the polar coordinate equation of tooth profile with thermal expansion can be expressed as: here, r b and  k1 denote the base radius and pressure angle, respectively. r k1 , θ k1 and are the involute curve parameters after thermal deformation.
The actual thermal deformation of tooth profile can be obtained by adding the deformation of tooth thickness and tooth height. The tooth profile with thermal expansion can be derived as: where s, r,  and  k are the tooth thickness, radius, pressure angle and pressure angle at point k of indexing circle before thermal deformation.
In order to quantitatively describe the difference between the actual tooth profile curve and the theoretical tooth profile curve along the meshing line, the involute error Δe is defined (Ref. [30,37]

Dynamic model of a gear-rotor-bearing system
The finite element model of the gear-rotor-bearing system is composed of gear rotor, meshing gear, flexible shaft and supporting bearing. In the finite element method, the number of nodes and the number of degrees of freedom are selected according to the actual situation, and the shaft is simplified to the Timoshenko beam element while the bearing is simplified to the support stiffness matrix and damping matrix. The following mainly introduces the four parts of the modeling and system model assembly process.

Gear rigid rotor
The gear rigid rotor model is shown in Fig 2. This model considers six DOFs. x and y are the translational displacements of the rotor in the X and Y direction, respectively. θ x and θ y are the torsion angles of the gear rotor around the Z-axis in the XY plane and the X-axis around the YZ plane, respectively. z is the axial displacement of the rotor in the Z direction, and θ z is the rotation angle of the gear rotor about the Z-axis. x y z q (19) The dynamics equations are expressed as the following through the Lagrange equation: (20) where M d denotes the mass matrix, G d represents the gyroscopic matrix and F d is force due to eccentricity, which can be found in Ref [23]. Fig 3 shows the meshing model of spur gear pair. Springs and damping elements are used to simulate the TVMS and meshing damping of the gears. The Cartesian inertial local coordinates are respectively established at each gear center Op and Og. As shown in Fig 3, x j , y j and z j are the translational displacements of gear j (j=p,g), and θ x , θ y and θ z are the rotational displacement around direction X、Y and Z, respectively. Considering the translation and torsional displacement, the dynamic deformation of the gear pair along the meshing line can be obtained as:

Gear pair meshing modeling
e t Vq (21) where δ is the dynamic transmission error (DTE), which is one of the main parameters that affect the dynamic characteristics of the gear system. V is the coefficient vector and q g is the gear pair displacement vector. The expressions of V and q g are as follows: yg zg x y z x y z q (23) where,  is pressure angle, r bp and r bg are the base circle radius of the driving wheel and driven wheel, respectively. e(t) is the static transmission error, and it is generally expressed in the form of Fourier series as: here, e 0 represents the amplitude of static transmission error.  represents the gear meshing frequency,  is the initial phase, The mesh force between the gears can be expressed as the sum of the elastic force generated by the spring element and the damping force generated by the damping element: where k * m and c m are the thermal TVMS and damping.  m1 and  m0 are the backlash functions, and the expressions can be expressed as: where b * denotes half of the thermal backlash, which can be derived as: 27) here, b is half of the backlash, Δe p and Δe g denote thermal profile error of the driving wheel and driven wheel.
where a i and b i denote the non-linear coefficient of the friction force and moment, which can be found in Ref. [37]. The stiffness matrix and damping matrix of the meshing gear pair can be written as The equation of motions for gear pairs can be expressed as:

Flexible shaft element
As shown in Fig 4, the shaft is segmented, and each axis segment has two nodes. Timoshenko beam theory is used to model the shaft element. The deformation of the beam element includes bending deformation, shear deformation and torsional deformation. Each node has six degrees of freedom, namely three translational degrees of freedom x(s,t), y(s,t), z(s,t) and three rotational degrees of freedom θ x (s,t), θ y (s,t), θ z (s,t). The displacement function can be expressed as: (34) here, the specific elements of the shape function N of the shaft element can be referred to [23]. q e is the displacement vector of the two-node Timoshenko beam element, which can be written as:  (36) where the mass, gyroscopic, stiffness and load matrix are M e , G e , K e and F e .

Support bearing
To make the researched bearing modeling universal, the stiffness and damping parameters of the axial  2  20  30  10-11  20  30  2-3  20  30  11-12  20  30  3-4  20  30  12-13  20  30  4-5  10  30  13-14  10  30  5-6  10  30  14-15  10  30  6-7  20  30  15-16  20  30  7-8  20  30  16-17  20  30  8-9  20  30 17-18 20 30 Ignore the coupling effect between the bearing shaft in all directions and the rotational stiffness and damping around the z-axis. Then the force of the supporting bearing can be written as 3.5 Assembly of system dynamic model According to the above formula, gear meshing pair, flexible shaft and bearing modeling, the model of the entire gear-rotor-bearing system can be assembled using the finite element method. Adding the stiffness, mass, gyroscopic, and damping matrices of each part to the corresponding nodes respectively, the system dynamics equation can be obtained as follows: here, M s , G s , C s and K s are the mass、the gyroscopic、the damping and the stiffness matrix respectively. F s is the generalized force excitation and q s is the node displacement matrix.

Gear-rotor-bearing system model and parameters
The finite element node model of the gear-rotor-bearing system analyzed in this paper is shown in Fig 5. Both the input shaft and the output shaft are 9 nodes, and each node has six DOFs. The bearings are located on nodes 1 and 9 of the input shaft and nodes 10 and 18 of the output shaft, respectively. The meshing gear pairs are symmetrically distributed and are fixed on nodes 5 and 14. The parameters of the meshing gear pair, the flexible shaft parameters and the supporting bearings parameters are shown in Table.3, Table.4 and  Table.5, respectively. Three types of bearings are listed in Table.5. Bearing A# is the deep groove ball bearing with small axial stiffness. B# and C# are the angular contact ball bearings with larger axial stiffness. And the radial stiffness of C# is the largest. The temperature analyzed in this article is a constant ambient temperature in the gearbox, regardless of the temperature change with speed.

Vibration response
The gear-rotor-bearing dynamic model analyzed in this paper is modeled by the finite element method, which has many degrees of freedom. To reduce the calculation time, the Newmark method is used to solve the dynamic equations. In this paper, the initial backlash is b=35m and static transmission error amplitude is e 0 =20m.
The influence of temperature on the gear system can be explained by comparing the dynamic response under different temperatures. Fig 6a and Fig 6b illustrate the bifurcation diagram of dynamic transmission error (DTE) of gear system under different steady-state temperature fields. The system model considers the temperature dependent material properties, thermal node load, TVMS and thermal backlash. At 25 ℃, with the increase of rotating speed (500-20000rpm), the system undergoes periodic-1, chaotic, periodic-2, chaotic, multi-period, chaotic and periodic-2 motions. At 225 ℃, the system experiences periodic-1, shortly chaotic, periodic-1, chaotic, periodic-2, chaotic, periodic-2 and periodic-3 motions, and its dynamic behavior is more complex. In the range of 500-7300 rpm under 25℃, the system response is periodic-1 motion. But at 2850rpm under 225℃, there is a transient chaotic motion. At 13500rpm, the system motion changes from multi-period to chaos as the temperature increases. And in the range of 15500-20000rpm, the system changes from periodic-2 motion to periodic-3 motion.
However, in part of the speed range, the increase in temperature causes the system to move more easily. For example, with the increase of temperature, the system changes from chaotic to periodic-2 motion at 12850rpm, and the system changes from chaotic to quasi periodic-2 motion at 14500rpm. To analyze the cause of this phenomenon, the case without the thermal backlash is analyzed, as shown in Fig 6c. In this case, the range of chaotic motion speed of the system is larger. At 7100 rpm, Fig 6a and Fig 6b show that the system is in periodic-1 motion, while Fig 6c shows the system is in chaotic motion. Compared with the chaotic speed range in the medium-and high-speed region, the system is in chaos at 13500-14700 rpm under 25℃, 13000-14200rpm under 225℃ with variable backlash, and 12100-14300 rpm under 225℃ with variable backlash. The results indicate that enough backlash for thermal expansion can reduce the unpredictability of system motion caused by thermal effect to a certain extent.

Bearing force and axis orbit
To further understand the influence of the change of temperature field on the gear system, the following mainly analyzes the bearing force and axial orbit. Fig 7 shows the RMS values of the total bearing force of the bearing 2 of the gear-rotor-bearing system vary with speed at different temperatures. The bearing type is A#. The total bearing force is the square root of the force in the X, Y, and Z directions. As the temperature increases, the bearing force at each rotation speed gradually increases. For example, when the rotation speed is 1000 rpm, the total bearing force is 1088N at 25℃. At 125℃, the force becomes 1855N and increases by 70.5%. At 225℃, it becomes 3434N and increases by 215.6%. The peak values in Fig 7 well corresponds to the chaos phenomenon in the bifurcation diagram. At the speed of 2800rpm and in the ranges of 6000-7500rpm, 8100-8600rpm, 13100-14600rpm, the system response is chaotic motion in Fig 6a. In general, the bearing force increases with increasing temperature, but there is also special case. At 14400rpm, the bearing force at 125℃ is 3785N, but the bearing force is 3498N at 225℃. This may be due to the excessive vibration response under the temperature excitation. Fig 8 shows the vibration displacement amplitude of the driving axle under different bearings and different temperature fields at 19000rpm. The vibration displacement amplitude is the square root value of vibration displacement in X, Y and Z directions, which represents the distance from the theoretical node position (0,0,0). It can be found from Fig 8 that the vibration displacement increases with the increase of temperature. Node 1 (or node 9) has the largest growth rate. As shown in Fig 8a, compared with 9.8μm at 25℃, the amplitude is 81.5μm at 125℃, increased by 7.32 times, while the amplitude is 172.7μm at 225℃, increased by 17.61 times. Under the same condition, the vibration displacement amplitude decreases with the increase of bearing stiffness. The displacement of bearing A at node 1 is 172.6μm at 225℃, while bearings B# and C# are both 158μm.
The causes of the increase of the total bearing force are further analyzed. Fig 9 shows the proportion of the force in each direction when the support bearings are A#, B# and C# respectively. When the temperature is 25℃ (no thermal node load), the axial load ratio is 0, which means that there is no axial displacement. With the increase of temperature, the proportion of radial bearing force decreases gradually, and the axial bearing force increases rapidly. Obviously, the reason for the rapid increase of the total bearing force is the increase of axial force. Comparing (a), (b) and (c), increasing the axial stiffness can significantly increase the proportion of axial bearing force, while increasing the radial stiffness almost does not change the proportion of bearing force in all directions. The results show that for the gear bearing system considering thermal effect, the influence of axial stiffness is obviously greater than that of radial stiffness. Fig 10 shows the axis locus of bearing 2 under different bearing supports and different temperature fields at the same speed. The results show that: 1) the more the temperature rises, the more the axis trajectory deviates from the ideal position. 2) The deviation of bearing C# is smaller than that of bearings A# and B#. Because the stiffness of bearing C# is the largest in X and Y directions, its vibration displacement is the smallest under the same load.

Conclusions
In this paper, the influence of temperature effect on the dynamic response of gear-shaft-bearing system is studied. The concept of thermal node load is proposed, and its calculation formula is derived based on the thermal constitutive equation. The calculation accuracy of the equation is verified by the finite element method. In addition, considering the material properties, meshing stiffness and backlash related to temperature, a multi-degree dynamic model is established based on the finite element node method. The influence of thermal effect on the dynamic response is illustrated by analyzing the bifurcation diagrams of different steady-state temperature fields with the change of rotational speed. Additionally, the effects of different types of bearings on the bearing force and axis orbit are also discussed. The main conclusions are summarized as follows: (1) Temperature rise will change the nonlinear characteristics of the system, such as from periodic motion to chaotic motion. Considering the thermal effect is necessary to the gear dynamic modeling.
(2) The gear system should be equipped with proper backlash. The increase of temperature will reduce the backlash and restrain the chaotic motion of the system at some speeds.
(3) The bearing force of the system will increase significantly considering the thermal effect. When the temperature rise reaches 200 ℃, the axial bearing force will vary from 0 % to over 80%.
(4) By choosing the right type of bearing, the axial displacement can be significantly reduced. Increasing the stiffness of the radial bearing can reduce the axis orbit, and increasing the stiffness of the axial bearing can reduce the total vibration displacement. The main purpose of this paper is to study the influence of thermal effect on the gear-rotor-bearing system. The temperature field used is a steady uniform temperature field, and the generation and transmission of temperature are not considered. This part can be referred to Ref [37].