Nonlinear vibration phenomenon of an aircraft rub-impact rotor system due to hovering flight

Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

Contents lists available at SciVerse ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Nonlinear vibration phenomenon of an aircraft rub-impact rotor system due to hovering flight Lei Hou ⇑, Yushu Chen, Qingjie Cao School of Astronautics, Harbin Institute of Technology, Harbin 150001, PR China

a r t i c l e

i n f o

Article history: Received 27 December 2012 Received in revised form 13 June 2013 Accepted 15 June 2013 Available online 27 June 2013 Keywords: Nonlinear vibration Rub-impact Hovering flight Maneuver load

a b s t r a c t This paper focuses on the nonlinear vibration phenomenon caused by aircraft hovering flight in a rub-impact rotor system supported by two general supports with cubic stiffness. The effect of aircraft hovering flight on the rotor system is considered as a maneuver load to formulate the equations of motion, which might result in periodic response instability to the rotor system even the eccentricity is small. The dynamic responses of the system under maneuver load are presented by bifurcation diagrams and the corresponding Lyapunov exponent spectrums. Numerical analyses are carried out to detect the periodic, sub-harmonic and quasi-periodic motions of the system, which are presented by orbit diagrams, phase trajectories, Poincare maps and amplitude power spectrums. The results obtained in this paper will contribute an understanding of the nonlinear dynamic behaviors of aircraft rotor systems in maneuvering flight. Published by Elsevier B.V.

1. Introduction The rub-impact of rotor-to-stator becomes one of the most common faults in aircraft rotor systems, which often occurs due to the trend of reducing the radial clearance between the rotor and the stator [1]. The nonlinear dynamic behaviors of rub-impact rotor systems have been widely investigated to detect the chaotic motions [2–5]. The thermal effects of rotor-tostator rub on the rotor vibration response were discussed [6,7]. The lateral-torsional coupling effects on the nonlinear dynamic behavior of rubbing rotor systems were investigated [8,9]. The effect of Coulomb damping caused by blade tips on the rubbing rotor system was analyzed [10]. These studies showed that mass imbalance, defective bearing, rotor misalignment, shaft bowing and compressor blade failure are the main reasons of the rub-impact when the eccentricity of rotor system increases significantly. Different methods are proposed to detect the existence of the periodic solutions and the stabilities of rub-impact systems [11–13]. Groll [14] described a numerical algorithm based on the harmonic balance method to calculate the periodic response of a non-linear system under periodic excitation and the stability of the found periodic solutions to determine the bifurcations following a solution branch over varying system parameters via arc-length continuation. Lu [15] presented a method consisting of analytical and numerical techniques to solve the existence problem of rub-impact periodic motions to find the grazing circle motions and the single-impact periodic motions. Han [16] investigated the stability of periodic motions of a dual-disk rotor system analytically, which is agreed with the numerical results obtained by the direct numerical integration. Diken [17] used a perturbation technique to obtain the approximate linear equation for the rotor system of a thin disk located on a flexible shaft, which shows there exist two sub-harmonic transient vibrations.

⇑ Corresponding author. Tel.: +86 18003667477. E-mail address: [email protected] (L. Hou). 1007-5704/$ - see front matter Published by Elsevier B.V. http://dx.doi.org/10.1016/j.cnsns.2013.06.023

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

287

Much attention has been paid on the rub-impact rotor systems with other coupled faults over the past five years, including oil film instability, crack faults, initial permanent bow and so on. Rub-impact rotor systems supported by oil film journal bearings were studied by Chang-Jian [18–21] and Wang [22] to demonstrate the periodic, quasi-periodic, sub-harmonic, chaotic and multi-periodic motions. Ma [23] analyzed the rub-impact fault coupled with crack and that coupled with oil-film instability respectively which shows that the increase of low frequency amplitude can be viewed as the most distinguishable characteristic for the seriousness of the bearing rub-impact. Experimental investigation was proposed by Patel [24,25] for vibration response of multi-fault rotor system to extract the distinctive directional features of the rotor faults. A general model of a rub-impact rotor-bearing system with initial permanent bow was set up and studied by Shen [26]. Cao [27] investigated the nonlinear dynamic characteristics of a rub-impact rotor system with fractional order damping which shows that the system exhibits rich nonlinear dynamic behaviors, including period doubling bifurcation, sudden transition and quasiperiodic motion. Based on a full-degree-of-freedom Jeffcott rotor, Yuan [28] studied an axial rub-impact problem and disclosed an internal bending-torsion coupling phenomenon which concludes that the dynamic behaviors of axial rub-impact are quite distinct from those of radial ones, and that the dynamic behaviors of axial vibration are essential to the diagnostics of axial rub-impact. The motivation of this paper is to detect the nonlinear vibration phenomenon of an aircraft rotor system caused by maneuvering flight, which has been rarely studied in the previous literatures. The equations of motion proposed here enable us to investigate the nonlinear dynamics of the system considering the effects of the maneuver load induced by hovering flight, which is the maneuvering in a horizontal plane. Numerical technique is employed to detect the complex nonlinear dynamic behaviors of the rub-impact occurring in the system with a small eccentricity, including sub-harmonic motion and quasi-periodic motion due to the maneuver load. In this paper, a hovering flight model is presented and simplified as a maneuver load in the differential equations of the rotor system by using Lagrange equation. The 4th order Runge–Kutta method is employed to calculate the nonlinear response of the system. Firstly, the bifurcation diagrams for different level of maneuver load are plotted to indicate its influence to the dynamics of the system. Then, in the condition of E ¼ 0:25, the Lyapunov exponent spectrums corresponding to the bifurcation diagrams for G2 ¼ 0:9 and X ¼ 2:2 are plotted respectively to determine whether the complex dynamic responses are chaotic or quasi periodic. Finally, the nonlinear dynamic behaviors under different maneuver loads for E ¼ 0:25 and X ¼ 2:2 are demonstrated by vibration responses of the rotor center, phase plane portraits, Poincare maps, and amplitude power spectrums.

2. Modeling the problem 2.1. Hovering flight model Consider the Jeffcott rotor model as shown in Fig. 1, where o is the left endpoint of the shaft, which is assumed to be at the gravity center of the aircraft, o0 is the geometric center of the disk and om is the mass center of the disk with the eccentricity e from o0 . The supports are assumed to be with a linear stiffness k and a nonlinear cubic stiffness a. The motion of the disk is modeled by four degrees of freedom: the vertical displacement y, the horizontal displacement z, and the corresponding rotation angles hy and hz . The maneuvering flight model here in the study is assumed to be under the hovering flight condition as shown in Fig. 2, where oxyz is the rotor coordinate system fixed on the fuselage of the plane, corresponding to Fig. 1. v is the flight speed of the plane and xB;y is the angular velocity of hovering flight. 2.2. Rub-impact force model Fig. 3 is a schematic diagram of the rub-impact force model with an initial clearance r 0 between the rotor and the stator, the rotor speed is x. The radial impact force PN and the tangential rub force PT are expressed as [29]

Fig. 1. Schematic diagram of the rub-impact rotor system.

288

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

ωB , y

v

Fig. 2. Maneuvering flight model.

y

PN

Rotor r0

z

o0

ϕ

PT

ω o′

Stator Fig. 3. Rub-impact force model.



PN

 ¼

PT

8 <0

   r : kc 1  rr0



r 6 r0 r > r0

lc r

ð1Þ

;

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where kc is the radial stiffness of the stator, lc is the coefficient of Coulomb friction, r ¼ y2 þ z2 is the radial displacement of the rotor. The rub-impact forces can be written in y  z coordinates as the following by taking the coordinate transformation of 3



Py

 ¼

Pz 

Py Pz



8 <0

: kc 1  

¼



r

r > r0

lc y  z

 sin /  cos / cos /

r 6 r0

   y þ lc z r0 

 sin /

PN PT

:

ð2Þ

 ;

ð3Þ

where sin / ¼ yr ; cos / ¼  zr. 2.3. Equations of motion The equations of motion can be obtained by using the Lagrange equation in the following

  d @ ðLÞ @ ðLÞ @ ðDÞ  þ ¼ Qj dt @ q_ j @qj @qj

ðj ¼ 1; 2; 3; 4Þ;

ð4Þ

T T where L ¼ T t þ T r  V; q ¼ ½ y z hy hz  ; Q ¼ ½ P y Pz 0 0  . T t and T r are the translational and the rotational kinetic energy, V and D are the elastic potential energy and the dissipated energy respectively, written as [30,31]:

2 1  2 2 m v þ xB;y ðz þ e sin xtÞ þ ðy_  xe sin xtÞ þ ðz_ þ xe cos xt Þ ; 2  

2 1 1 T r ¼ J p x2 þ J d h_ y þ xB;y þ h_ 2z  J p xh_ z hy ; 2 2 Tt ¼

V¼

 2  2 1  2 2 1 2 2 2 k ðy þ l1 hz Þ þ ðy  l2 hz Þ þ z  l1 hy þ z þ l2 hy þ a ðy þ l1 hz Þ þ z  l1 hy 2 4  2 2 1 2 þ a ðy  l2 hz Þ þ z þ l2 hy ; 4

ð5Þ

ð6Þ

ð7Þ

289

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

D¼



2

2

2

2  1 ; c y_ þ l1 h_ z þ y_  l2 h_ z þ z_  l1 h_ y þ z_ þ l2 h_ y 2

ð8Þ

in which, v and xB;y represent the motion of the aircraft. The equations of motion are obtained as follows



 

€ þ kð2y þ l1 hz  l2 hz Þ þ c 2y_ þ l1 h_ z  l2 h_ z þ a ðy þ l1 hz Þ2 þ z  l1 hy 2 ðy þ l1 hz Þ my  2

2 þ a ðy  l2 hz Þ þ z þ l2 hy ðy  l2 hz Þ ¼ Py þ mx2 e cos xt;

ð9aÞ



 2   2 z  l1 h y m€z þ k 2z  l1 hy þ l2 hy þ c 2z_  l1 h_ y þ l2 h_ y þ a ðy þ l1 hz Þ þ z  l1 hy  2   2 þ a ðy  l2 hz Þ þ z þ l2 hy z þ l2 h y 



¼ Pz þ mx2 e sin xt þ mxB;y v ;

ð9bÞ

       2   2 2 J d €hy þ k  z  l1 hy l1 þ z þ l2 hy l2  a ðy þ l1 hz Þ þ z  l1 hy z  l1 hy l1 þ a ðy  l2 hz Þ





 2   þ z þ l2 h y z þ l2 hy l2 þ c  z_  l1 h_ y l1 þ z_ þ l2 h_ y l2 þ J p xh_ z ¼ 0;

ð9cÞ

 2

2 J d €hz þ kððy þ l1 hz Þl1  ðy  l2 hz Þl2 Þ þ a ðy þ l1 hz Þ þ z  l1 hy ðy þ l1 hz Þl1





 2

2  a ðy  l2 hz Þ þ z þ l2 hy ðy  l2 hz Þl2 þ c y_ þ l1 h_ z l1  y_  l2 h_ z l2  J p xh_ y ¼ J p xxB;y ;

ð9dÞ

where mxB;y v in (9b) and J p xxB;y in (9d) that are induced by the motion of the aircraft represent the maneuver loads in the rotor system. 2.4. Dimensionless equations The dimensionless equations of (4) can be obtained as the following by letting e0 ¼ r0 ; q1 ¼ ey0 ; q2 ¼ ez0 ; q3 ¼ ðl2 þl1 Þhy e0

; q4 ¼ ðl2 þle01 Þhz ; s ¼ xt.

      q001 þ q1 þ a1 q4 þ a2 q01 þ a3 q04 þ a4 q31 þ q1 q22 þ a5 2q1 q2 q3 þ 3q21 q4 þ q22 q4 þ a6 q1 q23 þ 3q1 q24 þ 2q2 q3 q4   þ a7 q23 q4 þ q34 ¼ P1 þ EX2 cos Xs;

ð10aÞ

      q002 þ q2 þ a1 q3 þ a2 q02 þ a3 q03 þ a4 q21 q2 þ q32 þ a5 q21 q3 þ 2q1 q2 q4 þ 3q22 q3 þ a6 2q1 q3 q4 þ 3q2 q23 þ q2 q24  3  þ a7 q3 þ q3 q24 ¼ P2 þ EX2 sin Xs þ G2 ;

ð10bÞ

    q003 þ x22 q3 þ b1 q2 þ b2 q03 þ b3 Xq04 þ b4 q02 þ b6 q21 q3 þ 2q1 q2 q4 þ 3q22 q3 þ b5 q21 q2 þ q32     þ b7 2q1 q3 q4 þ 3q2 q23 þ q2 q24 þ b8 q33 þ q3 q24 ¼ 0;

ð10cÞ

    q004 þ x22 q4 þ b1 q1 þ b2 q04  b3 Xq03 þ b4 q01 þ b6 2q1 q2 q3 þ 3q21 q4 þ q22 q4 þ b5 q31 þ q1 q22     þ b7 q1 q23 þ 3q1 q24 þ 2q2 q3 q4 þ b8 q23 q4 þ q34 ¼ G4 X;

ð10dÞ Jp x1 ðl1 þl2 Þ G2 Jd v

which is a coupled nonlinear dynamic system with the parameters shown in Appendix A, where G2 and G4 ¼ represent the maneuver loads. The approximate solutions of the equations can be obtained by using numerical methods. The calculation parameters of the system are as follows [12,32]

m ¼ 20 kg; J d ¼ 0:072 kg m2 ;

J p ¼ 0:144 kg m2 ;

l1 ¼ 0:25 m;

l2 ¼ 0:5 m;

k ¼ 1:57  106 N m1 ;

12

a ¼ 3:8  10 N m3 ; c ¼ 1800 N s=m; v ¼ 100 m s1 ; r0 ¼ 0:4 mm; kc ¼ 3  106 N m1 ; lc ¼ 0:01 pa s: ð11Þ 3. Numerical results and discussions 3.1. Effect of maneuver load To get the effects of the maneuver load on the dynamic characteristics of the rotor system, the 4th order Runge–Kutta method is employed in the following computations. The bifurcation diagrams are plotted in Fig. 4, Fig. 5 and Fig. 6, respectively, for q1 and q2 versus X for different values of E to show the responses of the system under different maneuver loads,

290

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

Fig. 4. Bifurcation diagrams for different level of E in the condition of G2 ¼ 0. a q1 versus X. b q2 versus X.

Fig. 5. Bifurcation diagrams for different level of E in the condition of G2 ¼ 0:5. a q1 versus X. b q2 versus X.

Fig. 6. Bifurcation diagrams for different level of E in the condition of G2 ¼ 1. a q1 versus X. b q2 versus X.

Fig. 7. Bifurcation diagrams for E ¼ 0:25. a q1 versus X for different level of G2 . b q1 versus G2 for different level of X.

where q1 represents the vertical vibration response, q2 refers to the horizontal vibration response, E is the normalized mass eccentricity which is also the ratio of eccentricity to the clearance between rotor and stator, and X is the normalized frequency of the imbalance exciting force. Fig. 4 shows the bifurcation diagrams for G2 ¼ 0, in which the complex nonlinear dynamic behaviors of the system under the perturbation of large mass eccentricity and high frequency are observed, including period doubling bifurcation, grazing bifurcation, chaos and quasi periodic motion which are similar to the results given in the current literatures (see [2–6]).

291

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

(b) Maximum Lyapunov exponent

Maximum Lyapunov exponent

(a) 0.02 0 -0.02 -0.04 -0.06

2

2.2

2.4

2.6

2.8

0.02 0 -0.02 -0.04 -0.06

3

0

0.5

1

1.5

G2

Ω

(c)

(d)

Fig. 8. Bifurcation diagrams (a and b) and the corresponding Lyapunov exponent spectrums (c and d) for E ¼ 0:25. a q1 versus X and c maximum Lyapunov exponent versus X for G2 ¼ 0:9. b q2 versus G2 and d maximum Lyapunov exponent versus G2 for X ¼ 2:2.

1.5

rubbing boundary

0.5

0.5

0.5

0.5

0

0

0

-0.5

q1

1

q1

1

-0.5

0

-1

-1

-1

-1

-1.5

-1.5

-1.5

-1

0 q2

1

-1

1.5

0 q2

1

-1

(b) 1.5

rubbing boundary

0 q2

1

1.5

rubbing boundary

1.5

rubbing boundary

0.5

0.5

0.5

0

0

0

q1

0.5

q1

1

q1

1

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-1

(e) 1.5

0 q2

1

-1

(f) 1.5

rubbing boundary

0 q2

1

1.5

rubbing boundary

1.5

rubbing boundary

0.5

0.5

0.5

0

0

0

q1

0.5

q1

1

q1

1

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

(i)

-1

0 q2

(j)

1

rubbing boundary

0

-1 1

1

-0.5

-1.5

0 q2

0 q2

(h)

1

-1

-1

(g)

1

-0.5

rubbing boundary

0

-1 1

1

-0.5

-1.5

0 q2

0 q2

(d)

1

-1

-1

(c)

1

-0.5

rubbing boundary

-0.5

-1.5

(a)

q1

1.5

rubbing boundary

1

-0.5

q1

1.5

rubbing boundary

1

q1

q1

1.5

-1

0 q2

(k)

1

-1

0 q2

1

(l)

Fig. 9. Orbit diagrams of rotor center in the condition of E ¼ 0:25 and X ¼ 2:2. a for G2 ¼ 0:259. b for G2 ¼ 0:26. c for G2 ¼ 0:547. d for G2 ¼ 0:548. e for G2 ¼ 0:549. f for G2 ¼ 0:55. g for G2 ¼ 0:582. h for G2 ¼ 0:7. i for G2 ¼ 0:79. j for G2 ¼ 0:8. k for G2 ¼ 1. l for G2 ¼ 1:5.

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297 1

1

0.5

0.5

0

0

0

-0.5

-0.5

-1

-1

-0.5

1985

1990

1995

2000

-1.5 -1

-0.5

0 q1

τ

(a)

0.5

-1.5 -1

1

-0.5

(b)

1

0 q1

0.5

2

2

1

1 q'2

q'2

q2

0 -1

-0.5 1980

1985

1990

1995

0

(e)

1 2 3 4 Normalized frequency

0.5

0

-2 -1

1

5

3

-0.5

0 q2

q2

τ

0

(d)

-1

-2 -0.5

2000

1

(c)

0.5

0

2

0

1

Power spectrum of q2

-1 1980

3

q'1

q'1

q1

0.5

1

Power spectrum of q1

292

(f)

0.5

2

1

0

1

0

1 2 3 4 Normalized frequency

(g)

5

(h)

Fig. 10. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 0:5. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

2

0.5

1

1.4

3 Power spectrum of q1

1

1.2

0

0.8 -1

1500

-2 -1

2000

0.6 -0.5

0 q1

τ

(a)

0.5

0.4 -0.5

1

τ

(e)

2000

0.3 0.25

0

0.2

-0.5

0.15 0

0.5 q2

(f)

1

0.1 -0.5

1 2 3 4 Normalized frequency

5

2

q'2

1 0.5

-1 -0.5

0

(d)

0.35

q'2

2

1500

0

1

0

-0.5 1000

1

(c)

1.5

0.5

0.5

2

q1

(b)

1

0

Power spectrum of q2

-0.5 -1 1000

q'1

q'1

1

1 0

0 q2

(g)

0.5

1.5 1 0.5 0

0

1 2 3 4 Normalized frequency

5

(h)

Fig. 11. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 0:548. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

In the same way, Fig. 5 for G2 ¼ 0:5 and Fig. 6 for G2 ¼ 1 are, respectively, illustrite the effects of different maneuver loads on the vibration response of the rotor system, which show that a series of complex dynamic behaviors exist, and may occur even in the case of small mass eccentricity when the maneuver load is large enough, which occurs in a section of lower frequency. This phenomenon is singnificantly different from the nonlinear phenomenon that occurs only in the case of very high frequency without maneuver load. Fig. 7 shows the variation of bifurcation diagram due to different maneuver load and different frequency, respectively, for E ¼ 0:25. From Fig. 7(a), it can be seen that the frequency region of periodic-2T motion and quasi periodic motion is broaden when G2 is increasing. From Fig. 7(b), it can be seen that the start point and the end point of G2 of periodic-2T motion and quasi periodic motion are getting larger when the frequency is increasing from X ¼ 2 to X ¼ 3.

293

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297 2

0.5

1

1.2

2 Power spectrum of q1

1

1

0

0.6 -1

-0.5 -1 1000

1500

-2 -1

2000

0.4 -0.5

0 q1

τ

(a)

0.5

0.2 -0.5

1

0

(b)

1

q2

q'2

0.5

2

0.4

1

0.3

0 -1

-1 1000

1500

0

-2 -1

2000

0

1 2 3 4 Normalized frequency

1

0.2

0 -1

2

5

2

-0.5

0 q2

q2

(e)

0

(d)

0.1

τ

0.5

1

0 -0.5

1

(c)

q'2

1.5

0.5

1.5

q1

Power spectrum of q2

q1

0

q'1

q'1

0.8

(f)

0.5

1.5 1 0.5 0

1

0

1 2 3 4 Normalized frequency

(g)

5

(h)

Fig. 12. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 0:55. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

2

1.2 1

1

0.5

3 Power spectrum of q1

1

0

0.6 -1

-0.5

1500

-2 -1

2000

0.4 -0.5

0 q1

τ

(a)

0.5

0

q'2

q'2

q2

0

0.2 -1

2000

τ

(e)

-2 -1

0.1 0

1 q2

(f)

2

0 -1

1 2 3 4 Normalized frequency

5

(d)

0.4

0

1500

0

2

0.3

-0.5

0

1

0.5

1

0.5

1

(c)

2

1

0.5

2

q1

(b)

1.5

-1 1000

0.2 -0.5

1

Power spectrum of q2

-1 1000

q'1

q'1

q1

0.8 0

-0.5

0 q2

(g)

0.5

1

1.5 1 0.5 0

0

1 2 3 4 Normalized frequency

5

(h)

Fig. 13. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 0:551. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

3.2. Lyapunov exponent analysis Fig. 8(a) and Fig. 8(b) show the bifurcation diagrams for G2 ¼ 0:9 and X ¼ 2:2 in the condition of E ¼ 0:25, where the normalized exciting frequency and the maneuver load are, respectively, taken as the bifurcation parameter, Fig. 8(c) and Fig. 8(d) show the corresponding Lyapunov exponent spectrums. According to Fig. 8(a) and Fig. 8(c), in the case of E ¼ 0:25 and G2 ¼ 0:9, the dynamic behavior of the rotor system is periodic from X ¼ 2 to X ¼ 2:13. The period doubling bifurcation is observed at X ¼ 2:131 and maintained until X ¼ 2:196. The response of the rotor system jumps back to periodic motion at X ¼ 2:203 and keeps stable until X ¼ 2:26. In the region of X ¼ 2:261 to X ¼ 2:497, the motion is quasi periodic. Then the response jumps to period 2-T motion and last until X ¼ 2:999. Finally, the system returns back to periodic motion at X ¼ 3.

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297 1

2

0.5

1

1.2

2 Power spectrum of q1

294

1

0

0

0.6 -1

-0.5

1500

-2 -1

2000

0.4 -0.5

0 q1

τ

(a)

0.5

0.2 -0.5

1

0

(b)

1

q2

q'2

0.5

-1 1000

0.6

1

0.4

0 -1

1500

0

1

-0.2 -1

2

(e)

(f)

1 2 3 4 Normalized frequency

5

(d)

0.2

-0.5

0 q2

q2

τ

0

2

0

-2 -1

2000

0.5 0

1

0 -0.5

1

(c)

2

q'2

1.5

0.5

1.5

q1

Power spectrum of q2

-1 1000

q'1

q'1

q1

0.8

0.5

1.5 1 0.5 0

1

0

(g)

1 2 3 4 Normalized frequency

5

(h)

1

0

-1 1000

0.8

0

0.6

-1

-0.5

1500

-2 -1

2000

-0.5

τ

0 q1

0.5

0.4 -0.5

1

(b) 0.4

1

1

0.3

-0.5 1000

q'2

2

0 -1

0

1500

2000

τ

(e)

-2 -0.5

0 q1

0.5 q2

(f)

1 0.5 0

0.5

0

1

1.5

5

2

0.2

0 -0.5

1 2 3 4 Normalized frequency

(d)

0.1

0

1.5

(c)

1.5

q'2

q2

(a)

0.5

2 Power spectrum of q1

0.5

1

Power spectrum of q2

2

q'1

1

q'1

q1

Fig. 14. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 0:582. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

0

0.5 q2

(g)

1

1.5 1 0.5 0

0

1 2 3 4 Normalized frequency

5

(h)

Fig. 15. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 0:7. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

From Fig. 8(b) and Fig. 8(d), in the case of E ¼ 0:25 and X ¼ 2:2, the response of the system keeps periodic from G2 ¼ 0 to G2 ¼ 0:259, and jumps to period 2-T motion at G2 ¼ 0:26. When G2 ¼ 0:548, the motion becomes quasi periodic. After the quasi periodic motion in the region of G2 ¼ 0:549 to G2 ¼ 0:795, the response of the system returns back to periodic motion and keeps until G2 ¼ 0:85. Then the response jumps to period 2-T motion at G2 ¼ 0:853, and returns to periodic motion at G2 ¼ 1:38 after a reversed bifurcation. 3.3. Nonlinear dynamics analysis Fig. 9 shows the rotor center orbit evolving with the maneuver load in the condition of E ¼ 0:25 and X ¼ 2:2. At G2 ¼ 0:26, the response of the rotor center turns from period 1-T motion into period 2-T motion in a sudden and the rub-impact hap-

295

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

q1

q'1

0.2

0.8

0.5

0.75

0

2 Power spectrum of q1

0.4

1

q'1

0.6

0.7

0 -0.5

-0.2 1500

-1 -0.5

2000

0 q1

τ

(a)

0.5

0

0.05

(b)

1.5

0.1 q1

0.15

1

0.24

0.5

0.22 q'2

q'2

q2

0 -0.5 -1

0 1000

1500

2000

0.5

1

0.2

0.16 0.1

1.5

(e)

1 2 3 4 Normalized frequency

5

2

0.2

0.3 q2

q2

τ

0

(d)

0.18

0

0.5

(c)

1

0.5

1

0

0.2

Power spectrum of q2

-0.4 1000

0.65

1.5

(f)

0.4

1.5 1 0.5 0

0.5

0

(g)

1 2 3 4 Normalized frequency

5

(h)

Fig. 16. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 0:79. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

0.5

0.5

0

0

-0.5

1990

1995

-0.5

-0.5

-1

-1

-1.5 -1

2000

-0.5

0 q1

τ

(a)

0.5

-1.5 -1

1

(b)

1.5

-0.5

0 q1

0.5

1

0.5

0.5 q'2

q'2

q2 0 1980

0 -0.5 -1

1985

1990

τ

(e)

1995

2000

0.5

1 q2

(f)

0.5

0

1.5

5

(d)

0

-1

1 2 3 4 Normalized frequency

2

-0.5

0

1

(c)

1

1

0.5

1.5

0

1

Power spectrum of q2

1985

2

q'1

0

-1 1980

1

q'1

q1

0.5

1

Power spectrum of q1

1

0

0.5

1 q2

(g)

1.5

1.5 1 0.5 0

0

1 2 3 4 Normalized frequency

5

(h)

Fig. 17. Vibration responses (a and e), phase trajectories (b and f), Poincare maps (c and g), and amplitude power spectrums (d and h) of the rotor center in the condition of E ¼ 0:25; X ¼ 2:2 and G2 ¼ 1. a q1 versus s. b andc q01 versus q1 . d amplitude power spectrum of q1 . e q2 versus s. f and g q02 versus q2 . h amplitude power spectrum of q2 .

pens at the same time. At G2 ¼ 0:548, the period 2-T motion becomes quasi periodic. The quasi periodic motion changes dramatically from G2 ¼ 0:548 to G2 ¼ 0:79. At G2 ¼ 0:8, the response returns to period 1-T motion, however, a slight partial rubimpact exists in this case. At G2 ¼ 1, the response performs a period 2-T motion. At G2 ¼ 1:5, the response returns to period 1-T motion with a heavy rub-impact. The phase trajectories, the Poincare maps, and the amplitude power spectrums of the complex nonlinear motions are shown in Figs. 10–17. In the amplitude power spectrums of the period 2-T motions, a frequency of half exciting force frequency is observed so that the motions are sub-harmonic resonance responses. In the Poincare maps, it turns from one point to two single points when the response of the rotor center turns from period 1-T motion into period 2-T motion at G2 ¼ 0:26.

296

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

And the two single points become two circles when the response turns from period 2-T motion to quasi periodic motion at G2 ¼ 0:548. At G2 ¼ 0:551, the two circles combine into one single closed curve. Then the closed curve becomes smaller and smaller when the maneuver load increases, and turns into a small circle at G2 ¼ 0:79. At G2 ¼ 0:8, the circle shrinks into a point when the response returns to period 1-T motion. At G2 ¼ 1, there are two single points in the Poincare map. The phase trajectories, the Poincare maps, and the amplitude power spectrums coincide with each other. 4. Conclusions In this paper, we have investigated the nonlinear vibration phenomenon of an aircraft rub-impact rotor system caused by hovering flight which is modeled as a maneuver load. The 4th order Runge–Kutta method has been employed to obtain the bifurcation diagrams, Lyapunov exponent spectrums, orbit diagrams, phase trajectories, Poincare maps and amplitude power spectrums that show the nonlinear behaviors of the system. The results presented here show that the sub-harmonic resonance of the system can occur due to the maneuver load, which results in the rub-impact even in a rotor system with a small eccentricity. The quasi periodic motions have also been obtained due to the small damping considered in the simulating model, which is similar to the results of researches studying the effect of damping on the regular orbits and chaotic orbits of nonlinear systems (see [33–35]). The results will contribute an understanding of the nonlinear dynamic behaviors of the maneuvering rotor systems. The future studies on this subject are being carried out in two aspects: one is the analytically investigaton of the parameter region for the sub-harmonic resonance of the system and the other is the effect of the maneuvering flight in the three dimensional space. Acknowledgement The authors would like to acknowledge the financial supports from the Natural Science Foundation of China (Grant Nos. 10632040 and 11072065). Appendix A

kðl2  l1 Þ a1 ¼ ; mx21 ðl1 þ l2 Þ

3 3 e20 a l2  l1

a7 ¼

mx21 ðl1 þ l2 Þ 2

3

;

e E¼ ; e0

2

kðl þ l Þ x22 ¼ 1 2 2 ; J d x1

b6 ¼

2c a2 ¼ ; mx1



2 2 e20 a l1 þ l2 J d x21

b1 ¼

;

cðl2  l1 Þ a3 ¼ ; mx1 ðl1 þ l2 Þ

x kc X¼ ; g¼ ; x1 mx21



2 2 k l2  l1

b7 ¼

J d x21

;

b2 ¼



3 3 e20 a l2  l1 J d x21 ðl1 þ l2 Þ

;



2 2 c l1 þ l2 J d x1

b8 ¼

;

2e2 a a4 ¼ 0 2 ; mx1 

P1



P2

J d x21 ðl1 þ l2 Þ

2

;

a6 ¼



2 2 e20 a l1 þ l2

b4 ¼

G4 ¼



2 2 c l2  l1 J d x1

;

b5 ¼

2

mx21 ðl1 þ l2 Þ

8 R61 <0     q1 þ lq2 ¼ ; 1 R>1 :g 1  R lq1  q2

Jp b3 ¼ ; Jd



4 4 e20 a l1 þ l2

e2 aðl2  l1 Þ a5 ¼ 0 2 ; mx1 ðl1 þ l2 Þ

G2 ¼



2 2 e20 a l2  l1 J d x21

;

xBy v ; x21 e0

;

J p xBy ðl1 þ l2 Þ : J d x1 e0

References [1] Chen YS, Zhang HB. Review and prospect on the research of dynamics of the aero-engine system. Acta Aeronautica et Astronautica Sinica 2011;32:1–22. [2] Muszynska Agnes, Goldman Paul. Chaotic responses of unbalanced rotor/bearing/stator systems with looseness or rubs. Chaos Solitons Fract 1995;5:1683–704. [3] Chu F, Zhang Z. Bifurcation and chaos in a rub-impact Jeffcott rotor system. J Sound Vib 1998;210:1–18. [4] Sun ZC, Xu JX, Zhou T. Analysis on complicated characteristics of a high-speed rotor system with rub-impact. Mech Mach Theory 2002;37:659–72. [5] Qin WY, Chen GR, Meng G. Nonlinear responses of a rub-impact overhung rotor. Chaos Solitons Fract 2004;19:1161–72. [6] Goldman Paul, Muszynska Agnes. Rotor-to-stator, rub-related, thermal/mechanical effects in rotating machinery. Chaos Solitons Fract 1995;5:1579–601. [7] Bachschmid N, Pennacchi P, Vania A. Thermally induced vibrations due to rub in real rotors. J Sound Vib 2007;299:683–719. [8] Edwards S, Lees AW, Lees AW. M.I. friswell. The influence of torsion on rotor/stator contact in rotating machinery. J Sound Vib 1999;225:767–78. [9] Khanlo HM, Ghayour M, Ziaei-Rad S. The effects of lateral-torsional coupling on the nonlinear dynamic behavior of a rotating continuous flexible shaftdisk system with rub-impact. Commun Nonlinear Sci Numer Simul 2012:1–15. [10] Sinha SK. Dynamic characteristics of a flexible bladed-rotor with Coulomb damping due to tip-rub. J Sound Vib 2004;273:875–919. [11] Gra¯pis O, Tamuzˇs V, Ohlson NG, Andersons J. Overcritical high-speed rotor systems, full annular rub and accident. J Sound Vib 2006;290:910–27. [12] Zhang HB, Chen YS. Bifurcation analysis on full annular rub of a nonlinear rotor system. Sci China E 2011;54:1977–85. [13] Zhang GF, Xu WN, Xu B, Zhang W. Analytical study of nonlinear synchronous full annular rub motion of flexible rotor–stator system and its dynamic stability. Nonlinear Dyn 2009;57:579–92. [14] von Groll G, Ewins DJ. The harmonic balance method with arc-length continuation in rotor/stator contact problems. J Sound Vib 2001;241:223–33. [15] Lu QS, Li QH, Twizell EH. The existence of periodic motions in rub-impact rotor systems. J Sound Vib 2003;264:1127–37.

L. Hou et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 286–297

297

[16] Han QK, Zhang ZW, Liu CL, Wen BC. Periodic motion stability of a dual-disk rotor system with rub-impact at fixed limiter. Vibro-Impact Dyn Ocean Syst 2009;44:105–19. [17] Diken H. Non-linear vibration analysis and sub-harmonic whirl frequencies of the Jeffcott rotor model. J Sound Vib 2001;243:117–25. [18] Chang-Jian CW, Chen CK. Chaos and bifurcation of a flexible rub-impact rotor supported by oil film bearings with nonlinear suspension. Mech Mach Theory 2007;42:312–33. [19] Chang-Jian CW, Chen CK. Chaos of rub-impact rotor supported by bearings with nonlinear suspension. Tribol Int 2009;42:426–39. [20] Chang-Jian CW, Chen CK. Non-linear dynamic analysis of rub-impact rotor supported by turbulent journal bearings with non-linear suspension. Int J Mech Sci 2008;50:1090–113. [21] Chang-Jian CW, Chen CK. Nonlinear analysis of a rub-impact rotor supported by turbulent couple stress fluid film journal bearings under quadratic damping. Nonlinear Dyn 2009;56:297–314. [22] Wang JG, Zhou JZ, Dong DW, Yan B, Huang CR. Nonlinear dynamic analysis of a rub-impact rotor supported by oil film bearings. Arch Appl Mech 2012;83:413–30. [23] Ma H, Yu T, Han QK, Zhang YM, Wen BC. Time-frequency features of two types of coupled rub-impact faults in rotor systems. J Sound Vib 2009;321:1109–28. [24] Patel Tejas H, Darpe Ashish K. Vibration response of a cracked rotor in presence of rotor-stator rub. J Sound Vib 2008;317:841–65. [25] Patel Tejas H, Darpe Ashish K. Coupled bending-torsional vibration analysis of rotor with rub and crack. J Sound Vib 2009;326:740–52. [26] Shen XY, Jia JH, Zhao M. Nonlinear analysis of a rub-impact rotor-bearing system with initial permanent rotor bow. Arch Appl Mech 2008;78:225–40. [27] Cao JY, Ma CB, Jiang ZD, Liu SG. Nonlinear dynamic analysis of fractional order rub-impact rotor system. Commun Nonlinear Sci Numer Simul 2011;16. 1443–1163. [28] Yuan ZW, Chu FL, Hao RJ. Simulation of rotor’s axial rub-impact in full degrees of freedom. Mech Mach Theory 2007;42:763–75. [29] Choi YS. On the contact of partial rotor rub with experimental observations. KSME Int J 2001;15:1630–8. [30] Gao J, Zhu PS, Gao ZH. Advanced flight dynamics. Beijing: National Defence Industry Press; 2007. [31] Zhu CS, Chen YJ. General dynamic model of aero-engine’s rotor system during maneuvering flight. J Aerosp Power 2009;24:371–7. [32] Zhong YE, He YZ, Wang Z, Li FZ. Rotor dynamics. Beijing: Tsinghua University Press; 1984. [33] Lichtenberg AJ, Lieberman MA. Regular and stochastic dynamics. New York: Springer-Verlag; 1992. [34] Cao Q, Wiercigroch M, Pavlovskaia EE, Grebogi C, Thompson JMT. Archetypal oscillator for smooth and discontinuous dynamics. Phys Rev E 2006;74:046218. [35] Cao Q, Xiong Y, Wiercigroch M. Resonances of the SD oscillator due to the discontinuous phase. J Appl Anal Comput 2011;1:183–91.