Numerical and experimental analysis of the effect of eccentric phase difference in a rotor-bearing system with bolted-disk joint

Bolted joints are widely used in industrial rotating machines to fasten the adjacent disks together and affect system dynamic properties. Therefore, there is a strong need to study their influence on the response of such systems. This paper investigated the effect of the eccentric phase difference of the disk and bolted-disk joint on rotor dynamics, which takes into account the time-varying bending stiffness of the bolted joint. The bolted-disk joint is modeled as a two-node element based upon the energy theorem and Lagrange’s principle, where the relative lateral displacement stiffness, relative bending stiffness, and coupling stiffness between the adjacent disks are considered. And then the dynamic model of the rotor-bearing system is derived based on the proposed bolted joint element and lumped mass modeling method. Combining with the Newmark-β integration scheme, the established model allows the dynamic response characteristics of the rotor-bearing system with the bolted joint to be predicted, and the impact of eccentric phase difference in the rotor system on the response to be investigated. The validity of the simulation results was confirmed by experiment. Through the modeling method proposed in this paper and obtained results, the bifurcation characteristics of the bolted joint rotor system can be predicted.

considered. And then the dynamic model of the rotorbearing system is derived based on the proposed bolted joint element and lumped mass modeling method. Combining with the Newmark-b integration scheme, the established model allows the dynamic response characteristics of the rotor-bearing system with the bolted joint to be predicted, and the impact of eccentric phase difference in the rotor system on the response to be investigated. The validity of the simulation results was confirmed by experiment. Through the modeling method proposed in this paper and obtained results, the bifurcation characteristics of the bolted joint rotor system can be predicted.
Keywords Rotor-bearing system Á Bolted-disk joint Á Eccentric phase difference Á Instability characteristics Á Largest Lyapunov exponent

List of symbols c
Bearing clearance c blx Damping coefficient at left bearing in the x-direction c bly Damping coefficient at left bearing in the y-direction c brx Damping coefficient at right bearing in the x-direction c bry Damping coefficient at right bearing in the y-direction C Damping matrix of the rotor system 1 Introduction Rotating machines have numerous engineering applications, such as generators, turbomachinery, and compressors. Bolted joints are frequently used in these machineries, the purpose being connecting the multiple rotating components together for transmission of power, and provide sufficient strength to maintain the integrity of the rotor system [1]. Figure 1 shows the typical bolted-disk joint in a certain type of compressor, where the adjacent disks are connected by bolts distributed in the circumference. During normal operation of the rotor system, these bolted-disk joints rotate at high speeds, of which the dynamic properties have a significant influence on the rotor dynamics [2]. When investigating the dynamics of the rotors assembled by bolted joints, these joints are usually regarded as rigid connections [3][4][5][6][7]. However, the relative displacements between the adjacent disks of the bolted joint and eccentric phase difference might change the vibration characteristics of rotor systems. Moreover, the eccentric phase difference in the rotor system would contribute to unstable operation of the system, high-level vibration, and potential damage of the rotating machinery under a high rotational speed [8]. Therefore, the prediction of dynamic behaviors of the bolted joint rotors must be carried out more and more carefully at the design stage which has to take into account the relative displacements and eccentric phase difference of bolted-disk joint in order to fabricate rotating machines as reliable as possible.
Rotating equipment is supported on bearings of some form, which allow free rotation and energy transfer [9]. Bearing is one of the most significant sources of nonlinearity, caused by a periodical change in the contact position between rolling elements and races [10,11]. A lot of researches are available to solve the vibration analysis problems of a flexible shaft carrying rigid single or multiple disks supported by the bearings. Villa et al. [12] studied the dynamic stability of a flexible rotor-bearing system in the frequency domain. The results show the dependency between the critical speed and the unbalance level. Torsvik et al. [13] described the development of a generic model incorporating bearings and examined model fidelity requirements for bearings. Gupta [14] performed a detailed parametric study to investigate the effect of flexibility of the shaft, rotational speed, and radial clearance on the dynamic stiffness of the ball bearing.
Chen et al. [15] applied correlation functions of the responses to distinguish harmonic frequencies of rotating structures under operational conditions. Chen et al. [16] investigated the complex curvature mode shapes for a wind turbine rotor system through finite element simulation and experimental. In addition, some experimental studies are also demonstrated by scholars [17][18][19][20][21][22].
Throughout the years, several works have been published on the subject of nonlinear rotor dynamics considering the bolted joint. Hu et al. [23,24] developed the dynamic model of rod-fastening rotor system, where the interaction between two disks was equivalent to a contact layer with nonlinear flexural stiffness, and then analyzed the bifurcation and chaos phenomena of the system under rub-impact and initial deformation of the shaft. Hei et al. [25,26] established a 12 degrees-of-freedoms (DOFs) dynamic model of rod fastening rotor based on Lagrange's equations to investigate the nonlinear dynamic behaviors, where the disks are considered to be connected by the spring with cubic stiffness. Wang et al. [27] established a dynamic model of rod fastening rotor-bearing system based on the D'Alembert principle, in which the contact interface is regarded as a nonlinear stiffness spring. Liu et al. [28] established the finite element model of a rod-fastening rotor-bearing system to study the effect of pre-tightening force on motion stability; the research indicated that the nonlinear contact at the mating interface led to the dynamic properties of the rod-fastening rotor system quite different from the integrate rotor. As for the typical bolted joint in the aero-engine, Qin et al. [5] derived the analytical bending stiffness model of the bolted disk-drum joint, and then the difference of dynamic properties between the jointed rotor and integrate one is investigated through a simulation study. Moreover, bolt loosening effect on rotor dynamics is also conducted combined with nonlinear finite element simulation [29]; the results show that resonant peaks will take place at subcritical speeds due to bolt loosening. Sun et al. [30] built a dynamic model of the drum-disk-shaft system, where the mechanical model of bolted joints is established by characterizing the assembly mechanical relationship at the mating surface.
Up to now, many scholars have studied the dynamic behavior of the rotor system with bolted joint. However, the relative displacements between the adjacent disks (relative lateral movement, the relative angular movement, and the coupling relationship of lateral movement and angular movement between adjacent disks of the bolted-disk joint) of the bolted joint and eccentric phase difference were not fully considered. Therefore, it is very important to include this feature in a proper rotor model if more realistic results are required. This article is devoted to deriving an effective dynamic model of the bolted-disk joint rotor system, and analyze the effect of the bolted-disk joint on system response through the established model.
In rotor dynamics, continuous and discrete models are two commonly used modeling methods [31]. The derived dynamic equations are usually expressed in partial differential forms. The discretization modeling methods can be divided into three categories: the transfer matrix modeling method [32], the finite Fig. 1 Sketch of a compressor of a certain type of aero-engine: a simplified compressor, b cross-section diagram element modeling method [33][34][35][36], and one type of lumped mass modeling method [37][38][39]. Compared with the other two, the lumped mass modeling method could characterize the main features of complex rotating machines. Meanwhile, there are not too many DOFs in this kind of model, which is convenient for numerical calculation and dynamic analysis. Therefore, this paper aims at revealing the effects of the eccentric phase difference on the dynamic characteristics of the rotor-bearing system with bolted-disk joints based on the lumped mass modeling method and the proposed bolted joint element.
The rest of the present paper is outlined as follows: Firstly, the motion equation of the rotor-bearing system with a bolted-disk joint is established in Sect. 2. In Sect. 3, the validity of the established model is examined. The effects of the eccentric phase difference in the rotor-bearing system on rotor dynamics are discussed in Sect. 4. Section 5 validates some of the simulation results through a rotor-bearing test rig with multiple disks. Finally, some conclusions are summarized in Sect. 6.

Dynamic modeling
To analyze the effect of eccentric phase difference on the dynamic characteristics of the rotor system, a rotor-bearing system with a bolted joint structure is considered, as depicted in Fig. 2, where m i i ¼ 1; 2; 3; 4; 5; 6 ð Þare lumped mass points and l i i ¼ 1; 2; 3; 4 ð Þdenote the length of the shafts. In this section, the mathematical models of the rotor system are developed based on the method applied in the literature [8], and the bolted-disk joint is modeled by a proposed bolted joint element with a total of eight DOFs. The generalized coordinates of the rotor system

Rigid disk and bolted-disk joint model
This paper uses three coordinates of the second kind Euler angle to describe the movement of the rigid disk and bolted-disk joint [40,41]. The utilization of more than one frame of reference requires establishing the relationships between the different coordinate systems. The local coordinate systems P-X 1 Y 1 Z 1 , P-X 2 Y 2 Z 2 , and P-X 3 Y 3 Z 3 are established on the neutral surface of the disk, as shown in Fig. 3. P-xyz is the global coordinate system, in which O is coincident with the center of the rigid disk and the z-axis is collinear with the undeformed centerline of the shaft. P-X 3 Y 3 Z 3 is attached to the disk and rotates with the shaft. The disk rotates h Y around the y axis in the global coordinate system, and angular velocity corresponding to the local coordinate system P-X 1 Y 1 Z 1 can be written as: Then, the disk rotates h X1 around the O 0 X 1 axis; the relationship of angular velocity corresponding to the   Fig. 3 The local coordinate system for a disk local coordinate system P-X 1 Y 1 Z 1 and P-X 2 Y 2 Z 2 can be obtained as: And then, the disk rotates h Z2 around the O 0 Z 2 axis; the relationship of angular velocity corresponding to the local coordinate system P-X 2 Y 2 Z 2 and P-X 3 Y 3 Z 3 can be obtained as: According to the literature [41], the following relationship can be obtained: Finally, the angular velocity of P-X 3 Y 3 Z 3 corresponding to the global coordinate system P-xyz can be described as: where x is the rotational speed of the rotor system; h Y is the rotational angle from global coordinate system P-xyz to local coordinate system P-X 1 Y 1 Z 1 ; h X1 is the rotational angle from local coordinate system P-X 1 Y 1 Z 1 to local coordinate system P-X 2 Y 2 Z 2 ; h Z2 is the rotation angle from local coordinate system P-X 2 Y 2 Z 2 to local coordinate system P-X 3 Y 3 Z 3 ; A 0 , A 1 and A 2 are the transformation matrix from local coordinate system P-X 3 Y 3 Z 3 to the global coordinate system P-xyz; the specific expressions are given in ''Appendix A''. Equivalent mass, polar moments of inertia, and diametral moments of inertia for the disk are defined as m 3 , J p3 , and J d3 ; a lump mass point with four DOFs (q d ¼ ½x 3 ; y 3 ; h x3 ; h y3 T ) was used to represent it. Thus, according to the literature [41][42][43], the kinetic energy T D is composed of the sum of the translational kinetic energy and the rotational kinetic energy introduced by the rotating around the x, y, and z axes, i.e., According to the literature [41] and the global coordinate system depicted in Fig. 2 And then substituting Eq. (4) into Eq. (5) and without considering the high order terms, the kinetic energy of the disk can be rewritten as follows: The Lagrange's equations can be written as [44] d dt where q i represents the generalized coordinates; T and U are the kinetic energy and strain energy, respectively; F qi is the external force corresponding to q i ; the symbol '' •'' refers to the differentiation with respect to the time. The disk is considered as a rigid body in the present paper, and only the kinetic energy is taken into account. Without the strain energy and external excitation forces, Lagrange's equation can be written as [42,43]: After inserting Eq. (6) into Eq. (8), the application of Lagrange's equations to the resulting kinetic energy T D yields the differential equations of motion of the rigid disk described by: where M d , G d , Q u d and Q g d are mass matrix, gyroscopic matrix, unbalanced force vector, and gravity vector of the rigid disk, and the specific expressions are written in ''Appendix B''.
The bolted-disk joint is treated as a two-node element with a total of eight DOFs (see Fig. 4), where k hx , k y and c y are the rotational stiffness of the bolted-disk joint while rotation around x-axes, lateral stiffness along with the y-axes, damping coefficients introduced by the lateral movement of the disk. The generalized displacement vector of the proposed jointed element is given as: It can be found that the disks are symmetric with respect to the neutral plane of the bolted-disk joint, the relationships of stiffness and damping coefficients can then be written as and c x ¼ c y ¼ c s . So the kinetic energy and strain energy of the bolted joint element can be written as: where M J , G J , and K J are mass, damping, gyroscopic, and stiffness matrices of the jointed element, respectively; Q u J and Q g J are generalized unbalanced force vector and gravity vector, respectively. The matrices are expressed as shown in ''Appendix C''.

Rotor-bearing model
In order to study the nonlinear dynamic characteristics of the flexible rotor-bearing system with bolted-disk joint, a mathematical model of the rotor system depicted in Fig. 2 is derived according to the following assumptions: 1. The disks are treated as rigid ones, and the deformation is ignored. 2. The displacements of the rotor system in axial directions are negligible, and each lumped mass point has two lateral displacements (translation along the x-axis and y-axis) and two angular displacements (rotation around the x-axis and yaxis) DOFs, respectively. 3. The lumped mass points corresponding to the bolted-disk joint are regarded as the bolted joint element, and the others are connected by massless shaft sections.
In this section, the lumped mass modeling method [8,38] was used to derive the motion equations, where the stiffness of the shafts is determined based on the flexibility influence coefficient method in mechanics  Fig. 4 Sketch of the bolteddisk joint element in yzplane: a the bending stiffness, b the lateral stiffness of materials [37]. The overall rotor system with 24 DOFs can be derived as follows: where Q b is the bearing force vector of the rotor system; Q u and Q g are unbalanced force vector and gravity vector, respectively; C is the damping matrix of the global system; q represents the generalized displacement vector of the global system, which is given as follows: M, G, and K are mass matrix, gyroscopic matrix, and stiffness matrix of the rotor system, and The calculation of matrix elements of K x and K y is given as follows: In practical engineering, the damping of engineering structures is mostly assumed by the Rayleigh damping theory, where the damping matrix can be obtained by superposition of mass matrix and stiffness matrix [37]. This damping matrix construction method has a lot of numerical advantages and can meet the needs of rotor dynamic analysis.
In this paper, the Rayleigh damping form is adopted to determine the overall damping matrix (C), which can be expressed by the following [45]: where a and b are the Rayleigh damping parameters [46]. The external force vector Q b is formed by bearing force, which can be expressed as: As for the ball bearings, the nonlinear bearing force generated in the vertical and horizontal directions can be calculated by the following [47][48][49] where k B is the contact stiffness; N b is the number of balls; c is the radial clearance of the bearing; h j is the angle location of the jth rolling ball, which can be calculated by: where x c is the rotational speed of the bearing cage, x c ¼ x Â r=ðR þ rÞ; R and r are the outer and inner radius of the bearing, respectively. According to the literature [10], the so-called varying compliance (VC) frequency can be described by The VC vibration is the inherent characteristic of the bearings and Eq. (28) can be applied for the ; identification of the nonlinearity introduced by the bearing in the rotor dynamic analysis [17].

Non-dimensional model
For the convenience of the solution and avoid excessive truncation errors, dimensionless transformations of Eq. (14) are performed by the following: where c is bearing clearance. Substituting Eq. (29) into Eq. (14), the dimensionless form motion equations of the rotor system can be written as: whereQ b is the dimensionless bearing force vector, ; 0 T ;Q e is the dimensionless unbalanced force vector, which can be expressed as: And the dimensionless bearing force in x and ydirection can be calculated by the following: here,h

Model verification
In order to verify the established dynamic model of the rotor-bearing system and proposed bolted joint element, in this section, the simulation results of the two cases were compared with the available results reported in the literature [10,37]. In the present paper, the Newmark integration method is adopted due to the unconditionally stable and robustness in the numerical simulation [8,50].
3.1 Case 1: Comparing with the simple rotorbearing system in literature [10] To verify the modeling method, the steady-state vibration waveforms at low operating speeds are generated for a rotor system (see Fig. 5) and compared with the results available in the literature [10] for a similar system. The lumped mass model of the rotor system with a bolted-disk joint is derived as follows: whereM,C,G,K, andq are mass matrix, damping matrix, gyroscopic matrix, stiffness matrix, and generalized displacement vector of the global system;Q u , Q b , andQ g are unbalance force vector, bearing force vector, and gravity vector, respectively (detail in ''Appendix D'').
According to the available researches [3,51], the bending stiffness of the bolted joint structure presents piecewise linear properties, and the transition takes place when the external force equal to preload [52]. The relationship between bending stiffness k h and relative rotation angle at the transition point U 0 can be described by [46]: where k h1 is the bending stiffness at the first bending stage;k h2 is the bending stiffness at the second bending stage;U represents the relative angular of the bolteddisk joint and can be described by: The bending stiffnesses at two bending stages are defined to be the same in case 1, where k h1 ¼ k h2 ¼ 8 Â 10 6 Nm, and the lateral stiffness of the bolted joint is specified as k S ¼2 Â 10 9 N/m. The coupling stiffness of the bolted joint is set as:k 1 ¼ k 2 ¼ 2:17 Â 10 12 N, and k 0 h ¼ 4:34 Â 10 15 N. Since the dynamic responses obtained through the lumped mass model are to be compared with a rigid model in [10], the length of the shafts is reduced to 0.05 m to achieve an excellent matching of results. The other detailed parameters of the rotor and bearing are listed in Tables 1 and 2, respectively.
The horizontal and vertical displacements under 300 rpm are compared, as shown in Fig. 6, where a good consistent can be seen of the time-domain waveforms in both x and y directions. And then the introduced bearing force can be proved to be accurate, but the accuracy of the established rotor model with bolted joint still needs further verification.
3.2 Case 2: Comparing with the rotor system in the literature [37] To further validate the established model and proposed bolted joint element, the steady-state vibration waveforms corresponding to multiple speeds of the rotor system as shown in Fig. 7 are compared with the results in the literature [37]. The bearing in the left is simulated by a spring-damping model, and the verification needs to replace the nonlinear bearing force in Eq. (30) with the nonlinear oil film force. Then the non-dimensional bearing forceQ b should be replaced by the following: where g represents the viscosity of lubricating oil. L, D, and c are sliding bearing length, diameter, and radial clearance, respectively. And where the functions V, S, G, and a are given in Eq. (41): Physical dimensions are shown in Fig. 7, and the other model parameters are listed in Table 3.
The bending stiffness at two bending stages is defined as k h1 ¼ k h2 ¼ 8 Â 10 6 Nm, and the lateral stiffness of the joint k S ¼2 Â 10 15 N/m. The coupling stiffness of the bolted joint is set as:k 1 ¼ k 2 ¼ 6:185 Â 10 10 N, and k 0 h ¼ 7:93 Â 10 6 N. The vibration waveforms of the rotor system are obtained using the Newmark method and compared with the results in the literature [37], as shown in Figs. 8 and 9. It can be seen that the timedomain waveforms have a good consistent with the results in the literature [37], and then the accuracy of the established model can be proved. It needs to be emphasized that in numerical simulation case 1, the speed is low, which cannot fully validate the modeling method and the derived bolted joint element. But it can be proved that the calculation of bearing force is correct. And in case 2, the bolted-disk joint element and established model can be verified by comparing the steady-state responses of the system with the results in the literature [37] under multiple speeds. Moreover, the modeling method in literature [10] is different from that in this paper; therefore the length of the left and right shafts is reduced in steps until an excellent  Fig. 7 Sketch of the rotor system with a bolted-disk joint  matching of results is achieved. And the total mass of the rotor is maintained equal to that taken by Chen [10]. In numerical validation case 2, the modeling method of literature [37] and present paper are the same. Therefore, a good consistent of the time-domain waveforms can be seen. The above reference works did not consider the bolted joint, so the reported results only used for model validation in this Section.

Numerical simulation and discussion
A parametric study of the effect of the eccentric phase difference in the rotor-bearing system on motion stability is performed in this section. The dynamic characteristics of the rotor system are investigated through bifurcation diagrams, Largest Lyapunov exponent (LLE), and frequency spectrums. Moreover, according to the literature [29], the bending stiffness of the bolted joint will not change significantly while there are some loosened bolts, which can be regarded as an extreme uneven pre-tightening force situation. Therefore, the influence of uneven pre-tightening force is not considered in our work.

Bifurcation and nonlinear response analysis
According to the literature [37], the two unbalanced loading conditions corresponding to the in-phase and out-of-phase eccentricities are shown in Fig. 10. This Vibration waveforms of the rotor system in the present paper under multiple speeds: a 2800 rev/min, b 4400 rev/min, c 6000 rev/ min section presents results for the chaotic response of the bolted joint rotor system induced by eccentric phase difference corresponding to these two unbalanced loading conditions. By comparing and analyzing the difference between the two cases, the influence of the eccentric phase difference between the disk and the bolted-disk joint is investigated. Figure 2, Tables 2, and 4 summarize the parameters of the structural parameters of the rotor system and journal bearings used in this study. The elements in the stiffness matrix of the bolted-disk joint element are defined as:k h1 ¼ 2 Â 10 8 Nm, k h2 ¼ 2 Â 10 7 Nm, U 0 ¼ 1 Â 10 À6 rad, k S ¼1:97 Â 10 13 N/m, k 1 ¼ k 2 ¼ 9:89 Â 10 9 N, and k 0 h ¼ 3:29 Â 10 6 N. The bending stiffness is determined according to the relative rotational angle between the disks of bolted joint obtained at each numerical integration step. The stiffness matrix of the rotor system can then be obtained by introducing the bending stiffness into Eq. (18), and then the damping matrix of the rotor system can be calculated by Eq. (23). For the first loading condition as shown in Fig. 10a,  u 1 = 1.1838 9 10 -4 kgÁm, and u 2 = u 3 = u 1 /2 are inphase, for the second loading condition as shown in Fig. 10b, the eccentricities are the same as the first loading condition but out-of-phase. The dynamical simulation of Eq. (30) is carried out by using the Newmark method, where the time step is specified as 2p/512. The time response corresponding to the first 160 periods of the system was discarded from the sampled data to ensure that the responses reached steady-state conditions. The steady-state responses during the last 80 periods were employed for the plots shown below, and the bifurcation diagrams, Largest Lyapunov exponents, and spectrum cascade within x [ [600, 15000] rev/min corresponding to the inphase eccentricities are obtained, as shown in Fig. 11. Figure 11a and c shows bifurcation and LLE versus rotational speed within x [ [600, 15000] rev/min under in-phase eccentricities loading condition, which delivers rich nonlinear phenomena. The system is in a stable motion state before 2580 rev/min; the values of the LLE are always negative numbers. When the rotating speed is at the range of [2640 * 4500] rev/ min, the motions of the system go through the variation from chaotic motion state to stable motion state. The values of the LLE are 0.017 at x = 2640 rev/min, it ends at x = 2700 rev/min and turns to a stable motion state, and the corresponding values of the LLE are negative numbers within x [ [2760 * 4500] rev/min. The vibration response is presented in Fig. 12 when x = 2700 rev/min. It can be seen that the points in the Poincaré map are distributed in a disorderly state and the rotor orbit is complicated. In this paper, the ratio of VC frequency to fundamental frequency is 3.08, which can be obtained by Eq. (28) and Table 2. Therefore, the fundamental frequency, VC frequency, their combination frequency components, and the continuous frequency spectra are the main frequency components in the spectrum as shown in Fig. 12b. This is because, under the in-phase eccentricity loading condition, the bearing vibrates violently at the current speed, resulting in an unstable motion state. When rotating speed rises to  Fig. 11 Bifurcation diagram, three-dimensional spectrum, and Largest Lyapunov exponents at disk 2 in the rotor system under inphase eccentricity loading condition: a bifurcation diagram; b three-dimensional spectrum; c Largest Lyapunov exponents 4560 rad/s, the system enters into chaotic motion, and the values of the LLE are always positive numbers at 4560 rev/min \ x \ 5640 rev/min. During the speed range of [5700 * 6720] rev/min, the system enters a stable motion state after chaos motion. Vibration response at x = 5700 rev/min is presented in Fig. 13. As can be seen, the orbit is a loop; the fundamental frequency and its harmonics are the main frequency components. It can be found from the Poincaré map that the system undergoes the quasi-periodic state because the points formed a quasiperiodic loop. When rotating speed rises to 6780 rad/s, the system enters into chaotic motion, and the values of the LLE are positive at 6780 rev/min \ x \ 7260 rev/min. Then, the system changes into a long-term stable motion as the rotating speed increases. However, chaotic motion occurs in three separate rotating speed at x = 7560 rev/min, 7680 rev/min, and 12,120 rev/min. The  Fig. 11c. It can be seen from Fig. 14 that the axis orbit is complicated; the fundamental frequency, its harmonic, sub-harmonic, and complicate combination frequency components are illustrated in Fig. 14b. From the Poincaré map, it can be found that the system undergoes a chaotic motion state. It can be seen from Fig. 11 that with the rotating speed increasing, the segments of the rotating speeds of the stable and chaos motions appear by turn due to the severe vibration of the bearing. The corresponding values of the LLE also change positively or negatively as presented in Fig. 11c. Figure 11b shows the cascade spectrum of the system under in-phase eccentricity loading conditions.
Corresponding to the bifurcation diagram and Largest Lyapunov exponents (Fig. 11a, b), the whole vibration process of the system is illustrated. The results show that when the operating speed of the rotor system is lower than 5700 rev/min, the fundamental frequency f r and the harmonic 2f r in the three-dimensional spectrum are the main components. As the rotational speed increases, the harmonic 3f r appears. When the rotational speed continuously increases, the subharmonic frequency starts to emerge and the combined frequency components about 0.5f r and f r appear in the three-dimensional spectrum simultaneously.
The Bifurcation diagram, three-dimensional spectrum, and Largest Lyapunov exponents of the rotor system under out-of-phase eccentricity loading condition are shown in Fig. 15. The vibration response of the system presents more stable compared with the case of in-phase eccentricity loading condition (see Fig. 11). The system is in a stable motion state before 1680 rev/min, and the values of the LLE are always  Fig. 15 Bifurcation diagram, three-dimensional spectrum, and Largest Lyapunov exponents at disk 2 in the rotor system under out-ofphase eccentricity loading condition: a bifurcation diagram; b three-dimensional spectrum; c Largest Lyapunov exponents negative numbers as illustrated in Fig. 15c. The first chaotic motion state appears when x = 1740 rev/min, which is much earlier than that of Fig. 11. At the range of [1740 * 1800] rev/min, the motions of the system go through the variation from stable motion state to chaotic motion state, the values of the LLE are always positive numbers. As Fig. 16 shows, the axis orbit is complicated; the fundamental frequency, VC frequency, their combination frequency components, and the continuous frequency spectra are illustrated in Fig. 16b. From the Poincaré map, it can be found that the points in the Poincaré map are distributed in a disorderly state, and the system undergoes a chaotic motion state. When the rotating speed is at the range of [1860 * 15,000] rev/min, the system changes into a long-term stable motion state, the corresponding values of the LLE are always negative numbers. It can be observed from Fig. 17 that the dynamic behavior of the rotor is stable when the rotational speed is 6780 rev/min, which can be found by a closed circle in Poincaré maps. Moreover, the vibration amplitude of the rotor system is also smaller, which can be deduced by comparing Figs. 11b and 15b.
Although the chaos phenomenon of intermittency has also been observed under the out-of-phase eccentricity loading condition, the resultant force of the eccentric forces generated by the rotating of the disk and bolteddisk joint under the current operating conditions is reduced, resulting in a more stable vibration of the bearing.
As can be seen, the loading condition about the eccentric phase has an important effect on the motion stability of the system. The rotor system will shift into a chaotic state at a lower rotational speed under the condition of out-of-phase eccentricity loading condition; however, the motion state of the system at high speeds will be more stable compared with the in-phase eccentricity loading condition. The root cause lies in the fact that the resultant force of the eccentric force under the out-of-phase eccentricity loading condition is reduced, resulting in the reduction of the bearing force generated by the bearing so that the motion state is more stable.
It should be mentioned that the pre-tightening force only affects the transition point U 0 of the bending stiffness of the bolted joint [5,7]. Therefore, we define the bending stiffness at first and second bending stage (k h1 and k h2 ), as well as the transition point U 0 to simulate the bending behavior of the bolted joint.

Effect of eccentric phase difference in the bolted-disk joint
The effect of unbalanced loading conditions corresponding to the in-phase and out-of-phase eccentricities on the chaotic response of the bolted joint rotor system is investigated in Sect. 4.1. The effects of eccentric phase difference in the bolted-disk joint on the rotor system's chaotic response are discussed below. The model parameter values include the eccentric phase of disk 1 u 1 ¼ 0 , the eccentric phase of disk 3 u 3 ¼ 0 , and the eccentric phase of disk 2 was selected as the control parameter in this section. The time response corresponding to the first 160 periods of the system was discarded from the sampled data to ensure that the responses reached steady-state conditions. The steady-state responses during the last 80 periods were employed for the plots shown below. The LLEs of the rotor system were plotted to identify the rotating speed region of chaotic motion. The bifurcation diagrams, Largest Lyapunov exponents, and spectrum cascade within x [ [600, 15000] rev/ min corresponding to u 2 ¼ 90 and u 2 ¼ 180 are illustrated in Figs. 18 and 19, respectively. The LLEs can be used to find changes in chaotic behavior. As can be seen by comparing Figs. 11c and 18c and 19c, increasing the eccentric phase of disk 2 is seen to eliminate chaotic motion at high rotating speed region ([8000, 15000] rev/min), and the LLEs always show negative values as shown in Fig. 19c. For the chaotic behavior in the middle-speed range ([4000, 8000] rev/min), it can be seen that with the eccentric phase of disk 2 increasing, the segments of the rotating speeds of the unstable motions are more. Moreover, as  Fig. 18 Bifurcation diagram, three-dimensional spectrum, and Largest Lyapunov exponents of the rotor-bearing system (lumped mass point 2) at u 2 ¼ 90 : a bifurcation diagram, b three-dimensional spectrum; c Largest Lyapunov exponents the eccentric phase of disk 2 increases, the speed range of the unstable motion is gradually narrowing from a wider speed range. This will cause the stable rotation speed to become unstable in the middle-speed range, as the eccentricity of disk 2 increases. As for the low rotation speed region, the motion state of the system becomes stable as the eccentric phase of disk 2 increases, which can be confirmed by the corresponding values of LLEs are always negative when the eccentric phase of disk 2 is equal to 180°. This is  Fig. 20 The rotor-bearing test rig with multiple disks: a the test rig, b test system because the resultant force of eccentricity of the bolted-disk joint is gradually weakened as u 2 continuously increased, resulting in the resultant force of the eccentric force of the disk and the bolted joint in the rotor system is reduced so that the bearing vibration is more stable.
In the three-dimensional spectrum corresponding to u 2 ¼ 0 ; 90 ; 180 illustrated in Figs. 11b, 18b, and  19b, it can be seen that the 1/2 subharmonic frequency, 3/2 subharmonic frequency, and third harmonic frequency are disappeared with the increase in u 2 . Moreover, the second harmonic frequency component is gradually enhanced. This is because the eccentric effect of the bolted-disk joint is gradually weakened as u 2 continuously increased. From previous simulation results, it can be found that the rotor's motion state should be quite different with different eccentric phases of the disks in the bolted-disk joint.

Experimental study
In order to verify the simulation results, the experimental study is conducted through a rotor-bearing test rig with multiple disks, as shown in Fig. 20. As can be seen in Fig. 20a, two adjacent disks (Disk 1 and Disk 2) are used to simulate the bolted-disk joint, so that the eccentric phase can be controlled by screwing the bolts on the disks conveniently. The sketch of the rotor test rig is shown in Fig. 21, and the structure parameters of the rotor are shown in Table 5. The rotor test rig is supported by NSK7007 bearing, and its parameters are listed in Table 6. The rotor is driven by a frequencycontrolled motor, and the eddy current transducers are applied to measure the vibration responses in vertical, where the sampling frequency is 1613 Hz. Thirty 6 mm diameter screwed holes are drilled in disk 1, disk 2, and disk 3 along the circumferential. By screwing the bolts on the disks, the eccentric phase of the disks can be controlled conveniently.
In the experiments, the bolts were screwed at the same position in disk 1 and disk 3; the position of the   bolt screwed in disk 2 is selected as a control parameter. The experimental studies are divided into two cases, which are illustrated by Fig. 22. In experimental case 1, the bolt screwed in disk 2 is localized in the same position as disk 1 and disk 3 (see Fig. 22a). As for experimental case 2, the bolt screwed in disk 2 is localized in the opposite position as disk 1 and disk 3 (see Fig. 22b). The experiment cases are used to study the effect of in-phase and out-of-phase eccentricities of two disks in the bolted-disk joint on the rotor dynamics, where the first case (case 1) is referred to as the in-phase eccentricity loading condition and the second case (case 2) is referred to as the out-of-phase eccentricity loading condition. The three-dimensional spectrum of the measured responses corresponding to the operating speed of [30,1500] rev/min for case 1 is shown in Fig. 23a; the fundamental frequency, second harmonic frequency, third harmonic frequency, and fourth harmonic frequency can be observed as the main frequency components of the rotor system. The three-dimensional spectrum for case 2 is demonstrated by Fig. 23b; it can be found that the third harmonic frequency is significantly weakened, which is consistent with the simulation results observed in Sect. 4.2. In order to investigate the influence of the eccentric phase of the disk in the bolted-disk joint in detail, the measured vibration responses under 855 rev/min for case 1 and case 2 are shown in Figs. 24 and 25, respectively. It can be seen that the third harmonic frequency component is weakened with the increase in the eccentric phase of disk 2. Moreover, the quasi-periodic loop is broken in the Poincaré map for case 1, which means that the motion state is unstable at the current rotational speed. As for case 2, there are a large number of points lying on a closed circle in the Poincaré map; the motion is quasi-periodic. Figures 26 and 27 show the vibration responses at 1500 rev/min for case 1 and case 2. It can be seen that And the quasi-period loop is broken for case 2, which means that the motion state is unstable. From the vibration responses in experimental studies, it can be concluded that the motion state of the system under the in-phase and out-of-phase eccentricities conditions may be completely different at the same speed. All the results observed from experimental studies are consistent with the simulation results obtained in Sect. 4.2, and then the simulation results are verified.
It should be mentioned that the motivation of the theoretical simulation and experimental study in this paper is to obtain some general results, about the effect of the eccentric phase difference in a rotor-bearing system with the bolted-disk joint. That is to say, the simulation results will also be reflected in the experimental study, and the difference between the parameters of the simulation model and test rig will not affect this general phenomenon. Simulation and experimental study results show that some key harmonic frequencies will disappeared with the increase in eccentric phase difference in the rotor system, which can be regarded as a general result.

Conclusions
Developing a deep understanding of the nonlinear dynamic properties of a bolted-disk joint rotor system is of fundamental importance, due to their wide application in large-scale power-generation and aeroengines. The bolted-disk joint forms a fundamental part of such industrial rotating machines, potentially introducing a significant effect on the motion stability of the system due to the eccentric phase difference. A lumped mass model of the flexible bolted-disk joint rotor-bearing system was established, where the stiffness of the shafts is determined based on flexibility influence coefficient method in mechanics of materials, and time-varying bending stiffness of bolted joint was taken into account. In the modeling process of the bolted joint, a two-node element was developed based on Lagrange's equation; the lateral displacement stiffness, bending stiffness, and coupling stiffness between the adjacent disks were considered. This paper investigated the bifurcation behavior of a flexible rotor-bearing system with a bolted-disk joint under different eccentric force loading conditions, and the following conclusions were obtained. The chaos phenomenon of intermittency was identified by the corresponding bifurcation diagrams, the Largest Lyapunov exponent, and Poincaré maps. The simulation research shows that because the VC vibration of the system is weakened under the condition of out-of-phase eccentricity, the segments of the rotating speeds of the quasiperiodic motion under such loading conditions are much wider than those under the in-phase eccentric loading condition. Moreover, the rotor will enter into a chaotic state at a lower speed. Moreover, the vibration amplitude of the rotor system is also smaller under the condition of outof-phase eccentricity. The influence of the eccentric phase difference of the bolt-disk joint on motion stability was theoretically investigated and clarified. It is found out that increasing the eccentric phase difference of bolted joint is seen to eliminate chaotic motion at high rotating speed region, and the segments of the rotating speeds of the unstable motions are more. Furthermore, the motion state of the system becomes stable as the eccentric phase difference of the bolted joint increases at the low rotation speed. The third harmonic frequency will disappear with the increase in u 2 . Moreover, the second harmonic frequency component will gradually be enhanced. This is because the resultant force of eccentricity of the bolted-disk joint is gradually weakened as u 2 continuously increased, resulting in the resultant force of the eccentric force of the disk and the bolted joint in the rotor system is reduced. Finally, the influence of the eccentric phase difference of the bolt-disk joint on system responses was validated qualitatively by the experiment. The present study provides useful tools for the analysis and predict bolted-disk joint rotorbearing system vibrations and could be extended to consider the effect of other physical parameters of bolted joint on rotor dynamics.