Nonlinear Dynamic Characteristics of Multidisk Rod Fastening Rotor with Axial Unbalance Mass Distribution

During long-term operation, the blades of rod fastening rotor are very prone to fault and may cause uneven axial mass distribution. +e dynamic characteristics of a multidisk rod fastening rotor with axial unbalance mass distribution are studied. A four-disk rod fastening rotor model was established by the lumped mass method. +e mass unbalance on different disks was used to simulate blade faults in different parts, and nonlinear factors such as oil film force and disk contact effect were considered. +e fourth-order Runge-Kutta method was used to solve the model, and the spectrum and axis trajectory curves were obtained. +e influence of fault location and rotating speed on rotor dynamic characteristics was analysed when there were blade faults on two disks.+e results show that a large, unbalanced mass will bring strong nonlinear characteristics.+e rod fastening rotor may be in different motion states at different rotating speeds. When the fault position of two disks changes, the changes of bearing spectrum, axis trajectory, and phase difference have regularity. +e results can provide a theoretical basis for condition monitoring and fault diagnosis of rod fastening rotor.


Introduction
Blades are attached to the rotors of many mechanical equipment using multi-disk rod fastening rotor. Due to the large changes of temperature and pressure in airway, longterm operation has an impact on blade health. Statistics show that blade scaling, corrosion, wear, damage, and fracture account for a large proportion [1]. Blade fault may affect the efficiency of machinery at least, or it may develop into blade failure accident and cause serious secondary damage to the equipment. erefore, it is of great engineering significance to study the fault mechanism and dynamic characteristics of a multidisk rod fastening rotor.
Blade scaling, corrosion, and other faults cause the change of shafting mass distribution, resulting in the change of shafting vibration response characteristics. Gao et al. [2] studied the mechanism of vibration and impact of the rotor system when the eccentricity excitation of the two-stage centrifugal pump impeller changes suddenly. e results show that the sudden change of eccentricity is very destructive to the system. Chen et al. [3] studied the influence of unbalanced load caused by blade shedding on the dynamic characteristics of rotor system based on the finite element method and analysed the influence of rotating speed, unbalance, disk eccentricity, and support stiffness on the unbalanced vibration response of the system. Yang et al. [4] considered the dynamic behaviour of nonlinear rotor systems under two coupling faults: rotor misalignment fault and unbalance fault. Taghipour et al. [5] used three different rotor systems to study the nonlinear dynamics of a flexible rotor shaft with an unbalanced disk. Hong et al. [6] studied the lateral-torsional coupling vibration of the rotor system with large unbalance. e modal characteristics of an unbalanced rotor system and the relationship between modal and response of unbalanced rotor are obtained. Wang and Jiang [7] deduced the dynamic model of aeroengine rotor system in case of unbalanced misalignment coupling fault. e effects of different rotating speeds, mass eccentricity, angle eccentricity, and parallel eccentricity on complex vibration response are analysed by time domain and spectrum, and the accuracy of the model is verified by a double-rotor test rig. e research on the rod fastening rotor mainly focuses on the contact characteristics and fault dynamics. Xu et al. [8] proposed a new contact model based on Persson contact theory and optimized the solution method of contact stiffness of distributed rod fastening rotor. Miao et al. [9] studied the dynamic modelling and correction method of contact interface of rod fastening rotor based on thin-layer element with linear constitutive relationship and verified the model through experiment. Hu et al. [10][11][12] studied the nonlinear dynamics of two-disk rod fastening rotor with nonlinear contact stiffness under rub impact, crack, and initial bending faults. Li et al. [13] studied the nonlinear dynamic characteristics of bending-torsion coupling of three-disk rod fastening rotor under rub impact fault. e results show that the greater the rotating speed, the greater the influence of rub impact force on the rotor. Wang et al. [14,15] studied the influence of coupling faults such as uneven contact stress distribution between two disk and initial bending of the system due to the existence of cracks on rotor dynamic characteristics and verified it through experiment. Cui et al. [16] used finite element method to study the dynamic characteristics of rod fastening rotor under different preload and bearing assembly clearance. e results show that the critical speed of the rotor system increases nonlinearly with the increase of preload. When the rotor passes the critical speed, the amplitude of unbalanced response increases with the increase of preload and bearing assembly clearance. Wu et al. [17] established the mechanical model of circumferential distributed rod and the nonlinear contact stiffness matrix characterizing the contact effect between disks. e effects of preload, rod diameter, and contact surface roughness are analysed. Hei et al. [18] studied the dynamic characteristics of rod fastening rotor supported by finite length sliding bearing. An improved Newmark method is proposed. On this basis, the dynamic characteristics of rod fastening rotor under different speeds, disk eccentricity, bearing stiffness, and contact stiffness are analysed. e contact surface between the disks of multidisk rod fastening rotor has nonlinear stiffness characteristics. When the blades have scaling, corrosion, and other faults, the mass distribution along the axial direction may change, which may produce complex nonlinear vibration response characteristics. At present, there are few simulation studies on this kind of dynamic response characteristics of multidisk rod fastening rotor structure, and the number of disks considered is small, which cannot be in good agreement with the reality. In this paper, a multistage rotor is simplified into a four-disk rod fastening rotor structure, and its dynamic model is established by using the lumped mass method. At the same time, considering the nonlinear characteristics of oil film bearing stiffness and disk contact stiffness, the distributed mass imbalance is used to simulate blade faults on different disks, and the nonlinear dynamic response characteristics of the rotor are calculated and analysed.

Rotor Dynamics Modeling of Multidisk
Tie Rod 2.1. Rotor Structure. e structural diagram of multidisk rod fastening rotor system composed of four disks is shown in Figure 1. Four disks are fixed and connected together through circumferentially evenly distributed tie rods, and both ends of the rotor are supported by sliding bearings.

Oil Film Force of Sliding Bearing.
e nonlinear oil film force of bearing can be obtained by solving Reynolds equation. According to it, the dimensionless nonlinear oil film force distribution can be expressed as [19] e following is obtained: 2 Shock and Vibration (2)

Stiffness of Disk Contact Layer.
In order to facilitate, generally, there is a gap between the tie rod and the mating hole. When the relative displacement of the disk in x direction is less than the gap, the tie rod and the mating hole do not contact. Under the preload of the tie rod, the tangential stiffness of the friction shear layer between the disks is k 1 . When the relative displacement of the disks is larger than the clearance, the tie rod contacts the mating hole and produces an additional tangential stiffness k 2 , which is a piecewise linear function. e restoring force of the disk contact layer in the x direction generated by the additional tangential stiffness can be obtained from (3). Similarly, the restoring force in the y direction can be obtained: When ε is a small quantity, the restoring force in the x direction due to the additional stiffness can be expressed as [20] F � k 1 x 12 + k 2 x 3 12 . Figure 2 shows the schematic diagram of rotor fault location, O i is the centroid of each disk, and C i is the center of gravity of each disk (I � 1, 2, 3, 4). Simulate two states of rotor as follows.

Fault Simulation.
(1) Normal state: add initial unbalance eccentricities e i of 1e-5 m on each disk to make each disk meet G1 balance accuracy level at 1000 rpm. e included angle between the eccentricity direction and the x direction are 4 random numbers (227°, 35°, 100°, and 197°).
(2) Fault status: unbalance masses on different disks are used to characterize the fault such as scaling, corrosion, wear, damage, and fracture of the blade. e fault unbalance eccentricity e b1 and e b2 with the size of 1 × 10 − 4 m is added to disks 2 and 3. e total unbalance eccentricity e f1 and e f2 of disks 2 and 3 are the vector sum of initial unbalance eccentricity and x y z Contact layer Bearing 1 Bearing 2 Tie rod e included angle between the two fault unbalance eccentricities is α.

Motion Equation and Parameter of Pull Rod Rotor.
Take disk 1 as an example for force analysis, which is shown in Figure 3, where C 1 is the center of gravity of the disk and O 1 is the centroid of the rotor.
According to the theorem of mass center motion and moment of momentum, the motion differential equations of six lump masses and 12 degrees of freedom of the rod fastening rotor system can be obtained. Since the rotor is symmetrical structure, only the motion equations of bearing 1 and disks 1 and 2 are listed as follows: e parameters of the rod fastening rotor system are shown in Table 1.
ere are many ways to solve such a system of equations, such as Runge Kutta method, harmonic balance method, and HB-AFT method [21]. e fourthorder Runge Kutta method is used to solve the motion equation of the rotor system, and the vibration responses of different rotating speeds and different parts of the rotor under normal and fault conditions are calculated, respectively.
e rotor vibration at bearing 1 is analysed as an example in the following.  the rotor has been in single period motion state in the analysed speed range, while, under fault conditions, the system shows strong nonlinear characteristics, and the motion states are more complex and diverse.

Influence of Rotating Speed on Rotor Dynamic Characteristics
In order to analyze the motion state of the system at each speed, Figure 5 shows the phase trajectory and Poincare diagram under normal state and fault state at speeds of 2000 rpm, 5000 rpm, 8000 rpm, and 18000 rpm, igure 3: Stress analysis diagram of disk.  Shock and Vibration 5 corresponding to several motion states in the bifurcation diagram. Under normal conditions, at the speeds of 2000 rpm, 5000 rpm, and 8000 rpm, the phase trajectory diagram is circular or oval, and the Poincare diagram has only one point, and it means the system is in single period motion state. At 18000 rpm, the phase trajectory and Poincare diagram are not completely coincident. Although there is still one value corresponding to 18000 rpm on the bifurcation diagram, the system has changed to quasiperiodic motion state. erefore, in the normal state, the system has been in a single period motion state at low speed. With the increase of speed, the system will gradually enter quasiperiodic motion or even chaotic state under the influence of nonlinear oil film force. Under fault state, when the rotating speed is 2000 rpm, the system is in the single period motion state, like normal state. At 5000 rpm, the phase trajectory is concave inward, there are two points in the diagram, and the system is in double motion state. At 8000 rpm, the phase trajectory diagram does not coincide, there are multiple points in the Poincare diagram, and they form a circle, indicating that the system is in quasiperiodic motion state. At this time, the oil film whirl has a great impact on the system. At 18000 rpm, the phase trajectory diagram is chaotic and irregular, and the points of Poincare diagram are more open and can no longer form a closed ring, which means that the system is in chaotic state. e oil film force may have developed or gradually developed into oil film oscillation. Figure 6 shows the spectrum diagram under normal state and fault state at speeds of 2000 rpm, 5000 rpm, 8000 rpm, and 18000 rpm. Under normal state, because the system is in single period motion state at these speeds, the spectrum has only one-time rotating frequency component at the four rotational speeds. Under fault state, for similar reasons as previously mentioned, the spectrum has mainly one-time rotating frequency component at 2000 rpm. At 5000 rpm, the system is in double motion state, in addition to one-time rotating frequency component, there is obviously half rotating frequency component. At 8000 rpm, the system is in quasiperiodic motion state, many nonoctave components appear in the spectrum. At 18000 rpm, the system is in chaotic motion state, and many sideband components appear. Figures 7 and 8 are the vibration time domain diagram and axis trajectory diagram under normal state and fault state at speeds of 2000 rpm, 5000 rpm, 8000 rpm, and 18000 rpm, corresponding to several motion states in the bifurcation diagram. Under normal conditions, the time domain diagrams are sine curve, the axis trajectory is small and the shape is circular or oval. Under fault state, when the rotating speed is 2000 rpm, the system is in the single cycle motion state, and the time domain diagram is also sine curve. However, the amplitude is much greater than the normal state, and the axis trajectory is a quasiellipse; At 5000 rpm, the system is in II motion period, and the axis trajectory is like a small ellipse in a large ellipse. At 8000 rpm and 18000 rpm, the system is in quasiperiodic or chaotic motion period. e time domain diagram shows complex waveform, and the shape of axis trajectory is complex and confused, which is no longer stable to one shape.
In Figure 9, the spectrum waterfall diagram in normal state and blade fault state is shown. Under normal conditions, there are mainly one-time rotating frequency components in the spectrum, and two-time rotating frequency components are too small to see clearly. Compared with the normal state, the spectrum components in the fault state are relatively rich, the frequency doubling components increase significantly, but two-time rotating frequency components are also weak. e frequency spectrum of 0-7000 rpm and 18000-20000 rpm have half rotating frequency components, which are nonlinear components caused by oil film whirl. e frequency spectrum of 15000-20000 rpm has some components and sideband composition less than half rotating frequency, and the amplitudes increase with the increase of rotating speed, even bigger than of one-time frequency. e sliding bearing gradually changes from oil film whirl state to oil film oscillation state.    Shock and Vibration 7 In Figure 10, the vibration response curve of disk 1 under normal state and fault state is shown; the speed range is 0-20000 rpm. Under normal state, the disk has three order critical speed, which is 2000 rpm, 7000 rpm, and 12000 rpm, respectively. In the blade fault state, compared with the normal state, the critical speed increases to 3400 rpm, 7500 rpm, and 13000 rpm, respectively, and the amplitude of each critical speed increases greatly. By comparing this  Shock and Vibration figure with Figure 4, under normal state, the system has been in single period motion state, the response curve is very smooth. But under fault state, the response fluctuates around 8000 rpm and after 15000 rpm. e system is in unstable quasiperiodic or chaotic motion state at these speeds. It can be seen that the motion state of the system has a great influence on the shape of the response curve.

Influence of Blade Fault Location on Rotor Dynamic Characteristics
Select three rotating speeds of 2000, 5000, and 8000 rpm and change the included angle of fault unbalance eccentricity of disks 1 and 2 to study the influence of blade fault position on the dynamic characteristics of rod fastening rotor. Figure 11 shows the Bifurcation diagram of α. It can be seen that different blade fault location angles have obvious effects on the motion state of the system. When the rotating speed is 2000 rpm, the motion state of the system does not change with α and is always in single period motion state. At 5000 rpm, as α changes, the system may be in single period, double period, quasiperiodic, or chaotic motion state. And the intervals of each motion state are roughly symmetrically distributed. At 8000 rpm, the system is in single period motion state in most intervals, in quasiperiodic or chaotic motion state in the range near 0°and 360°. Figure 12 shows the spectrum waterfall diagram varying with α and the amplitude of each frequency doubling component. When the rotating speed is 2000 rpm, there are only one-time and two-time rotating frequency components in the spectrum. With α gradually increasing, the amplitudes of one-time and two-time rotating frequency components first decrease and then increase, both amplitudes reach the minimum value when α is 180°, and the minimum value is 0. When the rotating speed is 5000 rpm, there are obvious components of half rotating frequency, one and half rotating frequency, one-time rotating frequency, and two-time rotating frequency, and the frequency components are rich and complex. With α gradually increasing, the amplitudes of one-time and two-time rotating frequency components first decrease and then increase, the amplitudes of two-time rotating frequency components reach the minimum value when α is 200°, and the amplitudes of one-time rotating frequency components reach the minimum value when α is 230°. e amplitudes of half rotating frequency components and one and half rotating frequency components are affected strongly by the nonlinear oil film force, so that the variation trend is more complex, which is approximately along the 180°symmetrical distribution. When the rotating speed is 8000 rpm, under some α values, the spectrum components under the value are complex, and other frequency components other than integral frequency components appear but very small. With α gradually increasing, the amplitude of frequency doubling first decreases and then increases, and the amplitudes of one-time and two-time rotating frequency components reach the minimum value when α reach 295°a nd 280°, respectively. e amplitudes of half and one and half rotating frequency components are small to zero other than that when α is near 100°. Figure 13 shows the waterfall diagram of axis trajectory varying with α. It can be seen from the figure that when the rotating speed is 2000 rpm, the axis trajectory is the same as that in the fault state in Figure 8(a), which is similar to the shape of an ellipse. With α gradually increasing, the track shape remains unchanged, but the size changes. e trajectory is minimum at 180°. When the rotating speed is 5000 rpm, α is 0°to 50°and 320°to 360°, the trajectory shape is relatively regular, which is a large ellipse with a small ellipse, between 50°and 320°, and nonlinear factors have a great impact, so that the shape of axis trajectory is complex and diverse, the system may be in quasiperiodic or chaotic motion period in these speed ranges, and the trajectory gets minimum when α reach 200°. When the rotating speed is 8000 rpm, with α gradually increasing, the shape of the axis trajectory is multiple ellipses that do not coincide, and the trajectory gets to minimum when α reach 300°. Figure 14 shows the variation curve of bearing phase varying with α. ere are two kinds of bearing phase difference to study such as the phase difference in x direction and y direction of two bearings. e value range of phase difference is between − 180°and 180°. Positive value represents that the phase of bearing 1 is ahead of that of bearing 2, while negative value represents lag. In the figure, we can see that the phase differences in x direction and y direction have  same trend to change. When the rotating speed is 2000 rpm, with α gradually increasing, they first increase from 0°to 100°, then decrease to about − 50°, and finally return to 0°. When the rotating speed is 5000 rpm, the phase difference first decreases from 0°to − 180°, during which the phase is negative, and bearing 2 is in lag compared with bearing 1. When α is 180°, the two bearings are reversed in phase, and then the phase difference suddenly changes from − 180°to 180°and then decreases to 0°. When the rotating speed is 8000 rpm, with α gradually increasing, the phase difference first decreases from − 25°to − 175°and then changes slowly when α is in the range of 100°to 250°. At this range, the two bearings are close to reversed in phase. Finally, the phase difference increases back to − 25°. We can get that in different rotating speed, bearings have different phase differences, and their changing trend is also different varying with α.

Conclusion
In this paper, a rod fastening rotor dynamics model is established by using the lumped mass method. Considering the nonlinear oil film force and contact characteristics, the blade scaling, corrosion, wear, and other faults are simulated by unbalance. e effects of rotating speed and fault location on the rod fastening rotor dynamics are studied when two disks have blade faults at the same time, the following conclusions are drawn.
(1) Under normal conditions, the rotor is always in single cycle motion, and the time domain diagram is sine curve. e axis trajectory is standard circle or ellipse. Only one-time rotating frequency components and two-time rotating frequency components appear in the spectrum, and the amplitude of twotime rotating frequency components are too weak compared with one-time components. ree critical speeds appear in the response curve. In case of blade failure, the rotor has different motion states at different rotating speeds, the influence of nonlinear factors is more obvious, the time domain diagram curve and axis trajectory are more complex and diverse, the spectrum diagram has many components, the critical speed and the corresponding amplitude under each critical speed increase remarkably, and the response amplitude fluctuates at the rotating speed of quasiperiodic or chaotic state. (2) When the two disks of the rod fastening rotor have blade faults, the fault location will affect the dynamic characteristics of the rotor system. With the changes of the included angle between the two fault unbalance eccentricities, the motion state of the system, the frequency doubling component of the vibration spectrum, the amplitude of each frequency doubling, the axis trajectory, and the phase difference between the two bearings are different and show a certain regularity under different rotating speeds, which can provide a certain reference basis for the condition monitoring and fault diagnosis of the rod fastening rotor structure of rod fastening rotor.
Because the rod fastening rotor dynamics studied in this paper has strong nonlinearity, limited by test difficulty and conditions, it is difficult to carry out test verification, and the experimental research is not listed in this paper and will be taken in the next stage.

Nomenclature
x, y: Displacement in x and y directions f x , f y : Dimensionless nonlinear oil film force in x and y directions F x , F y : Components of nonlinear oil film force of sliding bearing in x and y directions F xb1 , F yb1 , F xb2 , F yb2 : Components of nonlinear oil film force of sliding bearings 1and 2 in x and y directions X, Y: Dimensionless displacement of the shaft segment at the sliding bearing position in x and y directions _ X, _ Y: Dimensionless velocity in x and y directions c: Sliding bearing clearance δ: Sommerfeld correction factor m P : Half of the disk mass F cx , F cy : Restoring force of the disk contact layer generated by the additional tangential stiffness in x and y directions k: Shaft segment stiffness k 1 : Linear stiffness coefficient of contact layer k 2 : Nonlinear stiffness coefficient of contact layer a: Relative displacement between disks ε: Clearance between tie rod and mating hole x 12 : Relative displacement between two adjacent plates c b1 , c b2 , c i (i � 1, 2, 3, 4): Damping of bearings 1 and 2 and disks 1, 2, 3, and 4 m b1 , m b2 , m i (i � 1, 2, 3, 4): Mass of bearings 1 and 2 and disks 1, 2, 3, and 4 e b2 , e b3 : Blade failure unbalance mass eccentricity θ b2 , θ b3 : Position angle of blade failure unbalance mass eccentricity e i (i � 1, 2, 3, 4): Initial unbalance mass eccentricity of disks 1, 2, 3, and 4 θ i (i � 1, 2, 3, 4): Position angle of initial unbalance mass eccentricity of disks 1, 2, 3, and 4 ω: Rotating speed t: Time of rotor rotation R: Radius of sliding bearing L: Length of sliding bearing μ: Viscosity of oil g: Gravitational acceleration α: Angle between two fault unbalance eccentricities f n : Rotating frequency of rotor.

Data Availability
Data are obtained by measurement, and simulation data are provided in the manuscript.

Conflicts of Interest
e authors declare that there are no conflicts of interest.