# Dynamical Behavior Analysis of Rubbing Rotor System under Asymmetric Oil Film Force

Authors: Youfu Tang; Feng Lin; Qian Zou
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Rubbing is one of the most common and significant faults that exist in the rotor system. It exhibits extremely complicated dynamic behavior. Existing dynamic models are primarily based on the symmetrical rotor structure, in which the oil film forces at both ends of support are considered identical. However, in practice, the oil film force may be differently affected by multiple factors (e.g., viscosity of lubricants, the thickness of oil film, and clearance of journal bearing). In this study, a novel dynamic model of the rubbing rotor system under asymmetric oil film force was developed. Furthermore, based on the nonlinear rotor dynamic theory, its dynamic behavior with different parameters was analyzed, and the corresponding chaotic features were extracted. The results indicated that the evolution law of chaotic motion was more complicated, and the chaotic region of system response was obviously wider under the asymmetric oil film force than the symmetrical oil film force.


---

## Body

## 1. Introduction

In high-speed rotating machinery, the rubbing is one of the most common and negative faults, easily leading to instability of the rotor system. It may cause the equipment to vibrate robustly or even bend and destroy the entire shaft. Thus, the mechanical mechanism and feature extracting of the rubbing rotor system should be studied, which is of great implication to the prevention and elimination of the fault. The rubbing fault of the rotor system is a robust nonlinear phenomenon with complex dynamic behavior [1]. As the rotor speed and load have been improved continuously, the nonlinear characteristics of the rotor system will be of higher implication. Then, in the system, multiple solutions, jumps, quasi-periods, periodic solution bifurcations, and chaotic motions will take place. Researches on rubbing are still not enough on rubbing is still not enough. On the one hand, the friction phenomenon often appears unexpectedly, and it is difficult to perform experiments in the practical field. On the other hand, the rotor rubbing involves the dynamic characteristics of the nonsmooth system, and the rubbing rotor system has different theoretical and numerical analysis methods from those of smooth systems. Though many research results have been achieved, the types of faults and their motion characteristics are still not comprehensively known. In particular, the bifurcation and chaos behavior in the rubbing process has not been fully understood. Accordingly, the nonlinear dynamic behavior of the rubbing rotor system needs in-depth study, which can provide a reliable design and lay a fault diagnosis basis for the safe and stable operation of the rotor system.At present, the rubbing dynamics model is generally split into two major mechanical models, described by constrained differential systems and piecewise smooth systems, respectively [2, 3]. The constrained differential system, focusing on the initial impact phase of the rubbing, is capable of objectively describing the rubbing process from the dynamic perspective. The first complete constrained differential system was developed by Muszynska [4]. He studied the rubbing mechanism in accordance with the law of conservation of momentum. In the modeling process, the additional elastic force and damping force, acting on the rotor during a collision and the displacement of the rubbing process, were neglected. Shaw elucidated the nonlinear dynamics of the single-degree-of-freedom collision system and compared the theoretical analysis with the experimental results to verify the validity of the model [5]. Wang, considering the effects of local surface deformation, contact stress, and elastic wave, optimized the model of the rotating static rubbing model [6]. Beatty proposed a mathematical model for the nonlinear frictional force in a piecewise linear form. Through numerical simulation and experimental verification, the viewpoint of diagnosing such faults has been obtained [7]. Choy and Padovan discussed the mutual rubbing action between the rotor and the stator in the rotor structure. The numerical simulation method was employed to reveal the internal relationship between friction, energy, rubbing process, and collision separation angle, and the effects of stator stiffness, friction coefficient, unbalanced load, and system damping on the rubbing force were analyzed [8]. Kim and Noah, in accordance with the bifurcation theory, used a modified Jeffcott model to determine the starting point of the aperiodic winding motion [9].In the study of the rotor rubbing problem, the stiffness of the shaft is usually assumed to be linear. This assumption can to some extent explain some phenomena of rotor rubbing. However, as new materials with nonlinear physical properties have been used increasingly, the stress-strain relationship of materials is no longer linear and required to be expressed by nonlinear constitutive relations. Shaw studied the vibration of an elongated shaft with nonlinear stiffness theoretically and thus demonstrated the existence of chaotic motion [10]. Cveticanin, based on the cubic polynomial function of displacement, explained the physical nonlinear factors of the material. The averaging method was used to analyze the major resonance response of the nonlinear elastic material rotor system under cyclic excitation, and the rotating shaft of different materials was tested to achieve the frequency response. The curve confirmed the vibration abrupt phenomenon of the rotor system with nonlinear stiffness [11]. Zavodney, using the multiscale method, studied the fixed point and stability of a single-degree-of-freedom system with cubic nonlinear stiffness. The nonlinear problem analysis method was employed to analyze the coexistence of chaotic response and periodic attractors under certain conditions. The phenomenon was verified through numerical integration [12]. Furthermore, when the relative rotational speed of the rotating stator is high, the friction coefficient between the rotating stators will be associated with the relative speed. Thus, the effect of nonlinear friction must be considered when establishing the rotor rub model.In the above studies on the dynamic model of the rubbing rotor, the symmetrical rotor structure has been primarily used; i.e. the oil film forces at both ends have been considered identical. According to the running condition and experimental data of the rotating machine in actual work, the oil film force significantly affects the normal operation of the rotating machine, and the different oil film support at both ends will complicate the dynamic characteristics of the rotor system. Accordingly, it is necessary to study the dynamics of rotor systems supported by different oil films at both ends. Literature [13–15] studied the effect of asymmetric lubricating oil condition on the dynamic behavior of the rotor-bearing system, whereas this model only considered the effect of the viscosity of the lubricating oil on the oil film force, instead of the effect of linear stiffness and nonlinear friction. In fact, there are many nonlinear factors in the rubbing rotor (e.g., nonlinear oil film force, nonlinear stiffness, and nonlinear friction). To more accurately describe the various nonlinear dynamic behaviors of the rotor system, an asymmetrical rotor system dynamics model with multiple nonlinear factors should be developed.When the rubbing rotor system fails, the dynamic process of the system will be complicated, and the vibration response will exhibit rich nonlinear features, especially chaotic motion characteristics. In recent years, scholars have focused on the segmental smooth rubbing model and studied the nonlinear dynamic behavior of the rubbing rotor system from different perspectives. To discuss the effect of rubbing on the rotor system's next harmonic vibration in the parameter excitation, Childs studied the Jeffcott model in polar form and used the approximation method [16]. Begg analyzed the effect of dry friction on rotor motion stability during full ring rubbing [17]. To analyze the ordered and chaotic response of the rotor structure with gap and rubbing, Goldman and Muszynska used the discontinuous variable transformation method and the averaging method [18]. Paidoussis built a model with gap and mass imbalance, considered the dry friction coefficient and rotor imbalance eccentricity, and demonstrated the effect of these two parameters through phase diagram, bifurcation diagram, Poincaré map, and Lyapunov exponent [19]. Adiletta conducted theoretical and experimental studies on rigid rotors, identified chaotic motion, and took measures to suppress chaotic occurrence [20]. Lahmar, by analyzing the dynamics of the eddy bearing, proposed the optimal short bearing theory and proved that, with the increase in the rotational speed, the effect of the eddy on the dynamic behavior of the rotor-bearing dynamic system will be increasingly obvious [21]. Andres explored the nonlinear dynamics of different oil film supported rotor systems [22]. Enrich investigated the chaotic vibration phenomenon of the rotor at high rotational speeds [23]. Goldman studied the chaotic behavior of the rubbing rotor system [24]. Sinha studied rotor whirl movement caused by shaft bending, rotor imbalance, or rotor misalignment, which resulted in radial rubbing [25]. From the above analysis, it is not difficult to see that the rubbing rotor has strong nonlinear dynamic behavior, and the focus of research on the rubbing rotor was primarily based on theoretical analysis and numerical simulation, which can provide an accurate theoretical basis for fault diagnosis of rubbing rotor system.The contributions of this paper are as follows:( 1 ) Based on the consideration of nonlinear oil film force, nonlinear stiffness, and nonlinear friction force, a dynamic model of rubbing rotor system under asymmetric oil film force is established.( 2 ) The nonlinear dynamical behavior of the rubbing rotor system is studied by some numerical analysis methods, such as bifurcation diagram, Poincare section diagram, time domain wave diagram, axle center trail, phase diagram, and frequency spectrum diagram.( 3 ) The effect of some system parameters on bifurcation and chaotic behavior is principally researched. The analysis results show that the effect of asymmetric oil film force on the dynamical behavior of the rotor system has serious influence with rotation speed.


## 2. Dynamic Model of Rubbing Rotor System under Asymmetric Oil Film Force

### 2.1. Failure Mechanism of Rubbing Rotor System

Fault mechanism research acts as an objective basis for condition monitoring and fault diagnosis. During the operation of the rotor, the anomalies or malfunctions are reflected by changes in their status signals. The study of the rotor failure mechanism primarily aims to build a reasonable physical or mathematical model from the perspective of kinematics and dynamics to study the fault causes and state effects of the rotor system. It attempts to master the fault formation and development process of the rotor, to understand its inherent nature and its characteristics, and to provide sensitive parameters and relevant features for fault diagnosis and condition monitoring. The factors inducing rotor rubbing can be generally split into two categories: external factors and internal factors. Under normal circumstances, the vortex amplitude of the rotor is small, and generally no rubbing occurs, but if it is affected by external disturbances (e.g., fundamental vibration, load change, and external damping), the vortex amplitude of the rotor may increase dramatically and lead to rubbing. In the meantime, the eccentricity of the rotor system is another trigger of the rubbing. The mass eccentricity makes the gap between the rotor and the stator too small, leading to the occurrence of the rubbing, which significantly affects the dynamic characteristics of the rotor. The oil film force is a type of nonlinear motion coupling force. Under certain speed conditions, the oil film force will also stimulate the sharp increase in the whirl, thus causing the rubbing to occur. Thus, in this paper, we conduct a detailed nonlinear analysis of the rubbing faults with nonlinear oil film forces.The rotor rubbing fault can cause the rotor to malfunction, and it even has catastrophic consequences for the entire equipment. When the rotor rubbing fault occurs, there will be many nonlinear factors (e.g., nonlinear oil film force, stiffness asymmetry, and gyro effect). The interaction between them complicates the rotor nonlinear coupling frictional fault. When the rotational speed is high, the nonlinear friction will also significantly affect the dynamic characteristics of the rotor system, making the dynamic response of the rotor system more abundant. These nonlinear factors significantly affect the rotor system. Accordingly, the rubbing rotor is a typical nonsmooth nonlinear dynamic system that exhibits complex dynamic characteristics, which is of high implication to its dynamics research.


### 2.2. Description of the Nonlinear Dynamic System

The rubbing rotor system refers to a nonlinear dynamic system, and the nonlinear dynamic system can be expressed by a multiparameter finite dimensional second order ordinary differential equation system. When described by a state vector, the differential equations are written as(1)dudt=Ft,λ,uu∈Rn,t>0where u=q1q2; λ denotes the system parameter.This paper focuses on the study of the stable periodic solution and the bifurcation of the rubbing rotor system, which involves the periodic solution of the system, the subharmonic solution, the almost periodic solution, and chaotic solution. Due to the complexity of the problem, in most cases, the study of nonlinear rotor dynamics problems cannot be solved by analytical methods but by means of numerical analysis only. Numerical methods for periodic boundary value problems include the shooting method, the harmonic balance method, and the difference method. However, as computer computing speed continues to rise, many scientists and technicians in recent years have prioritized the direct numerical integration method. In this study, the Runge-Kutta method was used to solve the motion differential equation.


### 2.3. Nonlinear Oil Film Force Model of Short-Bearing Unsteady

A key issue in the study of oil film instability is to build a model of oil film force being reasonable and convenient for practical application. The final research results are directly affected by the difference in the oil film force model. For the instability caused by the linear oil film force based on eight dynamic characteristic coefficients, the research has been quite mature and used widely. As the rotor is inclined to have a larger size, higher speed, and flexibility, the dynamic behavior of the rotor system will be more complex, and quasi-periodic and chaotic motions will occur, which cannot be explained by the linear oil film force model. Thus, it is inevitable to use the nonlinear oil film force model to study many complex dynamic phenomena in the rotor system.Figure1 gives a schematic view of the journal bearing structure, where O denotes the bearing center; O′ is the center of the journal; Oxy is a fixed coordinate system;O′ηξ is a rotating coordinate system, assuming the center of the journal O′at any moment The position is e1, ϕ, taking polar coordinate system O′rθ, O′ is the pole, ξ is represented as the polar coordinate line of O′O direction, and α is counterclockwise measured based on the minimum oil film thickness. The dimensionless unsteady nonlinear oil film force along the ξ and η directions of the circular axis in polar coordinates is expressed as(2)fξfη=-C1C2C2C3ε′ϕ′-0.5εwhere ε=e1/δ; e1=x2+y2 denotes the eccentricity of the journal center; ϕ is the center deviation angle of the journal; tan⁡ϕ=y/x; τ=Ωt; ε′=dε/dτ; ϕ′=dϕ/dτ; Cii=1,2,3 is a function of ε, ε′ and ϕ′, and the expression is given by equation (7). fξ=Fξ/SW; fη=Fη/SWfη=Fη/SW; Sommerfeld number S=ηΩRL/WR/δ2L/2R2; η is the viscosity coefficient of the lubricating oil; Ω is the angular velocity of the shaft diameter; L is the width of the bearing; δ is the average clearance of the bearing; R is the radius of the journal; W is the static load of the bearing.Figure 1
Schematic diagram of journal bearing.To study the problem conveniently, the expression of dimensionless unsteady nonlinear oil film force in the rotating coordinate system was transformed into a fixed coordinate system. The dimensionless unsteady nonlinear oil film force in the fixed coordinate system is expressed as(3)fXfY=-CX′Y′-KXYSubsequently, the unsteady nonlinear oil film force of the short bearing is expressed as(4)FXFY=ηΩRLRδ2L2R2-CX′Y′-KXYwhere fX=FX/SW; fY=FY/SW. X, Y denote the dimensionless coordinates of the axis; X=x/δ; Y=y/δ; τ=Ωt; X′=x˙/Ωδ; Y′=y˙/Ωδ. In the above and below equation, · indicates that the time t is derived, and (′) means τ is deriving, where c and k are the transient damping matrix and stiffness matrix of the unsteady nonlinear oil film force, and each element in the matrix is a function of the transient disturbance position and the transient disturbance velocity.(5)C=C11C12C21C22,K=12-C2C3-C3-C2(6)C11=C1cos2ϕ+C3sin2ϕ-2C2sin⁡ϕcos⁡ϕC12=C21=C2cos2ϕ-sin2ϕ+C1-C3sin⁡ϕcos⁡ϕC22=C1sin2ϕ+C3cos2ϕ+2C2sin⁡ϕcos⁡ϕFor short bearings,Ci is expressed as(7)C1=4εε′A3A2+2-5ε2ε2ϕ′-1/221-ε22A2-ε4ϕ′-1/222+2+4ε21-ε25/2ΔϕC2=8A4ϕ′-1/23A2-ε4ϕ′-1/222C3=4εε′AA2+ε2-2ε2ϕ′-1/221-ε22A2-ε4ϕ′-1/222+21-ε23/2Δϕwhere A=(ε′)2+(φ′-0.5)2ε2, Δϕ=π+2tan-1(ε′ε/A(1-ε2)1/2), ε=e1/δ=X2+Y2, tan⁡ϕ=Y/X.


### 2.4. Nonlinear Dynamic Model of Rubbing Rotor System

Since the rubbing process of the rotor is sophisticated, involving many factors, this study only considered the effect of nonlinear oil film force, nonlinear friction force, and nonlinear stiffness on the rubbing rotor system. The developed dynamic model is illustrated in Figure2.Figure 2
Dynamic model of rubbing rotor system.
(a)
Plan of rubbing rotor (b)
Model of local rubbingBoth ends of the rotor were supported by a symmetrical sliding bearing.O1 denotes the geometric center of the left bearing bush; O2 is the geometric center of the rotor; Oe is the center of mass of the rotor; O3 is the geometric center of the bearing bush at the right end. kc denotes the stator stiffness;k is the linear stiffness term coefficient of the elastic shaft material; k1 is the nonlinear stiffness term coefficient of the elastic shaft material. c1, c3 refer to the damping coefficient of rotor at the left and right bearings, respectively. c2 is the damping coefficient of the rotor disc. m1, m3 are the concentrated mass of the rotor at the left and right bearings, respectively, and m2 is the equivalent concentrated mass of the rotor at the disc, i.e., the eccentricity of the centroid of the disc. No mass elastic axis exists between the rubbing disc and the bearing, and Fx1, Fy1, Fx3, Fy3 refer to the left and right end sliding bearings, respectively. The nonlinear oil film force acts on the rotating shaft, and Px and Py are the components of the rubbing force in the x and y directions, respectively.The rotor system with rubbing fault is virtually a nonlinear vibration system with piecewise linear stiffness. For the convenience of research, the rubbing rotor was simplified in the study, regardless of the thermal effect of friction and the collision force between the rotors. Besides, the stator was assumed to be elastic deformation. Based on the above assumptions, the rotor local rub model is shown in Figure2(b), where Pn denotes the normal collision force; Pt is the tangential frictional force; θ is the angle between the radial direction of the rubbing point and the x axis.When studying the rubbing force model of the rubbing rotor, it has been usually assumed that when the rubbing occurs, the friction between the rotor and the stator follows the Coulomb friction law, and the friction coefficient is constant. However, when the relative speed of the rotor and the stator is high, the friction coefficient between the rotor and the stator will have a relation to the relative velocity. Then, the relationship between the friction coefficient and the relative velocity can be expressed by the following equation:(8)fn=f+bvnwhere f denotes the friction coefficient when the speed is not considered; b is the speed influence factor; n is the number of speed terms. For the convenience of calculation, n is taken as 1 in this study.Assume that the gap between the rotor and the stator ise0 at rest with no radial rubbing and no rubbing force. When the radial displacement e of the rotor core is larger than the gap e0 of the moving and moving parts, Radial rubbing will occur; its normal collision force Pn and tangential nonlinear friction force Pt are expressed as(9)Pn=kce-e0Pt=f+bvPnwhere kc denotes the stator radial stiffness; f is the friction coefficient when the speed is not considered; v=x˙2+y˙2. For the relative sliding speed between the moving parts, x˙=dx/dt, y˙=dy/dt; b denotes the speed influence factor; e=x2+y2 is the radial displacement of the rotor axis.When the rubbing occurs, the components of thex and y directions decomposed by the rubbing force on the xy coordinate system are, respectively, expressed as(10)PxPy=-cos⁡θsin⁡θ-sin⁡θ-cos⁡θPnPtSubstituting (9) into (10), the friction is transformed to(11)PxPy=-kce-e0e1-f+bvf+bv1xyFor the convenience of discussion, dimensionless displacement was introducedX=x/δ; Y=y/δ; X′=dX/dτ=x˙/δΩ; Y′=dY/dτ=y˙/δΩ;ω denotes the rotational angular velocity of the rotor; τ=Ωt; δ is the average clearance of the bearing; then the friction force is x-y. Relative and dimensionless displacement in the direction X, Y is expressed as(12)PXPY=-kce-e0δe1-f+bvf+bv1XYFor a rotor-bearing system with nonlinear stiffness, the nonlinear stiffness is expressed as a polynomial function of vibration displacement or velocity. In the analysis of the nonlinear vibration of such rotor-bearing systems, each of the polynomials in the vibration equations with polynomial nonlinear stiffness has a clear physical implication, representing different sources of nonlinearity, which can be demonstrated and actualized. The observation and experimental results of the rotor-bearing system are equally rich in nonlinear vibration phenomena, and they are easy to analyze and solve. According to the experimental results, the physical nonlinear factors of the shaft material derived from the neutral axis elongation and symmetry are expressed by the sum of the linear term and the cubic term, respectively, and the nonlinear stress-strain relationship is defined as(13)σ=E1ε+E2ε3Accordingly, the corresponding nonlinear elastic force of the rotating shaft is expressed as(14)Tx=kx+k1x2+y2xTy=ky+k1x2+y2ywhere k and k1 refer to the linear and nonlinear stiffness term coefficients of the shaft material, respectively.Assume that the radial displacement of the shaft section at the left end bearing isx1, y1; the radial displacement of the disc is x2, y2; the radial displacement of the shaft section at the right end bearing is x3, y3. In this paper, the physical nonlinear factors of the shaft material are expressed by the sum of the linear term and the cubic term, and the corresponding dimensioned nonlinear elastic force is defined as(15)Tx1=kx1-x2+k1x1-x22+y1-y22x1-x2Ty1=ky1-y2+k1x1-x22+y1-y22y1-y2Tx2=k2x2-x1-x3+k12x2-x1-x32+2y2-y1-y322x2-x1-x3Ty2=k2y2-y1-y3+k12x2-x1-x32+2y2-y1-y322y2-y1-y3Tx3=kx3-x2+k1x3-x22+y3-y22x3-x2Ty3=ky3-y2+k1x3-x22+y3-y22y3-y2For the convenience of qualitative analysis, the same dimensionless transformation was introduced as the rubbing force and the oil film force to generate the dimensionless nonlinear elastic force.(16)TX1=λ12X1-X2+a1λ12X1-X22+Y1-Y22X1-X2TY1=λ12Y1-Y2+a1λ12X1-X22+Y1-Y22Y1-Y2TX2=λ222X2-X1-X3+a1λ222X2-X1-X32+2Y2-Y1-Y322X2-X1-X3TY2=λ222Y2-Y1-Y3+a1λ222X2-X1-X32+2Y2-Y1-Yy322Y2-Y1-Y3TX3=λ32X3-X2+a1λ32X3-X22+Y3-Y22X3-X2TY3=λ32Y3-Y2+a1λ32X3-X22+Y3-Y22Y3-Y2where λi denotes frequency ratio, λi=1/Ωk/mi=1,2,3; a1 is nonlinear stiffness ratio; a1=k1δ2/k.The differential equation of motion of a rotor system with nonlinear stiffness under the action of asymmetric oil film force is written as(17)m1x¨1+c1x˙1+Tx1x1,x2,y1,y2=Fx1x1,y1,x˙1,y˙1m1y¨1+c1y˙1+Ty1x1,x2,y1,y2=Fy1x1,y1,x˙1,y˙1-m1gm2x¨2+c2x˙2+Tx2x1,x2,x3,y1,y2,y3=Pxx2,y2+m2EΩ2cos⁡Ωtm2y¨2+c2y˙2+Ty2x1,x2,x3,y1,y2,y3=Pyx2,y2+m2EΩ2sin⁡Ωt-m2gm3x¨3+c3x˙3+Tx3x2,x3,y2,y3=Fx3x3,y3,x˙3,y˙3m3y¨3+c3y˙3+Ty3x2,x3,y2,y3=Fy3x3,y3,x˙3,y˙3-m3gwhere Tx1, Ty1, Tx2, Ty2, Tx3 and Ty3 all refer to the dimensional elastic forces of the rotating shaft, which are specifically expressed in equation (15). Px, Py are the components of the rotor system rubbing force in the x, y direction, respectively, as expressed in equation (15). To make the equation universal, dimensionless parameters were introduced as follows:(18)Xi=xiδYi=yiδXi′=x˙iδΩYi′=y˙iδΩXi′′=x¨iδΩ2Yi′′=y¨iδΩ2i=1,2,3Substituting the above parameters into equation (17) for dimensionless transformation, the dimensionless form of the differential equation of motion of the rotor system with nonlinear stiffness under the action of asymmetric oil film force is(19)X1′′+2ζ1λ1X1′+TX1=w1fX1α1Y1′′+2ζ1λ1Y1′+TY1=w1fY1α1-GX2′′+2ζ2λ2X2′+TX2=PXm2δΩ2+pcos⁡τY2′′+2ζ2λ2Y2′+TY2=PYm2δΩ2+psin⁡τ-GX3′′+2ζ3λ3X3′+TX3=w3fX3α3Y3′′+2ζ3λ3Y3′+TY3=w3fY3α3-Gwhere λi denotes the frequency ratio, λi=1/Ωk/mi(i=1,2,3); p is a dimensionless eccentricity; p=E/δ. ζi is the damping ratio, ζi=ci/2mik/mi(i=1,2,3); a1 is a nonlinear stiffness ratio; a1=k1δ2/k. G is a dimensionless external load, G=g/δΩ2. Wi is the static load acting on the bearing. wi denotes the load ratio, wi=Wi/mig(i=1,3); αi=Ω2δ/Sig(i=1,3); Si is the Sommerfeld number. fx1, fy1, fx3 and fy3 refer to the dimensionless nonlinear oil film forces defined by equation (3) in the x and y directions, respectively. The components Px and Py are, respectively, the components of the dimensionless rubbing force in the x and y directions, as expressed in equation (12); Tx1, Ty1, Tx2, Ty2, Tx3 and Ty3 are, respectively, the components of the dimensionless nonlinear elastic force in the x and y directions, which are specifically expressed in equation (16).For the convenience of analysis and calculation, a comprehensive parameterσi and dimensionless speed ω are introduced. σi reflects the effect of various factors, e.g., lubricant viscosity, rotor weight, bearing clearance, and aspect ratio on the nonlinear oil film force; their expressions are (20)σi=αiSi=Wi2δηiRiLi2δg42RiLi4i=1,3(21)ω=Ωδg


## 2.1. Failure Mechanism of Rubbing Rotor System

Fault mechanism research acts as an objective basis for condition monitoring and fault diagnosis. During the operation of the rotor, the anomalies or malfunctions are reflected by changes in their status signals. The study of the rotor failure mechanism primarily aims to build a reasonable physical or mathematical model from the perspective of kinematics and dynamics to study the fault causes and state effects of the rotor system. It attempts to master the fault formation and development process of the rotor, to understand its inherent nature and its characteristics, and to provide sensitive parameters and relevant features for fault diagnosis and condition monitoring. The factors inducing rotor rubbing can be generally split into two categories: external factors and internal factors. Under normal circumstances, the vortex amplitude of the rotor is small, and generally no rubbing occurs, but if it is affected by external disturbances (e.g., fundamental vibration, load change, and external damping), the vortex amplitude of the rotor may increase dramatically and lead to rubbing. In the meantime, the eccentricity of the rotor system is another trigger of the rubbing. The mass eccentricity makes the gap between the rotor and the stator too small, leading to the occurrence of the rubbing, which significantly affects the dynamic characteristics of the rotor. The oil film force is a type of nonlinear motion coupling force. Under certain speed conditions, the oil film force will also stimulate the sharp increase in the whirl, thus causing the rubbing to occur. Thus, in this paper, we conduct a detailed nonlinear analysis of the rubbing faults with nonlinear oil film forces.The rotor rubbing fault can cause the rotor to malfunction, and it even has catastrophic consequences for the entire equipment. When the rotor rubbing fault occurs, there will be many nonlinear factors (e.g., nonlinear oil film force, stiffness asymmetry, and gyro effect). The interaction between them complicates the rotor nonlinear coupling frictional fault. When the rotational speed is high, the nonlinear friction will also significantly affect the dynamic characteristics of the rotor system, making the dynamic response of the rotor system more abundant. These nonlinear factors significantly affect the rotor system. Accordingly, the rubbing rotor is a typical nonsmooth nonlinear dynamic system that exhibits complex dynamic characteristics, which is of high implication to its dynamics research.


## 2.2. Description of the Nonlinear Dynamic System

The rubbing rotor system refers to a nonlinear dynamic system, and the nonlinear dynamic system can be expressed by a multiparameter finite dimensional second order ordinary differential equation system. When described by a state vector, the differential equations are written as(1)dudt=Ft,λ,uu∈Rn,t>0where u=q1q2; λ denotes the system parameter.This paper focuses on the study of the stable periodic solution and the bifurcation of the rubbing rotor system, which involves the periodic solution of the system, the subharmonic solution, the almost periodic solution, and chaotic solution. Due to the complexity of the problem, in most cases, the study of nonlinear rotor dynamics problems cannot be solved by analytical methods but by means of numerical analysis only. Numerical methods for periodic boundary value problems include the shooting method, the harmonic balance method, and the difference method. However, as computer computing speed continues to rise, many scientists and technicians in recent years have prioritized the direct numerical integration method. In this study, the Runge-Kutta method was used to solve the motion differential equation.


## 2.3. Nonlinear Oil Film Force Model of Short-Bearing Unsteady

A key issue in the study of oil film instability is to build a model of oil film force being reasonable and convenient for practical application. The final research results are directly affected by the difference in the oil film force model. For the instability caused by the linear oil film force based on eight dynamic characteristic coefficients, the research has been quite mature and used widely. As the rotor is inclined to have a larger size, higher speed, and flexibility, the dynamic behavior of the rotor system will be more complex, and quasi-periodic and chaotic motions will occur, which cannot be explained by the linear oil film force model. Thus, it is inevitable to use the nonlinear oil film force model to study many complex dynamic phenomena in the rotor system.Figure1 gives a schematic view of the journal bearing structure, where O denotes the bearing center; O′ is the center of the journal; Oxy is a fixed coordinate system;O′ηξ is a rotating coordinate system, assuming the center of the journal O′at any moment The position is e1, ϕ, taking polar coordinate system O′rθ, O′ is the pole, ξ is represented as the polar coordinate line of O′O direction, and α is counterclockwise measured based on the minimum oil film thickness. The dimensionless unsteady nonlinear oil film force along the ξ and η directions of the circular axis in polar coordinates is expressed as(2)fξfη=-C1C2C2C3ε′ϕ′-0.5εwhere ε=e1/δ; e1=x2+y2 denotes the eccentricity of the journal center; ϕ is the center deviation angle of the journal; tan⁡ϕ=y/x; τ=Ωt; ε′=dε/dτ; ϕ′=dϕ/dτ; Cii=1,2,3 is a function of ε, ε′ and ϕ′, and the expression is given by equation (7). fξ=Fξ/SW; fη=Fη/SWfη=Fη/SW; Sommerfeld number S=ηΩRL/WR/δ2L/2R2; η is the viscosity coefficient of the lubricating oil; Ω is the angular velocity of the shaft diameter; L is the width of the bearing; δ is the average clearance of the bearing; R is the radius of the journal; W is the static load of the bearing.Figure 1
Schematic diagram of journal bearing.To study the problem conveniently, the expression of dimensionless unsteady nonlinear oil film force in the rotating coordinate system was transformed into a fixed coordinate system. The dimensionless unsteady nonlinear oil film force in the fixed coordinate system is expressed as(3)fXfY=-CX′Y′-KXYSubsequently, the unsteady nonlinear oil film force of the short bearing is expressed as(4)FXFY=ηΩRLRδ2L2R2-CX′Y′-KXYwhere fX=FX/SW; fY=FY/SW. X, Y denote the dimensionless coordinates of the axis; X=x/δ; Y=y/δ; τ=Ωt; X′=x˙/Ωδ; Y′=y˙/Ωδ. In the above and below equation, · indicates that the time t is derived, and (′) means τ is deriving, where c and k are the transient damping matrix and stiffness matrix of the unsteady nonlinear oil film force, and each element in the matrix is a function of the transient disturbance position and the transient disturbance velocity.(5)C=C11C12C21C22,K=12-C2C3-C3-C2(6)C11=C1cos2ϕ+C3sin2ϕ-2C2sin⁡ϕcos⁡ϕC12=C21=C2cos2ϕ-sin2ϕ+C1-C3sin⁡ϕcos⁡ϕC22=C1sin2ϕ+C3cos2ϕ+2C2sin⁡ϕcos⁡ϕFor short bearings,Ci is expressed as(7)C1=4εε′A3A2+2-5ε2ε2ϕ′-1/221-ε22A2-ε4ϕ′-1/222+2+4ε21-ε25/2ΔϕC2=8A4ϕ′-1/23A2-ε4ϕ′-1/222C3=4εε′AA2+ε2-2ε2ϕ′-1/221-ε22A2-ε4ϕ′-1/222+21-ε23/2Δϕwhere A=(ε′)2+(φ′-0.5)2ε2, Δϕ=π+2tan-1(ε′ε/A(1-ε2)1/2), ε=e1/δ=X2+Y2, tan⁡ϕ=Y/X.


## 2.4. Nonlinear Dynamic Model of Rubbing Rotor System

Since the rubbing process of the rotor is sophisticated, involving many factors, this study only considered the effect of nonlinear oil film force, nonlinear friction force, and nonlinear stiffness on the rubbing rotor system. The developed dynamic model is illustrated in Figure2.Figure 2
Dynamic model of rubbing rotor system.
(a)
Plan of rubbing rotor (b)
Model of local rubbingBoth ends of the rotor were supported by a symmetrical sliding bearing.O1 denotes the geometric center of the left bearing bush; O2 is the geometric center of the rotor; Oe is the center of mass of the rotor; O3 is the geometric center of the bearing bush at the right end. kc denotes the stator stiffness;k is the linear stiffness term coefficient of the elastic shaft material; k1 is the nonlinear stiffness term coefficient of the elastic shaft material. c1, c3 refer to the damping coefficient of rotor at the left and right bearings, respectively. c2 is the damping coefficient of the rotor disc. m1, m3 are the concentrated mass of the rotor at the left and right bearings, respectively, and m2 is the equivalent concentrated mass of the rotor at the disc, i.e., the eccentricity of the centroid of the disc. No mass elastic axis exists between the rubbing disc and the bearing, and Fx1, Fy1, Fx3, Fy3 refer to the left and right end sliding bearings, respectively. The nonlinear oil film force acts on the rotating shaft, and Px and Py are the components of the rubbing force in the x and y directions, respectively.The rotor system with rubbing fault is virtually a nonlinear vibration system with piecewise linear stiffness. For the convenience of research, the rubbing rotor was simplified in the study, regardless of the thermal effect of friction and the collision force between the rotors. Besides, the stator was assumed to be elastic deformation. Based on the above assumptions, the rotor local rub model is shown in Figure2(b), where Pn denotes the normal collision force; Pt is the tangential frictional force; θ is the angle between the radial direction of the rubbing point and the x axis.When studying the rubbing force model of the rubbing rotor, it has been usually assumed that when the rubbing occurs, the friction between the rotor and the stator follows the Coulomb friction law, and the friction coefficient is constant. However, when the relative speed of the rotor and the stator is high, the friction coefficient between the rotor and the stator will have a relation to the relative velocity. Then, the relationship between the friction coefficient and the relative velocity can be expressed by the following equation:(8)fn=f+bvnwhere f denotes the friction coefficient when the speed is not considered; b is the speed influence factor; n is the number of speed terms. For the convenience of calculation, n is taken as 1 in this study.Assume that the gap between the rotor and the stator ise0 at rest with no radial rubbing and no rubbing force. When the radial displacement e of the rotor core is larger than the gap e0 of the moving and moving parts, Radial rubbing will occur; its normal collision force Pn and tangential nonlinear friction force Pt are expressed as(9)Pn=kce-e0Pt=f+bvPnwhere kc denotes the stator radial stiffness; f is the friction coefficient when the speed is not considered; v=x˙2+y˙2. For the relative sliding speed between the moving parts, x˙=dx/dt, y˙=dy/dt; b denotes the speed influence factor; e=x2+y2 is the radial displacement of the rotor axis.When the rubbing occurs, the components of thex and y directions decomposed by the rubbing force on the xy coordinate system are, respectively, expressed as(10)PxPy=-cos⁡θsin⁡θ-sin⁡θ-cos⁡θPnPtSubstituting (9) into (10), the friction is transformed to(11)PxPy=-kce-e0e1-f+bvf+bv1xyFor the convenience of discussion, dimensionless displacement was introducedX=x/δ; Y=y/δ; X′=dX/dτ=x˙/δΩ; Y′=dY/dτ=y˙/δΩ;ω denotes the rotational angular velocity of the rotor; τ=Ωt; δ is the average clearance of the bearing; then the friction force is x-y. Relative and dimensionless displacement in the direction X, Y is expressed as(12)PXPY=-kce-e0δe1-f+bvf+bv1XYFor a rotor-bearing system with nonlinear stiffness, the nonlinear stiffness is expressed as a polynomial function of vibration displacement or velocity. In the analysis of the nonlinear vibration of such rotor-bearing systems, each of the polynomials in the vibration equations with polynomial nonlinear stiffness has a clear physical implication, representing different sources of nonlinearity, which can be demonstrated and actualized. The observation and experimental results of the rotor-bearing system are equally rich in nonlinear vibration phenomena, and they are easy to analyze and solve. According to the experimental results, the physical nonlinear factors of the shaft material derived from the neutral axis elongation and symmetry are expressed by the sum of the linear term and the cubic term, respectively, and the nonlinear stress-strain relationship is defined as(13)σ=E1ε+E2ε3Accordingly, the corresponding nonlinear elastic force of the rotating shaft is expressed as(14)Tx=kx+k1x2+y2xTy=ky+k1x2+y2ywhere k and k1 refer to the linear and nonlinear stiffness term coefficients of the shaft material, respectively.Assume that the radial displacement of the shaft section at the left end bearing isx1, y1; the radial displacement of the disc is x2, y2; the radial displacement of the shaft section at the right end bearing is x3, y3. In this paper, the physical nonlinear factors of the shaft material are expressed by the sum of the linear term and the cubic term, and the corresponding dimensioned nonlinear elastic force is defined as(15)Tx1=kx1-x2+k1x1-x22+y1-y22x1-x2Ty1=ky1-y2+k1x1-x22+y1-y22y1-y2Tx2=k2x2-x1-x3+k12x2-x1-x32+2y2-y1-y322x2-x1-x3Ty2=k2y2-y1-y3+k12x2-x1-x32+2y2-y1-y322y2-y1-y3Tx3=kx3-x2+k1x3-x22+y3-y22x3-x2Ty3=ky3-y2+k1x3-x22+y3-y22y3-y2For the convenience of qualitative analysis, the same dimensionless transformation was introduced as the rubbing force and the oil film force to generate the dimensionless nonlinear elastic force.(16)TX1=λ12X1-X2+a1λ12X1-X22+Y1-Y22X1-X2TY1=λ12Y1-Y2+a1λ12X1-X22+Y1-Y22Y1-Y2TX2=λ222X2-X1-X3+a1λ222X2-X1-X32+2Y2-Y1-Y322X2-X1-X3TY2=λ222Y2-Y1-Y3+a1λ222X2-X1-X32+2Y2-Y1-Yy322Y2-Y1-Y3TX3=λ32X3-X2+a1λ32X3-X22+Y3-Y22X3-X2TY3=λ32Y3-Y2+a1λ32X3-X22+Y3-Y22Y3-Y2where λi denotes frequency ratio, λi=1/Ωk/mi=1,2,3; a1 is nonlinear stiffness ratio; a1=k1δ2/k.The differential equation of motion of a rotor system with nonlinear stiffness under the action of asymmetric oil film force is written as(17)m1x¨1+c1x˙1+Tx1x1,x2,y1,y2=Fx1x1,y1,x˙1,y˙1m1y¨1+c1y˙1+Ty1x1,x2,y1,y2=Fy1x1,y1,x˙1,y˙1-m1gm2x¨2+c2x˙2+Tx2x1,x2,x3,y1,y2,y3=Pxx2,y2+m2EΩ2cos⁡Ωtm2y¨2+c2y˙2+Ty2x1,x2,x3,y1,y2,y3=Pyx2,y2+m2EΩ2sin⁡Ωt-m2gm3x¨3+c3x˙3+Tx3x2,x3,y2,y3=Fx3x3,y3,x˙3,y˙3m3y¨3+c3y˙3+Ty3x2,x3,y2,y3=Fy3x3,y3,x˙3,y˙3-m3gwhere Tx1, Ty1, Tx2, Ty2, Tx3 and Ty3 all refer to the dimensional elastic forces of the rotating shaft, which are specifically expressed in equation (15). Px, Py are the components of the rotor system rubbing force in the x, y direction, respectively, as expressed in equation (15). To make the equation universal, dimensionless parameters were introduced as follows:(18)Xi=xiδYi=yiδXi′=x˙iδΩYi′=y˙iδΩXi′′=x¨iδΩ2Yi′′=y¨iδΩ2i=1,2,3Substituting the above parameters into equation (17) for dimensionless transformation, the dimensionless form of the differential equation of motion of the rotor system with nonlinear stiffness under the action of asymmetric oil film force is(19)X1′′+2ζ1λ1X1′+TX1=w1fX1α1Y1′′+2ζ1λ1Y1′+TY1=w1fY1α1-GX2′′+2ζ2λ2X2′+TX2=PXm2δΩ2+pcos⁡τY2′′+2ζ2λ2Y2′+TY2=PYm2δΩ2+psin⁡τ-GX3′′+2ζ3λ3X3′+TX3=w3fX3α3Y3′′+2ζ3λ3Y3′+TY3=w3fY3α3-Gwhere λi denotes the frequency ratio, λi=1/Ωk/mi(i=1,2,3); p is a dimensionless eccentricity; p=E/δ. ζi is the damping ratio, ζi=ci/2mik/mi(i=1,2,3); a1 is a nonlinear stiffness ratio; a1=k1δ2/k. G is a dimensionless external load, G=g/δΩ2. Wi is the static load acting on the bearing. wi denotes the load ratio, wi=Wi/mig(i=1,3); αi=Ω2δ/Sig(i=1,3); Si is the Sommerfeld number. fx1, fy1, fx3 and fy3 refer to the dimensionless nonlinear oil film forces defined by equation (3) in the x and y directions, respectively. The components Px and Py are, respectively, the components of the dimensionless rubbing force in the x and y directions, as expressed in equation (12); Tx1, Ty1, Tx2, Ty2, Tx3 and Ty3 are, respectively, the components of the dimensionless nonlinear elastic force in the x and y directions, which are specifically expressed in equation (16).For the convenience of analysis and calculation, a comprehensive parameterσi and dimensionless speed ω are introduced. σi reflects the effect of various factors, e.g., lubricant viscosity, rotor weight, bearing clearance, and aspect ratio on the nonlinear oil film force; their expressions are (20)σi=αiSi=Wi2δηiRiLi2δg42RiLi4i=1,3(21)ω=Ωδg


## 3. Dynamic Behavior of Rubbing Rotor System Caused by Different Integrated Parameters

### 3.1. The Overall Analysis

Due to the nonlinear nonlinearity, nonlinear rubbing force, and nonlinear force of the oil film force, the fourth-order Runge-Kutta method was employed in this study to numerically analyze equation (18) and obtains the rotor system under different parameters. The vibration response which explains the effect of parameter changes on the chaotic motion of the rubbing rotor system. Since this paper is based on the nonlinear stiffness, nonlinear rubbing force, and nonlinear oil film force, the mechanical model was developed. Thus, when selecting the bifurcation parameters, it is primarily the main influencing factors of the three; i.e., the selection is infinite. The rotational speed ω, the dimensionless eccentricity p, the nonlinear stiffness ratio a1, the stator stiffness kc, and the velocity influence factor b were taken as the bifurcation parameters. The integration step size for each cycle in the calculation was 0.01π. A total of 800 cycles were calculated, the first 500 cycles were discarded, and the next 300 cycles were taken. The parameters of the rotor system were the following: f=0.2, δ=0.1mm, c1=1050N·s/m, c2=2100N·s/m, c3=1050N·s/m, k=6.0×106 N/m, e0=0.05 mm, critical speed Ω0=774.6 rad/s, and dimensionless critical speed ratio ω0=2.47.The comprehensive parameters reflect the effect of various factors, e.g., lubricant viscosity, rotor weight, bearing clearance, and aspect ratio on the nonlinear oil film force. The effect of comprehensive parameters on the nonlinear dynamic behavior of the rotor system should be clarified. Assuming that the nonlinear stiffness ratio of the rotor systema1=0.3, the speed influence factor b=0.01, the dimensionless eccentricity p=0.325, and the stator stiffness kc=4.5×107N/m, other system parameters remain constant. In this section, the nonlinear dynamic behavior of the rubbing rotor system with different comprehensive parameters was analyzed from two aspects. On the one hand, the bifurcation diagram was used to analyze the evolution law of the system chaotic motion. On the other hand, the time domain waveform graph and phase were employed. Graphs, spectrograms, and Poincare sections were used to analyze the motion state of the system in detail.Figures3(a)–3(d) give the bifurcation diagrams of the displacement response of the Y-direction disc with the dimensionless rotation speed of the rubbing rotor system under different combined parameters. The bifurcation diagram reveals that when the rotational speed is lower than the critical speed, the effect of the oil film force is small, and the system motion is primarily periodic motion. As the rotational speed continues to rise, after the critical speed is reached, the oil film force acts to make the rotor move. The state becomes complicated, and chaotic motion and double-cycle motion begin to appear. When the dimensionless speed is greater than 3.72, i.e., the actual working speed reaches 1.5 times of the critical speed, the system appears oil film whirl under the action of oil film force. Figure 3(b) suggests that the amplitude has increased. With the further increase in the rotational speed, the oil film whirl makes the system rubbing fault more serious, and the system motion has a large range of chaotic motion and quasi-periodic motion. It can be preliminarily analyzed from the bifurcation diagram that the rotor system response leads to chaotic motion in the form of paroxysmal bifurcation. To further judge the evolution law of the chaotic motion of the system, the Poincare section diagram is incremented by 0.1 for each increment of 0~10, so that the evolution law of chaotic motion of different parameter combination systems changes within a given speed range. When the comprehensive parameters are the same as shown in Figure 3(a), the system motion change follows the law: period 1 motion → period 2 motion → chaotic motion → periodic 5 motion → chaotic motion → periodic motion 2 → chaotic motion → quasi-periodic motion → Cycle 3 motion → chaotic motion. When the combination of integrated parameters changes, i.e., the oil film forces at both ends of the rotor system are different, the chaotic region is more obvious, and the period of motion is included in the chaotic region, as shown in Figure 3(b). When the two ends of the rotor are supported by different supports, the road of the system moving into and out of chaos also changes. In Figure 3(c), when σ1=50,σ3=400, the system motion change follows the law: cycle 1 motion → cycle 2 motion → cycle 1 motion → chaotic motion (intermediate multicycle motion) → cycle 2 motion → cycle 5 motion → chaotic motion → quasi-periodic motion → periodic 1 motion → periodic 3 motion → chaotic motion; the middle contains multiperiod motion; in Figure 3(d), when σ1=100,σ3=400, the system motion change follows the law: period 1 motion → period 2 motion → chaotic motion (including 2 chaotic islands) → multiperiod motion → chaotic motion (intermediate multiperiod motion) → periodic motion 2 → chaotic motion (multicycle motion in the middle) → quasi-periodic motion → periodic 1 motion → periodic 3 motion.Figure 3
Bifurcation of rubbing rotor system at different integrated parameters.
(a)
σ 1 = 150, σ3=150 (b)
σ 1 = 150, σ3=300 (c)
σ 1 = 50, σ3=400 (d)
σ 1 = 100, σ3=400


### 3.2. The Motion State Analysis of the System

Figures4(a) and 5(a) show the chaotic motion form of the rubbing rotor in detail when the integrated parameters σ1=150 and σ3=150. On the axis trajectory map and phase diagram, the rotor motion is obviously disordered at this time, the continuum component is more obvious on the spectrogram, the frequency components are concentrated at the low frequency, and the corresponding Poincare cross-section is a distribution point with a certain geometric shape set. Figures 4(b) and 5(b) give the time domain waveform diagram, the axial center trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when ω=6.2, respectively, and the axial trajectory and the phase trajectory are in a space loop of one width. With the movement, the obvious “shooting vibration” phenomenon can be observed on the time domain waveform diagram, and the corresponding Poincare sectional diagram displays a closed curve, revealing that the rotor system performed quasi-periodic motion at this time. Figures 4(c) and 5(c) show the time domain waveform diagram, the axial center trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when ω=6.8, respectively. The axial trajectory diagram and the phase diagram are single loop lines and the spectrum, respectively. There is only one spectral line on the graph, and one fixed point is displayed on the corresponding Poincare cross-section, which suggests that the rotor system performed the period 1 motion. Figures 4(d) and 5(d) are the time domain waveform diagram, the axial center trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when ω=7.5, respectively, and the axial trajectory diagram and phase diagram are composed of 3 loop lines. Three discrete spectral lines are clear on the spectrogram, and three fixed points are displayed on the corresponding Poincare cross-section, implying that the rotor system performed periodic motion.Figure 4
System response at different nondimension rotation speeds whenσ1=σ3=150.
(a)
ω = 3.2 (b)
ω = 6.2  (c)
ω = 6.8 (d)
ω = 7.5Figure 5
Poincare section at different nondimension rotation speeds whenσ1=σ3=150.
(a)
ω = 3.2 (b)
ω = 6.2 (c)
ω = 6.8 (d)
ω = 7.5


## 3.1. The Overall Analysis

Due to the nonlinear nonlinearity, nonlinear rubbing force, and nonlinear force of the oil film force, the fourth-order Runge-Kutta method was employed in this study to numerically analyze equation (18) and obtains the rotor system under different parameters. The vibration response which explains the effect of parameter changes on the chaotic motion of the rubbing rotor system. Since this paper is based on the nonlinear stiffness, nonlinear rubbing force, and nonlinear oil film force, the mechanical model was developed. Thus, when selecting the bifurcation parameters, it is primarily the main influencing factors of the three; i.e., the selection is infinite. The rotational speed ω, the dimensionless eccentricity p, the nonlinear stiffness ratio a1, the stator stiffness kc, and the velocity influence factor b were taken as the bifurcation parameters. The integration step size for each cycle in the calculation was 0.01π. A total of 800 cycles were calculated, the first 500 cycles were discarded, and the next 300 cycles were taken. The parameters of the rotor system were the following: f=0.2, δ=0.1mm, c1=1050N·s/m, c2=2100N·s/m, c3=1050N·s/m, k=6.0×106 N/m, e0=0.05 mm, critical speed Ω0=774.6 rad/s, and dimensionless critical speed ratio ω0=2.47.The comprehensive parameters reflect the effect of various factors, e.g., lubricant viscosity, rotor weight, bearing clearance, and aspect ratio on the nonlinear oil film force. The effect of comprehensive parameters on the nonlinear dynamic behavior of the rotor system should be clarified. Assuming that the nonlinear stiffness ratio of the rotor systema1=0.3, the speed influence factor b=0.01, the dimensionless eccentricity p=0.325, and the stator stiffness kc=4.5×107N/m, other system parameters remain constant. In this section, the nonlinear dynamic behavior of the rubbing rotor system with different comprehensive parameters was analyzed from two aspects. On the one hand, the bifurcation diagram was used to analyze the evolution law of the system chaotic motion. On the other hand, the time domain waveform graph and phase were employed. Graphs, spectrograms, and Poincare sections were used to analyze the motion state of the system in detail.Figures3(a)–3(d) give the bifurcation diagrams of the displacement response of the Y-direction disc with the dimensionless rotation speed of the rubbing rotor system under different combined parameters. The bifurcation diagram reveals that when the rotational speed is lower than the critical speed, the effect of the oil film force is small, and the system motion is primarily periodic motion. As the rotational speed continues to rise, after the critical speed is reached, the oil film force acts to make the rotor move. The state becomes complicated, and chaotic motion and double-cycle motion begin to appear. When the dimensionless speed is greater than 3.72, i.e., the actual working speed reaches 1.5 times of the critical speed, the system appears oil film whirl under the action of oil film force. Figure 3(b) suggests that the amplitude has increased. With the further increase in the rotational speed, the oil film whirl makes the system rubbing fault more serious, and the system motion has a large range of chaotic motion and quasi-periodic motion. It can be preliminarily analyzed from the bifurcation diagram that the rotor system response leads to chaotic motion in the form of paroxysmal bifurcation. To further judge the evolution law of the chaotic motion of the system, the Poincare section diagram is incremented by 0.1 for each increment of 0~10, so that the evolution law of chaotic motion of different parameter combination systems changes within a given speed range. When the comprehensive parameters are the same as shown in Figure 3(a), the system motion change follows the law: period 1 motion → period 2 motion → chaotic motion → periodic 5 motion → chaotic motion → periodic motion 2 → chaotic motion → quasi-periodic motion → Cycle 3 motion → chaotic motion. When the combination of integrated parameters changes, i.e., the oil film forces at both ends of the rotor system are different, the chaotic region is more obvious, and the period of motion is included in the chaotic region, as shown in Figure 3(b). When the two ends of the rotor are supported by different supports, the road of the system moving into and out of chaos also changes. In Figure 3(c), when σ1=50,σ3=400, the system motion change follows the law: cycle 1 motion → cycle 2 motion → cycle 1 motion → chaotic motion (intermediate multicycle motion) → cycle 2 motion → cycle 5 motion → chaotic motion → quasi-periodic motion → periodic 1 motion → periodic 3 motion → chaotic motion; the middle contains multiperiod motion; in Figure 3(d), when σ1=100,σ3=400, the system motion change follows the law: period 1 motion → period 2 motion → chaotic motion (including 2 chaotic islands) → multiperiod motion → chaotic motion (intermediate multiperiod motion) → periodic motion 2 → chaotic motion (multicycle motion in the middle) → quasi-periodic motion → periodic 1 motion → periodic 3 motion.Figure 3
Bifurcation of rubbing rotor system at different integrated parameters.
(a)
σ 1 = 150, σ3=150 (b)
σ 1 = 150, σ3=300 (c)
σ 1 = 50, σ3=400 (d)
σ 1 = 100, σ3=400


## 3.2. The Motion State Analysis of the System

Figures4(a) and 5(a) show the chaotic motion form of the rubbing rotor in detail when the integrated parameters σ1=150 and σ3=150. On the axis trajectory map and phase diagram, the rotor motion is obviously disordered at this time, the continuum component is more obvious on the spectrogram, the frequency components are concentrated at the low frequency, and the corresponding Poincare cross-section is a distribution point with a certain geometric shape set. Figures 4(b) and 5(b) give the time domain waveform diagram, the axial center trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when ω=6.2, respectively, and the axial trajectory and the phase trajectory are in a space loop of one width. With the movement, the obvious “shooting vibration” phenomenon can be observed on the time domain waveform diagram, and the corresponding Poincare sectional diagram displays a closed curve, revealing that the rotor system performed quasi-periodic motion at this time. Figures 4(c) and 5(c) show the time domain waveform diagram, the axial center trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when ω=6.8, respectively. The axial trajectory diagram and the phase diagram are single loop lines and the spectrum, respectively. There is only one spectral line on the graph, and one fixed point is displayed on the corresponding Poincare cross-section, which suggests that the rotor system performed the period 1 motion. Figures 4(d) and 5(d) are the time domain waveform diagram, the axial center trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when ω=7.5, respectively, and the axial trajectory diagram and phase diagram are composed of 3 loop lines. Three discrete spectral lines are clear on the spectrogram, and three fixed points are displayed on the corresponding Poincare cross-section, implying that the rotor system performed periodic motion.Figure 4
System response at different nondimension rotation speeds whenσ1=σ3=150.
(a)
ω = 3.2 (b)
ω = 6.2  (c)
ω = 6.8 (d)
ω = 7.5Figure 5
Poincare section at different nondimension rotation speeds whenσ1=σ3=150.
(a)
ω = 3.2 (b)
ω = 6.2 (c)
ω = 6.8 (d)
ω = 7.5


## 4. Dynamic Behavior of Rubbing Rotor System Caused by Different Eccentric

### 4.1. The Overall Analysis

In rotating machinery, the rotor, due to processing and assembly, may be eccentric, causing the rotor system to vibrate and cause rubbing, especially for large or high-speed rotating machinery. When the vibration caused by eccentricity exceeds a certain limit, it will often lead to operational failure, even a catastrophic accident. Accordingly, the rotor eccentricity is a critical parameter directly related to the dynamic balance of the rotor and the stability of the operation, which is true especially for the frictional rotor. In the nonlinear rotor dynamics, the bifurcation and chaos of the rotor eccentricity response to the rotor system behavior have an important effect.The eccentricity is the main factor affecting the unbalanced force of the rotor system. It is important to understand the effect of the eccentricity on the rotor system response. The rotor system parameters areω=3.3, a1=0.3, b=0.01, and kc=4.5×107 N/m. This section will analyze the effect of eccentricity on the nonlinear dynamic behavior of the rubbing rotor system from two aspects. On the one hand, the combination of different comprehensive parameters, the bifurcation diagram and the local analysis diagram with the dimensionless eccentricity as the bifurcation parameter, analyzes the effect of the eccentricity change on the motion state of the system from the whole and the local; on the other hand, they make different eccentricity. The bifurcation diagram of the bifurcation parameter is used to analyze the effect of different eccentricity on the evolution law of chaotic motion from the whole motion process.Figures6(a)–6(c) give the bifurcation diagrams of the displacement response in the Y direction of the rubbing rotor system with the dimensionless eccentricity under different combined parameters, respectively. The bifurcation diagram suggests that the eccentricity was very small, and the rotor system performed a stable periodic motion. With the rise in the eccentricity, the area of ​​the rubbing part increased gradually, thereby aggravating the rubbing of the rotor, and thus the rotor system began to appear chaotic motion and period-doubling motion. To further verify the evolution law of chaotic motion of the system, in the dimensionless eccentricity between 0 and 0.5, the Poincare section diagram was added every 0.005. The bifurcation diagram reveals that the rotor system with different parameters entered and left the chaotic system. The path of motion is not exactly the same. The main difference is that the parameters of the transformation of the motion form were different, and the range of chaotic regions that appear is also different. When the comprehensive parameters are the same (see Figure 6(a)), the system motion change follows the law: period 1 motion→chaotic motion→cycle 4 motion→chaotic motion (middle-cycle motion in the middle)→cycle 4 motion→chaotic motion. When the combination of integrated parameters changes, i.e., the oil film forces at both ends of the rotor system are different, the range of the chaotic region becomes wider. See Figure 6(b). The system motion change follows the law: period 1 motion → chaotic motion → period 5 motion →Chaotic motion→Period 4 motion→Chaotic motion→Period 1 motion→Chaotic motion (multiperiod motion in the middle). When σ1=50,σ3=400 (see Figure 6(c)), the system motion change follows the law: period 1 motion→chaotic motion→period 5 motion→chaotic motion→cycle 4 motion→chaos Exercise → Cycle 1 movement.Figure 6
Bifurcation of rubbing rotor system at different integrated parameter combination.
(a)
σ 1 = 150 , σ 3 = 150 (b)
σ 1 = 100 , σ 3 = 150 (c)
σ 1 = 50 , σ 3 = 400


### 4.2. The Motion State Analysis of the System Based on Poincare Section

In this section, when the integrated parameters of the rubbing rotor systemσ1=100,σ3=150, the time domain waveform diagram, the axial trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram under different eccentricities are elucidated. The motion characteristics of the system chaos and double period were analyzed, as shown in Figures 7 and 8.Figure 7
System response at different nondimension rotation speeds whenσ1=100, σ3=150.
(a)
p = 0.3 (b)
p = 0.37  (c)
p = 0.4675Figure 8
Poincare section at different nondimension rotation speeds whenσ1=100, σ3=150.
(a)
p = 0.3 (b)
p = 0.37 (c)
p = 0.4675Figures7(a) and 8(a) show the period 5 motion pattern of the rubbing rotor in detail when the eccentricity p=0.3. Five intertwined loop lines are observed on the trajectory and phase diagram of the axis. Five major discrete spectral lines appear on the spectrogram, and five fixed points are displayed on the corresponding Poincare cross-section.Figures7(b) and 8(b) show the time domain waveform, axial trajectory, phase diagram, spectrogram, and Poincare section at p=0.37, respectively. The axial trajectory and phase trajectory are disordered. There exists a continuous line on the corresponding line, and the corresponding Poincare section is a set of distribution points with a certain geometry, suggesting that the rotor system performed chaotic motion.Figures7(c) and 8(c) give the time domain waveform diagram, the axial trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when p=0.4675, respectively. The axial trajectory diagram and the phase diagram are composed of single loop lines, the spectrum is only one spectral line on the graph, and one fixed point is displayed on the corresponding Poincare cross-section. It suggests that the rotor system performed the period 1 motion.


## 4.1. The Overall Analysis

In rotating machinery, the rotor, due to processing and assembly, may be eccentric, causing the rotor system to vibrate and cause rubbing, especially for large or high-speed rotating machinery. When the vibration caused by eccentricity exceeds a certain limit, it will often lead to operational failure, even a catastrophic accident. Accordingly, the rotor eccentricity is a critical parameter directly related to the dynamic balance of the rotor and the stability of the operation, which is true especially for the frictional rotor. In the nonlinear rotor dynamics, the bifurcation and chaos of the rotor eccentricity response to the rotor system behavior have an important effect.The eccentricity is the main factor affecting the unbalanced force of the rotor system. It is important to understand the effect of the eccentricity on the rotor system response. The rotor system parameters areω=3.3, a1=0.3, b=0.01, and kc=4.5×107 N/m. This section will analyze the effect of eccentricity on the nonlinear dynamic behavior of the rubbing rotor system from two aspects. On the one hand, the combination of different comprehensive parameters, the bifurcation diagram and the local analysis diagram with the dimensionless eccentricity as the bifurcation parameter, analyzes the effect of the eccentricity change on the motion state of the system from the whole and the local; on the other hand, they make different eccentricity. The bifurcation diagram of the bifurcation parameter is used to analyze the effect of different eccentricity on the evolution law of chaotic motion from the whole motion process.Figures6(a)–6(c) give the bifurcation diagrams of the displacement response in the Y direction of the rubbing rotor system with the dimensionless eccentricity under different combined parameters, respectively. The bifurcation diagram suggests that the eccentricity was very small, and the rotor system performed a stable periodic motion. With the rise in the eccentricity, the area of ​​the rubbing part increased gradually, thereby aggravating the rubbing of the rotor, and thus the rotor system began to appear chaotic motion and period-doubling motion. To further verify the evolution law of chaotic motion of the system, in the dimensionless eccentricity between 0 and 0.5, the Poincare section diagram was added every 0.005. The bifurcation diagram reveals that the rotor system with different parameters entered and left the chaotic system. The path of motion is not exactly the same. The main difference is that the parameters of the transformation of the motion form were different, and the range of chaotic regions that appear is also different. When the comprehensive parameters are the same (see Figure 6(a)), the system motion change follows the law: period 1 motion→chaotic motion→cycle 4 motion→chaotic motion (middle-cycle motion in the middle)→cycle 4 motion→chaotic motion. When the combination of integrated parameters changes, i.e., the oil film forces at both ends of the rotor system are different, the range of the chaotic region becomes wider. See Figure 6(b). The system motion change follows the law: period 1 motion → chaotic motion → period 5 motion →Chaotic motion→Period 4 motion→Chaotic motion→Period 1 motion→Chaotic motion (multiperiod motion in the middle). When σ1=50,σ3=400 (see Figure 6(c)), the system motion change follows the law: period 1 motion→chaotic motion→period 5 motion→chaotic motion→cycle 4 motion→chaos Exercise → Cycle 1 movement.Figure 6
Bifurcation of rubbing rotor system at different integrated parameter combination.
(a)
σ 1 = 150 , σ 3 = 150 (b)
σ 1 = 100 , σ 3 = 150 (c)
σ 1 = 50 , σ 3 = 400


## 4.2. The Motion State Analysis of the System Based on Poincare Section

In this section, when the integrated parameters of the rubbing rotor systemσ1=100,σ3=150, the time domain waveform diagram, the axial trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram under different eccentricities are elucidated. The motion characteristics of the system chaos and double period were analyzed, as shown in Figures 7 and 8.Figure 7
System response at different nondimension rotation speeds whenσ1=100, σ3=150.
(a)
p = 0.3 (b)
p = 0.37  (c)
p = 0.4675Figure 8
Poincare section at different nondimension rotation speeds whenσ1=100, σ3=150.
(a)
p = 0.3 (b)
p = 0.37 (c)
p = 0.4675Figures7(a) and 8(a) show the period 5 motion pattern of the rubbing rotor in detail when the eccentricity p=0.3. Five intertwined loop lines are observed on the trajectory and phase diagram of the axis. Five major discrete spectral lines appear on the spectrogram, and five fixed points are displayed on the corresponding Poincare cross-section.Figures7(b) and 8(b) show the time domain waveform, axial trajectory, phase diagram, spectrogram, and Poincare section at p=0.37, respectively. The axial trajectory and phase trajectory are disordered. There exists a continuous line on the corresponding line, and the corresponding Poincare section is a set of distribution points with a certain geometry, suggesting that the rotor system performed chaotic motion.Figures7(c) and 8(c) give the time domain waveform diagram, the axial trajectory diagram, the phase diagram, the spectrogram, and the Poincare section diagram when p=0.4675, respectively. The axial trajectory diagram and the phase diagram are composed of single loop lines, the spectrum is only one spectral line on the graph, and one fixed point is displayed on the corresponding Poincare cross-section. It suggests that the rotor system performed the period 1 motion.


## 5. Conclusions

Our results showed that, considering the nonlinear oil film force, the rubbing rotor system could appear chaotic motion state in a large parameter range, and the evolution law of chaotic motion would be more complicated. With the increase in the rotational speed and eccentricity, the chaotic region of the system response under the action of asymmetric oil film force was obviously wider, much larger than that of the system under the action of symmetrical oil film force.



---

*Source: 1253646-2019-07-14.xml*