Vibration response analysis of a gear-rotor-bearing system considering steady-state temperature

Thermal effect is a crucial factor leading to failure of gear system. However, the lack of comprehensive nonlinear dynamics model limits the further study of thermal effects. The constitutive relation of beam element considering steady-state temperature is reconstructed, and thermal node load is formulated. Considering the influences of thermal expansion and temperature on material properties, a more comprehensive dynamic model of gear-rotor-bearing system is established based on the finite element node method. Nonlinear friction, high-speed gyroscopic effect, thermal-related time-varying meshing stiffness (TVMS) and thermal backlash are included in the model. The effects of temperature and bearing type on the vibration response of gear system are analyzed. 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 could restrain the chaotic motion caused by temperature rise. Moreover, the temperature significantly increases the axial bearing force, and the appropriate bearing could reduce the axial displacement. This research can further understand the influence of temperature on the dynamic response of gear system and guide the design of gear system.

Keywords Steady-state temperature Á Thermal node load Á Thermal backlash Á Vibration response Á Bearing force Backlash functions k Thermal expansion coefficient r

List of symbols
Stress matrix X Spin speed of the gear rotor n i ði ¼ 1; 2Þ Correction coefficient of the filletfoundation stiffness of tooth pair i

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 the incremental harmonic balance method to establish a single-degreeof-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-rotor-bearing 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-rotorbearing model coupled with lateral torsion and analyzed the effect of the gyroscopic effect on the lateraltorsional 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°C. 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°C. Zhou et al. [29] regarded the gear as a three-dimensional 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, 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 lumpedparameter 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-degree-of-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-offreedom 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 showed 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-offreedom 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. In Section 2, the expressions of thermal node load, TVMS and thermal expansion error are proposed under the steady temperature field and verified by the finite element method. Section 3 considers the supporting bearing, flexible shaft and meshing gear, introduces time-varying meshing stiffness, backlash, transmission error excitation, rotor gyroscopic effect and other nonlinear factors and 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 Sect. 5.

Thermal node load
The gear systems under high-speed and heavy-load working conditions generate huge heat due to the frictions between the tooth surfaces and between bearings. 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 Duhamel-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 e is the total strain, e E and e 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. For the three-dimensional solid model element, the thermal strain can be expressed as follows: here, e i0 and c 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 DT is the amount of temperature change. When considering the temperature factor, the constitutive relation can be expressed as: where e and r are the strain matrix and stress matrix. D is the transformation matrix between e and r. When thermal strain is generated, if the threedimensional 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: where N is the total number of displacement degrees and q i is the displacement vector of the element. Taking Eq. (4) into Eq. (5), we get: where B denotes the transformation matrix between e and q and N denotes beam element shape function matrix.
Then we can write Eq. (6) as: where K is the element stiffness matrix and F T is defined as the thermal node load.
The thermal node load of Timoshenko beam element can be expressed as: where l, I s , l and A represent the length of the beam element, 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. The modeling is divided into four steps: the first step, a three-dimensional line model of B31 element with a length of 20 mm is established in the commercial finite element software. The temperature-dependent material properties as shown in Table 1 are defined, and the linear expansion coefficient is set to be 1.17 9 10 -5 1/°C. The second step creates a circular beam section with a radius of 15 mm and assign the section. The third step, constraints are added to the whole analysis model, and the degrees of freedom h z are released, while the degrees of freedom in other directions are locked to simulate the actual rotational pairs of the flexible axis. The fourth step adds the temperature field. Set the initial steady-state temperature field value as 25°C, and set the simulation temperature from 75 to 425°C. 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 n 1 and n 2 denote the correction coefficients of the fillet-foundation stiffness, k h , k a , k b , k s and k f are the nonlinear 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 DT, the polar coordinate equation of tooth profile with thermal expansion can be expressed as: here, r b and a k1 denote the base radius and pressure angle, respectively. r k1 and h k1 are the involute curve parameters after thermal deformation. k is the thermal expansion coefficient. 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, a and a 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 De is defined (Ref. [30,37]). For the convenience of calculation, suppose a k1 = a k , then combining Eqs. (16) and (17), we can obtain the following expansion: 3 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. h x is the angular displacement of the gear rotor in the YZ plane, which is considered as a small rotation around the X-axis. h y is the angular displacement in the XZ plane, which is considered as a small rotation around the Y-axis. z is the axial displacement of the rotor in the Z direction, and h z is the rotation angle of the gear rotor around the Z-axis. X is the spin speed of the gear rotor. Assume The dynamics equations are expressed as the following through the Lagrange equation: 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]. Figure 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 O p and O g . As shown in Fig. 3, x j , y j and z j are the translational displacements of gear j (j = p,g), and h x , h y and h 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
where d 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: V ¼ sin a cos a 0 0 0 r bp À sin a À cos h 0 0 0 r bg ½ ð22Þ where a is pressure angle, r bp and r bg are the base circle radius of the driving wheel and driven wheel, respectively. Fig. 2 Gear-rotor disc model 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. x 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. c m1 and c 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: here, b is half of the backlash, De p and De g denote thermal profile error of the driving wheel and driven wheel. The direction of the meshing force F m is along the line-of-action, and the direction of the friction force is along the off-line-of-action. When the friction force and friction moment are considered, the force produced by friction effect can be expressed as: V f ¼ À sin a þ a 1 cos a À cos a À a 1 sin a 0 0 0 Àr bp þ b 1 ½ sin a À a 2 cos a cos a þ a 2 sin a 0 0 0 where a i and b i denote the nonlinear 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: where: Fig. 3 Model of gear pair M g ¼ diag m p m p m p J dp J dp J pp m g m g m g J dg J dg J pg ½ ð32Þ

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 h x (s,t), h y (s,t), h z (s,t).
The displacement function can be expressed as: 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: Using Lagrange's equation to process the above formula, we can get the dynamic equations of the shaft element: 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 vibration and torsional vibration of the bearing should be considered. The stiffness matrix K b and damping matrix C b of the supporting bearing are as follows 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-rotorbearing system analyzed in this paper is shown in Fig. 5. Both the input shaft and the output shaft are nine 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 Tables 3, 4 and 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 paper 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 lm and static transmission error amplitude is e 0 = 20 lm. Calculate 800 cycles at each speed, take 200 points in each cycle, and finally take the results of the last 400 cycles for analysis.  The influence of temperature on the gear system can be explained by comparing the dynamic response under different temperatures. Figure 6a and b illustrates the bifurcation diagram of dynamic transmission error (DTE) of gear system under different steadystate temperature fields. The system model considers the temperature-dependent material properties, thermal node load, TVMS and thermal backlash. Above each bifurcation diagram are the phase diagram (blue curve) and Poincaré diagram (red points) at 2850 rpm, 7100 rpm, 9400 rpm, 12700 rpm, 13800 rpm, 14800 rpm and 19100 rpm. At 25°C (Fig. 6a), with the increase of rotating speed (500-16000 rpm), the system undergoes periodic-1, chaotic, periodic-1, chaotic, quasi periodic, chaotic, multi-period, chaotic, periodic-2 and periodic-1 motions. At 225°C (Fig. 6b), the system undergoes periodic-1, chaotic, periodic-1, chaotic, periodic-3, chaotic, multi-period, periodic-2, chaotic, periodic-2 and periodic-1 motions, and its dynamic behavior is more complex. At 2850 rpm under 25°C, the system response is periodic-1 motion. But at the same speed under 225°C , there is a transient chaotic motion. At 9400 rpm, the system motion changes from quasi-periodic (two red circles in Poincaré diagram) to chaos as the temperature increases. At 13800 rpm, the system motion changes from multi-period to chaos as the temperature increases.
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 12,700 rpm, and the system changes from chaotic to periodic-2 motion at 14,800 rpm. To analyze the cause of this phenomenon, the case 225°C 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 b shows 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 14,100-14,900 rpm under 25°C, 13,100-14500 rpm under 225°C with thermal backlash, and 12,200-14,300 rpm under 225°C without thermal backlash. The effects of temperature and backlash on vibration response are further analyzed in Fig. 7. Subplots a, b and c are the time history, phase and Poincaré diagram, and frequency spectrum of the system at 2850 rpm, respectively. When the temperature is 25°C , the Poincaré section have a single point, and the frequency components include one mesh frequency, which indicate the system is in periodic-1 motion. When the temperature is 225°C, there are multiple irregular points in the Poincare section diagram, and there are infinite frequency components in the frequency spectrum, which indicate the system is in chaotic motion. Subplots d, e and f in Fig. 7 are the time history, phase and Poincaré diagram, and   frequency spectrum of the system at 12,700 rpm under 225°C, respectively. When the temperature is 225°C with thermal backlash, the Poincaré section have two points, and the frequency component includes 1/2 mesh frequency, which indicate the system moves in period-2 motion. When the temperature is 225°C without thermal backlash, there are many irregular points in the Poincare section, and the wireless Fig. 7 System vibration responses. a The time history, b phase and Poincaré diagram, and c frequency spectrum at 2850 rpm; d the time history, e phase and Poincaré diagram, and f frequency spectrum at 12,700 rpm frequency component appears in the frequency spectrum, which indicate that the system is in chaotic motion. Although the increase of temperature leads to the unpredictable motion state at some speeds and complex vibration response, the vibration response at some speeds will be restrained due to the reduction of thermal backlash caused by thermal expansion.

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. Figure 8 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 1088 N at 25°C. At 125°C, the force becomes 1855 N and increases by 70.5%. At 225°C, it becomes 3434 N and increases by 215.6%. The peak values in Fig. 8 well corresponds to the chaos phenomenon in the bifurcation diagram. At the speed of 2800 rpm and in the ranges of 6000-7500 rpm, 8100-8600 rpm, 13,100-14600r pm, 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 14400 rpm, the bearing force at 125°C is 3785 N, but the bearing force is 3498 N at 225°C. This may be due to the excessive vibration response under the temperature excitation. Figure 9 shows the vibration displacement amplitude of the driving axle under different bearings and different temperature fields at 19000 rpm. 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. 9 that the vibration displacement increases with the increase of temperature. Node 1 (or node 9) has the largest growth rate. As shown in Fig. 9a, compared with 9.8 lm at 25°C , the amplitude is 81.5 lm at 125°C, increased by 7.32 times, while the amplitude is 172.7 lm at 225°C, 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 lm at 225°C, while bearings B# and C# are both 158 lm.
The causes of the increase of the total bearing force are further analyzed. Figure 10 shows the proportion of the force in each direction when the support bearings are A#, B# and C#, respectively. When the temperature is 25°C (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 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. Figure 11 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°C, 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].