Nonlinear Dynamic Characteristics of Rod Fastening Rotor with Preload Relaxation

Compared to ordinary rotor, rod fastening rotor has the advantages of lighter weight, higher strength and easier installation, so it is widely used in gas turbine. However, in the process of a long-term operation, the rod may be deformed due to the influence of alternating load, high temperature and other uncertain factors. In serious cases, it can even lead to a major accident. The discontinuous characteristic of rod fastening rotor leads to great differences in dynamic characteristics compared to ordinary rotor. Based on the Herz contact theory and GW contact model, the contact effect between two discs was studied, and the relationship among the contact load, the distance between two disks and the equivalent bending stiffness was obtained. Findings show the bending stiffness to decrease nonlinearly with the decrease in contact load. The lumped mass method was used to establish the rotor model. The contact effect was considered and the Runge–Kutta method was used to solve the model. Combined with the bifurcation diagram, time domain diagram and spectrum diagram, the influence of contact stiffness on rotor dynamic characteristics was analyzed. The results show that the dynamic characteristics of the rod fastening rotor are rich due to the influence of nonlinear factors. In the case of uniform relaxation, the contact stiffness has different effects on the response state and frequency doubling amplitude of the system at different speeds, which is mainly related to the motion state of the system. In the case of non-uniform relaxation, the degree of relaxation does not affect the motion state of the system, but only changes the amplitude of vibration. The results provide theoretical support for condition monitoring and fault diagnosis of rod fastening rotor.


Introduction
Gas turbine is widely used in power generation, aviation industry and other fields because of its high efficiency and high reliability. Rod fastening rotor has reliable strength and convenient installation, and has strength characteristics close to the integral rotor, so it is widely used in heavy-duty gas turbine. However, during a long-term operation, the tie rod will produce fatigue deformation or even fracture due to the action of alternating load, and the tie rod is prone to creep at high temperature. These problems will lead to the looseness of the tie rod and the preload between the discs will be reduced, resulting in the change of the dynamic characteristics of the rotor system, and there are also serious hidden dangers of failure. Therefore, it is of great engineering significance to study the rotor dynamic characteristics of rod fastening rotor for gas turbine design and manufacturing, safety monitoring and fault diagnosis.
Many scholars have made outstanding contributions to the design and dynamics of gas turbine rotor [1][2][3][4][5][6]. Rusin et al. [1] analyzed the stress state and performance of gas turbine rotor under unstable operating conditions. In order to solve the problem of gas turbine vibration, Thompson et al. [2] developed an analytical model to calculate engine rotor dynamics and structural response. Ehrich [3] described some of the nonlinear problems of gas turbine rotor. Amin [4] established a nonlinear model of a gas turbine rotor, fitted the nonlinear stiffness of the bearing, and studied the nonlinear characteristics of the rotor.
The difference between rod fastening rotor and integral rotor lies in the contact characteristics between discs, thus, there are many studies on the contact mechanism of contact plane [7][8][9][10][11][12]. Brizmer et al. [7] theoretically studied the elastic-plastic contact behavior between a deformable sphere and a rigid plane under full contact conditions. Kucharski et al. [8] studied the contact problem of rough and rigid planes by theoretical and experimental methods. Zhao et al. [11] introduced a virtual material layer at the contact surface of the disks to characterize the contact effect. Cui et al. [12] established a rough surface model between the discs of rod fastening rotor by using the autocorrelation function method, generated the geometric model by finite element software and carried out contact simulation calculation, then obtained the relationship between preload and contact stiffness of the rough unit area.
Some scholars [13][14][15][16][17][18] have studied the contact characteristics of a rod fastening rotor; mainly focusing their research on the stiffness characteristics and influence of the preload, roughness and other factors on the contact characteristics. For example, Li et al. [13] carried out finite element contact analysis on a rectangular micro element with a rough surface based on elastic-plastic theory, and obtained the normal and tangential interface of the contact stiffness under different loads according to the relationship between force and deformation. The tie rod rotor connection is similar to the flange structure; some scholars have carried out research on the modeling of flange structure. Berne et al. [19] modeled the stiffness characteristics of a flanged rotor joint and found that its moment stiffness is a nonlinear function of bending moment. The results show that the value of moment stiffness basically depends on the bolt tightening force, axial force and other factors. Schwingshackl et al. [20] found that the ability to model a flange depends primarily on two things: load distribution on the joint, the number and distribution of nonlinear elements.
Other scholars have studied the dynamic characteristics of rod fastening rotor. Cheng et al. [21] established the dynamics model of rod fastening rotor with squeeze film damper, studied the influence of system parameters on the rotor bistable characteristics and found that the rotor has a frequency doubling component in the bistable region. Liu et al. [22] established the dynamic model of a rod fastening rotor-bearing system, studied the system stability using Floquet theory and compared it to ordinary rotor. He et al. [23] proposed a new method to determine the normal stiffness of the contact surface. It is according to the finite element analysis results of the force and deformation of the micro body element on the elastic-plastic rough surface, which improved the calculation accuracy. Hu et al. [24][25][26] studied the nonlinear dynamic response of a two-disk rod fastening rotor with transverse crack, initial bending, rub impact and other faults, and analyzed the influence of system parameters such as speed, crack size, crack angle, contact surface damping coefficient, initial bending and rub impact stiffness on the dynamic characteristics. Li et al. [27] studied the nonlinear dynamic characteristics of a three-disk bending torsional coupling rod fastening rotor under rub impact fault, concluded that the influence of torsional vibration on bending vibration gradually increases with the increase in rotating speed. Wang et al. [28,29] studied the dynamic characteristics of rod fastening rotor under the coupling fault state, which are the non-uniform contact stressors of the disks, crack, initial bending and others, and carried out experimental verification. Hei et al. [30,31] studied the nonlinear dynamic characteristics of coupling faults, which are the rod fastening rotor with looseness of support and rub impact. It was found that the rotor behavior in the fault state was more complex, and the bifurcation point moved backward. Da et al. [32] compared the dynamic characteristics of the three-dimensional and bilinear nonlinear stiffness models of the rod fastening rotor, and found that the reasons for the sudden change of the amplitude of the two models were different, and the dynamic characteristics near the unstable response region were different. Xu et al. [33,34] deduced the vibration equation of rod fastening rotor system considering various factors of connecting rod in detail. Based on the fractal contact theory, the equivalent bending stiffness was calculated, and the effects of the preload, tie rod detuning rate and the number of detuning tie rods on the natural frequency were studied.
When the tie rod is loose, the preload decreases, and the contact characteristics between rotor disks change, which will affect the vibration response of rotor, which has hidden trouble. Most of the research on rod fastening rotor dynamics only considers the contact effect between rotor disks when building models; there is little research on how the contact effect affects the dynamic characteristics. In this paper, the lumped mass method is used to establish the two-disk rod fastening rotor dynamic model. The influence of preload between discs on contact stiffness is studied based on fractal contact theory, and the influence of preload relaxation on the dynamic characteristics of rotor is also studied. This study provides a theoretical basis for the design and the fault diagnosis of the rod fastening rotor.

Structure
The structure of a two-disk rod fastening rotor is shown in Figure 1. The two disks are fastened and connected together through six tie rods evenly distributed along the circumference. The disks are connected to the shaft, and both ends of the elastic shafts are supported by two sliding bearings.

Oil Film Force of Sliding Bearing
When the rotor is rotating, the journal of sliding bearing will be affected by oil film force, and the force on the journal is shown in Figure 2. The dimensionless nonlinear oil film force on the journal can be obtained by solving the Reynolds equation [35]: The parameters in Equation (1) are in Appendix A.

Bending Stiffness of Disc Contact Layer
As shown in Figure 3, the contact between the discs is equivalent to a structure of uniformly distributed spring and hinge, which can be further equivalent to a torsion spring [36]. The function of the hinge is to limit the linear displacement of the two discs, and the main function of the uniformly distributed spring is to provide the bending stiffness of the contact layer and limit the angular displacement of the two disks.

Contact layer
Uniform spring Hinge Torsion spring When the two discs are pressed together by the preload of tie rods, the first contact part is the micro convex body on the rough surface. With the increase in the preload, the number of micro convex bodies in contact increases, and the micro convex bodies produce elastic or even plastic deformation. In order to study the contact stiffness characteristics of a rough surface, the shape of micro convex body is usually regularized, and the spherical micro convex body model is conducive to analyze the contact characteristics of rough surface from the perspective of mechanics. According to Hertz contact theory, the contact between two elastic spheres can be simplified as the contact between an elastic sphere and a smooth rigid plane. According to the Persson contact theory, the height of micro convex body on the plane obeys Gaussian distribution. It can be concluded that the total load of the contact layer is: and the bending stiffness is: The detailed derivation of Equations (2) and (3) are in Appendix B. The surface contact model parameters [37] are as shown in the Table 1. The relationship between the contact load and the distance between the two disks is shown in Figure 4, and the relationship between the contact load and the bending stiffness is shown in Figure 5. The load W has a nonlinear relationship with the distance h. With the increase in h, the contact load W decreases, and the rate of reduction decreases. The relationship between bending stiffness and contact load is also weakly nonlinear, with the decrease in preload load, the bending stiffness decreases.

Dynamic Model
The physical model is simplified and equivalent to a lumped mass model composed of concentrated mass at the bearings, disks and massless elastic axis. The translational degrees of freedom in x and y directions and rotational degrees of freedom around x and y axes of the concentrated masses are considered, and the axial vibration and torsional vibration of the rotor are ignored. Based on the principle of Da Lambert, the dynamic model of rod fastening rotor, which has 16 degrees of freedom, can be obtained.
The parameters in Equation (4) are in Appendix C.
The system parameters in the dynamic model are shown in Table 2

Calculation Results
The fourth-order Runge-Kutta method and simulation software MATLAB are used to solve the dynamic equations. According to the above analysis of the contact effect between discs, the flexural stiffness Kc mainly affects the displacement angle of the disc θ and ψ, which are the displacement angle of the disc rotating around the x axis and y axis. Therefore, the following analysis is about these two angular displacements. When the preload is in a saturated state (50 kN), the dynamic equations are solved. Figure 6 shows the bifurcation diagram about the angle displacement of disk 1, expressed as θ1, varying with rotating speed ω. The speed range is 300~2000 rad/s. Figure 7 shows the time domain diagrams and spectrum diagrams about angular displacement θ1 at several speeds. The abscissa of the spectrum is the ratio of frequency f to rotating frequency fr. When the rotating speed is at 500 rad/s, the time domain curve of the angle displacement θ1 is sinusoidal, and the spectrum diagram has only a component at rotating frequency, which is mainly caused by mass unbalance, and the corresponding bifurcation diagram is in a period one motion state. At 800 rad/s, the time domain curve is no longer sinusoidal but is still periodic. Affected by the nonlinear oil film force, the spectrum has a component at half of rotating frequency, and this component is greater than the component at rotating frequency, and the corresponding bifurcation diagram is in period two motion state. Compared to 800 rad/s, the frequency component of 1100 rad/s at half of rotating frequency is greatly increased, and that at rotating speed is basically unchanged, the corresponding bifurcation diagram is thus in a quasi-periodic or chaotic motion state. At 1800 rad/s, the chaotic state of the system intensifies, the time domain curve is irregular, and there are rich frequency components in the spectrum diagram, such as half, one time and two times rotating frequency, and there are more sideband components around these components.

Influence of Uniform Relaxation of Preload on Dynamic Characteristics
Uniform relaxation means that the overall preload of the disks decreases but the tension of each tie rod is equal. The uniform relaxation of the preload will reduce the contact stiffness Kc, which will affect the rotor dynamic characteristics of the system. Figure 8 shows the bifurcation diagrams about angle displacement θ1 with the change of contact stiffness Kc, and the range of stiffness is 1 × 10 6~5 × 10 7 N/rad. When the rotating speed is at 500 rad/s, the bifurcation diagram is a horizontal straight line, and there are two horizontal straight lines at 800 rad/s, indicating that the system state is basically not affected by the change of stiffness at these two speeds. Corresponding to Figure 6, the system has been in period one and period two motion state. At 1100 rad/s, the contact stiffness Kc is between 1.1 × 10 7 N/rad and 5 × 10 7 N/rad, displacement response is affected by the bending stiffness. With the increase in stiffness, the displacement amplitude first increases and then decreases. The system state in other stiffness ranges is not affected by the change of stiffness. When the rotating speed is at 1800 rad/s, the system is in a strong chaotic state with the change of stiffness.  Figure 9 shows spectrum waterfall diagrams about the polarization angle displacement of disk 1 with the contact stiffness Kc at the above four speeds. When the rotating speed is at 500 rad/s, with the change of stiffness Kc, the spectrum diagram mainly has a component at rotating frequency, the component at two times frequency is very small, and the amplitude is not affected by the change of stiffness. At 800 rad/s, the amplitude is still not affected by the change of stiffness, but there is a larger component at half of the rotating frequency. This is because when the speed reaches a certain degree, the component at half of the rotating frequency generated by the oil film whirl of sliding bearing gradually occupies the main component. At 1100 rad/s, with the contact stiffness increases, the component at half of the rotating frequency first increases and then stabilizes, and the component at rotating frequency remains unchanged. At 1800 rad/s, there are more frequency components and sideband components because the system is in the strong chaotic state.
The amplitudes of components at different frequencies shown in Figure 9 can be extracted to obtain Figure 10. They are amplitudes at half, one time, two times and three times rotating frequency, and change with the varying contact stiffness. When the rotating speed is at 500 rad/s and 800 rad/s, the amplitude is not affected by the change of stiffness. At 1100 rad/s, the amplitudes at half, one time and three times rotating frequency first increase and then decrease with the increase in stiffness, and get the maximum value when Kc is about 1.5 × 10 7 N/rad. The amplitude at two times rotating frequency is not affected by the change of stiffness. At 1800 rad/s, the amplitude at each frequency fluctuates strongly, and the fluctuation period is irregular.

Additional Moment Model
Non-uniform relaxation means that the tension of all tie rods is not exactly the same. When the preload relaxation degree of each tie rod of the rod fastening rotor system is different, it will lead to uneven distribution of the preload between the discs, resulting in an additional bending moment at both ends of the disks [38]. The magnitude and direction of the additional bending moment are determined by the uneven distribution of the preload. Suppose this moment is Ms, as shown in Figure 11, the moment Ms can be decomposed in the rotating coordinate system ζO'η, the speed in this coordinate system is the rotation speed of the rotor. Mζ and Mη are the decomposition amount of Ms in the rotating coordinate system. Define s as the tie rod looseness rate, which is the degree of non-uniform relaxation of the preload. Suppose only one tie rod is loose, the calculation formula is as follows: Convert the additional bending moment in the rotating coordinate system into the additional bending moment in the fixed coordinate system. According to the angle conversion relationship: The additional bending moment should be added to the external force matrix. After that, the external force matrix F in Equation (4) is in Appendix C. Figure 12 shows the bifurcation diagrams about the polarization angle displacement of disk 1 varying with the tie rod looseness rate s at different speeds. When the rotating speed is at 500 rad/s, the system is in period one motion state, the bifurcation diagram is approximately a straight line. At 800 rad/s, the system is in period two motion state, the bifurcation diagram is approximately two parallel lines. At 1100 rad/s and 1800 rad/s, the system is in chaotic motion state, the degree of chaos is more serious at 1800 rad/s. There are also some common qualities. Regardless of the speed and the state of the system, with the increase in s, the vibration angular displacement always offsets in negative direction. This is due to the increase in s, resulting in the increase in the additional bending moment, which will bring a negative constant to the external force term of the system. In addition, with the increase in s, the motion state of the system will not change.  Figure 13 is a waterfall diagram about the polarization angle displacement of disk 1 varying with the tie rod looseness rate s at different speeds. When the rotating speed is at 500 rad/s, the spectrum diagram has components at one time and two times. At 800 rad/s and 1100 rad/s, the spectrum diagram has an additional component at half rotating frequency. At 1800 rad/s, in addition to the above components, the spectrum diagram has component at three times the rotating frequency and side-band components.

Influence of Tie Rod Looseness Rate on Dynamic Characteristics
The amplitudes of components at different frequencies shown in Figure 13 can be extracted to obtain Figure 14. When the rotating speed is at 500 rad/s, 800 rad/s and 1100 rad/s, the variation trend of spectrum amplitude at the four frequencies is the same. The difference is that when the rotating speed is at 500 rad/s, there are spectrums only at one time and two times rotating frequencies. And when the rotating speed is at 1100 rad/s, the system is in a slightly chaotic motion state, and the amplitude fluctuates with s. At 1800 rad/s, the system is in a strongly chaotic motion state, so the amplitude fluctuates sharply with s.  =500 rad/s =800 rad/s =1100 rad/s =1800 rad/s

Conclusions
Firstly, the dynamic model of the rod fastening rotor system is established. Based on the fractal contact theory of rough surface, the relationship between the contact stiffness of the two disks and the contact load was studied from the mechanism, and the influence of the contact stiffness on the displacement bifurcation characteristics and spectrum components of the rod fastening rotor system was explored. Then, the effects of uniform relaxation and non-uniform relaxation of preload on rotor system dynamics were studied. The following conclusions can be drawn from this study: (1) The contact load of the rod fastening rotor has a nonlinear relationship with the distance between the two disks. With the increase in distance, the preload between the disks decreases, and the decreasing rate decreases gradually. The relationship between contact stiffness and contact load is weakly nonlinear. With the decrease in the preload load, the contact stiffness decreases.
(2) Affected by nonlinear factors, the rod fastening rotor system has rich dynamic characteristics at different speeds. When the rotation speed is small (under 760 rad/s), it is in a period one motion state, the spectrum diagram has a main component at rotating frequency. As the rotational speed increases (760~1600 rad/s), the system will appear in a period two motion state, and the spectrum diagram has components at rotating frequency and half rotating frequency. At one of the rotational speeds (1000~1200 rad/s), the system is in a weak chaotic state, and the component at half rotating frequency increases. At higher rotational speeds (1600~2000 rad/s), the system is in a chaotic state, many sideband components appear in the frequency spectrum, and the time-domain vibration curve is no longer periodic.
(3) When the preload of the discs is uniformly relaxed, the influence of the contact stiffness of two discs on dynamic characteristics of the rod fastening rotor system is different at different speeds, which is mainly related to the motion state of the system. When the system is in the stable periodic motion state, such as period one (500 rad/s) and period two (800 rad/s), the vibration displacement, spectrum amplitude and other characteristics are basically not affected by the change of contact stiffness. When the system is in a quasiperiodic or chaotic motion state (1100 rad/s), the system responses show a certain variation law with the change of contact stiffness. When the system is in a strong chaotic state (1800 rad/s), the system responses fluctuate strongly with the change of contact stiffness, and the fluctuation period is not regular.
(4) When the preload of the discs is non-uniformly relaxed and the system is in a stable periodic motion state (500 rad/s and 800 rad/s), along with the increase in the tie rod looseness rate s, the amplitude increases linearly at two times rotating frequency, and remain unchanged at the other frequencies. When the system is in a chaotic motion state (1100 rad/s, 1800 rad/s), the amplitude fluctuates with s. The degree of fluctuation depends on the degree of chaos.
This study provides theoretical support for the design and manufacture of rod fastening rotor and its state monitoring and fault diagnosis. Due to the strong nonlinearity of the system limited by test conditions, it is difficult to carry out test verification. An experimental study will be carried out in the next stage of our research.

Conflicts of Interest:
The authors declare that there is no conflict of interest.

Ae
Actual contact area between elastomer and rigid plane  Figure A1. Schematic diagram of contact between elastic sphere plane and rigid plane.
According to the elastic mechanics analysis, the deformation δ, the radius a and the contact load of single micro convex body We are: The actual contact area between the elastic sphere and the rigid plane is: According to GW contact theory, the height of micro convex body on the machining plane obeys Gaussian distribution. Figure A2 is a micro schematic diagram of the contact of the two planes. h is defined as the equivalent spacing between the two planes. Only the micro convex body with height z > h on the elastic plane will contact the rigid plane. The probability of contact between two planes is calculated from the probability distribution:  According to the above formula, the actual contact area of the two planes is: According to the above formula, the relationship between the load Pc and h is nonlinear. If the contact load between two planes is known, the distance h between two planes can be calculated. The load and distance will be equivalent to the contact action between two contact planes to a distributed spring with an elastic coefficient of kc.