# Dynamic Analysis of a Hybrid Squeeze Film Damper Mounted Rub-Impact Rotor-Stator System

Authors: Cai-Wan Chang-Jian
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

An investigation is carried out on the systematic analysis of the dynamic behavior of the hybrid squeeze-film damper (HSFD) mounted a rotor-bearing system with strongly nonlinear oil-film force and nonlinear rub-impact force in the present study. The dynamic orbits of the system are observed using bifurcation diagrams plotted using the dimensionless rotating speed ratio as control parameters. The onset of chaotic motion is identified from the phase diagrams, power spectra, Poincaré maps, bifurcation diagrams, maximum Lyapunov exponents, and fractal dimension of the rotor-bearing system. The dynamic behaviors are unlike the usual ways into chaos (1T⇒2T⇒4T⇒8T⇒16T⇒32T⋯⇒ chaos or periodic ⇒ quasi-periodic ⇒ chaotic), it suddenly gets in chaos from the periodic motion without any transition. The results presented in this study provide some useful insights into the design and development of a rotor-bearing system for rotating machinery that operates in highly rotating speed and highly nonlinear regimes.


---

## Body

## 1. Introduction

Squeeze-film damper (SFD) bearing is actually a special type of journal bearing with its journal mechanically prevented from rotating but free to vibrate within the clearance space. The hybrid squeeze-film damper (HSFD) and the porous squeeze-film damper (PSFD) are the well-known applications of SFD and also useful for industry. Some literatures discussed dynamic behaviors in SFD bearings and also found many interesting and useful results. Holmes et al. [1] published a paper dealing with aperiodic behavior in journal bearings and what may very well have been the first paper about aperiodic behavior in journal bearing systems. Nikolajsent and Holmes [2] reported their observation of nonsynchronous vibrations in a test rig of a flexible, symmetric rotor on two identical plain journal bearings supported by centralized squeeze-film dampers. Sykes and Holmes [3] showed experimental observations of subharmonic motion in squeeze film bearings and linked this to possible precursors of chaotic motion. At the same time, Kim and Noah [4] analyzed the bifurcation of a modified Jeffcott rotor with bearing clearance. Ehrich [5] used a simple numerical model of a Jeffcott rotor mounted on a nonlinear spring. It was found that the vibratory response in the transition zone midway between adjacent zones of subharmonic response has all the characteristics of chaotic behavior. Zhao et al. [6] discussed the subharmonic and quasiperiodic motions of an eccentric squeeze film damper-mounted rigid rotor system. Brown [7] studied a simple model of a rigid and hydrodynamically supported journal bearing, using a short bearing theory. Theoretical and experimental investigations were reported by Adiletta et al. [8–10] in which a rigid rotor in short bearings would have subharmonic, quasiperiodic, and chaotic motion for suitable values of the system parameters. Sundararajan and Noah [11] proposed a simple shooting scheme along with an arc-length continuation algorithm with applications to periodically forced rotor systems. The occurrence of periodic, quasiperiodic and chaotic motions was predicted for various ranges of rotor speeds. Chang-Jian and Chen [12–16] presented a series of papers discussing about flexible rotor supported by journal bearings under nonlinear suspension and also combined with rub-impact effect, turbulent effect and micropolar lubricant into consideration. They found very bountiful nonperiodic responses occurring in rotor-bearing systems, and the studies would help engineers or scientists escape undesired motions in either designing or operating rotor-bearing systems.Although virtually all physical phenomena in the real world can be regarded as nonlinear, most of these phenomena can be simplified to a linear form given a sufficiently precise linearization technique. However, this simplification is inappropriate for high-power, high rotating speed system and its application during the design and analysis stage may result in a flawed or potentially dangerous operation. As a result, nonlinear analysis methods are generally preferred within engineering and academic circles. The current study performs a nonlinear analysis of the dynamic behavior of a rotor-bearing system equipped with hybrid squeeze-film damper under nonlinear rub-impact force effect. The nondimensional equation of the rotor-bearing system is then solved using the fourth-order Runge-Kutta method. The nonperiodic behavior of this system is characterized using phase diagrams, power spectra, Poincaré maps, bifurcation diagrams, Lyapunov exponents, and the fractal dimension of the system.


## 2. Mathematical Modeling

Figure1 shows a rotor supported on HSFDs in parallel with retaining springs. The bearing consists of four hydrostatic chambers and four hydrodynamic regions. The oil film supporting force is dependent on the integrated action of hydrodynamic pressure and hydrostatic pressure of HSFD. Figure 2(a) represents the cross-section of HSFD and rub-impact rotor-stator model. The structure of this kind bearing should be popularized to consist of 2N (N=2,3,4…) hydrostatic chambers and 2N hydrodynamic regions. In this study, oil pressure distribution model in the HSFD is proposed to integrate the pressure distribution of dynamic pressure region and static pressure region.Figure 1
Schematic illustration of hybrid squeeze-film damper mounted the rotor-bearing system.Cross-section of HSFD model and rub-impact rotor-stator model.
(a)
(b)


### 2.1. The Instant Oil Film Supporting Force for HSFD

To analyze the pressure distribution, the Reynolds equation for constant lubricant properties and noncompressibility should be assumed, then the Reynolds equation is introduced as follows [12]:
(2.1)1R2∂∂θ(h3∂p∂θ)+1R2∂∂z(h3∂p∂z)=-12μΩ∂h∂θ+12μ∂h∂t.
The supporting region of HSFD should be divided into three regions: static pressure region, rotating direction dynamic pressure region, and axial direction dynamic pressure region, as shown in Figure 2. In the part of HSFD with -a≤z≤a, the long bearing theory is assumed and Reynolds equation is solved with the boundary condition of static pressure region pc,i acquiring the pressure distribution p0(θ). In the part of HSFD with a≤|z|≤L/2, the short bearing theory is assumed and solves the Reynolds equation with the boundary condition of p(z,θ)|z=±a=p0(θ) and p(z,θ)|z=±L/2=0, yielding the pressure distribution in axis direction dynamic pressure region p(z,θ). Finally, a formula of pressure distribution in whole supporting region is obtained.According to the above conditions, the instant oil film pressure distribution is as follows. The instant pressure in rotating direction within the range of-a≤z≤a is
(2.2)p0(θ)={pc,i;π2(i-1)-β2-φb≤θ≤π2(i-1)+β2-φbpi(θ);π2(i-1)+β2-φb≤θ≤π2i-β2-φb,i=1,2,3,4,
where
(2.3)pi(θ)=pc,i+6μR2δ2ε˙ε[1(1+εcos⁡⁡θ)2-1(1+εcos⁡⁡θi1)2]+c1∫θi1θ1δ3(1+εcos⁡⁡θ)3dθ-∫θi1θ12φ˙bμδεR2cos⁡⁡θδ3(1+εcos⁡⁡θ)3dθi=1,2,3,4,(2.4)c1=pc,i+1-pc,i-(6μR2/δ2)(ε˙/ε)[1/(1+εcosθi2)2-1/(1+εcosθi1)2]∫θi1θi2(1/δ3(1+εcosθ)3)dθ+∫θi1θi2(12φ˙bμδεR2cosθ/δ(1+εcosθ)3)dθ∫θi1θi2(1/δ3(1+εcosθ)3)dθ,θi1=(i-1)π2+β2-φb,θi2=iπ2-β2-φb,i=1,2,3,4.The instant pressure in the axis direction within the range ofa≤|z|≤L/2 is
(2.5)p(θ,z)=(L2-|z|){[A1(θ)φ˙bε+A2(θ)ε˙](a-|z|)+p0(θ)1L/2-a},
where (2.6)A1(θ)=6μδsin⁡⁡θδ3(1+εcos⁡⁡θ)3;A2(θ)=6μδcos⁡⁡θδ3(1+εcos⁡⁡θ)3.The instant oil film forces of the different elements are determined by integrating (2.2) and (2.5) over the area of the journal sleeve. In the static pressure region, the forces are
(2.7)Frs=∑i=14pci2aR[sin(π2(i-1)+β2-φb)-sin(π2(i-1)-β2-φb)],Fτs=∑i=14pci2aR[cos⁡⁡(π2(i-1)-β2-φb)-cos⁡⁡(π2(i-1)+β2-φb)].
In the rotating direction dynamics pressure region, the forces are
(2.8)Frc=∑i=14∫θi1θi2pi(θ)R2acosθdθ,Fτc=∑i=14∫θi1θi2pi(θ)R2asin⁡⁡θdθ.
In the axial direction dynamic pressure region, the forces are
(2.9)Fra=∫-L/2-adz∫02πp(θ,z)Rcosθdθ+∫L/2adz∫02πp(θ,z)Rcosθdθ,Fτa=∫-L/2-adz∫02πp(θ,z)Rsin⁡⁡θdθ+∫L/2adz∫02πp(θ,z)Rsin⁡⁡θdθ.
The resulting damper forces in the radial and tangential directions are determined by summing the above supporting forces. It is as follows:
(2.10)Fr=Frs+Frc+Fra,Fτ=Fτs+Fτc+Fτa.


### 2.2. Rub-Impact Force

Figure2(b) shows the radial impact force f1 and the tangential rub force f2. f1 and f2 could be expressed as [17]
(2.11)f1=(e-δ)kcf2=(f+bv)f1,ife≥δ.
Then we could get the rub-impact forces in the horizontal and vertical directions as follows:
(2.12)Rx=-(e-δ)kce[X-(f+bv)Y]Ry=-(e-δ)kce[(f+bv)X+Y].


### 2.3. Dynamics Equation

The equations of rotor motion in the Cartesian coordinates can be written as(2.13)mx¨+dx˙+kx=mρω2cosωt+fx+kx0+Rx,my¨+dy˙+ky=mρω2sin⁡ωt+fy+ky0+Ry.
The origin of the o-xyz-coordinate system is taken to be the bearing center Ob. Dividing these two equations by mcω2 and defining a nondimensional time ϕ=ωt and a speed parameter s=ω/ωn, one obtains the following nondimensionalized equations of motion:
(2.14)X′′+DsX′+1s2X=Ucos⁡⁡ϕ+Bs(XFr-YFτε)+X0s2+Rxmcω2,(2.15)Y′′+DsY′+1s2Y=Usinϕ+Bs(YFr+XFτε)+Y0s2+Rymcω2.
Equations (2.14)~(2.15) describe a nonlinear dynamic system. In the current study, the approximate solutions of these coupled nonlinear differential equations are obtained using the fourth-order Runge-Kutta numerical scheme.


## 2.1. The Instant Oil Film Supporting Force for HSFD

To analyze the pressure distribution, the Reynolds equation for constant lubricant properties and noncompressibility should be assumed, then the Reynolds equation is introduced as follows [12]:
(2.1)1R2∂∂θ(h3∂p∂θ)+1R2∂∂z(h3∂p∂z)=-12μΩ∂h∂θ+12μ∂h∂t.
The supporting region of HSFD should be divided into three regions: static pressure region, rotating direction dynamic pressure region, and axial direction dynamic pressure region, as shown in Figure 2. In the part of HSFD with -a≤z≤a, the long bearing theory is assumed and Reynolds equation is solved with the boundary condition of static pressure region pc,i acquiring the pressure distribution p0(θ). In the part of HSFD with a≤|z|≤L/2, the short bearing theory is assumed and solves the Reynolds equation with the boundary condition of p(z,θ)|z=±a=p0(θ) and p(z,θ)|z=±L/2=0, yielding the pressure distribution in axis direction dynamic pressure region p(z,θ). Finally, a formula of pressure distribution in whole supporting region is obtained.According to the above conditions, the instant oil film pressure distribution is as follows. The instant pressure in rotating direction within the range of-a≤z≤a is
(2.2)p0(θ)={pc,i;π2(i-1)-β2-φb≤θ≤π2(i-1)+β2-φbpi(θ);π2(i-1)+β2-φb≤θ≤π2i-β2-φb,i=1,2,3,4,
where
(2.3)pi(θ)=pc,i+6μR2δ2ε˙ε[1(1+εcos⁡⁡θ)2-1(1+εcos⁡⁡θi1)2]+c1∫θi1θ1δ3(1+εcos⁡⁡θ)3dθ-∫θi1θ12φ˙bμδεR2cos⁡⁡θδ3(1+εcos⁡⁡θ)3dθi=1,2,3,4,(2.4)c1=pc,i+1-pc,i-(6μR2/δ2)(ε˙/ε)[1/(1+εcosθi2)2-1/(1+εcosθi1)2]∫θi1θi2(1/δ3(1+εcosθ)3)dθ+∫θi1θi2(12φ˙bμδεR2cosθ/δ(1+εcosθ)3)dθ∫θi1θi2(1/δ3(1+εcosθ)3)dθ,θi1=(i-1)π2+β2-φb,θi2=iπ2-β2-φb,i=1,2,3,4.The instant pressure in the axis direction within the range ofa≤|z|≤L/2 is
(2.5)p(θ,z)=(L2-|z|){[A1(θ)φ˙bε+A2(θ)ε˙](a-|z|)+p0(θ)1L/2-a},
where (2.6)A1(θ)=6μδsin⁡⁡θδ3(1+εcos⁡⁡θ)3;A2(θ)=6μδcos⁡⁡θδ3(1+εcos⁡⁡θ)3.The instant oil film forces of the different elements are determined by integrating (2.2) and (2.5) over the area of the journal sleeve. In the static pressure region, the forces are
(2.7)Frs=∑i=14pci2aR[sin(π2(i-1)+β2-φb)-sin(π2(i-1)-β2-φb)],Fτs=∑i=14pci2aR[cos⁡⁡(π2(i-1)-β2-φb)-cos⁡⁡(π2(i-1)+β2-φb)].
In the rotating direction dynamics pressure region, the forces are
(2.8)Frc=∑i=14∫θi1θi2pi(θ)R2acosθdθ,Fτc=∑i=14∫θi1θi2pi(θ)R2asin⁡⁡θdθ.
In the axial direction dynamic pressure region, the forces are
(2.9)Fra=∫-L/2-adz∫02πp(θ,z)Rcosθdθ+∫L/2adz∫02πp(θ,z)Rcosθdθ,Fτa=∫-L/2-adz∫02πp(θ,z)Rsin⁡⁡θdθ+∫L/2adz∫02πp(θ,z)Rsin⁡⁡θdθ.
The resulting damper forces in the radial and tangential directions are determined by summing the above supporting forces. It is as follows:
(2.10)Fr=Frs+Frc+Fra,Fτ=Fτs+Fτc+Fτa.


## 2.2. Rub-Impact Force

Figure2(b) shows the radial impact force f1 and the tangential rub force f2. f1 and f2 could be expressed as [17]
(2.11)f1=(e-δ)kcf2=(f+bv)f1,ife≥δ.
Then we could get the rub-impact forces in the horizontal and vertical directions as follows:
(2.12)Rx=-(e-δ)kce[X-(f+bv)Y]Ry=-(e-δ)kce[(f+bv)X+Y].


## 2.3. Dynamics Equation

The equations of rotor motion in the Cartesian coordinates can be written as(2.13)mx¨+dx˙+kx=mρω2cosωt+fx+kx0+Rx,my¨+dy˙+ky=mρω2sin⁡ωt+fy+ky0+Ry.
The origin of the o-xyz-coordinate system is taken to be the bearing center Ob. Dividing these two equations by mcω2 and defining a nondimensional time ϕ=ωt and a speed parameter s=ω/ωn, one obtains the following nondimensionalized equations of motion:
(2.14)X′′+DsX′+1s2X=Ucos⁡⁡ϕ+Bs(XFr-YFτε)+X0s2+Rxmcω2,(2.15)Y′′+DsY′+1s2Y=Usinϕ+Bs(YFr+XFτε)+Y0s2+Rymcω2.
Equations (2.14)~(2.15) describe a nonlinear dynamic system. In the current study, the approximate solutions of these coupled nonlinear differential equations are obtained using the fourth-order Runge-Kutta numerical scheme.


## 3.Analytical Tools for Observing Nonlinear Dynamics of Rotor-Bearing System

In the present study, the nonlinear dynamics of the rotor-bearing system equipped with HSFD shown in Figure1 are analyzed using Poincaré maps, bifurcation diagrams, the Lyapunov exponent and the fractal dimension. The basic principles of each analytical method are reviewed in the following subsections.


### 3.1. Dynamic Trajectories and Poincaré Maps

The dynamic trajectories of the rotor-bearing system provide a basic indication as to whether the system behavior is periodic or nonperiodic. However, they are unable to identify the onset of chaotic motion. Accordingly, some other form of analytical method is required. In the current study, the dynamics of the rotor-bearing system are analyzed using Poincaré maps derived from the Poincaré section of the rotor system. A Poincaré section is a hypersurface in the state-space transverse to the flow of the system of interest. In nonautonomous systems, points on the Poincaré section represent the return points of a time series corresponding to a constant intervalT, where T is the driving period of the excitation force. The projection of the Poincaré section on the y(nT) plane is referred to as the Poincaré map of the dynamic system. When the system performs quasiperiodic motion, the return points in the Poincaré map form a closed curve. For chaotic motion, the return points form a fractal structure comprising many irregularly distributed points. Finally, for nT-periodic motion, the return points have the form of n discrete points.


### 3.2. Power Spectrum

In this study, the spectrum components of the motion performed by the rotor-bearing system are analyzed by using the Fast Fourier Transformation method to derive the power spectrum of the displacement of the dimensionless dynamic transmission error. In the analysis, the frequency axis of the power spectrum plot is normalized using the rotating speed,ω.


### 3.3. Bifurcation Diagram

A bifurcation diagram summarizes the essential dynamics of a rotor-train system and is therefore a useful means of observing its nonlinear dynamic response. In the present analysis, the bifurcation diagrams are generated using two different control parameters, namely the dimensionless unbalance coefficient,β, and the dimensionless rotating speed ratio, s, respectively. In each case, the bifurcation control parameter is varied with a constant step, and the state variables at the end of one integration step are taken as the initial values for the next step. The corresponding variations of the y(nT) coordinates of the return points in the Poincaré map are then plotted to form the bifurcation diagram.


### 3.4. Lyapunov Exponent

The Lyapunov exponent of a dynamic system characterizes the rate of separation of infinitesimally close trajectories and provides a useful test for the presence of chaos. In a chaotic system, the points of nearby trajectories starting initially within a sphere of radiusε0 form after time t an approximately ellipsoidal distribution with semiaxes of length εj(t). The Lyapunov exponents of a dynamic system are defined by λj=limt→∞(1/t)log(εj(t)/ε0), where λj denotes the rate of divergence of the nearby trajectories. The exponents of a system are usually ordered into a Lyapunov spectrum, that is, λ1>λ2>⋯>λm. A positive value of the maximum Lyapunov exponent (λ1) is generally taken as an indication of chaotic motion [16].


### 3.5. Fractal Dimension

The presence of chaotic vibration in a system is generally detected using either the Lyapunov exponent or the fractal dimension property. The Lyapunov exponent test can be used for both dissipative systems and nondissipative (i.e. conservative) systems, but is not easily applied to the analysis of experimental data. Conversely, the fractal dimension test can only be used for dissipative systems but is easily applied to experimental data. In contrast to Fourier transform-based techniques and bifurcation diagrams, which provide only a general indication of the change from periodic motion to chaotic behavior, dimensional measures allow chaotic signals to be differentiated from random signals. Although many dimensional measures have been proposed, the most commonly applied measure is the correlation dimensiondG defined by Grassberger and Procaccia due to its computational speed and the consistency of its results. However, before the correlation dimension of a dynamic system flow can be evaluated, it is first necessary to generate a time series of one of the system variables using a time-delayed pseudo-phase-plane method. Assume an original time series of xi={x(iτ);i=1,2,3,…N}, where τ is the time delay (or sampling time). If the system is acted upon by an excitation force with a frequency ω, the sampling time, τ, is generally chosen such that it is much smaller than the driving period. The delay coordinates are then used to construct an n-dimensional vector X=(x(jτ),x[(j+1)τ],x[(j+2)τ],…,x[(j+n-1)τ]), where j=1,2,3,…(N-n+1). The resulting vector comprises a total of (N-n+1) vectors, which are then plotted in an n-dimensional embedding space. Importantly, the system flow in the reconstructed n-dimensional phase space retains the dynamic characteristics of the system in the original phase space. In other words, if the system flow has the form of a closed orbit in the original phase plane, it also forms a closed path in the n-dimensional embedding space. Similarly, if the system exhibits a chaotic behavior in the original phase plane, its path in the embedding space will also be chaotic. The characteristics of the attractor in the n-dimensional embedding space are generally tested using the function ∑i,j=1NH(r-|xi-xj|) to determine the number of pairs (i,j) lying within a distance |xi-xj|<rin{xi}i=1N, where H denotes the Heaviside step function, N represents the number of data points, and r is the radius of an n-dimensional hypersphere. For many attractors, this function exhibits a power law dependence on r as r→0, that is c(r)∝rdG. Therefore, the correlation dimension, dG, can be determined from the slope of a plot of [logc(r)] versus [logr]. Chen and Yau [18] showed that the correlation dimension represents the lower bound to the capacity or fractal dimension dc and approaches its value asymptotically when the attracting set is distributed more uniformly in the embedding phase space. A set of points in the embedding space is said to be fractal if its dimension has a finite noninteger value. Otherwise, the attractor is referred to as a “strange attractor.” To establish the nature of the attractor, the embedding dimension is progressively increased, causing the slope of the characteristic curve to approach a steady-state value. This value is then used to determine whether the system has a fractal structure or a strange attractor structure. If the dimension of the system flow is found to be fractal (i.e. to have a noninteger value), the system is judged to be chaotic.In the current study, the attractors in the embedding space were constructed using a total of 60000 data points taken from the time series corresponding to the displacement of the system. Via a process of trial and error, the optimum delay time when constructing the time series was found to correspond to one third of a revolution of the system. The reconstructed attractors were placed in embedding spaces with dimensions ofn = 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20, respectively, yielding 10 different [logc(r)] versus [logr] plots for each attractor. The number of data points chosen for embedding purposes (i.e., 60000) reflects the need for a compromise between the computation time and the accuracy of the results. In accordance with Grassberger and Procaccia [19], the number of points used to estimate the intrinsic dimension of the attracting set in the current analysis is less than 42M, where M is the greatest integer value less than the fractal dimension of the attracting set.


## 3.1. Dynamic Trajectories and Poincaré Maps

The dynamic trajectories of the rotor-bearing system provide a basic indication as to whether the system behavior is periodic or nonperiodic. However, they are unable to identify the onset of chaotic motion. Accordingly, some other form of analytical method is required. In the current study, the dynamics of the rotor-bearing system are analyzed using Poincaré maps derived from the Poincaré section of the rotor system. A Poincaré section is a hypersurface in the state-space transverse to the flow of the system of interest. In nonautonomous systems, points on the Poincaré section represent the return points of a time series corresponding to a constant intervalT, where T is the driving period of the excitation force. The projection of the Poincaré section on the y(nT) plane is referred to as the Poincaré map of the dynamic system. When the system performs quasiperiodic motion, the return points in the Poincaré map form a closed curve. For chaotic motion, the return points form a fractal structure comprising many irregularly distributed points. Finally, for nT-periodic motion, the return points have the form of n discrete points.


## 3.2. Power Spectrum

In this study, the spectrum components of the motion performed by the rotor-bearing system are analyzed by using the Fast Fourier Transformation method to derive the power spectrum of the displacement of the dimensionless dynamic transmission error. In the analysis, the frequency axis of the power spectrum plot is normalized using the rotating speed,ω.


## 3.3. Bifurcation Diagram

A bifurcation diagram summarizes the essential dynamics of a rotor-train system and is therefore a useful means of observing its nonlinear dynamic response. In the present analysis, the bifurcation diagrams are generated using two different control parameters, namely the dimensionless unbalance coefficient,β, and the dimensionless rotating speed ratio, s, respectively. In each case, the bifurcation control parameter is varied with a constant step, and the state variables at the end of one integration step are taken as the initial values for the next step. The corresponding variations of the y(nT) coordinates of the return points in the Poincaré map are then plotted to form the bifurcation diagram.


## 3.4. Lyapunov Exponent

The Lyapunov exponent of a dynamic system characterizes the rate of separation of infinitesimally close trajectories and provides a useful test for the presence of chaos. In a chaotic system, the points of nearby trajectories starting initially within a sphere of radiusε0 form after time t an approximately ellipsoidal distribution with semiaxes of length εj(t). The Lyapunov exponents of a dynamic system are defined by λj=limt→∞(1/t)log(εj(t)/ε0), where λj denotes the rate of divergence of the nearby trajectories. The exponents of a system are usually ordered into a Lyapunov spectrum, that is, λ1>λ2>⋯>λm. A positive value of the maximum Lyapunov exponent (λ1) is generally taken as an indication of chaotic motion [16].


## 3.5. Fractal Dimension

The presence of chaotic vibration in a system is generally detected using either the Lyapunov exponent or the fractal dimension property. The Lyapunov exponent test can be used for both dissipative systems and nondissipative (i.e. conservative) systems, but is not easily applied to the analysis of experimental data. Conversely, the fractal dimension test can only be used for dissipative systems but is easily applied to experimental data. In contrast to Fourier transform-based techniques and bifurcation diagrams, which provide only a general indication of the change from periodic motion to chaotic behavior, dimensional measures allow chaotic signals to be differentiated from random signals. Although many dimensional measures have been proposed, the most commonly applied measure is the correlation dimensiondG defined by Grassberger and Procaccia due to its computational speed and the consistency of its results. However, before the correlation dimension of a dynamic system flow can be evaluated, it is first necessary to generate a time series of one of the system variables using a time-delayed pseudo-phase-plane method. Assume an original time series of xi={x(iτ);i=1,2,3,…N}, where τ is the time delay (or sampling time). If the system is acted upon by an excitation force with a frequency ω, the sampling time, τ, is generally chosen such that it is much smaller than the driving period. The delay coordinates are then used to construct an n-dimensional vector X=(x(jτ),x[(j+1)τ],x[(j+2)τ],…,x[(j+n-1)τ]), where j=1,2,3,…(N-n+1). The resulting vector comprises a total of (N-n+1) vectors, which are then plotted in an n-dimensional embedding space. Importantly, the system flow in the reconstructed n-dimensional phase space retains the dynamic characteristics of the system in the original phase space. In other words, if the system flow has the form of a closed orbit in the original phase plane, it also forms a closed path in the n-dimensional embedding space. Similarly, if the system exhibits a chaotic behavior in the original phase plane, its path in the embedding space will also be chaotic. The characteristics of the attractor in the n-dimensional embedding space are generally tested using the function ∑i,j=1NH(r-|xi-xj|) to determine the number of pairs (i,j) lying within a distance |xi-xj|<rin{xi}i=1N, where H denotes the Heaviside step function, N represents the number of data points, and r is the radius of an n-dimensional hypersphere. For many attractors, this function exhibits a power law dependence on r as r→0, that is c(r)∝rdG. Therefore, the correlation dimension, dG, can be determined from the slope of a plot of [logc(r)] versus [logr]. Chen and Yau [18] showed that the correlation dimension represents the lower bound to the capacity or fractal dimension dc and approaches its value asymptotically when the attracting set is distributed more uniformly in the embedding phase space. A set of points in the embedding space is said to be fractal if its dimension has a finite noninteger value. Otherwise, the attractor is referred to as a “strange attractor.” To establish the nature of the attractor, the embedding dimension is progressively increased, causing the slope of the characteristic curve to approach a steady-state value. This value is then used to determine whether the system has a fractal structure or a strange attractor structure. If the dimension of the system flow is found to be fractal (i.e. to have a noninteger value), the system is judged to be chaotic.In the current study, the attractors in the embedding space were constructed using a total of 60000 data points taken from the time series corresponding to the displacement of the system. Via a process of trial and error, the optimum delay time when constructing the time series was found to correspond to one third of a revolution of the system. The reconstructed attractors were placed in embedding spaces with dimensions ofn = 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20, respectively, yielding 10 different [logc(r)] versus [logr] plots for each attractor. The number of data points chosen for embedding purposes (i.e., 60000) reflects the need for a compromise between the computation time and the accuracy of the results. In accordance with Grassberger and Procaccia [19], the number of points used to estimate the intrinsic dimension of the attracting set in the current analysis is less than 42M, where M is the greatest integer value less than the fractal dimension of the attracting set.


## 4. Numerical Results and Discussions

The nonlinear dynamic equations presented in (2.14) to (2.15) for the HSFD rotor-bearing system with strongly nonlinear oil-film force and nonlinear rub-impact force were solved using the fourth-order Runge-Kutta method. The time step in the iterative solution procedure was assigned a value of π/300, and the termination criterion was specified as an error tolerance of less than 0.0001. The time series data corresponding to the first 800 revolutions of the rotor was deliberately excluded from the dynamic analysis to ensure that the analyzed data related to steady-state conditions. The sampled data were used to generate the dynamic trajectories, Poincaré maps, and bifurcation diagrams of the spur rotor system in order to obtain a basic understanding of its dynamic behavior. The maximum Lyapunov exponent and the fractal dimension measure were then used to identify the onset of chaotic motion. The rotating speed ratio s is one of the most significant and commonly used as a control parameter in analyzing dynamic characteristics of bearing systems. Accordingly, the dynamic behavior of the current rotor-bearing system was examined using the dimensionless rotating speed ratio s as a bifurcation control parameter.The bifurcation diagram in Figure3 shows the long-term values of the rotational angle, plotted with rotor displacement against the dimensionless speed s  without rub-impact effect. Qualitatively different behavior was observed at values of s within the range of 0<s<5. It can be seen that the dynamic motion of rotor trajectory in low speed is T-periodic motion both in X and Y directions, and it drops to a lower spatial displacement mode at the speed s=2.3. As the speed is increased, the T-period motion loses its stability at s=2.52, and a 2T-periodic motion starts to build up. The jump phenomenon is also occurred under 2T-periodic motion at s=2.7. As the speed is further increased, the 2T-periodic motion loses its stability at s=2.82, and a T-periodic motion suddenly appears. The rotor trajectory, the Poincaré map, and the displacement power spectrum in the X and Y directions at s=2.6 are given in Figure 4, from which the 0.5-subharmonic motion is shown by the double loops of the rotor trajectory, two discrete points in the Poincaré map and peaks at 0.5 in the power spectrum. The pressure distributions in the four oil chambers are shown in Figure 5. It can be seen that the variations of pressure distributions are periodic, and the period is the same with the rotor trajectory.Bifurcation diagram ofX(nT) (a) and Y(nT) (b) versus rotor speed s (without rub-impact effect).
(a)
(b)Subharmonic motion ats=2.6 (case 1); (a) Rotor trajectory; (b) Poincaré map; (c) and (d) Displacement power in X and Y directions (without rub-impact effect).
(a)
(b)
(c)
(d)Pressure distribution in the static pressure chamber ats=2.6 (without rub-impact effect).
(a)
(b)
(c)
(d)Figures6(a) and 6(b) show the bifurcation diagrams for the rotor displacement against the dimensionless rotating speed ratio with rub-impact effect. Compared with bifurcation results without rub-impact effect, bifurcation results with rub-impact effect show that dynamic trajectories perform strongly nonperiodic at low rotating speeds, but it would escape nonperiodic motions to periodic motions. The bifurcation diagrams show that the geometric centers of rotor in the horizontal and vertical directions perform nonperiodic motion or the so-called chaotic motion at low values of the rotating speed ratio, that is, s<0.61. Figures 7, 8, and 9 represent phase diagrams, power spectra, Poincaré maps, Lyapunov exponents, and the fractal dimensions of pinion center with s=0.32,0.36,ands=0.42, respectively. The simulation results show that phase diagrams show disordered dynamic behaviors with s=0.32,0.36 and s=0.42; power spectra reveal numerous excitation frequencies; the return points in the Poincaré maps form some geometrically fractal structures, but the maximum Lyapunov exponent is positive with s=0.36 and maximum Lyapunov exponent is negative with s=0.32 and s=0.42. Thus, the results show that the dynamic trajectory performs chaotic motion with s=0.36, but they present no chaotic motions with s=0.32ands=0.42. Figures 10 and 11 are phase diagrams and Poincaré maps for the route of subharmonic motion into chaos, out of chaos to periodic response at different rotating speed ratios of s (with rub-impact effect). Unlike the usual ways into chaos (1T⇒2T⇒4T⇒8T⇒16T⇒32T⋯⇒chaos or periodic⇒quasi-periodic⇒chaotic), it suddenly gets in chaos from the periodic motion without any transition or suddenly escape from irregular motions into periodic motions in accordance with phase diagrams and Poincaré maps.Bifurcation diagram of HSFD rotor-bearing system using dimensionless rotating speed coefficient,s, as bifurcation parameter (with rub-impact effect).
(a)
(b)Simulation results obtained for rotor-bearing system withs=0.32 (with rub-impact effect).
(a)
Phase diagram
(b)
Power spectrum
(c)
Lyapunov exponent
(d)
Poincaré map
(e)
Fractal dimensionSimulation results obtained for rotor-bearing system withs = 0.36 (with rub-impact effect).
(a)
Phase diagram
(b)
Power spectrum
(c)
Lyapunov exponent
(d)
Poincaré map
(e)
Fractal dimensionSimulation results obtained for rotor-bearing system withs=0.42 (with rub-impact effect).
(a)
Phase diagram
(b)
Power spectrum
(c)
Lyapunov exponent
(d)
Poincaré map
(e)
Fractal dimensionPhase diagrams for the route of subharmonic motion into chaos, out of chaos to periodic response at different rotating speed ratios ofs (with rub-impact effect).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(q)
(r)
(s)
(t)
(u)Poincaré maps for the route of subharmonic motion into chaos, out of chaos to periodic response at different rotating speed ratios ofs (with rub-impact effect).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(q)
(r)
(s)
(t)
(u)


## 5. Conclusions

A hybrid squeeze-film damper mounted rotor-bearing system with nonlinear oil-film force and nonlinear rub-impact force has been presented and studied by a numerical analysis of the nonlinear dynamic response in this study. The dynamics of the system have been analyzed by reference to its dynamic trajectories, power spectra, Poincaré maps, bifurcation diagrams, maximum Lyapunov exponents, and fractal dimensions. The bifurcation results can be observed that HSFD may be used to improve dynamic irregularity. The system with rub-impact force effect may be a strongly nonlinear effect, and the bifurcation results show that HSFD mounted rotor-bearing system with rub-impact force effect present nonperiodic motions at low rotating speeds and perform periodic motions at high rotating speeds. The results will enable suitable values of the rotating speed ratio to be specified such that chaotic behavior can be avoided, thus reducing the amplitude of the vibration within the system and extending the system life.



---

*Source: 279827-2012-07-24.xml*