Nonlinear Analysis of Rod Fastened Rotor under Nonuniform Contact Stiffness

Rod fastened rotor is widely used in gas turbine, aero engine, and other occasions. The bending stiﬀness of the contact interface directly aﬀects the stable operation of the rotor. Dynamic model of the rod fastened rotor-bearing system has been established considering nonuniform stiﬀness of interface. The motion equation of this system has been deduced from Lagrange’s equations. The linear dynamic characteristics of this rotor has been investigated, such as Campbell diagram, critical speed, and formation, and the nonlinear characteristics of this system, such as chaos and bifurcation, has been investigated too. The result shows that “bistable state” characteristic appeared on the rod fastened rotor system; that is, there are two critical speeds for each order, and they are all positive precession critical speed, and the amplitude response to the lower critical speed is much larger than that its counterparts to the higher critical speed. In terms of nonlinear characteristics, the rod fastened bearing system has experienced periodic 1 motion, multiple periodic motion, quasi-periodic motion, periodic 1 motion, and chaotic motion successively.


Introduction
e rotor-bearing system is the core component of rotating machines. Its dynamic characteristics affect the stable operation of rotary machinery directly. In order to prevent mechanical unstable operation, it is necessary to study the stability and nonlinear characteristics of rotor-bearing system. Studies on integral rotors are presented in literature [1][2][3][4][5][6]. Ehrich [1] studied the bifurcation phenomenon of Jeffcott rotor-bearing system and investigated the subharmonic vibration phenomenon in the rotor system. Wang [2] established the dynamic model of rub-impact rotor system and studied the periodic response stability of the system by Floquet theory. Gardner et al. [3] analyzed the nonlinear motion of the rotor system under the long bearing after linear instability and studied the subcritical and supercritical bifurcation in the method of multiscale. Bonello et al. [4] studied the nonlinear dynamic response of an extruded oil film damper-rotor system by harmonic balance method. Choi and Naoh [5] established the rotor model of bearing rubbing and analyzed the subharmonic, harmonic, and superharmonic vibration caused by rubbing. Yu et al. [6] studied the friction phenomenon of integral rotor by experimental method and discussed the influence of eccentricity, damping, and other factors on the nonlinear characteristics of this rotor. e researches above focus on the integral rotor-bearing system. Rod fastened rotor has many advantages, such as high strength, light weight, and being easy to assemble and disassemble. e rod fastened rotor-bearing system is widely used in gas turbines and aeroengines (see Figure 1). Different from the integral rotor, the rod fastened rotor is not integral, and the disks at all levels are connected by the rod bolts. erefore, its dynamic characteristics are quite different from the integral rotor. We have the following literature on the research of rod fastened rotor [7][8][9][10][11][12]. Hu et al. [7] considered the initial bending of the rotor when calculating the dynamics of the rod fastened rotor and systematically investigated the influence of rotational speed and initial bending on the dynamic response of the system. Hei et al. [8] analyzed the dynamic characteristics of the fixed tilting pad bearing on the rotor-bearing system and studied the influence of the moment of inertia on the nonlinear characteristics of the system. Li et al. [9] built a model of rod fastened rotor by the three-dimensional finite element method and investigated the nonlinear dynamic characteristics of the rotor by Floquet method. Based on the harmonic balance method, the influence of system parameters on the nonlinearity of the rod fastened rotor was proposed by Da et al. [10]. Wang et al. [11] analyzed the influence of machining errors such as mass eccentricity on the nonlinear characteristics of rotor-bearing system of rod by the three-dimensional finite element method and Floquet theory. Yi et al. [12] proposed a calculation model considering mass effect for the flexible rotor-ball bearing system and analyzed the influence of preload and other factors on the nonlinear characteristics of the rotor subsystem. e stiffness of the contact interface has a great influence on the dynamic characteristics of the rod fastened rotor. Unsuitable or nonuniform stiffness will generate system instability. e stiffness of the rod fastened rotor has been studied as follows [13][14][15][16][17][18]. Based on the theory of double friction, Liu et al. [13] analyzed the influence of preload force, gravity bending moment and other factors on the bending stiffness of the interface and verified this model by modal test. Isa et al. [14] calculated the bending stiffness under the action of gravity bending moment based on bilinear model and verified the simulation by experimental research. Lu et al. [15] calculated the natural frequencies under different preload by using the three-dimensional finite element model and analyzed the relationship between the natural frequencies and the preload. Gao et al. [16] simplified the bending stiffness to the stiffness of bending and twisting spring and analyzed the influence of preload on the bending stiffness of contact interface. Based on Hertz contact theory and Greenwood model, Rao [17] established a contact stiffness model of the rod fastened rotor and verified the model with hammering test. Xia et al. [18] proposed a piecewise bilinear torsional deformation mode and studied the influence of the preloading force and the friction coefficient of the contact interface on the bending stiffness of the contact interface.
In this paper, the dynamic characteristics of rod fastened rotor under nonuniform stiffness have been studied, such as Campbell diagram, critical speed, and formation. e dynamic model of rod fastened rotor-bearing system has been established. e fourth-order Runge-Kutta method is used to get the solution of the nonlinear model. Bifurcation diagram, frequency spectrum, phase trajectory, and Poincare map have been applied to study the nonlinear dynamic phenomena of the rod fastened rotor.

Theoretical Analysis
e rod fastened rotor is complicated; in order to simplify this model, the following assumptions are proposed: (1) ere is no torsional vibration and axial vibration in vibration process (2) e disks at all levels do not separate, in vibration process (3) Bearings at both sides of the rod fastened rotor are identical

Equation of Motion.
As is shown in Figure 2, the rotors are closely connected through circumferential uniform rods, and the rotor is supported by sliding bearings. m b1 and m b2 are the lumped mass of the rotor at the sliding bearing, respectively, while m 1 , m 2 , and m 3 are the lumped mass of the rotor at the interface, where k is the stiffness of the shaft and k 2 and k 2 ′ are the maximum bending stiffness and the minimum bending stiffness of interface of rod fastened rotor, respectively. c 1 is the bearing damping and c 2 is the rotor damping. According to the literature [7], the directions of k 2 and k 2 ′ are perpendicular to each other, so the coordinate system is established with the direction of k 2 as the x-axis and the direction of k 2 ′ as the y-axis. e equation of the rod fastened rotor-bearing system can be deduced by Lagrange equation. Lagrange's equation can be described by where L is Lagrange function and L � V − U, V is the kinetic energy of the system, U is the potential energy of the system, q i and q i are the generalized coordinates and velocities of the system, D is the dissipation energy of the system, and f i is the generalized force in the direction of q i . e plane of the rotor center is zero potential energy reference surface of the gravitational potential energy. e total kinetic energy of the rod fastened rotor system can be expressed as follows: e total potential energy of the rod fastened rotor system can be expressed as follows:

Compressor
Rod Turbine Figure 1: Rod fastened rotor model.
e total dissipation energy of the rod fastened rotor system can be expressed as follows: represents the elastic potential energy of lumped mass of the rotor at the left sliding bearing, represents the elastic potential energy of disk3, and (1/2)k(x 2 b2 + y 2 b2 ) represents the dissipated energy the lumped mass of the rotor at the right sliding bearing.
In (4), represents the dissipated energy of lumped mass of the rotor at the left sliding bearing, represents the dissipated energy of the lumped mass of the rotor at the right sliding bearing. e generalized coordinates of the system are q i � (x 1 , y 1 , x 2 , y 2 ) T , and submitting (2)-(4) into (1), the equation of motion can be deduced as follows: where, F x1 , F y1 , F x2 , and F y2 are the components of the oil film forces on the left and right sides of the bearings in the xand y-directions, respectively. Based on Assumption 2, it can be considered that  Figure 3, where O b is the bearing center, O j is the journal center, θ is the attitude angle, F ε is the radial component of the nonlinear oil film force on the axis diameter, and F θ is the radial component of the nonlinear oil film force on the axis diameter. e oil film force of bearing can be obtained by solving the Reynolds equation.
Reynolds's equations are described by According to literature [19], when the oil film meets the following conditions, the oil is adiabatic, the flow of the oil is laminar, and the oil is incompressible; the oil film force can be expressed as follows: where μ is kinematic viscosity of oil film, L is the bearing length, R is the bearing diameter, ω is angular velocity of bearing, c is the bearing clearance, and ε is the eccentricity ratio.

Bending Stiffness.
Different from the integral rotor, rod fastened rotor can not be seen as a whole, the disks at all levels are connected by preload, and the bending stiffness of the contact interface plays an important role in the stability of the system. ere are two methods to calculate bending stiffness: one is based on the theory of contact mechanics, the Green-Woods model, to calculate the bending stiffness, and another is calculated with finite element method. In this paper, the finite element method has been used to calculate the stiffness. Figure 4 is a model of calculating the bending stiffness of a certain type of rod fastened rotor. e model consists of two disks and rods; there are four kinds of contact in this model: rough contact between the contact interface, no separation contact between the disk and rod, radial contact between two disks, plane contact between disk and rod nut, and boundary conditions and external load is shown in Figure 4. It should be noted that in order to adapt to the geometric characteristics of the disk, the rod nut has been modeled to be in the form of circle. Hexahedral elements have been used for mesh generation of this model. In order to improve the calculation accuracy, more precise grids have been arranged on the rod, as shown in Figure 5. e stiffness of disks is e stiffness of rod is e stiffness of rod fastened rotor is where M is the bending moment applied and α is the deflection angle of the left disk, E is the elastic modulus of the rod, A is the cross-sectional area of the rod, L 1 is the length of the rod, N is the number of rods, and r p is the radius at the center of the rod.

Dimensionless Equation of Motion.
In order to simplify the calculation, the dimensionless transformations are given as follows: According to (11), (5) can be simplified as where F x10 , F y10 , F x20 , and F y20 are nondimensional oil film force of the bearing.

Solving Method.
Considering the nonuniform bending stiffness of the rod fastened rotor, k 2 ≠ k 2 ′ , (12) has strong nonlinearity, and it is difficult to solve such problems by analytical method. Numerical integration methods are often used to solve the nonlinear problem. Because of its high precision, Runge-Kutta method can calculate the influence of system parameters on system response, and it is the main method to solve nonlinear differential equations. e normal form of the Runge-Kutta method is expressed as After determining the order number, the coefficient c i , λ i , μ ij can be determined through Taylor expansion and comparing coefficient.
e Fourth-order Runge-Kutta method is expressed as

Results and Discussion
e parameters of the rotor model have been shown in Table 1. Based on Assumption 2, the bearings on the left and right sides of the rotor are identical. Without considering nonuniform stiffness, the first-to third-order critical speed of this system is 470.63 rad/s, 1223.06 rad/s, and 2022.56 rad/s, respectively, as shown in Figure 6. e nonlinear dynamic behaviors of the rod fastened rotor-bearing system are performed by using the fourth-order Runge-Kutta method and implemented in MATLAB. Bifurcation diagram, time domain diagram, phase trajectory diagram, frequency spectrum, and Poincare map are used to illustrate the nonlinear characteristics. e nonuniform bending stiffness will generate the circumferential anisotropy of the rod fastened rotor in the process of motion. In the analysis of the rotor fastened rotorbearing system, the bending stiffness is an important factor. Figure 6 shows the Campbell diagram of the rod fastened rotor under different degrees of nonuniform. It can be seen from Figure 6 that when the bending stiffness of the contact interface is nonuniform, the rotor appears to be "bistable state" phenomenon; that is, there are two critical speed of each order, and all of them are positive precession critical speed, and the higher critical speed is almost not affected by the nonuniform stiffness, while the smaller critical speed is greatly affected by k 2 ′. Taking k 2 � 5k 2 ′ as an example, the firstto third-order formation of the rotor have been shown in Figure 7. As can be seen from Figure 7, the amplitude response to a smaller critical speed is much larger than its counterpart to a higher critical speed. Figures 8-10 are the bifurcation diagram of dimensionless horizontal displacement X1 at the right bearing location for three different types of stiffness k 2 � 5k 2 ′ , k 2 � 10k 2 ′ , k 2 � 20k 2 ′ . It can be seen from Figure 8 that, at lower rotation speed, ω < 2135 rad/s, and the rod fastened rotor system maintains periodic 1 motion, which is shown as Shock and Vibration 5 a closed loop in Figure 11(b) and one isolated point in Figure 11(d). e system turns into a quasi-periodic motion, at a speed of ω � 2135 rad/s. As shown in Figure 12(b), the phase trajectory presents a lot of closed loops that do not overlap each other, and the Poincare map of the quasiperiodic motion presents a closed loop in Figure 12(d). Experienced a short quasi-periodic motion, 2135 < ω < 2605 rad/s, the system returns into periodic-1 motion again. With the increasing speed, ω > 3745 rad/s, the system response finally presents a chaos motion, which can be seen in Figure 13(d), Poincare map appears as messy, and irregular scatter. Figure 9 is the bifurcation diagram of dimensionless horizontal displacement X1 at the right bearing location for k 2 � 10k 2 ′ . It can be seen from Figure 9 that the bifurcation diagrams of the system have a little difference.

Shock and Vibration
System keeps period 1 motion, at low rotational speed, ω < 1769 rad/s. With the increasing of rotational speed, 1769 < ω < 1987 rad/s, system comes into quasi-periodic motion. 1875 < ω < 1987 rad/s; the system returns to periodic motion in the form of a period 3 inverse bifurcation; as shown in Figure 14(b), the phase trajectory presents three closed loops, and the Poincare map of the quasi-periodic motion presents three isolate points in Figure 14(d). With the increasing of speed, the system enters chaotic motion after a short period 1 motion. Figure 10 is the bifurcation diagram of dimensionless horizontal displacement X 1 at the right bearing location   for k 2 � 20k 2 ′ . e system response becomes more complicated under this condition. Compared with Figures 8  and 9, the quasi-periodic motion happens earlier and the quasi-periodic region narrows with 1027 < ω < 1125 rad/ s. With the increase in the speed, system enters an unstable periodic motion, 1125 < ω < 1193 rad/s. At ω � 1193 rad/s, the system bifurcates again, and the system causes periodic 2 motion under certain conditions. Figure 13 is the response of the system at ω � 1250 rad/s, and the system presents periodic 2 motion, which is shown as two closed loops in Figure 15(b), two obvious frequency components in Figure 12(c), and two isolated points in Figure 15(d). As the rotational speed increases, the system enters chaotic motion without jump phenomenon.

Conclusion
Based on Lagrange's equation, the model of rod fastened rotor-bearing system under the condition of nonuniform bending stiffness of interface has been established in this paper. e linear dynamics behaviors have been explored, such as Campbell diagram and critical speed. e nonlinear dynamics behaviors have been studied by using the fourthorder Runge-Kutta method. e bifurcation diagram, vibration waveform, spectrum, phase trajectory, and Poincare map are given to illustrate the nonlinear dynamic phenomena of the system. e following conclusions can be drawn from the above research.
(1) With the increase of the rotating speed, the rod fastened rotor system has experienced periodic 1 motion, multiple periodic motion, and quasi-periodic motion successively. en, the system returns to periodic 1 motion by an inverted bifurcation. With the speed increasing, the system enters chaotic motion in the end. (2) With the increasing of k 2 /k 2 ′ , the critical value of the system entering the chaotic phases is smaller, and the interval of quasi-periodic motion gets smaller. (3) Because of the difference between k 2 and k 2 ′ , the rod fastened rotor system appears to be "bistable state" phenomenon. ere are two critical speeds of each order, and both of them are positive precession critical speed. e smaller critical speed is greatly influenced by k 2 ′ , and the amplitude response to the smaller critical speed is much larger than the counterparts to the higher critical speed.
is study can provide guidance for the failure caused by nonuniform bending stiffness of the interface of the rod fastened rotor, and, at the same time, it is helpful to further study the nonlinear dynamic characteristics of the rod fastened rotor. e further research work will focus on the modeling of the stability of the rod fastened rotor caused by the nonuniform stiffness.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.