Abstract

The impact-rubbing dynamic characteristics of the power turbine rotor with the hollow shaft and offset discs for aircraft engine are investigated, and the impact-rubbing analytical method for the complex rotor based on MDOP Timoshenko beam theory is proposed in this paper. Compared with the traditional approach, the novel method can obtain more data to satisfy the need of engineering. The Lagrange equation is adopted to derive the equations of motion for the rotor system, and the Newmark-β method is applied to solve the equations. The diagrams such as the bifurcation, axis trajectory, spectrum, and Poincaré map are obtained to research on the effect of the rotating speed, gap, and eccentricity on the vibration response. The finite element analysis was carried out to validate the correctness of the theoretical modeling method. The research results indicate that the power turbine rotor with the hollow shaft on operation shows the various nonlinear dynamic behaviors including the multiperiod, quasi-period, jumping phenomenon, and chaotic motions; there exists an optimal gap between the rotor and the stator from the perspective of the efficiency and the dynamics; the optimal gap should make system avoid the resulting chaos or the quasi-period motion for the stability and safety of the machinery.

1. Introduction

Impact-rubbing fault is an ordinary fault in rotation machinery. It is usually an indirect fault caused by misalignment, mass imbalance, or other factors. The fault is very harmful to a rotor system, and it not only affects the safe and stable operation of the system but also leads to more serious accidents such as the blade and the shaft fracture. Therefore, it is crucial to investigate the impact-rubbing mechanism of a rotor system. The gap between the stator and the rotor is closely related to the efficiency of the turboshaft engine. Too large gaps will make the working efficiency of the rotor system decrease, while too small gaps increase the risk of the impact-rubbing because of the vibration. Therefore, in order to pursue the higher efficiency and the lower risk of rubbing, researchers have been exploring an optimal gap value between the rotors and the stators for the complex working conditions. The influences of the clearance, rotational speed, rub-impact stiffness, and imbalance on the system dynamics have become a hot spot in the field of the rotor dynamics.

Many researchers have carried out lots of work on the impact-rubbing dynamics of the rotor system. Zhao et al. [1] developed a drop impact-rubbing equation for a bearing system considering four failure models to simulate the influence rules of the related parameters on the impact-rubbing behaviors of the falling rotor system. Xiang et al. [2] built a dynamic model of the offset rotor system with two discs under consideration of the impact-rubbing and the nonlinear oil film load and investigated the coupled dynamics of the system, and the research is meaningful for the stability of the rotor system. Von Groll and Ewins [3] obtained the dynamic responses of a nonlinear system based on the harmonic balance method. Sun et al. [4] calculated the accurate responses of the steady-state vibration in accordance with MHB-AFT and studied the stability of the impact-rubbing of a dual rotor system. Tai et al. [5] explored the stabilities of a single-point impact-rubbing rotor system according to Hill’s approach; the research proves that the smaller gap can improve the stability of the system. Mokhtar et al. [6] investigated the dynamics of the interaction between the rotor and the stator by applying the Lagrange multiplier approach and obtained the characteristic frequency component of rubbing. Saeed et al. also investigated the nonlinear dynamics of the rub-impact between the stator and the rotor [7, 8], complicated bearing system [911], and self-excited system [12]. Lu et al. [13] studied the dynamic characteristics of a rotor-bearing system with misalignment and impact-rubbing. Pan et al. [14] developed the complicated rotor system model consisting of the rolling bearings, casing, work condition load, rotor imbalance, and squeeze film damper in order to study the effect of the impact-rubbing stiffness, oil film gap, and bearing gap on the vibration responses of the rotor system. The nonlinear dynamic algorithm [15] and the bifurcation characteristics of a multipole magnetic bearing system [16] were focused on by investigators nowadays. Prabith and Krishna [17] determined the effect of the stiffness, friction coefficient, and gap between the stator and the rotor on the stabilities for an impact-rubbing rotor system. Zhang et al. [18] studied the random dynamics of a Jeffcott rotor with parametric uncertainty and impact-rubbing. Fu et al. [19] investigated the nonlinear vibration response of a dual rotor system with the indeterminacy factors. Pan et al. [20] constructed a novel forecasting model of the rotor-bearing system where exists the impact-rubbing between the blade and the stator. Fan et al. [21] presented a novel method of resolving invariable set in the high dimension to study impact-rubbing between the rotor and the stator. Xu et al. [22] studied the dynamic behaviors of a spindle-bearing-housing-belt system with rub by theory and FEM. Ma et al. [23, 24] investigated the impact-rubbing between a rotor and the different shape stators based on the point-point contact theory, and they also studied the fixed point impact-rubbing faults of the rotor with two discs based on the FEA and experiments. Wang et al. [25] explored the sudden imbalance impact-rubbing of a rotor with blade abscission by the theory and experiments. Song et al. [26] researched on the nonlinear dynamics of aircraft engines, especially focused on the impact-rubbing between the blades and the stators by the theory and experiments. Lu et al. [27] developed a kinetic model by using a finite element method based on a dual-rotor experimental platform. Jin et al. studied the nonlinear vibration characteristics of the dual-rotor system through the finite element method [28] and the experimental method [29, 30] and have made much progress.

Researchers have conducted much work on the model building and analyses of the impact-rubbing dynamics for the rotor and made a lot of achievements. Nevertheless, the impact-rubbing dynamics of the rotor with the hollow shaft and two offset discs are often studied by building the model of the less degree of freedom (DOF), and more data from the modeling of multiple degrees of freedom (MDOF) are strongly demanded. The vibration analysis method and technique based on multiple degrees of freedom (MDOF) proposed in the paper can be applied to fault diagnosis, vibration analysis, structure design and optimization, and vibration measurement. Hence, it is necessary to build the MDOF model of the impact-rubbing for the complex rotor. In this paper, the power turbine rotor from the dual rotor system of the aircraft engine is taken as the object of study, seen in Figure 1. The model of the hollow shaft with the two offset discs is developed to be equivalent to the power turbine rotor based on the MDOF Timoshenko beam theory. The effects of the rotational speed, stiffness, and eccentricity on the vibration response of the system are investigated. The research has an important reference value for the vibration analysis and structural design of the same type of rotating machinery.

2. Dynamic Modeling of Power Turbine Rotor

The power turbine rotor includes one hollow shaft and two discs. The shaft is supported by the bearings’ both ends, and it will be modeled as the beam element model; the unbalanced mass of the disc and the bending-torsion coupling effect of the rotor are taken into account in this model. The motion of the disc can be decomposed into the translation in a certain direction and the rotation around the rotation axis. Therefore, the coordinate system OXYZ is at rest, is set up by rotating a certain angle around the X axis, is built by rotating a certain angle around the axis, and is formed by rotating a certain angle around the axis, as shown in Figure 2.

The relation between the stationary coordinate system and the rotating coordinate system can be expressed as

According to the relation between the two coordinate systems, the sum of the translational kinetic energy and the rotational kinetic energy for the unbiased disc can be written aswhere and are the moment of inertia of the disc diameter and the polar moment of inertia of the disc, respectively, and is the mass of the disc. Equation (2) is substituted into the Lagrange equation, and equation (3) is obtained by ignoring the second order and above terms:where is the node vector of the noneccentric disc and and are the mass matrix and the gyro matrix of the disc element, respectively:

Figure 3 shows the disc model considering eccentricity. OXYZ is the stationary coordinate system. is the rotating coordinate system. is the eccentric mass. α is the initial phase angle. e is the eccentricity, and the kinetic energy caused by the eccentric mass can be written aswhere is the vector of the eccentricity in the stationary coordinate system. It can be expressed aswhere and are the OX3 and OY3 components of the eccentricity e in the rotational coordinate system , respectively. The Lagrange equation iswhere is the eccentric element node vector and is the generalized force vector of the eccentric mass. According to equation (7), one obtainswhere , , and are the mass matrix, the damping matrix, and the stiffness matrix of the eccentric mass element, respectively. Based on the above analysis, the formula of the element combined disc is as follows:

2.1. Beam Element Model for Hollow Shaft

The Timoshenko beam is used to establish the element matrix of the shaft considering the rotational inertia and the shear deformation of the beam element. Figure 4 shows a schematic diagram of the Timoshenko beam element. A beam element contains two nodes, and each node contains five degrees of freedom (DOF): the X-direction translation, the Y-direction translation, the rotation around the X-axis, the Y-axis rotation, and the Z-axis rotation. The coordinate system oxyz is established; and are the displacements of the first node of the beam element in the x and y directions, respectively. Similarly, and are the displacements of the second node of the beam element in the x and y directions, respectively. , and are the rotational angles of the first node in x, y, and z directions and the rotational angles of the second node in x, y, and z directions, respectively. Writing the above DOF as a vector, then the node displacement is expressed as

Using the interpolation function method, the displacement of any section in the element is represented by the product of the shape function and the node displacement:where (i = 1, 2, 3, 4), (i = 1, 2, 3, 4), and (i = 1, 2) are all shape functions.where , , E is the elastic modulus of the shaft material, is the moment of inertial for the section, is the shear modulus of the shaft, and A is the cross-sectional area of the shaft. Since the shaft is hollow, the shear influence factor of the beam element κ = 2 (1 + μ)/(4 + 3μ), and μ is Poisson’s ratio of the material. Therefore, the strain energy of the shaft element can be written aswhere J is the polar moment of inertial. Substituting the expression, one can getwhere is the element stiffness matrix in the form of the ten-order symmetric matrix.where

The kinetic energy of the shaft beam element iswhere Ω is rotating speed of the shaft and is the diameter moment of inertia of the hollow shaft. Substituting the formula, one can getwhere is the mass matrix of the element in the form of ten-order symmetric matrix.where and are the outer diameter and the inner diameter of the hollow shaft section, respectively. Considering the gyroscopic effect of the beam element, any infinitesimal segment ds is taken from the hollow shaft section, and the cross-section diagram of the shaft is shown in Figure 5.

The rotational inertia of diameter can be written as

The polar moment of inertia is

The rotation of the microelement around the diameter will produce the kinetic energy increment dT:where and are the arbitrary displacement expressed in the form of the product of the interpolation function and the node displacement vector, respectively. By integrating along the axial section of the entire beam element, one can obtain

According to the Lagrange equation, the gyro matrix of the element is obtained. The matrix is a ten-order antisymmetric matrix as follows:

2.2. Bearing Modeling

Bearing is simplified as a linear spring, and the stiffness matrix of support and the damping matrix are derived by the Lagrange equation:

The Rayleigh damping is used to describe in the total damping matrix of the system. The specific expression is as follows:where α and β are the mass matrix and the stiffness matrix scale coefficients, respectively:where and are the first- and second-order critical rotational speed values, respectively, and are the first- and second-order modal damping ratios, respectively, and α = 0.02 and β = 0.0005.

2.3. Impact-Rubbing Modeling

The research on the impact-rubbing of the power turbine rotor needs to establish a rubbing model which is close to the reality. The rub-impact fault will occur when the radial displacement of the rotor exceeds the gap between the stator and the rotor. The local rubbing diagram between the stator and the rotor is shown in Figure 6, where is the radial impact force of the rotor, is the tangential friction force of the rotor, ω is the rotational speed of the rotor system, and φ is the angle between the normal direction and the x axis at the rubbing point.

Considering that the radial deformation of the stator is linear, is the radial stiffness of the stator, is the radial displacement of the rotor, and the coefficient of friction between the rotor and the stator is f. The friction law conforms to Coulomb’s Law, and friction is proportional to the radial pressure acting on the friction point. Assuming that the clearance between the rotor and the stator is when the system is stationary, the impact-rubbing force iswhere . Decomposing the impact-rubbing force into the OXY coordinate system produces the following equations:where is the friction torque generated by the friction force on the disc.

3. Equations and Solutions

The equation of motion for the power turbine rotor system can be derived based on the beam element model, bearing model, and impact-rubbing force model. If the rotor system has N nodes, the system consists of N − 1 beam elements. Letting q be the displacement column vector of the power turbine rotor, the equation of motion for the power turbine rotor system is as follows:where

[M], [C], [G], [K], and [F] are the mass matrix, damping matrix, gyro matrix, stiffness matrix, and force vector of the overall system, respectively.

According to above modeling, the FE model of the power turbine rotor was established as shown in Figure 7. This model consists of 16 nodes and 15 beam elements. Disc 1 is loaded at node 11 with a diameter of 0.2 m, and disc 2 is loaded at node 12 with a diameter of 0.3 m. The left and right bearings are set at nodes 2 and 14, respectively. The material used for the discs and shaft is the same, and the related parameters of the power turbine rotor system are shown in Table 1.

The Newmark-β method is used to calculate the dynamic response of the power turbine rotor system. One establishes at t-moment and at -moment, and the first-order Taylor expansion of the velocity function is obtained:where is the acceleration of a point in the interval (), and it can be assumed

Substituting equation (33) into equation (32) leads to

The first-order Taylor expansion of the displacement function is obtained:

Supposing

According to the above relationship, and can be expressed as follows in the form of and :

Substituting equations (37) and (38) into equation (30) produceswhere

4. Validation

In order to validate the correctness of the modeling method proposed in this research, the finite element model of the rotor system was established according to the same parameters and constraints as theoretical modeling, as shown in Figure 8, and the time history curve and frequency-domain diagram of the vibrational displacement at the center of mass for the disc under unbalanced load are drawn; the comparison of the finite element method (FEM) result and the theoretical modeling method (TMM) result is conducted, as shown in Figures 9(a) and 9(b). It can be inferred from the figures that the theoretical and finite element calculations agree well, which verifies the correctness of the theoretical modeling method.

5. Results and Analyses

5.1. Rotating Speed

According to the equations in the previous section, the working speed, eccentricity, impact-rubbing, gap, and friction coefficient affect the vibration characteristics of the whole rotor system. The working rotational speed is a fundamental parameter for rotational machinery. Therefore, the working speed of the rotor is selected as the controlled variable, and the related parameters of the system are shown in Table 2. The working speed is selected at 1000 rad/sec∼2500 rad/sec for analysis, the eccentricity is 5 × 10−4 m, the impact-rubbing stiffness is 2 × 107 m, the gap is 6 × 10−4 m, and the coefficient of friction is 0.1. The time-domain data of the system vibration can be obtained after calculating, the transient vibration caused by the initial conditions in the data is abandoned, and the steady-state vibration is used for analysis. The bifurcation diagram of the vibration with change in the rotational speed can reflect the dynamic characteristics of the system at different speeds, such as the periodic vibration, multiple periods, quasi-period, chaos, and jump. The system response is calculated once for every 10 rad/sec with increasing of the rotating speed (after a lot of calculations and analyses, the step length is more appropriate for research focus than others), and the system bifurcation diagram of the displacement with the change in the working speed is drawn as shown in Figure 10. It is shown from Figure 10 that the system exhibits a variety of nonlinear dynamic behaviors such as the period, quasi-period, and chaotic motion; the bifurcation diagram is a single curve (period-1) when the speed is from 1000 rad/sec to 1240 rad/sec; when the speed is in the range of 1240 rad/sec∼1931 rad/sec, the system is under the state of the quasi-period, multiperiod, and chaotic motions; when the speed is from 1931 rad/sec to 2026 rad/sec, the system is always in the stable period-1 motion; when the speed is from 2026 rad/sec to 2329 rad/sec, the system is also in the quasi-periodic motion, multiperiod motion, and chaotic motion. The rotating speeds 1300 rad/sec, 1700 rad/sec, 2000 rad/sec, and 2200 rad/sec are taken to analyze the specific nonlinear dynamic behaviors.

Figure 11 shows the response diagram of the system at the rotating speed of 1300 rad/sec, including the trajectory of shaft center, Poincaré map, and spectrum. It can be seen from Figure 11(a) that the blue trajectory exceeds the red boundary curve of the system and is characterized by many irregular curves, and the system is in quasi-periodic motion. It can be observed from Figure 11(b) that there are many irregularly distributed points. Figure 11(c) shows more frequency divisions. Figure 12 presents the trajectory of the shaft center, Poincaré map, and spectrum at the rotating speed of 1700 rad/sec. Figure 12(a) shows that local impact-rubbing has occurred between the rotor and the stator. In Figure 12(b), there are three points from which it can be inferred that the system enters period-3 motions from the chaotic region, and Figure 12(c) confirms this point. It can be seen from Figure 13(a) that the system does not intersect with the boundary and moves within the boundary. There is only a single isolated point in Figure 13(b), and it is inferred that the system is undergoing period-1 motion. Figure 13(c) shows only one frequency of the rotating speed. Figure 14(a) presents the regular petal shape that is complex but not chaotic. Figure 14(b) shows an approximately closed triangle. There are X/4 and X/2 of the rotating frequency in Figure 14(c). In conclusion, the vibration characteristics of the rotor are sensitive to the rotational speed, the rotor shows complicated motions such as period-1, period-3, and the quasi-period, and the order of these motions is almost irregular. The working rotational speed of the rotor shall make the system in periodic motion to avoid the quasi-periodic or the chaotic motion.

5.2. Gap between Rotor and Stator

The gap between the rotor and the stator is an essential parameter for the efficiency and dynamics. From the perspective of safe operation, the radial gap should not be less than the largest displacement of the vibration. The gap interval of 0.2 mm∼0.8 mm from Table 3 is chosen to draw the bifurcation diagram of the displacement, as shown in Figure 15. It can be shown from Figure 15 that the bifurcation diagram is a single curve (period-1) when the gap is from 0.2 mm to 0.4 mm; the jump occurs when the gap is 0.4 mm; when the gap increases from 0.4 mm to 0.49 mm, the mutual transformation of period-1 and the multiperiod occurs; when the gap is between 0.49 mm and 0.57 mm, the rotor system is always in the stable period-1 motion; when the gap interval is between 0.57 mm and 0.775 mm, the system is in the quasi-period motion, multiperiod, and chaotic state. The gaps of 0.41 mm, 0.48 mm, 0.54 mm, 0.62 mm, 0.70 mm, and 0.77 mm are selected for analysis.

Figure 16 includes the trajectory of the shaft center, Poincaré, and spectrum diagrams when the gap is 0.41 mm. Figure 16(a) shows that the axis trajectory does not exceed the boundary curve, but the rotor and the stator collide at a certain point. There is only one point in Figure 16(b), indicating that the system is in the period-1 motion which is consistent with Figure 16(c). Figure 17(a) shows that the blue trajectory exceeds the red boundary curve of the system and many stacked ellipses appear when the gap is 0.48 mm. Figure 17(b) presents the quasi-periodic motion of the system. In Figure 17(c), there are many messy and small frequency components which are in accordance with Figure 17(b). It can be inferred from Figure 18(a) that there is no impact-rubbing when the gap is increasing up to 0.54 mm. Figures 18(b) and 18(c) show the same conclusion. In Figure 19(a), there are many circles stacked together like a petal shape. Figure 19(b) shows a distorted closed triangle. It can be observed from Figure 19(c) that there is 3X/10 of the speed frequency except for the speed frequency and the system is in quasi-period motions. When the gap is increasing up to 0.7 mm, the three diagrams nearly remain unchanged, as shown in Figure 20. Figure 21 presents no impact-rubbing of the system when the gap is 0.8 mm. In summary, the gap between the rotor and the stator has a great influence on the dynamic characteristics of the whole system. The dynamics and efficiency of the system should be comprehensively considered for the gap design. The system is in the period-1 motion as much as possible to avoid the multiperiod motion or chaos.

5.3. Impact-Rubbing Stiffness

The impact-rubbing stiffness has a certain correlation with the vibrational response of the rotor. Table 4 shows the parameters of the impact-rubbing for the rotor system in this calculation example. Figure 22 shows the bifurcation diagram of the displacement for the rotor with the change in impact-rubbing stiffness. It can be shown from Figure 22 that the jump occurs when the impact-rubbing stiffness is 0.5 × 108 N/m. The system is in the period-1 motion when the impact-rubbing stiffness is from 0∼1.4 × 108 N/m; when the impact-rubbing stiffness is greater than 1.4 × 108 N/m, the system is in the mutual transformation of the quasi-period motion and the chaotic motion. The impact-rubbing stiffnesses 0.5 × 108 N/m, 1.4 × 108 N/m, 3.3 × 108 N/m, and 4.3 × 108 N/m are selected to draw the trajectory of the shaft center, Poincaré, and spectrum diagrams, as shown in Figures 2325. Figures 23(b) and 24(b) show only one isolated point which can be identified by Figures 23(c) and 24(c) presenting only one frequency of the rotating speed. Figures 23(a) and 24(a) indicate that the system is in a stable period-1 motion. Figures 25(b) and 26(b) show many points which form a closed irregular shape. It can be seen from Figures 25(a) and 26(a) that the system is in full-cycle rubbing. It is inferred that the rotor vibration hardly changes within a certain range of the small stiffness which is called the ideal stiffness range, and when the stiffness exceeds a critical value, the rotor vibration gradually increases. This critical stiffness value is of great significance to the structural design in engineering. Generally, the impact-rubbing stiffness of the rotor is less than this critical value and in the ideal stiffness range so as to avoid working for a long time in nonperiodic motion conditions. This work has some guiding significance for the design of rotor and stator stiffness, which can be achieved by selecting material and surface treatment.

5.4. Eccentricity

The eccentricity has a significant effect on the vibration of the system. The influence of the eccentricity on the dynamics will be studied in this section. The related parameters are shown in Table 5 where the eccentricity is in the range of 0.1 mm∼1 mm, and the bifurcation diagram of the eccentricity is drawn, as shown in Figure 27 which shows a straight line; i.e., the system is in the state of period-1 when the eccentricity is 0.1 mm∼0.53 mm; when the eccentricity is 0.53 mm∼0.84 mm, the motion of the system is more complicated; when the eccentricity is in the range of 0.84 mm∼0.94 mm, the bifurcation diagram presents a straight line again; when the eccentricity is 0.94 mm∼1 mm, the system is chaotic and complex. The eccentricities 0.5 mm, 0.6 mm, 0.7 mm, 0.84 mm, 0.91 mm, and 0.98 mm are chosen to draw the corresponding trajectory, the Poincaré map, and the spectrum diagrams, as shown in Figures 2833.

Figure 28 shows that the system is in a stable period-1 motion and the rotor does not rub and impact with the stator. Comparing with Figure 28(a), Figure 29(a) is more complex; Figure 29(b) shows the system is approximately in a period-3 motion; with the increase of the eccentricity up to 0.91 mm, the amplitude of the 1X frequency increases gradually and exceeds the amplitude of the 1/3 frequency in the spectrum diagram, while the system enters into the period-1 motion from the period-3 motion in Poincaré diagrams, as shown in Figure 32(b); with the increase of eccentricity to 0.98 mm, the axis trajectory enters into the disorderly state, as shown in Figure 33(a). It can be inferred from the analysis above that the dynamic characteristics of the system are sensitive to the change of eccentricity, and the system mainly presents period-1, period-3, the quasi-period, and other motion types; the order of occurrence has no obvious rule. In engineering, the dynamic balance of the rotor system can usually be conducted to adjust the eccentricity, and the ideal value of the eccentricity is determined according to the vibrational level of the rotor. The specific dynamic behaviors should be analyzed according to the actual situation.

6. Conclusions

The impact-rubbing dynamic characteristics of the power turbine rotor system with the hollow shaft and two offset discs for aircraft engines are investigated in this paper. Timoshenko beam theory considering the MDOF is applied to establish the rotor model, and the Newmark-β method is used to solve the equations of motion for the system. The influence of the rotating speed, gap, and eccentricity on the vibration response of the system is studied; the finite element analysis was conducted to verify the correctness of the theoretical modeling method. The main conclusions are drawn as follows:(1)The power turbine rotor with the hollow shaft on operation shows the various nonlinear dynamic behaviors such as the multiperiod motion, quasi-periodic motion, jumping phenomenon, and the chaotic motion. The impact-rubbing stiffness and the eccentricity also have the significant effect on the bifurcation characteristics of the system.(2)The rotating speed has a fundamental influence on the nonlinear dynamics of the system. With the increase of the speed, the system switches among period-1, period-3, the multiple period, and the chaos motions. The working rotational speed should be designed in period-1 for stability; otherwise, the vibration of the system is so unpredictable that a potential risk may occur.(3)There exists an optimal gap between the stator and the rotor from the perspective of the mechanical efficiency and the dynamics of impact-rubbing. Too large gaps will decrease the efficiency of the rotational machinery, while too small gaps will increase the risk of the impact-rubbing. The optimal gap should make system avoid the chaos or the quasi-periodic motion.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This project was supported by the National Natural Science Foundation of China (Grant no. 52275118).