Bifurcation and Stability Analysis of a Bolted Joint Rotor System Contains Multi-Discs Subjected to Rub-Impact Effect

: In aero-engines, the rotor systems are frequently designed with multistage discs, in which the discs are fastened together through bolted joints. During operation, rotating machines are susceptible to rotor–stator rubbing faults. Those bolted joints are subjected to friction and impact forces during a rubbing event, leading to a dramatic change in mechanical properties at the con-tacting interfaces, inﬂuencing the rotor dynamics, which have attracted the attention of scholars. In the present work, a mathematical model, which considers the unbalance force, rotor dimensional properties, nonlinear oil-ﬁlm force and rub-impact effect, is developed to study the bifurcation and stability characteristics of the bolted joint rotor system containing multi-discs subjected to the rub-impact effect. The time-domain waveforms of the system are obtained numerically by using the Runge–Kutta method, and a bifurcation diagram, time domain waveforms, spectrum plots, shaft orbits and Poincar é maps are adopted to reveal the rotor dynamics under the effect of the rub-impact. Additionally, the inﬂuences of rubbing position at the multi-discs on rotor dynamic properties are also examined through bifurcation diagrams. The numerical simulation results show that the segments of the rotating speeds for rubbing are wider and more numerous, and the middle disc is subjected to the rub-impact. When the rub-impact position is far away from disc 1, the rubbing force has little effect on the response of disc 1. The corresponding results can help to understand the bifurcation characteristics of a bolted joint rotor system containing multi-discs subjected to the rub-impact effect.


Introduction
In aero-engines, the rotor systems are frequently designed with several substructures containing multiple adjacent discs and drums because of the assembly requirements [1][2][3][4]. The bolted joints are used to fasten adjacent stage discs together, and the purpose is to bring the substructures together for the transmission of power and to achieve appropriate stiffness [5][6][7][8]. The fastened discs are distributed on the circumference of the aero-engine compressor, referred to as the multi-disc bolted joint, as shown in Figure 1. The local stiffness at the contact interface is nonlinear due to the nonlinear contact at the mating interface of the bolted joint. Therefore, the dynamic analysis of the bolted joint rotor system should consider the effect of the bolted joint. Over the past decade, aero-engine manufacturers have focused their efforts on reducing rotor-stator clearance to improve efficiency [9,10]. Consequently, rotor-stator rub-impact has become one of the most common failures which may cause a serious accident. To ensure the safe operating of aero-engines, it is essential to investigate the stability characteristics of the bolted joint rotor system containing multi-discs subjected to the rub-impact effect in advance. However, according to the authors' knowledge, few reports can be found on rotor dynamic analyses under rub-impact considering the local nonlinear contact interface of the multi-disc bolted joint. Detailed studies on nonlinear behaviors of bolted joints within rotating shafts and the effect of these structures on rotor dynamics have been started in recent years [11]. Hu et al. [12] developed a rod fastening rotor system model considering the nonlinear oil-film force and the initial deflection caused by unbalanced preload. Then, the effect of the unbalanced preload on the parametric instability of the rotor system was studied by numerical simulation. Hei et al. [13] established the motion equations of the rod fastening rotor system to investigate the effect of speed, eccentricity and bending stiffness on the nonlinear dynamics of the rotor system. The gyroscopic effect and bending stiffness at the contact surface was considered in the above works to improve simulation accuracy. Wang et al. [14,15] developed a mathematical model to obtain the effect of the unbalance of multiple discs and the effect of damping on rotor dynamics, where the contact effects at the mating interface were simulated by a spring with cubic stiffness. The influence of the preload, structure parameters, contact surface and uneven preload on the nonlinear dynamic characteristics of the bolted joint rotor system was numerically studied by Wu et al. [16]. Sun et al. [17] determined the influence of joint interface assembly parameters on the natural frequency of the high-pressure rotors of aero-engines through finite element simulation and experimental studies. Jin et al. [18] proposed a data-driven modeling method referred to as the frequency sweep modeling method for bolted joint rotor systems, and a Weighted Extended Forward Orthogonal Regression algorithm was also proposed for the bolted joint structure [19]. Li et al. [20] used a mathematical model of a rod fastening rotor system deduced from Lagrange's equations, which consider the non-uniform nonlinear stiffness at the mating interface, to investigate the nonlinear dynamic properties of the rotor system. The problem of dynamic modeling and response analysis on discontinuous shafts connected by the bolted disc-drum joint was solved by Qin et al. [21,22] for the well-connected case and bolt-loosening scenario. The finite element method was adopted to gain insight into the piecewise linear stiffness of the bolted disc-drum joint, and its effect on rotor dynamics was studied in terms of a simple rotor with a bolted disc-drum joint.
In recent years, many researchers have focused on the dynamic characteristics of rotor systems under the rub-impact effect. Ma et al. [23,24] established the dynamic model of a simple rotor system considering rotor-stator rubbing fault based on the finite element modeling method to investigate the rubbing vibration features. A dedicated experimental test rig to simulate rotor-stator frictional characteristics was used for the validation of the numerical simulation results. It was observed that the period-1 motion, period-2 motion, and period-3 motion appear with increasing operating speeds from numerical and experimental study results. Zhao et al. [25,26] conducted a numerical study on the dynamic behaviors of a magnetic liquid double suspension bearing system with a clearance rubbing fault. In [27], the radial displacement responses were measured from the rotating shaft, and the results showed that the rubbing frequency was equal to the product of the fundamental frequency and the number of blades. Yang et al. [28] investigated the effects of the fixed point rubbing on dual-rotor systems numerically and experimentally. Several combination frequencies were observed, which can be thought of as the particular fault frequencies.
Based on a similar dual-rotor structure, Prabith et al. [29] established a dynamic model to investigate the bifurcation behaviors on rotor partial rub considering the synchronous forward and backward precession. Krishna et al. [30] numerically examined the nonlinear dynamics of the Jeffcott rotor, having physical contact between the rotor and the stationary parts. The harmonic balance method was adopted to investigate the nonlinear dynamic properties, and an experimental study was conducted through a rotor-stator rub-impact test rig. Phadatare et al. [31] investigated the complex bifurcation behavior of a large deflection rotor-bearing system subjected to rubbing fault by elucidating the route to chaos phenomena. Mokhtar et al. [32] studied the dynamic response features of a simple rotorcasing system undergoing rubbing faults. It was found that the sub-harmonic frequency components and the harmonics of fundamental frequency are present in the spectrums during the rotor-stator rubbing occurrence.
It can be seen from the above research that, although many scholars have begun to study the dynamic properties of the bolted joint rotor system and the nonlinear vibrations of a rotor-bearing system undergoing rubbing fault, parametric studies about bolted joint rotor systems under rub-impact are still not enough. Moreover, the multi-disc bolted joint structure in the aero-engine compressors has also not been considered in the existing studies. The existing works about the rubbing characteristics of bolted joint rotor systems mainly focus on the effect of contact stiffness and radial clearance on the dynamic properties of the rotor systems [33,34]. The present work is focused on the structural characteristics of the compressor in the aero-engine and aims to reveal the effects of rubbing position on the rotor dynamics of the bolted joint rotor system containing multi-discs. In addition, the results can help to understand the distinguishable characteristics for the position of rubbing fault that occurs in a multi-disc bolted joint rotor system. Therefore, to provide further insight into the bifurcation phenomenon of a multidisc bolted joint rotor system under the rub-impact effect, a dynamic model with twelve degrees of freedom is established in the present work. For this purpose, the dynamic model is created by employing the lumped parameter based on Lagrange's equations, where the contact effect of the mating interface in the multi-disc bolted joint is included in the modeling and analysis. The impact of the rubbing position on the dynamic behaviors of the multi-disc bolted joint rotating system is also demonstrated by bifurcation diagrams, time-domain vibration waveforms, frequency spectrums, shaft orbits and Poincaré maps. This paper is focused on the structural characteristics of the compressor in the aero-engine, and the research is based on a practical background. The novelties of the present work are given as follows: (i) The dynamic model of a bolted joint rotor system containing multi-discs subjected to rubbing fault is developed. (ii) The influence of rubbing position on the dynamic characteristics of a bolted joint rotor system with multi-discs is studied. (iii) The research can provide a reference for the recognition of the rubbing position in the rotor system of an aero-engine.
The structure of the rest of the paper is organized as follows. The governing equations of motion of the rotor-bearing system with multi-disc bolted disc joints are established in Section 2. The rub-impact model and non-dimensional model are also derived. Next, the established model is validated by comparing the simulation results with the other literature in Section 3. Moreover, a parametric analysis of the rubbing position in multi-disc bolted joints is performed in Section 4. Finally, some main conclusions are summarized in Section 5.

Dynamic Modeling of a Multi-Disc Bolted Joint Rotor-Bearing System
A schematic of a rub-impact multi-disc bolted joint rotor system supported by oillubricated bearings is shown in Figure 2, where the discs are fastened together by several long bolts. The lumped mass at two bearings is defined as m A and m B , respectively. m O1 , m O2 , m O3 and m O4 are the lumped mass of the four discs, and the eccentricities of the discs are defined as e 1 , e 2 , e 3 and e 4 , respectively. To study the dynamic properties of the Processes 2022, 10, 1763 4 of 17 rotor system efficiently, the dynamic model of the multi-disc bolted joint rotor system is simplified according to the assumptions shown below: (a) Rigid discs and the bearings are modeled through lumped mass points, and each point has two lateral degrees of freedom along with x and y axes. The torsional and axial vibration of the rotor is neglected; (b) The bearings at both sides are identical and modeled based on the Capone nonlinear oil-film force model, and the shafts are massless; (c) The rub-impact between the discs and stator is elastic, and the tangential force generated by the rotor-stator rub-impact meets the Coulomb friction law.

Nonlinear Oil-Film Force Model
The Capone model is widely used for calculating the nonlinear oil-film force, which has been proved as an effective model [33,[35][36][37]. Based upon the short bearing theory and Reynolds equation, the dimensionless oil-film pressure can be described by [38] where X and Y are dimensionless displacements at the bearing; Z denotes the axial displacement of the bearing; p is the dimensionless pressure; h is the thickness of the oil film; and R and L are the radius and width of the bearing, respectively. The dimensionless pressure p can then be derived through integrating Equation (1), which can be expressed by: Then, the dimensionless oil-film force model can be derived by integrating Equation (2) along with the lubricated arc of bearing: where V(X, Y, α), S(X, Y, α), G(X, Y, α) and α can be calculated by the following equations: S(X, Y, α) = X cos α + Y sin α Processes 2022, 10, 1763 Finally, the dimensional oil film force can be written as follows [35]: where σ is the modified Sommerfeld number, which can be calculated by [39,40]: where η represents the viscosity of lubricating oil; ω denotes the rotating speed of the rotor system; and L, R and c represent the bearing length, radius and radial clearance, respectively. Figure 3 shows the schematic diagram of the disc-stator rub-impact phenomenon with clearance. In this paper, the rotor-stator rubbing fault is a single-point rub-impact, where the tangential friction force meets the Coulomb friction law, and the radial impact meets the linear elastic deformation theory. Therefore, the radial impact force P N and the tangential rub force P T can be represented by [41]:

Rotor-Stator Rub-Impact Force Model
where k c represents the radial stiffness of rotor-stator contact; η denotes the friction coefficient; and r represents the radial displacement of the disc, which can be calculated by: According to Figure 3, the rubbing force in the x and y direction (P x and P y ) decomposed by P N and P T can then be expressed as: Processes 2022, 10, 1763 6 of 17 By substituting Equation (10) into Equation (12), the rub-impact forces in the x and y direction can be rewritten as follows:

Motion Equations of the Rotor System
According to Figure 2, the kinetic energy of the multi-disc bolted joint rotor system can be expressed by: where m i (i = A, B) represents the lumped mass at the left and right bearings; m Oi (i = 1, 2, 3, 4) represents the lumped mass at the disc centers; x i , y i (i = A, B) represent the displacements at the left and right bearings in x and y directions; and x Oi , y Oi (i = 1, 2, 3, 4) denote the lateral displacements of the discs in the x and y directions. The potential energy of the rotor system can be written as: where k is the stiffness of the shafts, and g represents the gravitational acceleration. The expression of the dissipation function can be obtained as: where c A and c B represent the linearized damping coefficients of bearings; c O1 and c O4 represent the damping coefficients introduced by the left and right shafts, respectively; and c f denotes the damping coefficient at the mating interface. Based upon the kinetic energy, elastic potential energy, gravitational potential energy and the dissipation function expressed by Equations (14)- (17), the motion equations of the rotor system can be obtained based on Lagrange's equations. According to the literature [42,43], Lagrange's equation is expressed as: where L = V − U − U g ; q i represents the generalized coordinate of the rotor system; the symbol "·" refers to the differentiation with respect to the time t; and F i is the external force. By substituting Equations (14)- (17) into Equation (18), the motion equations of the multi-disc bolted joint rotor system are deduced as follows: ..
Processes 2022, 10, 1763 y A represent the nonlinear oil film force in the x and y directions, respectively; P x and P y are the rub-impact force in the x and y directions; and F uxi , F uyi (i = 1, 2, 3, 4) denote the unbalanced forces of the discs in the x and y directions, which can be expressed by: where ω denotes the rotational speed of the rotor system. The restoring forces at the contact layer between the discs are:

Non-Dimensional Model
In order to avoid excessive truncation errors, the dimensionless transformations of Equation (19) are given as follows: By substituting Equation (22) into Equation (19), the dimensionless motion equations of the rotor system can be written as: where F u xi , F u yi (i = 1, 2, 3, 4) denote the dimensionless unbalanced forces of disc i in the x and y-direction, which can be calculated by the following equation: Furthermore, the rub-impact force relative to the dimensionless displacement X 1 , Y 1 can be expressed as: Processes 2022, 10, 1763 8 of 17

Model Verification
To prove the accuracy of the established model in Section 2.4, the steady-state responses of disc 1 in the x direction at ω = 1750 rad/s are obtained by using the Runge-Kutta method by neglecting the first 150 cycles of the simulation results. The system parameters used for model verification are given in Table 1, where the quality and the stiffness are defined to be the same as those in the literature [34]. The numerical integration result is compared with the result given in the literature by Hu et al. [34], which is shown in Figure 4. Table 1. Structure parameters of the multi-disc bolted joint rotor system for model verification.

Physical Description Parameters Values
Lumped masses (kg)  It can be found from Figure 4 that the time domain waveform, Poincaré map, shaft orbit and frequency spectrum obtained in this section have good consistency with the literature by Hu, thus proving that the dynamic model of the multi-disc bolted joint rotor system established in this paper is correct. Therefore, further numerical analysis can be conducted in the following sections.
Moreover, it should be mentioned that the rotor-bearing model studied in the present paper is slightly different from that in the literature [34], so the obtained results are slightly different. However, as the parameters of the rotor system are defined as consistently as possible with those in the literature [34], the obtained results can still prove the accuracy of the established model.

Numerical Simulation Results and Discussion
In the present section, the numerical simulation of a multi-disc bolted joint rotor system is carried out considering the disc casing rubbing fault. The numerical integration is conducted by using the fourth-order Runge-Kutta method, and the time step for numerical integration is defined as 2π/ω/512. Time-domain waveforms of disc 1 in the x direction are collected by neglecting the first 150 cycles of the transient responses. Bifurcation and stability response behaviors of the rotor system under rubbing fault are illustrated by the bifurcation diagram, time-domain response, Poincaré map, shaft orbit as well as frequency spectrums. Moreover, the effect of the rubbing position in the multi-disc bolted joint is also studied. The parameters of the rotor system are given in Table 2. Table 2. Structure parameters of the multi-disc bolted joint rotor-bearing system for numerical simulation.

Physical Description Parameters Values
Lumped

Bifurcation and Stability Analysis of the Multi-Disc Bolted Joint Rotor System
A bifurcation diagram is a useful tool used to study the dynamic properties of nonlinear dynamic systems. To demonstrate the dynamic characteristics of the multi-disc bolted joint rotor system, a bifurcation diagram with rotational speed varying within ω = [400~2000] rad/s, as the control parameter, is obtained, as shown in Figure 5. As the bifurcation diagram shows, there is only one point in the bifurcation diagram with a rotational speed below 544 rad/s, which means that the motion state under the current working condition has a synchronous vibration with period-1. To demonstrate the evolution process of the rotor system under rub-impact in-depth, the time-domain responses, shaft orbit, Poincaré map and spectrums with certain operating speeds are shown in Figures 6-13. In Figure 5, the system turns into period-2 motion at ω = 546 rad/s, and only two points are illustrated. In Figure 6, there are two points in the Poincaré map, and the fundamental frequency and the 1/2 subharmonic are the main components. It should be mentioned that the 1/2 subharmonic is introduced by the oil whirl, which results in the motion state of the rotor system entering period-2 motion. With the increase in rotational speed, the rotor turns into period-1 motion from period-2 motion at ω = 680 rad/s (shown in Figure 7). In Figure 7a, there is only one point in the Poincaré map and one peak amplitude in the frequency spectrum. Moreover, the shaft orbit is oval at the current rotational speed. As the rotational speed increases, the rotor system enters the chaotic motion from period-1 motion at ω = 754 rad/s. This is because the oil film oscillates violently at the current speed, resulting in an instable motion state, and the system enters a chaotic motion state. It can be found in Figures 5 and 8 that the rotor system motion is chaotic due to pieces of dispersed points coming forth in the Poincaré maps (see Figure 8d). As the rotational speed reaches 1046 rad/s, the motion state becomes period-2, as shown in Figure 9, where two points appear in the Poincaré map, and two points are correspondingly present in the bifurcation diagram (see Figure 5). The root cause lies in the fact that the oil-film oscillations weaken, and the oil whirl occurs at the current speed. Furthermore, it should be noted that the shaft orbits do not exceed rubbing boundary, which means that the rubbing fault does not occur corresponding to the above rotational speeds.       Processes 2022, 10, x FOR PEER REVIEW 13 of 18 Figure 10 shows the vibration waveform, rotor orbit, frequency spectrum and Poincaré map at ω = 1184 rad/s, and it shows that the system turns into period-6 motion, which can be verified by six isolated points in the Poincaré map and by the fundamental frequency and some subharmonics in the frequency spectrum, as shown in Figure 10c. This is because the oil whip occurs at the current speed, which leads to the subharmonics appearing in the frequency spectrum. Moreover, the rotor orbit exceeds the rubbing boundary, as shown by a red circle, which means that the rubbing fault occurs at the current speed. As speed continues to increase, the motion state turns into period-3 motion at ω = 1190 rad/s. There are three points illustrated in the Poincaré map (Figure 11d), and the fundamental frequency and its 1/2 and 1/3 subharmonics can be observed in the frequency spectrum because the oil whip phenomenon still occurs at the current speed. (Figure 11b). In Figure 12d, there are many points of disordered distribution in the Poincare map at ω = 1308 rad/s, indicating that the system enters chaotic motion. Therefore, the rotor system presents the period-3 motion state within the interval of 1190 rad/s < ω < 1306 rad/s. When the rotational speed reaches 1394 rad/s, inverse period-doubling bifurcation occurs, and the system leaves chaotic motion and turns into period-1 motion. Moreover, the bifurcation diagram depicted in Figure 5 shows that only one point is illustrated within 1394 rad/s < ω < 1660 rad/s. Finally, the system enters chaotic motion from period-1 motion at ω = 1662 rad/s, which is demonstrated by numerous disordered points in the Poincaré map and continuous frequency bands in the spectrum, as shown in Figure 13. This is because of the severe vibration introduced by the rubbing fault and the continuous frequency bands that appear in the frequency spectrum, resulting in the system entering a chaotic motion state.
From the results, we can see that the routes to the chaos of the rotor system under the rub-impact effect can be summarized as follows: period-1 motion, period-2 motion, period-1 motion, chaotic motion, period-1 motion, period-2 motion, period-6 motion, period-3 motion, chaotic motion, period-1 motion, period-2 motion and chaotic motion within 400 rad/s < ω < 2000 rad/s. Furthermore, it can be known that the rotor system undergoes a rub-impact effect within 1184 ≤ ω < 1392 and 1662 ≤ ω < 2000 from    Figure 10 shows the vibration waveform, rotor orbit, frequency spectrum and Poincaré map at ω = 1184 rad/s, and it shows that the system turns into period-6 motion, which can be verified by six isolated points in the Poincaré map and by the fundamental frequency and some subharmonics in the frequency spectrum, as shown in Figure 10c. This is because the oil whip occurs at the current speed, which leads to the subharmonics appearing in the frequency spectrum. Moreover, the rotor orbit exceeds the rubbing boundary, as shown by a red circle, which means that the rubbing fault occurs at the current speed. As speed continues to increase, the motion state turns into period-3 motion at ω = 1190 rad/s. There are three points illustrated in the Poincaré map (Figure 11d), and the fundamental frequency and its 1/2 and 1/3 subharmonics can be observed in the frequency spectrum because the oil whip phenomenon still occurs at the current speed. (Figure 11b). In Figure 12d, there are many points of disordered distribution in the Poincare map at ω = 1308 rad/s, indicating that the system enters chaotic motion. Therefore, the rotor system presents the period-3 motion state within the interval of 1190 rad/s < ω < 1306 rad/s. When the rotational speed reaches 1394 rad/s, inverse period-doubling bifurcation occurs, and the system leaves chaotic motion and turns into period-1 motion. Moreover, the bifurcation diagram depicted in Figure 5 shows that only one point is illustrated within 1394 rad/s < ω < 1660 rad/s. Finally, the system enters chaotic motion from period-1 motion at ω = 1662 rad/s, which is demonstrated by numerous disordered points in the Poincaré map and continuous frequency bands in the spectrum, as shown in Figure 13. This is because of the severe vibration introduced by the rubbing fault and the continuous frequency bands that appear in the frequency spectrum, resulting in the system entering a chaotic motion state.
From the results, we can see that the routes to the chaos of the rotor system under the rub-impact effect can be summarized as follows: period-1 motion, period-2 motion, period-1 motion, chaotic motion, period-1 motion, period-2 motion, period-6 motion, period-3 motion, chaotic motion, period-1 motion, period-2 motion and chaotic motion within 400 rad/s < ω < 2000 rad/s. Furthermore, it can be known that the rotor system undergoes a rub-impact effect within 1184 ≤ ω < 1392 and 1662 ≤ ω < 2000 from

Influence of Rubbing Position
To study the effect of the rubbing position on the system response characteristics of the rotor system, in this section, the bifurcation diagrams of disc 1 for the working condition of disc 2 subjected to rub-impact are obtained, which are shown in Figures 14-16. By comparing Figures 14 and 15, with numerical integration results illustrated in Figure 5, it can be observed that the dynamic behavior evolution process of the rotor system is different from that of the case of disc 1 subjected to rub-impact, which indicates that the rubbing position in the multi-disc bolted joint has a significant effect on the dynamic characteristics and rubbing speed range of the rotor system. To clearly show the impact of the rubbing position on the system response characteristics, the speed ranges of the system where rubimpact occurs and the speed ranges of chaotic motion are listed in Table 3. It can be found that the rubbing force acts on the rotor system at 1256 rad/s when disc 2 and disc 3 are subjected to rub-impact, which is higher than the case of disc 1 subjected to rubbing impact. This is because the nonlinear contact force introduced by the contact effect at mating interfaces suppresses the vibration of the intermediate disc in the multi-disc bolted joint, leading to a higher speed of rubbing occurrence. Moreover, the segments of the speed range for rubbing are more and wider for the case of disc 2 and disc 3 subjected to rubbing fault than that of disc 1 subjected to rub-impact, and this is caused by the coupling effect of the rubbing force and the contact force at the mating interfaces. Regarding the effect of rubbing position on the motion stability of the rotor system, it can be found in Table 3 that the segments of the chaotic motion speed range for the case of disc 2 and disc 3 subjected to rub-impact are more than the scenario of disc 1 subjected to rub-impact. Furthermore, the system enters chaotic motion at an earlier speed for the case of disc 2 and disc 3 subjected to rub-impact, which is caused by the coupling effect of the rubbing force and nonlinear contact force at mating interfaces.

Influence of Rubbing Position
To study the effect of the rubbing position on the system response characteristics of the rotor system, in this section, the bifurcation diagrams of disc 1 for the working condition of disc 2 subjected to rub-impact are obtained, which are shown in Figures 14-16. By comparing Figures 14 and 15, with numerical integration results illustrated in Figure 5, it can be observed that the dynamic behavior evolution process of the rotor system is different from that of the case of disc 1 subjected to rub-impact, which indicates that the rubbing position in the multi-disc bolted joint has a significant effect on the dynamic characteristics and rubbing speed range of the rotor system. To clearly show the impact of the rubbing position on the system response characteristics, the speed ranges of the system where rubimpact occurs and the speed ranges of chaotic motion are listed in Table 3. It can be found that the rubbing force acts on the rotor system at 1256 rad/s when disc 2 and disc 3 are subjected to rub-impact, which is higher than the case of disc 1 subjected to rubbing impact. This is because the nonlinear contact force introduced by the contact effect at mating interfaces suppresses the vibration of the intermediate disc in the multi-disc bolted joint, leading to a higher speed of rubbing occurrence. Moreover, the segments of the speed range for rubbing are more and wider for the case of disc 2 and disc 3 subjected to rubbing fault than that of disc 1 subjected to rub-impact, and this is caused by the coupling effect of the rubbing force and the contact force at the mating interfaces. Regarding the effect of rubbing position on the motion stability of the rotor system, it can be found in Table 3 that the segments of the chaotic motion speed range for the case of disc 2 and disc 3 subjected to rub-impact are more than the scenario of disc 1 subjected to rub-impact. Furthermore, the system enters chaotic motion at an earlier speed for the case of disc 2 and disc 3 subjected to rub-impact, which is caused by the coupling effect of the rubbing force and nonlinear contact force at mating interfaces.  Figure 10 shows the vibration waveform, rotor orbit, frequency spectrum and Poincaré map at ω = 1184 rad/s, and it shows that the system turns into period-6 motion, which can be verified by six isolated points in the Poincaré map and by the fundamental frequency and some subharmonics in the frequency spectrum, as shown in Figure 10c. This is because the oil whip occurs at the current speed, which leads to the subharmonics appearing in the frequency spectrum. Moreover, the rotor orbit exceeds the rubbing boundary, as shown by a red circle, which means that the rubbing fault occurs at the current speed. As speed continues to increase, the motion state turns into period-3 motion at ω = 1190 rad/s. There are three points illustrated in the Poincaré map (Figure 11d), and the fundamental frequency and its 1/2 and 1/3 subharmonics can be observed in the frequency spectrum because the oil whip phenomenon still occurs at the current speed. (Figure 11b). In Figure 12d, there are many points of disordered distribution in the Poincare map at ω = 1308 rad/s, indicating that the system enters chaotic motion. Therefore, the rotor system presents the period-3 motion state within the interval of 1190 rad/s < ω < 1306 rad/s. When the rotational speed reaches 1394 rad/s, inverse period-doubling bifurcation occurs, and the system leaves chaotic motion and turns into period-1 motion. Moreover, the bifurcation diagram depicted in Figure 5 shows that only one point is illustrated within 1394 rad/s < ω < 1660 rad/s. Finally, the system enters chaotic motion from period-1 motion at ω = 1662 rad/s, which is demonstrated by numerous disordered points in the Poincaré map and continuous frequency bands in the spectrum, as shown in Figure 13. This is because of the severe vibration introduced by the rubbing fault and the continuous frequency bands that appear in the frequency spectrum, resulting in the system entering a chaotic motion state.
From the results, we can see that the routes to the chaos of the rotor system under the rub-impact effect can be summarized as follows: period-1 motion, period-2 motion, period-1 motion, chaotic motion, period-1 motion, period-2 motion, period-6 motion, period-3 motion, chaotic motion, period-1 motion, period-2 motion and chaotic motion within 400 rad/s < ω < 2000 rad/s. Furthermore, it can be known that the rotor system undergoes a rub-impact effect within 1184 ≤ ω < 1392 and 1662 ≤ ω < 2000 from Figures 5,10, and 11-13.

Influence of Rubbing Position
To study the effect of the rubbing position on the system response characteristics of the rotor system, in this section, the bifurcation diagrams of disc 1 for the working condition of disc 2 subjected to rub-impact are obtained, which are shown in Figures 14-16. By comparing Figures 14 and 15, with numerical integration results illustrated in Figure 5, it can be observed that the dynamic behavior evolution process of the rotor system is different from that of the case of disc 1 subjected to rub-impact, which indicates that the rubbing position in the multi-disc bolted joint has a significant effect on the dynamic characteristics and rubbing speed range of the rotor system. To clearly show the impact of the rubbing position on the system response characteristics, the speed ranges of the system where rub-impact occurs and the speed ranges of chaotic motion are listed in Table 3. It can be found that the rubbing force acts on the rotor system at 1256 rad/s when disc 2 and disc 3 are subjected to rub-impact, which is higher than the case of disc 1 subjected to rubbing impact. This is because the nonlinear contact force introduced by the contact effect at mating interfaces suppresses the vibration of the intermediate disc in the multi-disc bolted joint, leading to a higher speed of rubbing occurrence. Moreover, the segments of the speed range for rubbing are more and wider for the case of disc 2 and disc 3 subjected to rubbing fault than that of disc 1 subjected to rub-impact, and this is caused by the coupling effect of the rubbing force and the contact force at the mating interfaces. Regarding the effect of rubbing position on the motion stability of the rotor system, it can be found in Table 3 that the segments of the chaotic motion speed range for the case of disc 2 and disc 3 subjected to rub-impact are more than the scenario of disc 1 subjected to rub-impact. Furthermore, the system enters chaotic motion at an earlier speed for the case of disc 2 and disc 3 subjected to rub-impact, which is caused by the coupling effect of the rubbing force and nonlinear contact force at mating interfaces.    The bifurcation diagram for the working condition of disc 4 subjected to rub-impact is shown in Figure 16. From Figures 5 and 16, it can be seen that, when the rub-impact position in the multi-disc bolted joint is far away from disc 1, the influence on the system response characteristics of disc 1 is not obvious. The root cause lies in the fact that, when the rubbing position is far away from disc 1, the rubbing force has little effect on the response of disc 1. In addition, the nonlinear contact force introduced by the contact effect at the mating interface between disc 1 and disc 2 is not affected by the rubbing effect, so the system response characteristics of disc 1 do not change significantly.

Conclusions
In this work, a rotor system involving the multi-disc bolted joint is modeled based on the energy theorem and Lagrange's principle, with the rubbing force acting on the disc. In addition, the accuracy of the established model is proved by comparing the simulation results of this work with the available results reported in the literature. Based on the fourth-order Runge-Kutta method, the nonlinear dynamic and bifurcation characteristics of the rotor system for the case of disc 1 subjected to rub-impact are investigated. The intermittent chaos phenomenon is observed through bifurcation diagrams, Poincare maps and spectrum plots. Moreover, by comparing the bifurcation diagrams under different rubbing positions in the multi-disc bolted joint, the dynamic properties of the rotor system and rubbing fault speed range are also examined. Some conclusions are summarized according to the parametric study as follows: 1.
The nonlinear contact force introduced by the contact effect at mating interfaces suppresses the vibration of the intermediate disc in the multi-disc bolted joint, leading to a higher speed of rubbing occurrence. Moreover, the segments of the rotating speed range for rubbing are more and wider for the case of the middle disc subjected to rub-impact.

2.
The speed segments of the chaotic motions for the case of the intermediate disc subjected to rub-impact are more than the scenario of disc 1 subjected to rub-impact. Furthermore, the system enters chaotic motion at an earlier speed.

3.
When the rub-impact position is far away from disc 1, the rubbing force has little effect on the system response characteristics of disc 1.
The corresponding results can help to understand the bifurcation characteristics of a bolted joint rotor system containing multi-discs subjected to the rub-impact effect.

Institutional Review Board Statement: Not applicable.
Data Availability Statement: Data will be made available on reasonable request.