Investigation on dynamic characteristics of a rod fastening rotor-bearing coupling system with fixed-point rubbing

The aim of this paper is to gain insight into the nonlinear vibration feature of a dynamic model of a gas turbine. First, a rod fastening rotor-bearing coupling model with fixed-point rubbing is proposed, where the fractal theory and the finite element method are utilized. For contact analysis, a novel contact force model is introduced in this paper. Meanwhile, the Coulomb model is adopted to expound the friction characteristics. Second, the governing equations of motion of the rotor system are numerically solved, and the nonlinear dynamic characteristics are analyzed in terms of the bifurcation diagram, Poincaré map, and time history. Third, the potential effects provided by contact degree of joint interface, distribution position, and amount of contact layer are discussed in detail. Finally, the contrast analysis between the integral rotor and the rod fastening rotor is conducted under the condition of fixed-point rubbing.


Introduction
As a typical high-tech large-scale military and civil power equipment, gas turbine has been widely used in the fields such as energy, transportation, and aviation. To satisfy the requirements of performance, processing, maintenance, and transportation, the combined structure is designed to replace the traditional integral structure and gradually adopted in the manufacturing process. Under this circumstance, rod fastening rotor has been one of the generally recognized structure design of heavy-duty gas turbine because of its advantages of easy assembly and maintenance [1] . Rod fastening rotor system is mainly composed of pull rods, discs, and blades. Several discs are integrated by the pressure and the friction force provided by the pull rods [2] . Since the global stiffness of the rotor system is closely related to the contact stiffness of the joint interface, it is significant for scholars to master its contact mechanism in advance. At present, three main methods for determining the contact stiffness can be concluded as equivalent stiffness method, stiffness matrix correcting method, and virtual material method. In the first one, the equivalent spring element and equivalent damper element are defined to describe the contact mechanism of joint interface. For a bolted beam, Hartwigsen et al. [3] identified the interface mechanical parameters from a transient excitation test. Utilizing Lagrange's equation, Hu et al. [4] proposed a dynamic model of rod fastening rotor-bearing system, in which the contact effects between the discs of the bolted joint were described by the cubic nonlinear stiffness. Qin et al. [5] developed an analytical model for the bending stiffness of the bolted disc-drum joints and implemented the joint stiffness into the dynamic model of a simple rotor connected through the bolted discdrum joint. Qin et al. [6][7] investigated the time-varying stiffness at the joint interface with bolt loosening by means of three-dimensional nonlinear finite element models and provided a general approach for the vibration analysis of a rotating cylindrical shell coupled with an annular plate. What is more, many researchers have made great endeavors on the stiffness matrix correcting method and virtual material method. Wu et al. [8] built a contact stiffness matrix capable of describing the contact effect between discs, where the lateral stiffness, shear stiffness, bending stiffness, and torsional stiffness of the joint interface were considered. Jalali et al. [9] used the thin-layer element theory to investigate the surface-to-surface contact interface of a compound structure.
The contact stiffness of the joint interface plays a crucial role in the vibration behaviors of a rod fastening rotor system. The most favorable model for this purpose is the Greenwood-Williamson (G-W) statistical model [10] . According to the boundary element method, Pohrt et al. [11] calculated the contact stiffness of the elastomer with fractal rough surface. Taking a two-dimensional fractal rough surface profile as the object, the load-compliance relation for the contact of this profile was estimated in Ref. [12]. Liu et al. [13] discussed the influence of friction on the normal contact stiffness of isotropic fractal interface when the elastic-plastic contact happened. When the surface fractal feature and normal loading conditions were considered, a new modeling approach was proposed to represent the tangential frictional stick-slip behaviors of contact interfaces in Ref. [14]. Jiang et al. [15] carried out research on interface stiffness based on the fractal-curve description of the topology. Moreover, the dynamic characteristics of the rod fastening rotor system have attracted extensive attention. By numerical simulation, Rimpel et al. [16] discussed the change law of natural frequency of a real tie bolt rotor system subject to different preload forces. From the perspective of nonlinear dynamics, Hei et al. [17] observed rich nonlinear phenomena in the rod fastening rotor system, including quasi period, bifurcation, and chaos.
The performance and efficiency of rotating machines can be enhanced by reducing the clearance between components in relative motion. However, this approach greatly increases the risk of rub-impact fault. Rub-impact can accelerate erosion, thermal fatigue, and even other catastrophic failures. In the past few decades, the full annular rub and partial rub have attracted a huge amount of interest already. Yu [18] studied the reverse full annular rub occurring in a two degrees-of-freedom rotor/seal model, in which the rubbing location was simulated away from the lumped rotor mass. When a modified Jeffcott rotor underwent either forward synchronous whirling or self-excited backward whirling motions with continuous stator contact, Vlajic et al. [19] investigated its complicated vibration behaviors. Jiang [20] gained deep insights into the global response characteristics of a piecewise smooth dynamical system with full annular rub and revealed the effects of parameters, including the coefficient of friction and impact stiffness. Yang et al. [21] simulated the dynamic behaviors of a dual-rotor system with partial rub and obtained the fault characteristics, such as combination forms in the frequency domain. Moreover, Yang et al. [22] investigated the fault characteristics of the rotor system with parametric uncertainties based on the polynomial chaos expansion.
Up to now, the contact mechanism of bolted joint interface and dynamic characteristics of rod fastening rotor system have been studied in depth. Nevertheless, the effects of multiple disc contacts on the vibration behaviors of the rotor system have been rarely studied. Moreover, compared with the integral rotor system, the pioneering contribution to the dynamic prediction of rod fastening rotor system with rub-impact fault is not enough. It is more fulfilling to include this feature in a rod fastening rotor model if more realistic and useful results are required.
This paper aims at studying the nonlinear dynamic characteristics of a rod fastening rotorbearing coupling system with rub-impact fault. The law of effects of joint interface in disc-disc contact is revealed in depth. On this basis, the comparison of fault responses between integral rotor and rod fastening rotor is conducted as well.

Mathematical model
In this section, the dynamic model of a circumferential rod fastening rotor system with rub-impact fault is established based on the rough fractal theory and finite element method. Meanwhile, a novel contact force model considering coating effects is adopted to solve the rub-impact fault. Figure 1(a) shows the actual structure of a circumferential rod fastening rotor system. By extracting the main structure features, the simplified physical model is given in Fig. 1(b). The model is formed by a flexible shaft, a disc group composed of eight discs with blades, eight long pull rods, and two pedestals. In the process of modeling, some basic assumptions are given as follows.
(i) The rotating shaft is considered as the Rayleigh beam, in which the shear effects and torsion deformation are neglected.
(ii) All the discs with blades are assumed to be rigid, and their deformations are ignored.
(iii) The balls of the bearing model are only rolling, and the bearing deformation is mainly considered as contact deformation between rolling balls and raceway.
(iv) All the long pull rods are evenly distributed on the discs along the circumference. (v) The contact phenomenon exists on the joint interface between any two discs.
(vi) In the process of rub-impact, the thermal effect and friction torque are ignored. In addition, the thermal effect caused by the joint interface is not considered.

Contact mechanism of joint interface
As shown in Fig. 2, two rough surfaces of joint interface contact each other, which can be regarded as an equivalent rough surface in contact with a rigid flat surface. By utilizing the rough fractal theory given in Ref. [23], the contour curve of the nth asperity before elastic deformation can be expressed as where G denotes the feature parameter of the rough surface profile, D denotes the fractal dimension of the two dimensional fractal rough surface, and l n is the sampling length of the nth asperity which obeys l n = 1/γ n . After obtaining Eq. (1), the curvature radius of the nth asperity shown in Fig. 2 can be further derived as Then, the expression of the height of the nth asperity is given by Correspondingly, the actual deformation of the nth asperity is approximately equal to According to the Hertz contact theory, the relation between the contact area and the contact force of the nth asperity takes the following forms: In the elastic deformation process of the nth asperity, the maximum value under the critical condition can be expressed as where E denotes the equivalent elastic modulus, E 1 and E 2 denote the elastic moduli of hard and soft materials, respectively, and κ denotes the coefficient relating to Poisson's ratio. The above derivations show the deformation mechanism of the nth asperity in the contact process. For a joint surface, the different asperities are distributed in the contact surface by referring to a certain probability density. For this case where the rough surface is in contact with the rigid plane, the area distribution function of the asperity conforms to the law of ocean area distribution [24] , namely, where a l denotes the maximum contact area of the asperity, and a denotes the contact area of the asperity. By integrating the above expression, the actual contact area of the whole rough surface can be obtained as Due to the multi-scale characteristics of the rough surface, it is assumed that the area distribution probability of different asperities in the rough plane is n n (a), and its value depends on the scale ordinal number n. When n n (a) = Cn(a), the actual contact area of the whole rough surface is also expressed as By comparing Eq. (8) with Eq. (9), the expression of the variable C is Therefore, the total elastic contact force and contact stiffness of the joint surface can be expressed respectively as in which the critical elastic contact area and critical elastic frequency index are given by After analyzing the normal contact mechanism of the joint interface, the bending stiffness of the rough contact layer shown in Fig. 3(a) is written as where I d denotes the second moment of the contact area.
According to the description given in Ref. [25], the tangential contact stiffness of the contact layer can also be obtained as where ν d denotes Poisson's ratio of the disc. It is noted that Eq. (15) can be used to modify the stiffness characteristic of the disc-disc contact segment. Rough contact layer Contact between two discs due to looseness Because the long rods are distributed in different positions of the disc set, their distances from the neutral layer of the rotating shaft are not the same. The distance between the neutral layer and the ith rod is assumed to be d i .
For the rod fastening rotor system, the angle between two discs in the process of whirling motion is set to θ (see Fig. 3(b)). Therefore, the moment acting on the ith rod can be expressed as in which r r denotes the radius of cross section of the ith rod, Δl i denotes the initial axial deformation caused by the preloaded, l i denotes the initial length of the ith rod, and E r denotes the elastic modulus of the rod. After that, the total moment of all the rods can be calculated by summation. The bending stiffness provided by the rods can be further given by Due to the parallel connection between the rough contact layer and the rods, the total bending stiffness of the disc-contact segment is modified as

Finite element discretization of rotor system
The flexible shaft is divided into several segments, and the disc set with long rods is seen as three parts, including the equivalent element of left discs, the disc-disc contact element, and the equivalent element of right discs.
Combining with the whirling motion of the rotor system, the flexible shaft is modeled as the Timoshenko beam. According to the vibration theory, there are four degrees of freedom at each node, including two translations and two rotations. Therefore, the generalized displacement vectors of the beam element in the coordinate planes of Oxz and Oyz can be expressed respectively as where x A , x B , y A , and y B represent four translational degrees of freedom at node A and node B respectively. Meanwhile, θ yA , θ yB , θ xA , and θ xB represent four rotational degrees of freedom at the two nodes.
On this basis, the translational inertial matrix of the beam element is given by where ρ s denotes the density of beam element, l s denotes the length of beam element, and A s denotes the area of cross section. The expressions of coefficients used in Eq. (20) are where μ s denotes the shear coefficient of the cross section, G s = E s /(2(1+ν s )) denotes the shear modulus, E s denotes the elastic modulus of beam element, and I s denotes the area moment of inertia.
Similarly, the rotational inertial matrix of the beam element can also be derived as in which the inertial coefficients take the following forms: The gyroscopic effect is a special feature of rotating machinery which can directly affect the dynamic characteristics of the system. Therefore, attention should be paid to the process of dynamic modeling. The gyroscopic effect is not only related to the structure size, but also to the rotational speed of the system. In the following part, the gyroscopic matrix of the beam element can be expressed as where ω denotes the rotational speed of the rotor system. Moreover, according to the deformation characteristics of the Timoshenko beam, the stiffness matrix of the beam element obeys where the coefficients of the stiffness matrix are When dealing with the disc set shown in Fig. 1(b), it can be treated as two kinds of types, including the disc-disc contact element and the equivalent element of left/right discs. From Fig. 4(a), it can be seen that the rough contact layer is sandwiched between two single discs.
Under this condition, the force transmission and deformation relationship of the series structure should be considered. Then, the stiffness matrix of the disc-disc contact element can be expressed as where the coefficients used in Eq. (27) take the following forms: where μ d denotes the shear coefficient of the cross section of the single disc, denotes the shear modulus of the single disc, E d denotes the elastic modulus of the single disc, and I d denotes the area moment of inertia of the single disc. Except for the disc-disc contact element, both the left discs and the right discs are equivalent to the lumped mass model, as shown in Fig. 4(b). Because the elastic deformation of the structure is ignored, each lumped mass element has only one node. Thus the mass matrix and gyroscopic matrix of the equivalent element of left/right discs are respectively written in a 2 × 2 form in the coordinate plane of Oxz, i.e., where m d denotes the mass of the lumped mass model, and J dd and J pd denote the diametric inertia and the polar inertia, respectively. According to the design of element division, the global mass, stiffness, and gyroscopic matrices of the rod fastening rotor system can be integrated. It should be noted that due to the existence of elastic support boundary, the global stiffness matrix needs to be modified. Consequently, the governing equations of motion of the rod fastening rotor system can be expressed as where u denotes the generalized displacement vector, M denotes the mass matrix, C denotes the damping matrix, K denotes the stiffness matrix, and Q denotes the generalized force vector.

Nonlinear constraint boundary
The nonlinear force provided by the constraint condition is studied in this section. The rolling bearing model is mainly composed of rolling balls, inner ring, outer ring, and bearing housing. The bearing housing is fixed on the foundation.
To describe the point-to-point contact mechanism between rolling balls and raceway, the Hertz contact theory is adopted here. Obviously, the contact forces are nonlinear, and their expressions in the two directions of Ox and Oy are where H represents the Heaviside function, k h denotes the point-point contact stiffness, and c b denotes the support damping of bearing. The contact displacement is determined by the initial gap of bearing model and the whirling motion of the rotor system. Therefore, the contact displacement of the ith rolling ball is written as where δ b0 denotes the initial gap of the bearing model, and x b and y b denote the vertical and lateral displacements of shaft node at the bearing position, respectively. According to the motion trajectory, it is easy to calculate the specific position of the ith rolling ball, namely, where R b and r b denote the outer and inner radii of the bearing model, respectively, and N b denotes the number of rolling balls.

Mechanical mechanism of rub-impact fault
Due to the complicated external suspension, such as fuel tank, cooler bin, and pipeline, the convex point is more likely to appear in the stator. Fixed-point rubbing means that the rotor contacts at a fixed point of the stator once in one cycle. When eccentricity exists in the equivalent element of left discs, to a remarkable extent delivered, fixed-point rubbing may be induced. More specifically, if the whirling amplitude of the disc is larger than the initial gap, the rub-impact fault happens. Otherwise, the fault never happens.
According to Ref. [26], a contact force model is adopted to reveal the mechanical mechanism of rub-impact under the influence coatings. The effectiveness of the contact model is verified in an impact experiment at different rotational speeds on a rotor test rig. The rotor-stator impact force can be expressed as where k s denotes the structure stiffness of the stator, and k l denotes the local contact stiffness of coatings.
In the above expression, the local deformation of coatings δ l and the global impact deformation δ N obey the following relation: At the same time, the Coulomb law is adopted to describe the friction characteristics between the rotor and the stator, so that where μ denotes the coefficient of Coulomb friction, θ s denotes the angular position of the stator, andẋ d andẏ d denote the vertical and lateral velocities of the disc, respectively.

Simulation and discussion
For a multiple degrees-of-freedom system with strong nonlinear factors, it is relatively difficult to analytically study its nonlinear characteristics. In this paper, the integral rotor system has 12 degrees of freedom, and the rod fastening rotor system has 16 degrees of freedom. The sources of nonlinearity are mainly reflected in two aspects. On the one hand, when the radial displacement of the rotor is greater than the initial gap, the rub-impact fault happens. Obviously, this leads to a clearance-type nonlinearity. On the other hand, the rolling bearing can bring stiffness nonlinearity to the system. Then the combination of the Runge-Kutta method and linear interpolation method is chosen. The structure parameters of rod fastening rotor-bearing coupling system are given in Table 1. The flow chart of dynamic prediction of the rod fastening rotor system is given in Fig. 5, in which there are four aspects, i.e., dynamic modeling, numerical simulation, dynamic characteristic analysis, and rub-impact fault investigation.
The motion state is checked to find if there has been a change of contact condition from no-rubbing to rubbing or from rubbing to no-rubbing. If this change has happened during a time step, linear interpolation will be used repeatedly to determine the time instant when this change happens, so that the tolerance value of rub-impact region is met. A case of 'no-rubbing to rubbing' is analyzed, in which δ 0 denotes the initial gap, Δt denotes the original time step, and γ is defined to describe the rub tolerance zone. Then, the case of 'no-rubbing to rubbing' can be expressed as where δ k and δ k+1 represent the radial displacements of the disc at the previous and following moments, respectively. If the condition of |δ k+1 − δ 0 | γ is satisfied, the time step remains the same. Otherwise, the new time step should be adopted by the linear interpolation, so that Thus, the radial displacement of the disc δ * k+1 is recalculated by utilizing the new time step Δt * and the initial condition of δ k . If the radial displacement δ * k+1 satisfies the rub tolerance condition, the iteration is finished.

Nonlinear responses of rotor system
In this section, fixed-point rubbing is not accounted provisionally, and the contact effect of rough joint interface between two discs is paid more attention. As one of useful ways of observing nonlinear vibration behavior, the bifurcation diagram is used to provide a summary of the essential dynamics.
For the case of no disc-disc contact and long rods, the rotor system can be seen as a common integral model to a great extent. At this time, the response bifurcation diagram of the system is depicted in Fig. 6(a). From this figure, it can be observed that the system responses are exhibited as 1 T -periodic motion in the range of rotational speed 300 rad · s −1 to 2 000 rad · s −1 .
At ω = 1 430 rad · s −1 , the first-order natural frequency of the rotor system is equal to the imbalance excitation frequency. After that, due to the nonlinear constraint boundary provided by the rolling bearings, the nonlinear jump phenomenon happens. With the increase in the rotational speed, the integral rotor system always keeps the regular periodic motion.
Under the same structure parameters and operating conditions, the disc-disc contact behavior and long rods are introduced into the system. The bifurcation diagram of the rod fastening rotor is shown in Fig. 6(b). Through the comparative analysis, the vibration response of the rotor system has changed significantly. The first-order resonance frequency of the system decreases to ω = 1 340 rad · s −1 . The main reason for this phenomenon is that the rough joint interface between discs can weaken the global stiffness of the system. Moreover, the region of complicated vibration forms, such as quasi period and even chaos, is observed in the region of high rotational speed. This just proves that for the rod fastening rotor system, the vibration response is obviously determined by the contact condition of discs. To fully understand the vibration characteristic of the rod fastening rotor system at different rotational speeds, the time history and Poincaré maps are used, respectively. As shown in Figs. 7(a) and 7(b), the 1 T -periodic motion happens around a position which is almost zero, and the vibration amplitude is about 9.38 × 10 −5 m. Meanwhile, only an isolated point is found in the Poincaré map. When the rotational speed mounts up to ω = 1 600 rad · s −1 , the original regular periodic motion will be broken and evolve into quasi-periodic motion. At this moment, the Poincaré map includes a closed petal curve constructed by discrete points, as shown in Figs. 7(c) and 7(d). According to the theory of nonlinear dynamics, the phenomenon in the Poincaré map illustrates the occurrence of quasi-periodic motion.
Additionally, the more complex nonlinear vibration behaviors can be observed in the response of rod fastening rotor system. When the rotational speed is ω = 1 700 rad · s −1 , the rotor system enters the window of chaotic motion. Under this circumstance, the Poincaré map is composed of a large number of irregular scattered points.
Next, the effects of distribution position on the nonlinear dynamic behavior of rod fastening rotor system are discussed. As shown in Fig. 8, there are four working conditions for describing the distribution position of the rough contact layer, including '2-6' distribution, '4-4' distribution, '5-3' distribution, and '6-2' distribution. In this paper, the '2-6' working condition is different from the '6-2' working condition. The main reason is that the disc group is not installed in the middle part of the rotating shaft, and thus the system is not an asymmetric structure with respect to the mid-point of the shaft. Then the comparative analysis of rod fastening rotor system under the above four working conditions is conducted. From Fig. 9, it appears clearly that the first-order resonance frequency of the system is obviously affected by the distribution position of the rough contact layer. For example, the firstorder resonance frequency of the system under '2-6' distribution is about ω = 1 300 rad · s −1 . While keeping the structure parameters unchanged, the resonance phenomenon happens at ω = 1 410 rad · s −1 under '6-2' distribution. Not only that, affected by the distribution state, the rotational speed range of complicated vibration response turns to be narrow, which means that the motion stability of the system is improved.
Different from Figs. 9(a) and 9(b), the system responses given in Figs. 9(c) and 9(d) are exhibited as the alternate forms of periodic motion, quasi-periodic motion, and chaos at the interval from 1 650 rad·s −1 to 2 000 rad·s −1 . Generally speaking, for rod fastening rotor system, the motion state and resonance characteristic are closely related to the distribution position of rough contact layer.

Effect of rough joint interface
This part focuses on the study of the contact characteristics of the rough joint interface by discussing the variation of the first-order resonance frequency. For the convenience of analysis, the stiffness influence coefficient is defined as e k , where its value is set to (1/3, 1/2, 2, 3).
Different from the integral rotor system, the inherent characteristics of the rod fastening rotor system are not only related to the rotating shaft and constraint boundary, but also affected by the contact state of joint interface between discs.
Through the forward sweep frequency analysis, the first-order resonance frequency of the rotor system with different stiffness coefficients is depicted in Fig. 10. The case of e k = 1/3 means that the contact degree between two discs is relatively light and the contact stiffness of the rough contact layer is reduced. Therefore, the first-order resonance frequency of the system decreases from ω = 1 340 rad · s −1 to ω = 1 230 rad · s −1 , as shown in Fig. 10(a). When the stiffness influence coefficient is e k = 1/2, the first-order resonance frequency is modified as  ω = 1 300 rad· s −1 (see Fig. 10(b)). On this basis, the natural characteristics of the rotor system are further investigated under the conditions of e k = 2 and e k = 3. It can be seen that the first-order resonance frequency could increase from ω = 1 300 rad · s −1 to ω = 1 360 rad · s −1 .

Effect of multiple disc contacts
In the previous works, for the rod fastening rotor system, the influence mechanism of a contact layer between two discs is mainly discussed. However, the contact analysis of multiple discs is rare in the existing research. Actually, multiple disc contacts may occur during the operation of rod fastening rotor. It is worth studying in order to improve the operation safety of the system. In the following work, the nonlinear characteristics of the rotor system affected by multiple disc contacts will be discussed by means of bifurcation diagram. Figure 11 shows the different contact conditions of multiple disc contacts, including two contact layers and three contact layers. When the two contact layers are considered in the dynamic model, the bifurcation diagram of the rotor system is given in Fig. 12(a). By comparing Fig. 12(a) with Fig. 6(b), it is clear that the first-order resonance frequency decreases from ω = 1 340 rad · s −1 to ω = 1 320 rad · s −1 . The system enters the window of quasi-periodic motion in the rotational speed range (380, 395) rad · s −1 . Moreover, the nonlinear vibration zone is shifted to the left as a whole.  When the rod fastening rotor system contains three contact layers, the bifurcation diagram is shown in Fig. 12(b). Under this condition, the first-order resonance frequency is about ω = 1 150 rad · s −1 . The complicated vibration region composed of periodic, quasi-periodic, and chaotic motion will occur in the lower speed region. Based on the above analysis, it can be concluded that multiple contact layers could weaken the global stiffness of the rotor system, reduce the natural frequency, and lead to nonlinear vibration of the system at low rotational speed, which is not conducive to smooth operation of the system to a certain extent.

Effect of fixed-point rubbing
In the above work, fixed-point rubbing considering coatings effects is temporarily ignored. When the amplitude of whirling motion of the rotor system is larger than the rotor-stator initial gap, fixed-point rubbing happens synchronously. Thus, the rotor-stator rub-impact is taken into consideration, and the rub-impact response differences between integral rotor and rod fastening rotor at different speeds are under contrastive analysis. Figure 13 shows the rub-impact state of integral rotor and rod fastening rotor at different rotational speeds. At the interval of (300, 580) rad · s −1 , the rub-impact fault never occurs due to the smaller whirling motion. With the increase in rotational speed, the whirling amplitude of disc becomes larger than the rotor-stator initial gap and fixed-point rubbing happens. When the value of the rotational speed is smaller than the first-order resonance frequency, the existence of the rough contact layer can strengthen the rub-impact degree. However, as the rotational speed continues to increase, more serious rub-impact occurs in the integral rotor system. Actually, these phenomena are closely related to the resonance characteristics of the rotor system, considering the effects of the contact layer in a disc group. Due to the disc-disc contact, the first-order critical speed of the rotor system tends to decrease, that is, the peak in the amplitude-frequency response of the system shifts to the left and the amplitude increases. Under this circumstance, the rod fastening rotor has the most serious rub-impact when the rotational speed is relatively low. A partial enlarged figure for describing the critical condition of 'no-rubbing to rubbing' is given in Fig. 13(b). It proves that due to the rough contact layer, the rub-impact is more likely to happen at the lower rotational speed.

Conclusions
In this paper, a dynamic model of rod fastening rotor-bearing coupling system is established, in which the contact mechanism of the rough joint interface between any two discs is analytically described through the fractal theory. After that, the nonlinear dynamic characteristics of the rod fastening rotor system are analyzed in terms of bifurcation diagram, time history, and Poincaré map. After that, the parameter analysis is carried out around the contact degree of the joint interface, the distribution position of contact layer, and multiple disc contacts. Finally, aiming at the fixed-point rubbing, the fault response characteristics of rod fastening rotor and integral rotor are analyzed. According to the calculation results, the following conclusions can be obtained.
(I) Different from the integral rotor system, the resonance characteristics of the rod fastening rotor system are affected not only by the structural design and support boundary, but also by the contact condition of the joint interface in the disc set. Due to the contact behavior of the joint interface, the first-order resonance frequency decreases from ω = 1 430 rad · s −1 to ω = 1 340 rad · s −1 .
(II) The existence of contact layer could lead to the occurrence of complicated nonlinear vibration, and the rotational speed range of nonlinear vibration is closely related to the distribution position and amount of contact layer.
(III) The contact behavior of the joint interface could promote the occurrence of fixed-point rubbing at a lower rotational speed, and aggravate the degree of rotor-stator rub in the range of (300, 1 300) rad · s −1 .