Study on the E ﬀ ect Mechanism of MRD to Rotor System

: Magnetorheological damper (MRD), a semi-active control device, has the stability of the passive control device and the adaptive ability of the active control device. It is considered to be an ideal replacement for the squeeze ﬁlm damper (SFD) and has been extensively investigated. However, the e ﬀ ect mechanism of MRD to a rotor system is still unclear, which will inﬂuence the structural design and control strategy of MRD. This paper aims to propose a general MRD excitation model for a rotor system and discover the e ﬀ ect mechanism. First, the e ﬀ ects of current and eccentricity on mechanical properties of MRD are studied in detail and approaches suppressing the vibration amplitude of a rotor system are put forward. Then, a general model considering MRD, ball bearing and squirrel cage is proposed. Based on the ﬁnite element method, the dynamic equations of a rotor-bearing-squirrel cage system with MRD and mass imbalance are presented. Finally, the inﬂuences of MRD on the dynamic characteristics of the rotor system are investigated, and the e ﬀ ect mechanism is discovered. In addition, the nonlinearity which is not desirable as it causes ﬂuctuating stresses in the main components of the rotor system and eventually leads to their failure due to fatigue, are observed. The changes of damping and sti ﬀ ness of MRD can alter the damping performance of MRD and the suitable range of current forms the zone of e ﬀ ective damping, which lays helpful guidance for the design and application of the MRD.


Introduction
Rotating machinery, especially the aero-engine, is developing towards high flexibility and high reliability, which puts forward new requirements for the vibration of the rotor system. When a rotor system passes through the critical speed, a large deflection of it is produced, which brings safety hazard to the normal operation of the machine. Therefore, the effective measurements to control the vibration are of great significance to the stability and safety of the rotating machine.
The supporting system is the only way to transmit vibration. So setting a damper between the bearing and the supporting structure is an effective way to reduce the load and vibration amplitude transmitted by the bearing. Squeeze film damper (SFD) is widely used in rotating machinery because of its simple structure and obvious effect. In 1963, Cooper [1] studied the damping mechanism of SFD to the rotor system with unbalanced force and verified that SFD could effectively reduce the vibration by experiment, which was helpful for rotor operation. Wang [2] studied the effect of SFD on the rotor system with misalignment rubbing coupling faults and found that SFD with appropriate film clearance could significantly suppress the nonsynchronous response and reduce the vibration amplitude. However, if the design parameters are inappropriate, SFD will deteriorate the vibration of the rotor system. Hahn [3] found that there was a limit to the unbalance of the rotor system. When the unbalance exceeded the limit, the dynamic characteristics of the rotor system supported on SFD were worse than those supported on the rigid supporting. Also, Gunter [4] found that the supporting damping and stiffness had optimum ranges, in which SFD could improve the rotor performance. When they were not in these optimum ranges, the rotor performance supported on SFD was significantly worse than that supported on the rigid supporting. So, SFD, a passive damper, has limitations in that it cannot adapt to system uncertainty and external interference.
Active control (e.g., controllable squeeze film damper [5,6], electromagnetic bearings [7,8], shape memory alloy actuator [9], piezoelectric actuator [10][11][12]) can overcome the limitations of the passive control through adjusting its parameters to adapt to the change of the external conditions. As a kind of intelligent material, magnetorheological fluid (MRF) can rapidly change from liquid state to quasi-solid state under the action of an external magnetic field, and the change process is continuous, controllable and reversible. Magnetorheological damper (MRD) with MRF as lubricating oil is a kind of controllable SFD. It can actively change stiffness and damping under the action of an external magnetic field. Due to these features, MRD has attracted more and more attention in the field of vibration control of rotating machinery. Based on the Reynolds equation, Wang [13,14] established the model of film pressure of MRD, in which MRF was modelled as Bingham constitutive law. The mechanical characteristics of MRD and the unbalanced response characteristics of the rotor system with MRD were analyzed by numerical simulation and verified by experiment. Zhu [15,16] proposed a disc-type MRD operating in shear mode and studied the effectiveness of MRD for suppressing the vibration of the rotor system under AC sinusoidal magnetic fields by experiment. Carmignani [17][18][19][20] presented a feasibility study of MRD and compared the two solutions. After that, an experimental table of MRD-rotor was developed and a series of experimental and theoretical researches were carried out. Considering the thermal effect, Ghaednia [21] established a dynamic model of MRD based on Bingham constitutive law and analyzed the influence of temperature and current on the steady-state response of the rotor system. Zapoměl [22] presented an enhanced mathematical model based on a bilinear constitutive law and studied the effect of MRD on the force transmission and vibration attenuation of a rotor system. Although the dynamic characteristics of rotor system with MRD has been extensively investigated and the theory is mature, the squirrel cage supporting coupled with rotor system with MRD and the effect mechanism of MRD to a rotor system are rarely involved.
In this work, the dynamic model of MRD is established in which MRF is modelled as a bilinear constitutive law. The mechanical properties of MRD are studied. Then MRD is introduced to a Jeffcott rotor system with squirrel cage supporting. The dynamic behaviors of the rotor system are analyzed, and the effect mechanism of MRD to it is discovered. Assuming that the fluid is incompressible and satisfies the short bearing theory, the balance equation and continuous equation of MRF can be expressed as:

Theory and Model of MRD
where p is the film pressure. v and w are the velocity of the fluid along the Y and Z directions, respectively. The bilinear constitutive law of MRF can be expressed as:  Under the action of a magnetic field, MRF behaves as a non-newtonian fluid which contains yield zone and unyielding zone as shown in Figure 2. In the unyielding zone, MRF behaves as a solid. L is the axial length of MRD. Zc is the point in the Z-direction at which the core begins to fill the entire clearance of MRD. Assuming that the fluid is incompressible and satisfies the short bearing theory, the balance equation and continuous equation of MRF can be expressed as: where p is the film pressure. v and w are the velocity of the fluid along the Y and Z directions, respectively. The bilinear constitutive law of MRF can be expressed as: where τ, τ y and τ c = α η α η −1 τ y are shear stress, yield stress and stress at the core boundary, respectively. η is the viscosity of MRF without external magnetic field. α η = η c /η is the viscosity ratio in which η c is the viscosity of core.
Under the action of a magnetic field, MRF behaves as a non-newtonian fluid which contains yield zone and unyielding zone as shown in Figure 2. In the unyielding zone, MRF behaves as a solid. L is the axial length of MRD. Zc is the point in the Z-direction at which the core begins to fill the entire clearance of MRD. Based on the bilinear constitutive law, a general mathematical model of MRD was proposed by Zapoměl et al. [22][23][24]. Using their method, the implicit Reynolds equation of MRD can be given as: Based on the bilinear constitutive law, a general mathematical model of MRD was proposed by Zapoměl et al. [22][23][24]. Using their method, the implicit Reynolds equation of MRD can be given as: h is the first derivative of film thickness h with respect to time t. p is the first derivative of film pressure p respect to Z. By solving Equation (5), p can be obtained. Then integrating p , film pressure p can be given as: where p A is the atmospheric pressure. When 0 ≤ Z ≤ Z c , the core fills the entire clearance of MRD and MRF behaves as a Newtonian fluid with viscosity η c . Based on the classical Reynolds equation, the film pressure can be obtained: The radial and circumferential film force of MRD can be given as: where R is the radius of middle clearance of MRD. The horizontal and vertical components of the oil force can be given as: where x o and y o are the displacements of bearing outer ring in the x and y directions, respectively.
The SG-MRF2035 type of MRF is selected. The relationship between the yielding shear stress and magnetic field strength by fitting can be expressed as: where a 1 = 0.0297, a 2 = 19.75, a 3 = 1.102 × 10 4 , a 4 = −2.482 × 10 4 , c 1 = −22.29, c 2 = 2.601 × 10 4 . Assuming that the leakage of magnetic flux can be ignored, the magnetic field strength can be given as H = IN/(2h), where I is the current intensity and N is the electric coil turn. divided into ten elements. Each node contains two translational and two rotational degrees of freedom along x and y directions.

The Investigated Rotor System
Assuming that the leakage of magnetic flux can be ignored, the magnetic field strength can be given as where I is the current intensity and N is the electric coil turn. Figure 3 shows a schematic diagram of a general coupling system model of MRD-bearing-rotor with supporting of squirrel cage spring. The rotor shaft is considered to be a Timoshenko beam and divided into ten elements. Each node contains two translational and two rotational degrees of freedom along x and y directions. According to the Lagrange equation, the dynamic equation of the rotor system [25] can be given as:

The Investigated Rotor System
where According to the Lagrange equation, the dynamic equation of the rotor system [25] can be given as: where . ω is the rotating speed of the rotor shaft. The mass matrix M, gyroscopic matrix J, and stiffness matrix K of the rotor system are given in the Appendix A. E is the elastic modulus. G is the shear modulus. ν is the Poisson ratio. I is the moment of inertia. ρ is the material density. A is the cross-sectional area. l is the length of the beam unit. φ is the shear coefficient. Q 1 and Q 2 are the generalized forces which contain unbalance excitations F ux = m d eω 2 cos(ωt) and F uy = m d eω 2 sin(ωt). m d is the equivalent lump mass at disc of the rotor system. . md is the equivalent lump mass at disc of the rotor system. Figure 4 shows a schematic diagram of the rotor supporting system. The outer ring of MRD is fixed with housing. The outer ring of the bearing is connected with the squirrel cage and limited to lateral bending vibration. The inner ring of the bearing is connected with the shaft and rotates together with it. Assuming that the balls are equidistantly spaced on the cage and there is pure rolling between the ball and the raceway, the bearing forces in x and y directions [26] can be obtained based on the Hertzian contact theory:  According to Newton's second law, the dynamic equation of the rotor supporting system can be given as Assuming that the balls are equidistantly spaced on the cage and there is pure rolling between the ball and the raceway, the bearing forces in x and y directions [26] can be obtained based on the Hertzian contact theory:

The Model of the Rotor Supporting System
where θ i = 2π n (i − 1) + ω cage t and ω cage = ( r i r i +r o )ω. C b is the contact stiffness. n is the number of balls. α is the contact angle. θ i is the angle position of the ith ball. ω cage is the angular velocity of the cage. r i and r o are the radiuses of the inner raceway and outer ring raceway, respectively. µ is the bearing clearance. H is the Heaviside function.
According to Newton's second law, the dynamic equation of the rotor supporting system can be given as where m b and c e are the equivalent lump mass and the equivalent damping coefficients at the journal of the rotor shaft, respectively. k a is the supporting stiffness of the squirrel cage.

Mechanical Properties of MRD
In this section, the mathematical model of MRD is calculated and the mechanical properties of it are analyzed. It is set that the rotor journal center performs a circular orbit at a constant radius of e B , and the line of the centers moves at a constant angular velocity Ω. The calculation parameters are given as e B = 6 × 10 −4 m, Ω = 350 rad/s, L = 0.05 m, R = 0.075 m, c = 1 × 10 −3 m, N = 100, η = 0.3 Pa · s, and α η = 500. Figure 5 shows the curve of the core volume of MRD varied with the applied current in the range of 0-1.5 A. It can be seen that when I = 0-0.1 A, the core volume is sensitive to the change of the current. With the increase of the current, the yield stress increases, which makes the core volume increase approximately linearly. When I > 0.1 A, the core volume changes slowly and starts to saturate. Figure 6 shows the positions of yield surface in MRD under several currents. When I = 0.5 A, the position of the yield surface moves significantly compared with that at I = 0.03 A, and the unyielding area is drastically enlarged. When I = 1 A and 1.5 A, the position of the yield surface moves slowly compared with that at I = 0.5 A. This phenomenon indicates that when the yield stress reaches a certain value, the position of the yield surface will be fixed at a certain position and not move anymore. In this case, the flow field is saturated with yield stress.

Mechanical Properties of MRD
In this section, the mathematical model of MRD is calculated and the mechanical properties of it are analyzed. It is set that the rotor journal center performs a circular orbit at a constant radius of . Figure 5 shows the curve of the core volume of MRD varied with the applied current in the range of 0-1.5 A. It can be seen that when I = 0-0.1 A, the core volume is sensitive to the change of the current. With the increase of the current, the yield stress increases, which makes the core volume increase approximately linearly. When I > 0.1 A, the core volume changes slowly and starts to saturate. Figure 6 shows the positions of yield surface in MRD under several currents. When I = 0.5 A, the position of the yield surface moves significantly compared with that at I = 0.03 A, and the unyielding area is drastically enlarged. When I = 1 A and 1.5 A, the position of the yield surface moves slowly compared with that at I = 0.5 A. This phenomenon indicates that when the yield stress reaches a certain value, the position of the yield surface will be fixed at a certain position and not move anymore. In this case, the flow field is saturated with yield stress.   Figure 7 shows the maximum film pressure and yield stress varied with the current in the range of 0-3 A. When I = 0-0.4 A, the maximum film pressure approximately linearly increases with the increase of the current. After that, it increases slowly and reaches a peak at I = 2 A, then remains basically unchanged. The change trend of the maximum yield stress keeps consistent with that of the maximum film pressure. Compared with Figure 5, it can be found that the peak point of the core volume lags far behind that of the maximum film pressure. That is to say, when the core fulfills the clearance of MRD, the maximum film pressure still increases with the increase of the current until the yield stress of MRF reaches a peak. Figure 8 shows the pressure distribution in the middle plane of MRD under several currents. It can be seen that when I = 0 A, MRF is in a state of Newtonian fluid. The curve of film pressure is continuous and smooth. The maximum film pressure Pmax is 0.92 MPa. When I = 0.5 A, the curve of film pressure has mutations at the positions of film formation and  behind that of the maximum film pressure. That is to say, when the core fulfills the clearance of MRD, the maximum film pressure still increases with the increase of the current until the yield stress of MRF reaches a peak. Figure 8 shows the pressure distribution in the middle plane of MRD under several currents. It can be seen that when I = 0 A, MRF is in a state of Newtonian fluid. The curve of film pressure is continuous and smooth. The maximum film pressure P max is 0.92 MPa. When I = 0.5 A, the curve of film pressure has mutations at the positions of film formation and fracture. Compared with I = 0 A, the film pressure increases significantly to 4.54 MPa and the increment is about 79.7%. In addition, the position of the maximum film pressure changes from 0.862 rad to 0.63 rad. When the current increases further to 2 A, the curve of film pressure is similar to that when I = 0.5 A, but the film pressure is further increased. The maximum film pressure P max at I = 2 A is 8.27 MPa. Compared with I = 0.5 A, the maximum film pressure increases by about 45.1%. When I = 3 A, the curve is basically consistent with that when I = 2 A, and the maximum film pressure remains unchanged. It can be seen from the above phenomena that the film pressure of MRD increases with the increase of the current and reaches saturation at a certain value of the current. The main reason is that with the increase of the current, the yield stress of MRF increases and the film becomes thinner, resulting in an increase of film pressure. When the two factors reach a peak, the film pressure reaches saturation. Figure 9 shows the corresponding three-dimensional graphs and contour maps of the film pressure distribution. Through comparison, it can be found that the shape of the contour map without the current is different from that with the current, which means that the nature of MRF changes. The main reason is that under the action of the current, MRF changes from Newtonian fluid to non-Newtonian fluid and the flow properties of it change.       When 0 ε is smaller than 0.75, the equivalent stiffness and damping coefficients increase slowly.
When 0 ε is bigger than 0.75, the film clearance is very thin and the equivalent stiffness and damping coefficients have a dramatical increase which will have a significant effect on the dynamic characteristics of a rotor system. The bigger the 0 ε , the smaller is the film clearance and the larger are the increments of the equivalent stiffness and damping coefficients.  Figure 10 shows the equivalent stiffness and damping coefficients of MRD varied with eccentricity ratio ε 0 in the range of [0.4, 0.9] when I = 0.5 A. It can be found that with the increase of ε 0 , the equivalent stiffness and damping coefficients increase. The change trends are similar. When ε 0 is smaller than 0.75, the equivalent stiffness and damping coefficients increase slowly. When ε 0 is bigger than 0.75, the film clearance is very thin and the equivalent stiffness and damping coefficients have a dramatical increase which will have a significant effect on the dynamic characteristics of a rotor system. The bigger the ε 0 , the smaller is the film clearance and the larger are the increments of the equivalent stiffness and damping coefficients. of 0 ε , the equivalent stiffness and damping coefficients increase. The change trends are similar.
When 0 ε is smaller than 0.75, the equivalent stiffness and damping coefficients increase slowly.
When 0 ε is bigger than 0.75, the film clearance is very thin and the equivalent stiffness and damping coefficients have a dramatical increase which will have a significant effect on the dynamic characteristics of a rotor system. The bigger the 0 ε , the smaller is the film clearance and the larger are the increments of the equivalent stiffness and damping coefficients.  Figure 11 shows the equivalent stiffness and damping coefficients of MRD varied with the current in the range of [0, 3] A when ε0 = 0.6. It can be seen that the equivalent stiffness coefficient increases approximately linearly at the beginning, raises slowly and reaches saturation at last. The change trend of the equivalent damping coefficient is similar to that of the equivalent stiffness coefficient, but the saturation point relatively lags behind. This is different from the eccentricity ratio  Figure 11 shows the equivalent stiffness and damping coefficients of MRD varied with the current in the range of [0, 3] A when ε 0 = 0.6. It can be seen that the equivalent stiffness coefficient increases approximately linearly at the beginning, raises slowly and reaches saturation at last. The change trend of the equivalent damping coefficient is similar to that of the equivalent stiffness coefficient, but the saturation point relatively lags behind. This is different from the eccentricity ratio effects that the equivalent stiffness and damping coefficients increase indefinitely with the decrease of the eccentricity ratio.

The Effects of MRD on the Rotor System
From the above analysis, it can be obtained that the current and eccentricity ratio have significant influences on the mechanical properties of MRD (e.g., core volume, oil film pressure, the equivalent stiffness and damping coefficients) which will affect the dynamic characteristics of the rotor system. In order to further discover the effect mechanism of MRD to a rotor system, MRD is introduced to a Jeffcott rotor, and a general coupling system model of MRD-bearing-rotor with supporting of squirrel cage is established. The numerical simulation is carried out, and the dynamic responses are analyzed. The calculation parameters of the rotor system are given as  Figure 12 shows the critical speed c ω and the vibration amplitude at this speed c A P varied with the applied current in the range of 0-0.6 A. As can be seen from the figure, when I = 0-0.15 A, MRD has almost no influence on the critical speed, and the value of it is 500 rad/s. The reason may be that the supporting stiffness of the rotor system is equal to the series stiffness of the squirrel cage, Figure 11. Equivalent stiffness and damping coefficients of MRD varied with the current.

The Effects of MRD on the Rotor System
From the above analysis, it can be obtained that the current and eccentricity ratio have significant influences on the mechanical properties of MRD (e.g., core volume, oil film pressure, the equivalent stiffness and damping coefficients) which will affect the dynamic characteristics of the rotor system. In order to further discover the effect mechanism of MRD to a rotor system, MRD is introduced to a Jeffcott rotor, and a general coupling system model of MRD-bearing-rotor with supporting of squirrel cage is established. The numerical simulation is carried out, and the dynamic responses are analyzed.  Figure 12 shows the critical speed ω c and the vibration amplitude at this speed AP c varied with the applied current in the range of 0-0.6 A. As can be seen from the figure, when I = 0-0.15 A, MRD has almost no influence on the critical speed, and the value of it is 500 rad/s. The reason may be that the supporting stiffness of the rotor system is equal to the series stiffness of the squirrel cage, bearing and MRD (as can be seen from Figure 4), and the stiffness of MRD in this current range is far less than the sum of the stiffness of the squirrel and bearing. When I ≥ 0.15 A, the stiffness of MRD reaches a certain level due to the increase of the yield stress of MRF, which has an influence on the critical speed of the rotor system. The higher the current, the larger the critical speed. When I = 0.5 A, MRD approaches quasi-rigid supporting, and the critical speed of the rotor system reaches the maximum value which is 690 rad/s. Compared with the critical speed at I = 0 A, the critical speed at I = 0.5 A increases by about 27.5%. When I ≥ 0.5 A, the critical speed remains relatively constant.
Appl. Sci. 2019, 9, x 12 of 16 and the squeeze effect is the source of MRD damping. Comparing the dynamic responses in Figure  13, it can be found that c A P at I = 0.5 A (0.31 mm) is larger than that at I = 0.3 A (0.25 mm). The increment is about 19.3%. The orbit size is larger and presents an unstable ellipse. In addition to the fundamental frequency, there appear irrational frequencies in the spectrogram and a closed loop in the Poincaré map, which indicates that the rotor system takes nonlinearity and changes from 1Tperiodic motion to quasi-periodic motion. The rotor system is unsteady. The reason may be that the large core volume and large deflection at critical speed lead to a too thin film so as to take nonlinearity, resulting in instability of the rotor system.  It can also be seen from Figure 12 that MRD has a complex influence on the vibration amplitude of the rotor system at critical speed. When I = 0-0.3 A, MRD mainly affects the vibration of the rotor system through its damping effect, which has a better suppression on AP c . The higher the current, the more obvious the damping effect. As shown in Figure 13, when I = 0 A and 0.3 A, there is an isolated point in the Poincaré map, which indicates that the rotor system presents a stable 1T-periodic motion. But when I = 0.3 A, AP c is 0.25 mm which is far less than 0.516 mm when I = 0 A. The decrement is about 51.6% and the size of the orbit is much smaller. After I = 0.3 A, with the increase of the current, the stiffness of MRD plays a major role, which weakens the damping effect of MRD and AP c gradually increases. These phenomena are consistent with the experiment results of Figure 13 in Ref. [27]. The reason may be that the increase of MRD stiffness limits the oil film to be squeezed and the squeeze effect is the source of MRD damping. Comparing the dynamic responses in Figure 13, it can be found that AP c at I = 0.5 A (0.31 mm) is larger than that at I = 0.3 A (0.25 mm). The increment is about 19.3%. The orbit size is larger and presents an unstable ellipse. In addition to the fundamental frequency, there appear irrational frequencies in the spectrogram and a closed loop in the Poincaré map, which indicates that the rotor system takes nonlinearity and changes from 1T-periodic motion to quasi-periodic motion. The rotor system is unsteady. The reason may be that the large core volume and large deflection at critical speed lead to a too thin film so as to take nonlinearity, resulting in instability of the rotor system.

Conclusions
To study the effect mechanism of MRD to a rotor system, a model of MRD based on bilinear constitutive law is built and introduced to a Jeffcott rotor. A general coupling system model of MRD-bearing-rotor with supporting of squirrel cage is developed using the finite element method. The mechanical properties of MRD are analyzed in detail, and the effect mechanism is discovered. Some conclusions are summarized as follows: 1.
The change of the yield surface position that will cause the change of flowing space is an important mechanical response of MRD to the applied current. With the increase of the yield stress, the position of the yield surface moves and the region of unyielding is gradually enlarged. When the yield stress hits a certain value, the flow field is saturated with it, and the yield surface is fixed at a certain position; 2.
With or without applying current, the mechanical properties of MRD are different. When the current is applied, the azimuth of the maximum film pressure moves and the film pressure changes abruptly at the positions of film formation and fracture. With the increase of the current, the film pressure increases obviously and the azimuth of the maximum film pressure remains constant; 3.
Eccentricity ratio and current have significant influences on the equivalent stiffness and damping coefficients of MRD. With the increase of the eccentricity or current, the equivalent stiffness and damping coefficients increase. The bigger the eccentricity ratio, the larger the increments of the equivalent stiffness and damping coefficients; 4.
MRD mainly affects the dynamic characteristics of the rotor system through its damping and stiffness effects, and the two effects are coupled. With the increase of the current, MRD first affects the dynamic characteristics of the rotor system mainly through its damping effect, which has an obvious suppression on the vibration of the rotor system at the critical speed. When the current increases to a certain value, the stiffness effect of MRD plays a major role, which limits the oil film to be squeezed and the damping effect is weakened. In addition, in the case of the large core volume and large deflection, the oil film is too thin and takes nonlinearities, resulting in instability of the rotor system. eccentricity Ω angular velocity of the line of the journal centers I external current applied on the clearance of MRD X, Y, Z coordinates fixed to the outer ring along circumferential, radial and axial directions x, y, z fixed coordinate system h oil film thickness p film pressure v, w velocity of the fluid along the Y and Z directions τ, τ y , τ c , shear stress, yield stress and stress at the core boundary of MRF ε eccentric ratio h 0 film thickness of the untextured bearing η viscosity of MRF without external magnetic field α η viscosity ratio η c viscosity of core L axial length of MRD Z c point in the Z-direction at which the core begins to fill the entire clearance of MRD . h first derivative of film thickness with respect to time p first derivative of film pressure p respect to Z p A atmospheric pressure R radius of middle clearance of MRD x o , y o displacements of bearing outer ring in the x and y directions, respectively H magnetic field strength N electric coil turn ω rotating speed of the rotor shaft M, J, K mass matrix, gyroscopic matrix and stiffness matrix of the rotor system E, G, ν elastic modulus, shear modulus and Poisson ratio ρ, A, l material density, cross-sectional area and length of the beam unit φ shear coefficient Q 1 , Q 2 generalized forces C b contact stiffness α contact angle θ i angle position of the ith ball ω cage angular velocity of the cage r i , r o raceway radiuses of the inner and outer ring, respectively µ bearing clearance m b , c e equivalent lump mass and the damping coefficients at the journal of the rotor shaft k a supporting stiffness of the squirrel cage spring  M is the mass matrix of which the assembly form is shown as Figure A1. The assembly form of J and K are the same as that of M.