Nonlinear Dynamic Analysis of Rotor-Bearing System with Cubic Nonlinearity

Nonlinear dynamic characteristics of a rotor-bearing system with cubic nonlinearity are investigated. +e 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.+e expression with the linear and the cubic nonlinear terms is adopted to characterize the synthetical nonlinearity of the rotor-bearing system.+e 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. +e 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. +e 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.+e research will contribute to a comprehensive understanding of the nonlinear dynamics of the rotor-bearing system.


Introduction
In rotational machinery, the rotor is commonly supported by the rolling element bearing. e 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. erefore, 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][6][7], ball size variation [8][9][10], surface waviness [11][12][13][14][15], the number of balls [16], imbalance [17,18], and internal radial clearance [19][20][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 rotorbearing 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. e nonlinear bearing force is replaced with an equivalent describing function gain. e 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. e 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. e research objective of this paper is to analyze the dynamic response of the rotor-bearing system under consideration of the cubic nonlinear factors.

Rolling Element Bearing.
e nonlinear dynamic in this research focuses on the global dynamics of the rotor-bearing system instead of the local characteristics. erefore, 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. e first step is to build the model of the rolling element bearing. Figure 1 shows the schematic diagram of the rolling bearing model. e rolling bearing includes the inner race, the outer race, the middle rolling elements, and the cage. e ball in the track has both the revolution around the center of the bearing (O 1 ) and the rotation around the center of the ball. e initial position of the mass center for the rotor at the bearing is O 1 . During the operation of the rotor, the position of the mass center for the rotor is O 2 . As shown in Figure 1, the inner race radius is R 1 and the outer race radius is R 2 . e 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 v 2 � 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 v 1 � ω 1 R 1 . 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 by e angular velocity of the ball revolution is given by erefore, the position angle of the nth (n � 1, 2, . . ., n) element is as follows: erefore, the radial contact deformation between the nth element and the inner race is given by where x j is the horizontal displacement of the rotor at the bearing at the left or the right end of the shaft, y j is the vertical displacement of the rotor at the bearing at the left or the right end of the shaft, and c 0 is the clearance between the bearing track and the element. Since the contact force exists  Shock and Vibration only when the contact deformation is positive, the contact force between each element and the inner race can be expressed as follows: where F n is the bearing force of the nth element on the rotor and K S is the Hertzian contact stiffness of the bearing. e components of the nth element bearing force on the rotor in the horizontal direction x and the vertical direction y are, respectively, as follows: e resultant force of the contact force between each element and the inner race is the force of the rolling bearing on the shaft, consequently, where F nX is the component of the bearing force for the nth element on the rotor in the horizontal direction and F nY 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.
ere 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]. e team has done a lot of research studies on rotor cubic nonlinearity [31]. e 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 In this model, the rotor is simplified into three lumped masses. e positions of points P, 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. m P , m R , and m L represent the equivalent lumped mass of the rotor at points P, R, and L, respectively; c p is the damping coefficients of the rotor at point P; c b is the damping coefficients of the rotor at both ends of the bearing. e horizontal and the vertical relative displacements at the points P, R, and L are x P , y P , x R , y R , x L , and y L , respectively. e horizontal and the vertical relative displacement differences between the point P and the point R are, respectively, as follows: e horizontal and the vertical relative displacement differences between the point P and the point L are, respectively, as follows: According to Newton's second theorem, the dynamic equation of the system can be obtained as follows: where e is the eccentricity of the disc; ω 1 is the angular velocity of the rotor; F XR and F YR 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 F XL and F YL 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.

Shock and Vibration
For the convenience of the research, equation (11) is dealt with dimensionless form by introducing the following terms: 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. e dimensionless governing equations are as follows:

Results and Analysis
e 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. e relevant calculation parameters of the rotor system are shown in Table 1. e dimension parameter δ 2 is 4 × 10 −5 mm. e system is analyzed by means of the bifurcation diagram, spectrum diagram, axis track diagram, Poincaré map, and the waterfall chart.

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. e model of the rolling element bearing selected in this paper is JIS6306. e relevant parameters are as follows: the inner radius R 1 � 40.1 mm, the outer radius R 2 � 63.9 mm, the number of balls N � 8, and the Hertzian contact stiffness K S � 13.34 × 10 9 N/m 3/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. e 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 X 1 versus the rotational speed Ω 1 rad/s when the bearing clearance c 0 is 10 μm and 20 μm, respectively. e nonlinear stiffness coefficient is 0 N/m 3 , 3 × 10 15 N/m 3 , 3 × 10 16 N/m 3 . and 3 × 10 17 N/m 3 , respectively.
It can be seen from Figures 3(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. e 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 × 10 15 N/m 3 , 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 k 1 � 3 × 10 17 N/m 3 , 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). e jumping phenomenon should be avoided to make the system run safely and stably in engineering practice. e 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 clearance c 0 � 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. e 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 Shock and Vibration 5 nonlinear characteristics are. erefore, it is of great significance to reduce the bearing clearance and design the suitable stiffness for the stable operation of the system.

Influence of Rotating Speed.
e 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. Figure 5(a) shows the bifurcation diagram of the nondimensional displacement in x-direction X 1 versus the rotational speed Ω 1 rad/s for e � 0.01 mm, c 0 � 5 μm, and k 1 � 3 × 10 16 N/m 3 ; Figure 5(b) shows the waterfall plot of the vibrational amplitude A versus the frequency ratio (f/f n ) with the changing rotating speed Ω 1 under the same parameters e � 0.01 mm, c 0 � 5 μm, and k 1 � 3 × 10 16 N/m 3 . 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.

Influence of Eccentricity.
Because of the uneven distribution of rotor mass, the rotor will produce eccentric force during the operation of the system. e existence of eccentric force leads to a big risk for the safety of the system. e 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. e imbalance is represented by the eccentricity of the rotor e. Due to the existence of nonlinear stiffness of the shaft, the displacement and elastic recovery force at the three focused positions are nonlinear. e influence of the mass eccentricity on the vibration characteristics of the nonlinear stiffness system is different from that of the linear stiffness system. Figure 7 (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. e bearing clearance is 10 μm, k 1 � 3 × 10 16 N/m 3 . 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 2: Schematic diagram of the rotor-bearing system.

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. e 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. e nonlinear characteristics of the system are significantly enhanced with the increase of the nonlinear stiffness. erefore, the material with lower nonlinear stiffness is more conducive to the stable operation of the system.
(2) e 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. e larger the imbalance in the low-speed area and the highspeed area, the more unstable the system response. (3) e 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) e 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

Data Availability
All data, models, and code generated or used during the study appear in the submitted article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.