# Nonlinear Dynamic Analysis of Rotor-Bearing System with Cubic Nonlinearity

Authors: Guofang Nan; Yujie Zhu; Yang Zhang; Wei Guo
Journal: Shock and Vibration
 (2021)

---

## Abstract

Nonlinear dynamic characteristics of a rotor-bearing system with cubic nonlinearity are investigated. The comprehensive effects of the unbalanced excitation, the internal clearance, the nonlinear Hertzian contact force, the varying compliance vibration, and the nonlinear stiffness of support material are considered. The expression with the linear and the cubic nonlinear terms is adopted to characterize the synthetical nonlinearity of the rotor-bearing system. The effects of nonlinear stiffness, rotating speed, and mass eccentricity on the dynamic behaviors of the system are studied using the rotor trajectory diagrams, bifurcation diagrams, and Poincaré map. The complicated dynamic behaviors and types of routes to chaos are found, including the periodic doubling bifurcation, sudden transition, and quasiperiodic from periodic motion to chaos. The research results show that the system has complex nonlinear dynamic behaviors such as multiple period, paroxysmal bifurcation, inverse bifurcation, jumping phenomena, and chaos; the nonlinear characteristics of the system are significantly enhanced with the increase of the nonlinear stiffness, and the material with lower nonlinear stiffness is more conducive to the stable operation of the system. The research will contribute to a comprehensive understanding of the nonlinear dynamics of the rotor-bearing system.


---

## Body

## 1. Introduction

In rotational machinery, the rotor is commonly supported by the rolling element bearing. The inner race of the bearing is affixed to a shaft and housing. Accordingly, there exists a coupling effect between the bearings and the rotor. It is essential to consider the rotor and bearings as an integrated part. Dynamic analysis of the rotor-bearing system is of great significance to exactly diagnose the fault of a rotor system. Therefore, much research on this topic has been carried out to analyze the dynamic characteristics of the rotor-bearing system.Jeffcott rotor has been widely applied to study the fundamental problems of rotor dynamics [1]. Afterwards, researchers carried out a more in-depth study on this basis, including the Alford force [2], the squeeze film dampers [3, 4], and the radial clearance between the elements and the races. A widely used model is that the rotor-bearing assembly is regarded as a spring-mass system where the bearings are modeled to be the massless springs and the rotor is modeled as a mass. For this model, investigators have studied the nonlinear dynamics of the rotor-bearing system caused by rotating speed effects [5–7], ball size variation [8–10], surface waviness [11–15], the number of balls [16], imbalance [17, 18], and internal radial clearance [19–21]. In addition to the commonly used two or three degrees of freedom rotor-bearing models, more complex models were presented to research on the dynamic behaviors of the rotor-bearing system. Aini et al. [22, 23] proposed a model with five DOFs, including the axial direction, two radial directions, and the yawing and the pitching modes about the two radial directions, respectively. Maraini and Nataraj [24] presented a technique for modeling the nonlinear behavior of a rotor-bearing system with the Hertzian contact, the internal clearance, and the imbalance. The nonlinear bearing force is replaced with an equivalent describing function gain. The research results show the presence of the resonance in bearings with both the clearance and the preload. And the numerical integration was conducted which is in agreement with the theoretical solution. Cesar et al. [25] studied the influence of unbalance levels on the nonlinear dynamics of a rotor-backup rolling bearing system. Cao et al. [26] researched on the vibration and the stability analysis of a rotor-bearing pedestal system due to the clearance fit, and experiments were carried out on a rotor test rig and good agreement verified the accuracy of the proposed dynamic model.Furthermore, other scientists have carried out research on the dynamics of the cubic nonlinear rotors. In the second section of mathematical modeling in [27], the motion formula of bearing center is proposed, the asynchronous vibration of a flexible rotor supported by two bearings with nonlinear spring is studied, and the chaotic characteristics of the flexible rotor system supported by nonlinear suspension bearing are determined. In [28], the dynamic equation of the rotor system is established by using the nonlinear stiffness model provided and ignoring the influence of gravity. A rotor dynamic model with both the linear and the cubic nonlinear stiffness terms is established. By numerical analysis, it is found that the chaotic motion of the system corresponds to large parameter space. In [29], the physical model and control equation are formulated, and the modeling of spring support is equivalent to isotropic spring with nonlinearity.It is concluded that the cubic nonlinear stiffness of the elastic support will cause the vibration and the hysteresis of the resonance curve when the rotor is running in the critical speed region, and the excessive amplitude of the rotor will also aggravate its unstable behavior. In summary, the research of this subject is of great significance and is a concrete problem in practical engineering, which is worthy of further and systematic research.In the rotor-bearing system, there are many nonlinear factors such as bearing clearance, the Hertzian contact force, imbalance, and support materials. The nonlinear factors from the rolling bearing and the rotor may affect the operation of the system together. At present, there are rarely articles to study the synthetical effect of these two types of nonlinear factors. The research objective of this paper is to analyze the dynamic response of the rotor-bearing system under consideration of the cubic nonlinear factors.


## 2. Modeling of the Rotor-Bearing System

### 2.1. Rolling Element Bearing

The nonlinear dynamic in this research focuses on the global dynamics of the rotor-bearing system instead of the local characteristics. Therefore, a simplified approximate model is developed to research the complicated system of the rotor bearing. It is important to establish the analytical model of the rotor-bearing system. The first step is to build the model of the rolling element bearing. Figure1 shows the schematic diagram of the rolling bearing model. The rolling bearing includes the inner race, the outer race, the middle rolling elements, and the cage. The ball in the track has both the revolution around the center of the bearing (O1) and the rotation around the center of the ball. The initial position of the mass center for the rotor at the bearing is O1. During the operation of the rotor, the position of the mass center for the rotor is O2. As shown in Figure 1, the inner race radius is R1 and the outer race radius is R2. The outer race is affixed to the bearing pedestal and its vibration is ignored, and the inner race is connected to the shaft; then, the speed of the outer race for the bearing is v2=0; if the speed of the inner race for the bearing is equal to the speed of the rotor, then the speed of the inner race is v1=ω1R1. A hypothesis for the study is proposed that the rolling elements are evenly arranged in the cage and there is no sliding friction between the ball and the inner and the outer raceways, and the linear velocity of the ball rotation is given byFigure 1
The schematic diagram of the rolling element bearing.(1)v3=v1+v22=12ω1R1.The angular velocity of the ball revolution is given by(2)ω3=v3R3=ω1R1/2R1+R2/2=ω1R1R1+R2.Therefore, the position angle of thenth (n = 1, 2, …, n) element is as follows:(3)θn=ω3t+2πN.Therefore, the radial contact deformation between thenth element and the inner race is given by(4)Δn=xjcosθn+yjsinθn−c0,where xj is the horizontal displacement of the rotor at the bearing at the left or the right end of the shaft, yj is the vertical displacement of the rotor at the bearing at the left or the right end of the shaft, and c0 is the clearance between the bearing track and the element. Since the contact force exists only when the contact deformation is positive, the contact force between each element and the inner race can be expressed as follows:(5)Fn=KSΔn3/2,Δn>0,0,Δn≤0,where Fn is the bearing force of the nth element on the rotor and KS is the Hertzian contact stiffness of the bearing. The components of the nth element bearing force on the rotor in the horizontal direction x and the vertical direction y are, respectively, as follows:(6)FnX=Fncosθn,FnY=Fnsinθn.The resultant force of the contact force between each element and the inner race is the force of the rolling bearing on the shaft, consequently,(7)FXj=∑n=1NFnX,FYj=∑n=1NFnY,where FnX is the component of the bearing force for the nth element on the rotor in the horizontal direction and FnY is the component of the bearing force for the nth element on the rotor in the vertical direction.Another important model for this study is the rotor system. Based on the rolling element bearing, the rotor model will be built in this section.Among the nonlinear sources of the rotor-bearing system, one is caused by materials of the bearing seat, such as the rubber, the copper, and the aluminum alloys and the other is the radial clearance of rolling bearing, the VC vibration, the ball surface ripple, and the defect, etc. Some scholars have carried out the related experiments and fitted the elastic force and the bearing deformation. There is a mathematical relationship of the cubic nonlinearity, and the fitting results are compared with the traditional analytical model. It is proved that the cubic nonlinear fitting method is more suitable [30]. The team has done a lot of research studies on rotor cubic nonlinearity [31]. The mathematical relationship between the elastic restoring force and the deformation in the rotor-bearing system is adopted in this paper. According to the research results of Cvetanin, the stiffness of the shaft increases with the increase of displacement, and its relationship between the stiffness and the displacement can be expressed by the linear equation with the cubic term. Figure 2 shows the schematic diagram of the rotor-bearing system, and the elastic restoring force of the system is given by(8)Fkx=kx+k1x2+y2x,Fky=ky+k1x2+y2y.Figure 2
Schematic diagram of the rotor-bearing system.In this model, the rotor is simplified into three lumped masses. The positions of pointsP, R, and L are, respectively, located at the mass center of the disk, the mass center of the rotor at the right bearing, and the mass center of the rotor at the left bearing. mP, mR, and mLrepresent the equivalent lumped mass of the rotor at points P, R, and L, respectively; cp is the damping coefficients of the rotor at point P; cb is the damping coefficients of the rotor at both ends of the bearing. The horizontal and the vertical relative displacements at the points P, R, and L are xP, yP, xR, yR, xL, and yL, respectively. The horizontal and the vertical relative displacement differences between the point P and the point R are, respectively, as follows:(9)xPR=xP−xR,yPR=yP−yR.The horizontal and the vertical relative displacement differences between the pointP and the point L are, respectively, as follows:(10)xPL=xP−xL,yPL=yP−yL.According to Newton’s second theorem, the dynamic equation of the system can be obtained as follows:(11)mPx¨P+cPx˙P+kxPR+k1xPR2+yPR2xPR+kxPL+k1xPL2+yPL2xPL=mPeω12cosω1t,mPy¨P+cPy˙P+kyPR+k1xPR2+yPR2yPR+kyPL+k1xPL2+yPL2yPL=mPeω12sinω1t−mPg,mRx¨R+cbx˙R+k−xPR+k1xPR2+yPR2−xPR=FXR,mRy¨R+cby˙R+k−yPR+k1xPR2+yPR2−yPR=FYR−mRg,mLx¨L+cbx˙L+k−xPL+k1xPL2+yPL2−xPL=FXL,mLy¨L+cby˙L+k−yPL+k1xPL2+yPL2−yPL=FYL−mLg,where e is the eccentricity of the disc; ω1 is the angular velocity of the rotor; FXR and FYR are the horizontal and the vertical components of the reaction force for the rolling bearing on the rotor at the right end of the rotor, respectively; and FXL and FYL are the horizontal and the vertical components of the reaction force of the rolling bearing on the rotor at the left end of the rotor, respectively.For the convenience of the research, equation (11) is dealt with dimensionless form by introducing the following terms:(12)τ=ω1t,X1=x1δ2,Y1=y1δ2,X2=x2δ2,Y2=y2δ2,X3=x3δ2,Y3=y3δ2,X˙1=x˙1ω1δ2,Y˙1=y˙1ω1δ2,X˙2=x˙2ω1δ2,Y˙2=y˙2ω1δ2,X˙3=x˙3ω1δ2,Y˙3=y˙3ω1δ2,X¨1=x¨1ω12δ2,Y¨1=y¨1ω12δ2,X¨2=x¨2ω12δ2,Y¨2=y¨2ω12δ2,X¨3=x¨3ω12δ2,Y¨3=y¨3ω12δ2,ρ=eδ2,ξ1=c1m1ω1,ξ2=c2m2ω1,ξ3=c1m3ω1,η1=km1ω12,η2=k1δ22m1ω12,η3=km2ω12,η4=k1δ22m2ω12,η5=km3ω12,η6=k1δ22m3ω12,G=gδ2ω12,M2=m2δ2ω12,M3=m3δ2ω12,where 1 corresponds to P of equation (11); 2 corresponds to R of equation (11); 3 corresponds to L of equation (11); τ is the nondimensional time; δ2 is a dimension parameter; and X and Y are the nondimensional displacements.The dimensionless governing equations are as follows:(13)X¨1=−ξ1X˙1−η0x1−x2−η1x1−x2x1−x22+y1−y22−η0x1−x3−η1x1−x3x1−x32+y1−y32+ρcosτ,Y¨1=−ξ1Y˙1−η1Y1−Y2−η2Y1−Y2Y1−Y22+Y1−Y22−η1Y1−Y3−η2Y1−Y3X1−X32+Y1−Y32+ρsinτ−G,X¨2=−ξ2X˙2+η3X1−X2+η4X1−X2X1−X22+Y1−Y22+Fx2M2,Y¨2=−ξ2Y˙2+η3Y1−Y2+η4Y1−Y2X1−X22+Y1−Y22+Fy2M2,X¨3=−ξ3X˙3+η5X1−X3+η6X1−X3X1−X32+Y1−Y32+Fx3M3,Y¨3=−ξ3Y˙3+η5Y1−Y3+η6Y1−Y3X1−X32+Y1−Y32+Fy3M3..


## 2.1. Rolling Element Bearing

The nonlinear dynamic in this research focuses on the global dynamics of the rotor-bearing system instead of the local characteristics. Therefore, a simplified approximate model is developed to research the complicated system of the rotor bearing. It is important to establish the analytical model of the rotor-bearing system. The first step is to build the model of the rolling element bearing. Figure1 shows the schematic diagram of the rolling bearing model. The rolling bearing includes the inner race, the outer race, the middle rolling elements, and the cage. The ball in the track has both the revolution around the center of the bearing (O1) and the rotation around the center of the ball. The initial position of the mass center for the rotor at the bearing is O1. During the operation of the rotor, the position of the mass center for the rotor is O2. As shown in Figure 1, the inner race radius is R1 and the outer race radius is R2. The outer race is affixed to the bearing pedestal and its vibration is ignored, and the inner race is connected to the shaft; then, the speed of the outer race for the bearing is v2=0; if the speed of the inner race for the bearing is equal to the speed of the rotor, then the speed of the inner race is v1=ω1R1. A hypothesis for the study is proposed that the rolling elements are evenly arranged in the cage and there is no sliding friction between the ball and the inner and the outer raceways, and the linear velocity of the ball rotation is given byFigure 1
The schematic diagram of the rolling element bearing.(1)v3=v1+v22=12ω1R1.The angular velocity of the ball revolution is given by(2)ω3=v3R3=ω1R1/2R1+R2/2=ω1R1R1+R2.Therefore, the position angle of thenth (n = 1, 2, …, n) element is as follows:(3)θn=ω3t+2πN.Therefore, the radial contact deformation between thenth element and the inner race is given by(4)Δn=xjcosθn+yjsinθn−c0,where xj is the horizontal displacement of the rotor at the bearing at the left or the right end of the shaft, yj is the vertical displacement of the rotor at the bearing at the left or the right end of the shaft, and c0 is the clearance between the bearing track and the element. Since the contact force exists only when the contact deformation is positive, the contact force between each element and the inner race can be expressed as follows:(5)Fn=KSΔn3/2,Δn>0,0,Δn≤0,where Fn is the bearing force of the nth element on the rotor and KS is the Hertzian contact stiffness of the bearing. The components of the nth element bearing force on the rotor in the horizontal direction x and the vertical direction y are, respectively, as follows:(6)FnX=Fncosθn,FnY=Fnsinθn.The resultant force of the contact force between each element and the inner race is the force of the rolling bearing on the shaft, consequently,(7)FXj=∑n=1NFnX,FYj=∑n=1NFnY,where FnX is the component of the bearing force for the nth element on the rotor in the horizontal direction and FnY is the component of the bearing force for the nth element on the rotor in the vertical direction.Another important model for this study is the rotor system. Based on the rolling element bearing, the rotor model will be built in this section.Among the nonlinear sources of the rotor-bearing system, one is caused by materials of the bearing seat, such as the rubber, the copper, and the aluminum alloys and the other is the radial clearance of rolling bearing, the VC vibration, the ball surface ripple, and the defect, etc. Some scholars have carried out the related experiments and fitted the elastic force and the bearing deformation. There is a mathematical relationship of the cubic nonlinearity, and the fitting results are compared with the traditional analytical model. It is proved that the cubic nonlinear fitting method is more suitable [30]. The team has done a lot of research studies on rotor cubic nonlinearity [31]. The mathematical relationship between the elastic restoring force and the deformation in the rotor-bearing system is adopted in this paper. According to the research results of Cvetanin, the stiffness of the shaft increases with the increase of displacement, and its relationship between the stiffness and the displacement can be expressed by the linear equation with the cubic term. Figure 2 shows the schematic diagram of the rotor-bearing system, and the elastic restoring force of the system is given by(8)Fkx=kx+k1x2+y2x,Fky=ky+k1x2+y2y.Figure 2
Schematic diagram of the rotor-bearing system.In this model, the rotor is simplified into three lumped masses. The positions of pointsP, R, and L are, respectively, located at the mass center of the disk, the mass center of the rotor at the right bearing, and the mass center of the rotor at the left bearing. mP, mR, and mLrepresent the equivalent lumped mass of the rotor at points P, R, and L, respectively; cp is the damping coefficients of the rotor at point P; cb is the damping coefficients of the rotor at both ends of the bearing. The horizontal and the vertical relative displacements at the points P, R, and L are xP, yP, xR, yR, xL, and yL, respectively. The horizontal and the vertical relative displacement differences between the point P and the point R are, respectively, as follows:(9)xPR=xP−xR,yPR=yP−yR.The horizontal and the vertical relative displacement differences between the pointP and the point L are, respectively, as follows:(10)xPL=xP−xL,yPL=yP−yL.According to Newton’s second theorem, the dynamic equation of the system can be obtained as follows:(11)mPx¨P+cPx˙P+kxPR+k1xPR2+yPR2xPR+kxPL+k1xPL2+yPL2xPL=mPeω12cosω1t,mPy¨P+cPy˙P+kyPR+k1xPR2+yPR2yPR+kyPL+k1xPL2+yPL2yPL=mPeω12sinω1t−mPg,mRx¨R+cbx˙R+k−xPR+k1xPR2+yPR2−xPR=FXR,mRy¨R+cby˙R+k−yPR+k1xPR2+yPR2−yPR=FYR−mRg,mLx¨L+cbx˙L+k−xPL+k1xPL2+yPL2−xPL=FXL,mLy¨L+cby˙L+k−yPL+k1xPL2+yPL2−yPL=FYL−mLg,where e is the eccentricity of the disc; ω1 is the angular velocity of the rotor; FXR and FYR are the horizontal and the vertical components of the reaction force for the rolling bearing on the rotor at the right end of the rotor, respectively; and FXL and FYL are the horizontal and the vertical components of the reaction force of the rolling bearing on the rotor at the left end of the rotor, respectively.For the convenience of the research, equation (11) is dealt with dimensionless form by introducing the following terms:(12)τ=ω1t,X1=x1δ2,Y1=y1δ2,X2=x2δ2,Y2=y2δ2,X3=x3δ2,Y3=y3δ2,X˙1=x˙1ω1δ2,Y˙1=y˙1ω1δ2,X˙2=x˙2ω1δ2,Y˙2=y˙2ω1δ2,X˙3=x˙3ω1δ2,Y˙3=y˙3ω1δ2,X¨1=x¨1ω12δ2,Y¨1=y¨1ω12δ2,X¨2=x¨2ω12δ2,Y¨2=y¨2ω12δ2,X¨3=x¨3ω12δ2,Y¨3=y¨3ω12δ2,ρ=eδ2,ξ1=c1m1ω1,ξ2=c2m2ω1,ξ3=c1m3ω1,η1=km1ω12,η2=k1δ22m1ω12,η3=km2ω12,η4=k1δ22m2ω12,η5=km3ω12,η6=k1δ22m3ω12,G=gδ2ω12,M2=m2δ2ω12,M3=m3δ2ω12,where 1 corresponds to P of equation (11); 2 corresponds to R of equation (11); 3 corresponds to L of equation (11); τ is the nondimensional time; δ2 is a dimension parameter; and X and Y are the nondimensional displacements.The dimensionless governing equations are as follows:(13)X¨1=−ξ1X˙1−η0x1−x2−η1x1−x2x1−x22+y1−y22−η0x1−x3−η1x1−x3x1−x32+y1−y32+ρcosτ,Y¨1=−ξ1Y˙1−η1Y1−Y2−η2Y1−Y2Y1−Y22+Y1−Y22−η1Y1−Y3−η2Y1−Y3X1−X32+Y1−Y32+ρsinτ−G,X¨2=−ξ2X˙2+η3X1−X2+η4X1−X2X1−X22+Y1−Y22+Fx2M2,Y¨2=−ξ2Y˙2+η3Y1−Y2+η4Y1−Y2X1−X22+Y1−Y22+Fy2M2,X¨3=−ξ3X˙3+η5X1−X3+η6X1−X3X1−X32+Y1−Y32+Fx3M3,Y¨3=−ξ3Y˙3+η5Y1−Y3+η6Y1−Y3X1−X32+Y1−Y32+Fy3M3..


## 3. Results and Analysis

The dynamics of the rotor-bearing system with the cubic nonlinearity are digitally simulated in Matlab. In order to calculate the nonlinear response of the system, the fourth-order Runge–Kutta method is applied to solve the steady-state response of the governing equations. The relevant calculation parameters of the rotor system are shown in Table1. The dimension parameter δ2 is 4 × 10−5 mm. The system is analyzed by means of the bifurcation diagram, spectrum diagram, axis track diagram, Poincaré map, and the waterfall chart.Table 1
Calculation parameters of the rotor system.
ParameterValueEquivalent mass at disk,mP (kg)32.1Equivalent mass at the left bearing,mL (kg)4.0Equivalent mass at the right bearing,mR (kg)4.0Damping coefficient at the disk,cP (Ns/m)2100Damping coefficient at the bearingcb (Ns/m)1050Linear stiffness coefficient of shaft,k (N/m)2.5 × 107


### 3.1. Influence of Cubic Nonlinear Stiffness

A novelty of this study lies in consideration of the cubic nonlinearity of the shaft. Firstly, the influence of the cubic nonlinear stiffness on the dynamic behaviors will be focused on in this section. The model of the rolling element bearing selected in this paper is JIS6306. The relevant parameters are as follows: the inner radiusR1 = 40.1 mm, the outer radius R2 = 63.9 mm, the number of balls N = 8, and the Hertzian contact stiffness KS=13.34×109N/m3/2. To study the influence of the nonlinear stiffness on the response of the system, the bifurcation diagrams of the system with the different nonlinear stiffness are drawn. The longitudinal coordinate of the bifurcation diagram is the vibration response including the negative response and the positive response for getting the proportional relationship of the two directions (the negative direction and the positive direction). Figures 3 and 4 are the bifurcation diagrams of the nondimensional displacement in x-direction X1 versus the rotational speed Ω1 rad/s when the bearing clearance c0 is 10 μm and 20 μm, respectively. The nonlinear stiffness coefficient is 0 N/m3, 3 × 1015 N/m3, 3 × 1016 N/m3. and 3 × 1017 N/m3, respectively.Figure 3
Bifurcation diagrams of the nondimensional displacement inx-direction X1 versus the rotational speed Ω1 rad/s for the clearance c0 = 10 μm. The jumping phenomenon occurs in the orange circles; chaotic motion occurs in the green circles. (a) k1 = 0 N/m3; (b) k1 = 3 × 1015 N/m3; (c) k1 = 3 × 1016 N/m3; (d) k1 = 3 × 1017 N/m3.
(a)(b)(c)(d)Figure 4
Bifurcation diagrams of the nondimensional displacement inx-direction X1 versus the rotational speed Ω1 rad/s for the clearance c0 = 20 μm. Note that P-1 means the periodic one motion, P-2 means the periodic two motion, and A means the aperiodic motion. (a) k1 = 0 N/m3; (b) k1 = 3 × 1015 N/m3; (c) k1 = 3 × 1016 N/m3; (d) k1 = 3 × 1017 N/m3.
(a)(b)(c)(d)It can be seen from Figures3(a)–3(d) that the system shows complex dynamic behaviors such as the periodic motion, the jump phenomena, quasiperiodic motion, and chaos, and the amplitude of the chaotic region increases and the range of the chaotic region changes with the enhancement of the nonlinear stiffness; in addition, the critical speed of the rotor increases and amplitude of the responses decreases with the increase of the stiffness. From Figure 3(a), it can be shown that there exists an obvious jumping phenomenon at the rotating speed of 760 rad/s. The response begins to change to chaotic from periodic at 1130 rad/s, and the transition path from the periodic one motion (P-1) to chaos is the periodic two motion (P-2) and then the system enters into the periodic one motion (P-1) once again, then enters into the periodic two motion (P-2) once more, and then into chaos, as illustrated in the local enlarged diagram in Figure 3(a). When the cubic nonlinear stiffness coefficient increases to 3 × 1015 N/m3, the speed scope of the jumping phenomenon expands to 775–885 rad/s, as shown in Figure 3(b); with the increase of the stiffness, the speed scope of jumping further expands to 780–1190 rad/s, as shown in Figure 3(c); When the stiffness increases to k1 = 3 × 1017 N/m3, only one jumping phenomenon exists at the beginning at the rotating speed of 792 rad/s and the chaos is more chaotic than that of smaller stiffness, as shown in Figure 3(d). It is crucial for rotational machinery to calculate the transition from the periodic motion to the aperiodic motion (A). The jumping phenomenon should be avoided to make the system run safely and stably in engineering practice.The clearance between the element and races has a certain effect on the nonlinear behaviors of the responses. In order to study the influence of the clearance, the bifurcation diagram for the clearancec0 = 20 μm is drawn, as illustrated in Figure 4. It can be seen from a comparison of Figures 3 and 4 that the clearance has an obvious influence on the dynamic characteristics, and the amplitude and range of the chaos are larger with the increase of the clearance. The periodic one motion (P-1) exists only in a very small speed range, and the aperiodic motion (A) occupies a larger percentage of the rotating speed range comparing with that of Figure 3; it is inferred that the larger the nonlinear stiffness coefficient is, the more obvious the nonlinear characteristics are. Therefore, it is of great significance to reduce the bearing clearance and design the suitable stiffness for the stable operation of the system.


### 3.2. Influence of Rotating Speed

The rotating speed is one of the key parameters affecting the dynamic behavior of the rotor system. In this section, the bifurcation diagrams, the waterfall plots, the time domain diagrams, the phase diagrams, and the Poincaré maps of the rotor system with nonlinear stiffness are drawn to exactly analyze the dynamics.Figure5(a) shows the bifurcation diagram of the nondimensional displacement in x-direction X1 versus the rotational speed Ω1 rad/s for e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3; Figure 5(b) shows the waterfall plot of the vibrational amplitude A versus the frequency ratio (f/fn) with the changing rotating speed Ω1 under the same parameters e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3. Figure 6 shows the time domain diagram, the phase trajectory, and the Poincaré section of the system at each speed where the dynamic characteristics of the system can be obtained more intuitively. It can be seen that the system has the periodic one motion, the periodic two motion, multiple periodic motion, and chaotic motion in the rotating speed span of 0–2000 rad/s, as demonstrated in Figure 5(a). To study the dynamic characteristic of the typical vibration, the typical speeds of the system are selected as 516 rad/s, 1490 rad/s, 1570 rad/s, and 1675 rad/s, and the corresponding diagrams are shown in Figures 6(a)–6(d). It can be observed from Figure 5(b) that the spectrum shows only 1X frequency instead of the 0.5X frequency or combined frequency in the lower Ω1 (0–1380 rad/s) because the varying compliance vibration is weaker when the rotating speed is lower, and the vibration induced by the eccentric force dominates the vibration; when the rotating speed increases, the spectrum diagram shows the combined frequencies including the 0.5X and 1X frequency, as shown in Figure 5(b) and spectrum diagrams of Figure 6.Figure 5
(a) Bifurcation diagram of the nondimensional displacement inx-direction X1 versus the rotational speed Ω1 rad/s for e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3. (b) Waterfall plot of displacement responses under the same parameters e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3.
(a)(b)Figure 6
Dynamic behaviors at the typical rotating speed. (a) P-1 motion: (a1) time history, (a2) spectrum, chart of axis track, (a3) phase trajectory, and (a4) Poincaré section whenω = 516 rad/s. (b) Multiple periodic motion: (b1) time history, (b2) spectrum, chart of axis track, (b3) phase trajectory, and (b4) Poincaré section when ω = 1490 rad/s. (c) P-2 motion: (c1) time history, (c2) spectrum, chart of axis track, (c3) phase trajectory, and (c4) Poincaré section when ω = 1570 rad/s. (d) Chaotic motion: (d1) time history, (d2) spectrum, chart of axis track, (d3) phase trajectory, and (d4) Poincaré section when ω = 1675 rad/s.
(a)(b)(c)(d)


### 3.3. Influence of Eccentricity

Because of the uneven distribution of rotor mass, the rotor will produce eccentric force during the operation of the system. The existence of eccentric force leads to a big risk for the safety of the system. The exciting forces induced by the mass imbalance are the major source of rotor vibrations in the field, so it is necessary to analyze the effect of eccentricity on the motion characteristics of the rotor-bearing system. In order to study the influence of mass eccentricity on the rotor-bearing system, the nonlinear stiffness system is analyzed in this section. The imbalance is represented by the eccentricity of the rotore. Due to the existence of nonlinear stiffness of the shaft, the displacement and elastic recovery force at the three focused positions are nonlinear. The influence of the mass eccentricity on the vibration characteristics of the nonlinear stiffness system is different from that of the linear stiffness system.Figure7 (a1) and (b1) are the bifurcation diagrams using Ω1 as the control parameter. Figure 7 (a2) and (b2) are the waterfall diagrams of the rotor system when the imbalance e = 0.002 mm and 0.01 mm, respectively. The bearing clearance is 10 μm, k1 = 3 × 1016 N/m3. Applying the dynamic analytical method similar to the previous section, it can be known that the orange circles in Figure 7 are in the chaotic state. Comparing Figure 7 (a1) and (b1), it can be deduced that as the mass eccentricity increases, the peak response of the system increases; the shape of the chaos is becoming more typical and the range of the chaos is changing wider. Comparing Figure 7 (a2) and (b2), it is found that with the increase in eccentricity, the peak response of the 1X frequency has a sharp increase (from 0.4 to 1.6) and the peak response of the 0.5X frequency has also a large increase (from 0.23 to 0.41); the resonance of the 1X frequency occurs at a larger rotating speed (from 650 rad/s to 1150 rad/s) and the resonance of the 0.5X frequency occurs at higher speed (from 900 rad/s to 1400 rad/s); however, the amplitude and the rotating speed of resonance for the 3X frequency remained stable (0.13 and 200 rad/s). To sum up, the mass eccentricity has an obvious effect on the dynamic behavior, with the increase of the eccentricity, the range and extent of the chaos has significant change, the amplitude and the rotational speed for occurrence of the 1X and 0.5X resonance increases, while the amplitude and the speed for 3X remain unchanged.Figure 7
The bifurcation and waterfall chart when (a)e = 0.002 mm, c0 = 10 μm, and k1 = 3 × 1016 N/m3 and (b) e = 0.01 mm, c0 = 10 μm, and k1 = 3 × 1016 N/m3.
(a)(b)


## 3.1. Influence of Cubic Nonlinear Stiffness

A novelty of this study lies in consideration of the cubic nonlinearity of the shaft. Firstly, the influence of the cubic nonlinear stiffness on the dynamic behaviors will be focused on in this section. The model of the rolling element bearing selected in this paper is JIS6306. The relevant parameters are as follows: the inner radiusR1 = 40.1 mm, the outer radius R2 = 63.9 mm, the number of balls N = 8, and the Hertzian contact stiffness KS=13.34×109N/m3/2. To study the influence of the nonlinear stiffness on the response of the system, the bifurcation diagrams of the system with the different nonlinear stiffness are drawn. The longitudinal coordinate of the bifurcation diagram is the vibration response including the negative response and the positive response for getting the proportional relationship of the two directions (the negative direction and the positive direction). Figures 3 and 4 are the bifurcation diagrams of the nondimensional displacement in x-direction X1 versus the rotational speed Ω1 rad/s when the bearing clearance c0 is 10 μm and 20 μm, respectively. The nonlinear stiffness coefficient is 0 N/m3, 3 × 1015 N/m3, 3 × 1016 N/m3. and 3 × 1017 N/m3, respectively.Figure 3
Bifurcation diagrams of the nondimensional displacement inx-direction X1 versus the rotational speed Ω1 rad/s for the clearance c0 = 10 μm. The jumping phenomenon occurs in the orange circles; chaotic motion occurs in the green circles. (a) k1 = 0 N/m3; (b) k1 = 3 × 1015 N/m3; (c) k1 = 3 × 1016 N/m3; (d) k1 = 3 × 1017 N/m3.
(a)(b)(c)(d)Figure 4
Bifurcation diagrams of the nondimensional displacement inx-direction X1 versus the rotational speed Ω1 rad/s for the clearance c0 = 20 μm. Note that P-1 means the periodic one motion, P-2 means the periodic two motion, and A means the aperiodic motion. (a) k1 = 0 N/m3; (b) k1 = 3 × 1015 N/m3; (c) k1 = 3 × 1016 N/m3; (d) k1 = 3 × 1017 N/m3.
(a)(b)(c)(d)It can be seen from Figures3(a)–3(d) that the system shows complex dynamic behaviors such as the periodic motion, the jump phenomena, quasiperiodic motion, and chaos, and the amplitude of the chaotic region increases and the range of the chaotic region changes with the enhancement of the nonlinear stiffness; in addition, the critical speed of the rotor increases and amplitude of the responses decreases with the increase of the stiffness. From Figure 3(a), it can be shown that there exists an obvious jumping phenomenon at the rotating speed of 760 rad/s. The response begins to change to chaotic from periodic at 1130 rad/s, and the transition path from the periodic one motion (P-1) to chaos is the periodic two motion (P-2) and then the system enters into the periodic one motion (P-1) once again, then enters into the periodic two motion (P-2) once more, and then into chaos, as illustrated in the local enlarged diagram in Figure 3(a). When the cubic nonlinear stiffness coefficient increases to 3 × 1015 N/m3, the speed scope of the jumping phenomenon expands to 775–885 rad/s, as shown in Figure 3(b); with the increase of the stiffness, the speed scope of jumping further expands to 780–1190 rad/s, as shown in Figure 3(c); When the stiffness increases to k1 = 3 × 1017 N/m3, only one jumping phenomenon exists at the beginning at the rotating speed of 792 rad/s and the chaos is more chaotic than that of smaller stiffness, as shown in Figure 3(d). It is crucial for rotational machinery to calculate the transition from the periodic motion to the aperiodic motion (A). The jumping phenomenon should be avoided to make the system run safely and stably in engineering practice.The clearance between the element and races has a certain effect on the nonlinear behaviors of the responses. In order to study the influence of the clearance, the bifurcation diagram for the clearancec0 = 20 μm is drawn, as illustrated in Figure 4. It can be seen from a comparison of Figures 3 and 4 that the clearance has an obvious influence on the dynamic characteristics, and the amplitude and range of the chaos are larger with the increase of the clearance. The periodic one motion (P-1) exists only in a very small speed range, and the aperiodic motion (A) occupies a larger percentage of the rotating speed range comparing with that of Figure 3; it is inferred that the larger the nonlinear stiffness coefficient is, the more obvious the nonlinear characteristics are. Therefore, it is of great significance to reduce the bearing clearance and design the suitable stiffness for the stable operation of the system.


## 3.2. Influence of Rotating Speed

The rotating speed is one of the key parameters affecting the dynamic behavior of the rotor system. In this section, the bifurcation diagrams, the waterfall plots, the time domain diagrams, the phase diagrams, and the Poincaré maps of the rotor system with nonlinear stiffness are drawn to exactly analyze the dynamics.Figure5(a) shows the bifurcation diagram of the nondimensional displacement in x-direction X1 versus the rotational speed Ω1 rad/s for e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3; Figure 5(b) shows the waterfall plot of the vibrational amplitude A versus the frequency ratio (f/fn) with the changing rotating speed Ω1 under the same parameters e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3. Figure 6 shows the time domain diagram, the phase trajectory, and the Poincaré section of the system at each speed where the dynamic characteristics of the system can be obtained more intuitively. It can be seen that the system has the periodic one motion, the periodic two motion, multiple periodic motion, and chaotic motion in the rotating speed span of 0–2000 rad/s, as demonstrated in Figure 5(a). To study the dynamic characteristic of the typical vibration, the typical speeds of the system are selected as 516 rad/s, 1490 rad/s, 1570 rad/s, and 1675 rad/s, and the corresponding diagrams are shown in Figures 6(a)–6(d). It can be observed from Figure 5(b) that the spectrum shows only 1X frequency instead of the 0.5X frequency or combined frequency in the lower Ω1 (0–1380 rad/s) because the varying compliance vibration is weaker when the rotating speed is lower, and the vibration induced by the eccentric force dominates the vibration; when the rotating speed increases, the spectrum diagram shows the combined frequencies including the 0.5X and 1X frequency, as shown in Figure 5(b) and spectrum diagrams of Figure 6.Figure 5
(a) Bifurcation diagram of the nondimensional displacement inx-direction X1 versus the rotational speed Ω1 rad/s for e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3. (b) Waterfall plot of displacement responses under the same parameters e = 0.01 mm, c0 = 5 μm, and k1 = 3 × 1016 N/m3.
(a)(b)Figure 6
Dynamic behaviors at the typical rotating speed. (a) P-1 motion: (a1) time history, (a2) spectrum, chart of axis track, (a3) phase trajectory, and (a4) Poincaré section whenω = 516 rad/s. (b) Multiple periodic motion: (b1) time history, (b2) spectrum, chart of axis track, (b3) phase trajectory, and (b4) Poincaré section when ω = 1490 rad/s. (c) P-2 motion: (c1) time history, (c2) spectrum, chart of axis track, (c3) phase trajectory, and (c4) Poincaré section when ω = 1570 rad/s. (d) Chaotic motion: (d1) time history, (d2) spectrum, chart of axis track, (d3) phase trajectory, and (d4) Poincaré section when ω = 1675 rad/s.
(a)(b)(c)(d)


## 3.3. Influence of Eccentricity

Because of the uneven distribution of rotor mass, the rotor will produce eccentric force during the operation of the system. The existence of eccentric force leads to a big risk for the safety of the system. The exciting forces induced by the mass imbalance are the major source of rotor vibrations in the field, so it is necessary to analyze the effect of eccentricity on the motion characteristics of the rotor-bearing system. In order to study the influence of mass eccentricity on the rotor-bearing system, the nonlinear stiffness system is analyzed in this section. The imbalance is represented by the eccentricity of the rotore. Due to the existence of nonlinear stiffness of the shaft, the displacement and elastic recovery force at the three focused positions are nonlinear. The influence of the mass eccentricity on the vibration characteristics of the nonlinear stiffness system is different from that of the linear stiffness system.Figure7 (a1) and (b1) are the bifurcation diagrams using Ω1 as the control parameter. Figure 7 (a2) and (b2) are the waterfall diagrams of the rotor system when the imbalance e = 0.002 mm and 0.01 mm, respectively. The bearing clearance is 10 μm, k1 = 3 × 1016 N/m3. Applying the dynamic analytical method similar to the previous section, it can be known that the orange circles in Figure 7 are in the chaotic state. Comparing Figure 7 (a1) and (b1), it can be deduced that as the mass eccentricity increases, the peak response of the system increases; the shape of the chaos is becoming more typical and the range of the chaos is changing wider. Comparing Figure 7 (a2) and (b2), it is found that with the increase in eccentricity, the peak response of the 1X frequency has a sharp increase (from 0.4 to 1.6) and the peak response of the 0.5X frequency has also a large increase (from 0.23 to 0.41); the resonance of the 1X frequency occurs at a larger rotating speed (from 650 rad/s to 1150 rad/s) and the resonance of the 0.5X frequency occurs at higher speed (from 900 rad/s to 1400 rad/s); however, the amplitude and the rotating speed of resonance for the 3X frequency remained stable (0.13 and 200 rad/s). To sum up, the mass eccentricity has an obvious effect on the dynamic behavior, with the increase of the eccentricity, the range and extent of the chaos has significant change, the amplitude and the rotational speed for occurrence of the 1X and 0.5X resonance increases, while the amplitude and the speed for 3X remain unchanged.Figure 7
The bifurcation and waterfall chart when (a)e = 0.002 mm, c0 = 10 μm, and k1 = 3 × 1016 N/m3 and (b) e = 0.01 mm, c0 = 10 μm, and k1 = 3 × 1016 N/m3.
(a)(b)


## 4. Conclusion

A rolling rotor-bearing system with cubic nonlinear stiffness is established; in the model, the effects of the unbalanced excitation, the bearing clearance, the nonlinear Hertzian contact force, the varying compliance vibration, and the nonlinear stiffness of shaft material on the system are fully considered. The influences of the bearing clearance, the shaft nonlinear stiffness, and the rotating speed on the dynamic behaviors of the system are studied, and the main conclusions are as follows:(1)
In engineering applications, the bearing clearance should be decreased as much as possible to reduce the risk of the system in chaos and jumping; the critical speed and the response peak of the system can be controlled by adjusting the stiffness of the shaft under specific scopes. The nonlinear characteristics of the system are significantly enhanced with the increase of the nonlinear stiffness. Therefore, the material with lower nonlinear stiffness is more conducive to the stable operation of the system.(2)
The rotor system supported by the rolling bearing not only has the frequency components of varying compliance vibration and the working rotating speed, but also has the frequency components of the frequency division, the frequency multiplication, and frequency combination. When the speed increases, the frequency components corresponding to the speed gradually become the dominant frequency components. The larger the imbalance in the low-speed area and the high-speed area, the more unstable the system response.(3)
The motion of the system will be in different motion states with the increase of the rotating speed under the condition of the same eccentricity. Under the condition of changing the eccentricity, the system usually maintains the same motion state, and only a few motion states will be deviated, usually at the rotating speed of the junction range for the chaotic motion and the periodic motion.(4)
The mass eccentricity has an obvious effect on the dynamic behavior, with the increase of the eccentricity, the range and extent of the chaos has significant change, the amplitude and the rotational speed for the occurrence of the 1X and 0.5X resonance increase, while the amplitude and the speed for 3X remain unchanged.ecam



---

*Source: 8878319-2021-05-25.xml*