Stability analysis and response of nonlinear rotor-seal system

M. Sayed1 , Y. S. Hamed2

1, 2Department of Mathematics and Statistics, Faculty of Science, Taif University, El-Taif, El-Haweiah, P.O. Box 888, Zip Code 21974, Kingdom of Saudi Arabia

1, 2Department of Engineering Mathematics, Faculty of Electronic Engineering, Menoufia University, Menouf 32952, Egypt

1Corresponding author

Journal of Vibroengineering, Vol. 16, Issue 8, 2014, p. 4152-4170.
Received 13 January 2014; received in revised form 28 November 2014; accepted 2 December 2014; published 30 December 2014

Copyright © 2014 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

In this paper, we study the stability of the nonlinear rotor-seal system using Liapunov’s first method. The mathematical solutions using multiple scales up to and including second order approximations is investigated. We extract all resonance cases from analytical solution and investigated. It is quite clear that some of the simultaneous resonance cases are undesirable in the design of such system as they represent some of the worst behavior of the system. The effects of various parameters on the behavior of the system and stability of the system are investigated numerically by response curve. Poincaré maps are used to determine stability and plot bifurcation diagrams.

Keywords: multiple time scale, resonance cases, stability.

1. Introduction

Rotating machinery is widely used in engineering. In order to increase the efficiency of rotating machinery, the clearance between rotor and stator has been designed to be smaller and smaller, however, which leads to a larger possibility of impact-rub. Rubbing may induce the rotor’s dynamic instability, blade break, and other serious accidents. Therefore, how to avoid rub and how to reduce the harm caused by rubs become an important problem in rotor dynamics. With the fast development of society and economy, there is a higher demand of improvement in speed, capacity, efficiency and safety reliability in the rotating machinery, such as great power turbo-generator units, shipboard steam turbines, high speed centrifugal compressors, aeroengine and high precision machine. Under some operating conditions, the nonlinear gas exciting force of labyrinth seals in the small gap can cause rotating machines to exhibit whirl-whip instability. Particular, the gas exciting force in ultra supercritical steam turbine unit shows significant influence on the system stability and sometimes can cause disastrous consequences [1].

Up to now, lots of researches have been done. However, the experimental research of gas exciting vibration in steam turbines is hard to proceed. Many researchers have designed the rotor-seal rigs to test the gas exciting force using the similarity theory [2-3]. However, the test accuracy is seriously affected by the testing conditions. In these years, computational fluid dynamics (CFD) and the computer technology have been developed. Recently, many researchers started to use the CFD method to simulate the gas exciting force of the labyrinth seals, and have obtained better results of the dynamic characteristics [4-5]. On the other hand, in 1965, Alford [6] first derived the formula of the gas exciting force when he researched on the stability of aeroengine. Ding and Zhang [7] applied both standard Galerkin method (SGM) and nonlinear Galerkin method (NGM) to numerically investigate the lateral vibration of a continuous rotor system with actions of oil-film force and seal fluid force. Improvement of the NGM procedure over the SGM procedure is discussed and the dynamics of nonlinear rotor system is revealed. Chen et al. [8] utilized the incremental harmonic balance method to obtain the solutions of a rotor-seal system, where the complicated nonlinearities are handled by the expansion method. Periodic, double-periodic and triple-periodic solutions are obtained in excellent agreement with numerical results, which shows the validity and efficacy of the proposed solution procedures.

Li et al. [9] used the Runge–Kutta method to solve the motion equation of the rotor/seal/bearing system. The dynamic characteristics of the rotor/bearing/seal system were analyzed with bifurcation diagrams, time-history diagrams, trajectory diagrams, Poincare maps and frequency spectrums. The numerical analysis indicates that the seal force and the oil-film force influence the nonlinear dynamic characteristics of the rotor system. Zhang and Chen [10] investigated analytically and numerically the full annular rub of a typical nonlinear Jeffcott rotor. Averaging method is used to find the synchronous response of governing equations. Routh-Hurwitz criteria are employed to analyze the stability of synchronous full annular rub solution and determine the boundaries of static and Hopf bifurcations. Finally, the response and onset condition of reverse dry whip are investigated numerically and at the same time, the influences of rotor parameters and rotation speed on the characteristics of the rotor response are investigated. Cheng et al. [11] investigated the non-linear dynamic behaviors of a rotor-bearing-seal coupled system by using Muszynska’s non-linear seal fluid dynamic force model and non-linear oil film force, and the result from the numerical analysis is in agreement with the one from the experiment. The bifurcation of the coupled system is analyzed under different operating conditions. Ding et al. [12] studied Hopf bifurcation and the self-excited vibration with parametric excitation of a single span rotor-seal system, applying the central manifold theorem and normal formal form theory. The results showed that non-synchronized whirl of the imbalanced rotor can either be a quasi-periodic motion resulting from the Hopf bifurcation or a half-frequency whirl from the period doubling bifurcation. Li and Chen [13] investigated the 1:2 sub-harmonic resonance of the labyrinth seals-rotor system, where the low-frequency vibration of steam turbines can be caused by the gas exciting force. The empirical parameters of gas exciting force of the Muszynska model are obtained by using the results of computational fluid dynamics (CFD).

Sayed and Mousa [14] investigated the influence of the quadratic and cubic terms on non-linear dynamic characteristics of the angle-ply composite laminated rectangular plate with parametric and external excitations. The method of multiple time scale perturbation is applied to solve the non-linear differential equations describing the system up to and including the second-order approximation. Two cases of the sub-harmonic resonances cases in the presence of 1:2 internal resonance are considered. The stability of the system is investigated using both frequency response equations and phase-plane method. It is quite clear that some of the simultaneous resonance cases are undesirable in the design of such system as they represent some of the worst behavior of the system. Such cases should be avoided as working conditions for the system. Eissa and Sayed [15-17] and Sayed [18] studied the effects of different active controllers on simple and spring pendulum at the primary resonance via negative velocity feedback or its square or cubic. Amer et al. [19] studied the dynamical system of a twin-tail aircraft, which is described by two coupled second order nonlinear differential equations having both quadratic and cubic nonlinearities, solved and controlled. The system is subjected to both multi-parametric and multi-external excitations. The method of multiple time scale perturbation is applied to solve the nonlinear differential equations up to the two order approximations. The stability of the system is investigated applying both frequency response equations and phase plane method. Two simple active control laws based on the linear negative velocity and acceleration feedback are used. Sayed and Hamed [20] studied the response of a two-degree-of-freedom system with quadratic coupling under parametric and harmonic excitations. The method of multiple scale perturbation technique is applied to solve the non-linear differential equations and obtain approximate solutions up to and including the second-order approximations.

Hamed et al. [21-23] studied USM model subject to multi-external or both multi-external and multi-parametric and both multi-external and tuned excitation forces. The model consists of multi-degree-of-freedom system consisting of the tool holder and absorbers (tools) simulating ultrasonic machining process. The advantages of using multi-tools are to machine different materials and different shapes at the same time. This leads to time saving and higher machining efficiency. Hamed et al. [24] presented the behavior of the nonlinear string beam coupled system subjected to external, parametric and tuned excitations for case 1:1 internal resonance. The stability of the system studied using frequency response equations and phase-plane method. It is found from numerical simulations that there are obvious jumping phenomena in the frequency response curves. Sayed and Kamel [25, 26] investigated the effect of different controllers on the vibrating system and the saturation control of a linear absorber to reduce vibrations due to rotor blade flapping motion. The stability of the obtained numerical solution is investigated using both phase plane methods and frequency response equations. Sayed et al. [27] investigated the non-linear dynamics of a two-degree-of freedom vibration system including quadratic and cubic non-linearities subjected to external and parametric excitation forces. There exist multi-valued solutions which increase or decrease by the variation of some parameters. The numerical simulations show the system exhibits periodic motions and chaotic motions. Sayed and Mousa [28] studied an analytical investigation of the nonlinear vibration of a symmetric cross-ply composite laminated piezoelectric rectangular plate under parametric and external excitations. Their study focused on the case of 1:1:3 primary resonance and internal resonance, and they verified the analytical results calculated by the method of multiple time scale by comparing them with the numerical results of the modal equations. The obtained results were verified by comparing the results of the finite difference method (FDM) and Runge-Kutta (RKM) method. When plotting the solutions to some nonlinear problems, the phase space can become overcrowded and the underlying structure may become obscured. To overcome these difficulties, a basic tool of Poincaré maps was proposed by Henri Poincaré [29] at the end of the nineteenth century.

2. Model of rotor-seal system

Assume that such an anisotropic rotor mounts at both rigid ends and rotates in the gas seals, and its rotation causes gas rotation and subsequently generates the gas exciting force. Suppose that the external force has a rotating character with synchronous frequencies due to the rotor unbalance. The mathematical model of the rotor-seal system in the stationary reference coordinates (x, y) as shown in Fig. 1 can be deduced as [13]:

(1)
M + m 0 0 M + m x ¨ y ¨ + D x e + D 2 m λ Ω - 2 m λ Ω D y e + D x ˙ y ˙
            + K x e + K - m λ 2 Ω 2 λ Ω D - λ Ω D K y e + K - m λ 2 Ω 2 x y = ρ cos Ω t G + ρ sin Ω t ,

where M is the rotor mass, m is the gas inertia effects, Kxe, Kye, Dxe, Dye are the stiffness and external damping on the x and y directions respectively and G is the gravitation effects of the rotor, ρ is the amplitude of the excitation force.

Fig. 1. Model of rotor-seal system

 Model of rotor-seal system

a)

 Model of rotor-seal system

b)

3. Mathematical analysis

The motion equation of a rotor-seal system (1) can be obtained by:

(2a)
x ¨ + ω 1 2 x = f x x , y , x ˙ , y ˙ + ρ c o s Ω t ,
(2b)
y ¨ + ω 2 2 y = f y x , y , x ˙ , y ˙ + ρ s i n Ω t ,

with initial condition x(0)=x˙(0)=y(0)=y˙(0)=0. Where x and y are the displacements of the shaft center, ω1 and ω2 the linear natural frequencies in x, y directions, and Ω the angular speed of the rotor, ρ is the amplitude of the excitation force, fx(x,y,x˙,y˙) and fy(x,y,x˙,y˙) are power series function of the nonlinear force on x and y directions given by:

f x x , y , x ˙ , y ˙ = α 1 x + α 2 x ˙ + α 3 y + α 4 y ˙ + α 5 x 2 + α 6 y 2 + α 7 x x ˙ + α 8 x y + α 9 x y ˙
              + α 10 x ˙ y + α 11 y y ˙ + α 12 x 3   + α 13 y 3 + α 14 x 2 x ˙ + α 15 x 2 y + α 16 x 2 y ˙ + α 17 x y 2
              + α 18 y 2 x ˙ + α 19 y 2 y ˙ + α 20 x y x ˙ + α 21 x y y ˙ ,
f y x , y , x ˙ , y ˙ = β 1 x + β 2 x ˙ + β 3 y + β 4 y ˙ + β 5 x 2 + β 6 y 2 + β 7 x x ˙ + β 8 x y + β 9 x y ˙
              + β 10 x ˙ y + β 11 y y ˙ + β 12 x 3   + β 13 y 3 + β 14 x 2 x ˙ + β 15 x 2 y + β 16 x 2 y ˙ + β 17 x y 2
              + β 18 y 2 x ˙ + β 19 y 2 y ˙ + β 20 x y x ˙ + β 21 x y y ˙ .

All coefficient of αi, βi are defined in [13]. The linear, non-linear coefficients and excitation amplitude are assumed to be:

(3)
α n = ε α ^ n ,       β n = ε β ^ n ,       ρ = ε ρ ^ ,       n = 1 , 2 , . . . , 21 ,

where ε is a small perturbation parameter and 0<ε<<1.

To consider the influence of the quadratic and cubic terms on non-linear dynamic characteristics of the non-linear rotor-seal system, we need to obtain the second-order approximate solution of Eqs. (2). We determine a second-order approximation for the system response by the method of multiple scales [30-31], which is a powerful tool in determining periodic solutions of small amplitude. For the second-order approximation, introducing three time scales:

(4)
T n = ε n t ,       n = 0 , 1 , 2 .

We seek an asymptotic solution of Eqs. (2) for x and y in the form of:

(5)
x t , ε = x 0 T 0 , T 1 , T 2 + ε x 1 T 0 , T 1 , T 2 + ε 2 x 2 T 0 , T 1 , T 2 + O ε 3 ,
(6)
y t , ε = y 0 T 0 , T 1 , T 2 + ε y 1 T 0 , T 1 , T 2 + ε 2 y 2 T 0 , T 1 , T 2 + O ε 3 .

Derivatives with respect to t are then transformed into:

(7)
d d t = T 0 T 0 t + T 1 T 1 t + T 2 T 2 t = D 0 + ε D 1 + ε 2 D 2 ,
(8)
d 2 d t 2 = D 0 2 + 2 ε D 0 D 1 + ε 2 ( D 1 2 + 2 D 0 D 2 ) ,

where Dn=/Tn, n= 0, 1, 2. Terms of O(ε3) and higher orders are neglected. Substituting Eqs. (5)-(6) and (7)-(8) into Eqs. (2) and equating the coefficients of similar powers of ε, one obtains the following set of ordinary differential equations:

Order ε0:

(9)
( D 0 2 + ω 1 2 ) x 0 = 0 ,
(10)
( D 0 2 + ω 2 2 ) y 0 = 0 .

Order ε1:

(11)
D 0 2 + ω 1 2 x 1 = - 2 D 0 D 1 x 0 + α ^ 1 x 0 + α ^ 2 D 0 x 0 + α ^ 3 y 0 + α ^ 4 D 0 y 0 + α ^ 5 x 0 2 + α ^ 6 y 0 2
              + α ^ 7 x 0 D 0 x 0 + α ^ 8 x 0 y 0 + α ^ 9 x 0 D 0 y 0 + α ^ 10 y 0 D 0 x 0 + α ^ 11 y 0 D 0 y 0 + α ^ 12 x 0 3 + α ^ 13 y 0 3
              + α ^ 14 x 0 2 D 0 x 0 + α ^ 15 x 0 2 y 0 + α ^ 16 x 0 2 D 0 y 0 + α ^ 17 x 0 y 0 2 + α ^ 18 y 0 2 D 0 x 0 + α ^ 19 y 0 2 D 0 y 0
              + α ^ 20 x 0 y 0 D 0 x 0 + α ^ 21 x 0 y 0 D 0 y 0 + ρ ^ c o s Ω t ,
(12)
D 0 2 + ω 2 2 y 1 = - 2 D 0 D 1 y 0 + β ^ 1 x 0 + β ^ 2 D 0 x 0 + β ^ 3 y 0 + β ^ 4 D 0 y 0 + β ^ 5 x 0 2 + β ^ 6 y 0 2
              + β ^ 7 x 0 D 0 x 0 + β ^ 8 x 0 y 0 + β ^ 9 x 0 D 0 y 0 + β ^ 10 y 0 D 0 x 0 + β ^ 11 y 0 D 0 y 0 + β ^ 12 x 0 3 + β ^ 13 y 0 3
              + β ^ 14 x 0 2 D 0 x 0 + β ^ 15 x 0 2 y 0 + β ^ 16 x 0 2 D 0 y 0 + β ^ 17 x 0 y 0 2 + β ^ 18 y 0 2 D 0 x 0 + β ^ 19 y 0 2 D 0 y 0
              + β ^ 20 x 0 y 0 D 0 x 0 + β ^ 21 x 0 y 0 D 0 y 0 + ρ ^ s i n Ω t .

Order ε2:

(13)
D 0 2 + ω 1 2 x 2 = - D 1 2 x 1 - 2 D 0 D 2 x 0 - 2 D 0 D 1 x 1 + α ^ 1 x 1 + α ^ 2 D 0 x 1 + D 1 x 0 + α ^ 3 y 1
              + α ^ 4 D 0 y 1 + D 1 y 0 + 2 α ^ 5 x 0 x 1 + 2 α ^ 6 y 0 y 1 + α ^ 7 x 1 D 0 x 0 + α ^ 7 x 0 D 0 x 1 + D 1 x 0
              + α ^ 8 x 0 y 1 + x 1 y 0 + α ^ 9 x 1 D 0 y 0 + α ^ 9 x 0 D 0 y 1 + D 1 y 0 + α ^ 10 y 1 D 0 x 0
              + α ^ 10 y 0 D 0 x 1 + D 1 x 0 + α ^ 11 y 1 D 0 y 0 + α ^ 11 y 0 D 0 y 1 + D 1 y 0 + 3 α ^ 12 x 0 2 x 1
              + 3 α ^ 13 y 0 2 y 1 + 2 α ^ 14 x 0 x 1 D 0 x 0 + α ^ 14 x 0 2 D 0 x 1 + D 1 x 0 + α ^ 15 x 0 2 y 1 + 2 α ^ 15 x 0 y 0 x 1
              + 2 α ^ 16 x 0 x 1 D 0 y 0 + α ^ 16 x 0 2 D 0 y 1 + D 1 y 0 + 2 α ^ 17 y 0 y 1 x 0 + α ^ 17 y 0 2 x 1 + 2 α ^ 18 y 0 y 1 D 0 x 0
              + α ^ 18 y 0 2 D 0 x 1 + D 1 x 0 + α ^ 19 y 0 2 D 0 y 1 + D 1 y 0 + 2 α ^ 19 y 0 y 1 D 0 y 0
              + α ^ 20 x 0 y 0 D 0 x 1 + D 1 x 0 + α ^ 20 x 0 y 1 + α ^ 20 x 1 y 0 D 0 x 0 + α ^ 21 x 0 y 0 D 0 y 1 + D 1 y 0
              + α ^ 21 x 0 y 1 + α ^ 21 x 1 y 0 D 0 y 0 ,
(14)
D 0 2 + ω 2 2 y 2 = - D 1 2 y 1 - 2 D 0 D 2 y 0 - 2 D 0 D 1 y 1 + β ^ 1 x 1 + β ^ 2 D 0 x 1 + D 1 x 0 + β ^ 3 y 1
              + β ^ 4 D 0 y 1 + D 1 y 0 + 2 β ^ 5 x 0 x 1 + 2 β ^ 6 y 0 y 1 + β ^ 7 x 1 D 0 x 0 + β ^ 7 x 0 D 0 x 1 + D 1 x 0
              + β ^ 8 x 0 y 1 + x 1 y 0 + β ^ 9 x 1 D 0 y 0 + β ^ 9 x 0 D 0 y 1 + D 1 y 0 + β ^ 10 y 1 D 0 x 0
              + β ^ 10 y 0 D 0 x 1 + D 1 x 0 + β ^ 11 y 1 D 0 y 0 + β ^ 11 y 0 D 0 y 1 + D 1 y 0 + 3 β ^ 12 x 0 2 x 1
              + 3 β ^ 13 y 0 2 y 1 + 2 β ^ 14 x 0 x 1 D 0 x 0 + β ^ 14 x 0 2 D 0 x 1 + D 1 x 0 + β ^ 15 x 0 2 y 1 + 2 β ^ 15 x 0 y 0 x 1
              + 2 β ^ 16 x 0 x 1 D 0 y 0 + β ^ 16 x 0 2 D 0 y 1 + D 1 y 0 + 2 β ^ 17 y 0 y 1 x 0 + β ^ 17 y 0 2 x 1 + 2 β ^ 18 y 0 y 1 D 0 x 0
              + β ^ 18 y 0 2 D 0 x 1 + D 1 x 0 + β ^ 19 y 0 2 D 0 y 1 + D 1 y 0 + 2 β ^ 19 y 0 y 1 D 0 y 0
              + β ^ 20 x 0 y 0 D 0 x 1 + D 1 x 0 + β ^ 20 x 0 y 1 + β ^ 20 x 1 y 0 D 0 x 0 + β ^ 21 x 0 y 0 D 0 y 1 + D 1 y 0
              + β ^ 21 x 0 y 1 + β ^ 21 x 1 y 0 D 0 y 0 .

The general solutions of Eqs. (9)-(10), can be written in the form:

(15)
x 0 = A 1 T 1 , T 2 e x p i ω 1 T 0 + c c ,
(16)
y 0 = A 2 T 1 , T 2 e x p i ω 2 T 0 + c c ,

where A1 and A2 are a complex function in T1, T2 and cc indicates the complex conjugates of the preceding terms. Substituting Eqs. (15)-(16) into Eqs. (11)-(12) and after eliminating the secular terms, the non-homogeneous solutions of Eqs. (11)-(12) are:

(17)
x 1 = E 1 e x p 2 i ω 1 T 0 + E 2 e x p 3 i ω 1 T 0 + E 3 e x p i ω 2 T 0 + E 4 e x p 2 i ω 2 T 0
            + E 5 e x p 3 i ω 2 T 0 + E 6 e x p i ω 1 ± ω 2 T 0 + E 7 e x p i ω 1 ± 2 ω 2 T 0
            + E 8 e x p i 2 ω 1 ± ω 2 T 0 + E 9 e x p i Ω T 0 + E 10 + c c ,
(18)
y 1 = E 11 e x p i ω 1 T 0 + E 12 e x p 2 i ω 1 T 0 + E 13 e x p 3 i ω 1 T 0 + E 14 e x p 2 i ω 2 T 0
            + E 15 e x p 3 i ω 2 T 0 + E 16 e x p i ω 1 ± ω 2 T 0 + E 17 e x p i ω 1 ± 2 ω 2 T 0
            + E 18 e x p i 2 ω 1 ± ω 2 T 0 + E 19 e x p i Ω T 0 + E 20 + c c ,

where Ei (i=1, 2, ..., 20) are complex functions in T1 and T2. Substituting Eqs. (15)-(18) into Eqs. (13)-(14) and after eliminating the secular terms, the non-homogeneous solutions of Eqs. (13)-(14) are:

(19)
x 2 = H 1 e x p ( 2 i ω 1 T 0 ) + H 2 e x p ( 3 i ω 1 T 0 ) + H 3 e x p ( 4 i ω 1 T 0 ) + H 4 e x p ( 5 i ω 1 T 0 )  
            + H 5 e x p ( i ω 2 T 0 ) + H 6 e x p ( 2 i ω 2 T 0 ) + H 7 e x p ( 3 i ω 2 T 0 ) + H 8 e x p ( 4 i ω 2 T 0 )
            + H 9 e x p ( 5 i ω 2 T 0 ) + H 10 e x p ( i ( ω 1 ± ω 2 ) T 0 ) + H 11 e x p ( i ( ω 1 ± 2 ω 2 ) T 0 )
            + H 12 e x p ( i ( 2 ω 1 ± ω 2 ) T 0 ) + H 13 e x p ( i ( 2 ω 1 ± 2 ω 2 ) T 0 ) + H 14 e x p ( i ( ω 1 ± 3 ω 2 ) T 0 )
            + H 15 e x p ( i ( ω 1 ± 4 ω 2 ) T 0 ) + H 16 e x p ( i ( 2 ω 1 ± 3 ω 2 ) T 0 ) + H 17 e x p ( i ( 3 ω 1 ± 2 ω 2 ) T 0 )
            + H 18 e x p ( i ( 3 ω 1 ± ω 2 ) T 0 ) + H 19 e x p ( i ( 4 ω 1 ± ω 2 ) T 0 ) + H 20 e x p ( i Ω T 0 )
            + H 21 e x p ( 2 i Ω T 0 ) + H 22 e x p ( 3 i Ω T 0 ) + H 23 e x p ( 4 i Ω T 0 ) + H 24 e x p ( 5 i Ω T 0 )
            + H 25 e x p ( i ( ω 1 ± Ω ) T 0 ) + H 26 e x p ( i ( 2 ω 1 ± Ω ) T 0 ) + H 27 e x p ( i ( ω 1 ± 2 Ω ) T 0 )
            + H 28 e x p ( i ( ω 1 ± 3 Ω ) T 0 ) + H 29 e x p ( i ( ω 1 ± 4 Ω ) T 0 ) + H 30 e x p ( i ( 2 ω 1 ± 2 Ω ) T 0 )
            + H 31 e x p ( i ( 2 ω 1 ± 3 Ω ) T 0 ) + H 32 e x p ( i ( 3 ω 1 ± Ω ) T 0 ) + H 33 e x p ( i ( 3 ω 1 ± 2 Ω ) T 0 )
            + H 34 e x p ( i ( 4 ω 1 ± Ω ) T 0 ) + H 35 e x p ( i ( ω 2 ± Ω ) T 0 ) + H 36 e x p ( i ( 2 ω 2 ± Ω ) T 0 )
            + H 37 e x p ( i ( ω 2 ± 2 Ω ) T 0 ) + H 38 e x p ( i ( 2 ω 2 ± 2 Ω ) T 0 ) + H 39 e x p ( i ( 2 ω 2 ± 3 Ω ) T 0 )
            + H 40 e x p ( i ( ω 2 ± 3 Ω ) T 0 ) + H 41 e x p ( i ( ω 2 ± 4 Ω ) T 0 ) + H 42 e x p ( i ( 3 ω 2 ± 2 Ω ) T 0 )
            + H 43 e x p ( i ( 3 ω 2 ± Ω ) T 0 ) + H 44 e x p ( i ( 4 ω 2 ± Ω ) T 0 ) + H 45 e x p ( i ( 2 ω 1 ± ω 2 ± Ω ) T 0 )
            + H 46 e x p ( i ( ω 1 ± ω 2 ± Ω ) T 0 ) + H 47 e x p ( i ( ω 1 ± ω 2 ± 2 Ω ) T 0 )
            + H 48 e x p ( i ( 3 ω 1 ± ω 2 ± Ω ) T 0 ) + H 49 e x p ( i ( ω 1 ± 3 ω 2 ± Ω ) T 0 )
            + H 50 e x p ( i ( ω 1 ± ω 2 ± 3 Ω ) T 0 ) Ω ) T 0 ) + H 51 e x p ( i ( 2 ω 1 ± 2 ω 2 ± Ω ) T 0 )
            + H 52 e x p ( i ( 2 ω 1 ± ω 2 ± 2 Ω ) T 0 ) + H 53 e x p ( i ( ω 1 ± 2 ω 2 ± 2 Ω ) T 0 ) + H 54 + c c ,
(20)
y 2 = G 1 e x p ( i ω 1 T 0 ) + G 2 e x p ( 2 i ω 1 T 0 ) + G 3 e x p ( 3 i ω 1 T 0 ) + G 4 e x p ( 4 i ω 1 T 0 )
            + G 5 e x p ( 5 i ω 1 T 0 ) + G 6 e x p ( 2 i ω 2 T 0 ) + G 7 e x p ( 3 i ω 2 T 0 ) + G 8 e x p ( 4 i ω 2 T 0 )
            + G 9 e x p ( 5 i ω 2 T 0 ) + G 10 e x p ( i ( ω 1 ± ω 2 ) T 0 ) + G 11 e x p ( i ( ω 1 ± 2 ω 2 ) T 0 )
            + G 12 e x p ( i ( 2 ω 1 ± ω 2 ) T 0 ) + G 13 e x p ( i ( 2 ω 1 ± 2 ω 2 ) T 0 ) + G 14 e x p ( i ( ω 1 ± 3 ω 2 ) T 0 )
            + G 15 e x p ( i ( ω 1 ± 4 ω 2 ) T 0 ) + G 16 e x p ( i ( 2 ω 1 ± 3 ω 2 ) T 0 ) + G 17 e x p ( i ( 3 ω 1 ± 2 ω 2 ) T 0 )
            + G 18 e x p ( i ( 3 ω 1 ± ω 2 ) T 0 ) + G 19 e x p ( i ( 4 ω 1 ± ω 2 ) T 0 ) + G 20 e x p ( i Ω T 0 )
            + G 21 e x p ( 2 i Ω T 0 )   + G 22 e x p ( 3 i Ω T 0 ) + G 23 e x p ( 4 i Ω T 0 ) + G 24 e x p ( 5 i Ω T 0 )
            + G 25 e x p ( i ( ω 1 ± Ω ) T 0 ) + G 26 e x p ( i ( 2 ω 1 ± Ω ) T 0 ) + G 27 e x p ( i ( ω 1 ± 2 Ω ) T 0 )
            + G 28 e x p ( i ( ω 1 ± 3 Ω ) T 0 ) + G 29 e x p ( i ( ω 1 ± 4 Ω ) T 0 ) + G 30 e x p ( i ( 2 ω 1 ± 2 Ω ) T 0 )
            + G 31 e x p ( i ( 2 ω 1 ± 3 Ω ) T 0 ) + G 32 e x p ( i ( 3 ω 1 ± Ω ) T 0 ) + G 33 e x p ( i ( 3 ω 1 ± 2 Ω ) T 0 )
            + G 34 e x p ( i ( 4 ω 1 ± Ω ) T 0 ) + G 35 e x p ( i ( ω 2 ± Ω ) T 0 ) + G 36 e x p ( i ( 2 ω 2 ± Ω ) T 0 )
            + G 37 e x p ( i ( ω 2 ± 2 Ω ) T 0 ) + G 38 e x p ( i ( 2 ω 2 ± 2 Ω ) T 0 ) + G 39 e x p ( i ( 2 ω 2 ± 3 Ω ) T 0 )
            + G 40 e x p ( i ( ω 2 ± 3 Ω ) T 0 ) + G 41 e x p ( i ( ω 2 ± 4 Ω ) T 0 ) + G 42 e x p ( i ( 3 ω 2 ± 2 Ω ) T 0 )
            + G 43 e x p ( i ( 3 ω 2 ± Ω ) T 0 ) + G 44 e x p ( i ( 4 ω 2 ± Ω ) T 0 ) + G 45 e x p ( i ( 2 ω 1 ± ω 2 ± Ω ) T 0 )
            + G 46 e x p ( i ( ω 1 ± ω 2 ± Ω ) T 0 ) + G 47 e x p ( i ( ω 1 ± ω 2 ± 2 Ω ) T 0 )
            + G 48 e x p ( i ( 3 ω 1 ± ω 2 ± Ω ) T 0 ) + G 49 e x p ( i ( ω 1 ± 3 ω 2 ± Ω ) T 0 )
            + G 50 e x p ( i ( ω 1 ± ω 2 ± 3 Ω ) T 0 ) Ω ) T 0 ) + G 51 e x p ( i ( 2 ω 1 ± 2 ω 2 ± Ω ) T 0 )
            + G 52 e x p ( i ( 2 ω 1 ± ω 2 ± 2 Ω ) T 0 ) + G 53 e x p ( i ( ω 1 ± 2 ω 2 ± 2 Ω ) T 0 ) + G 54 + c c ,

where Hi, Gi (i=1, 2, ..., 54) are complex functions in T1 and T2. From the above derived solutions, the reported resonance cases are:

• Primary resonance: Ωωn, n=1, 2.

• Sub-harmonic resonance: Ωkωn, k= 2, 3, 4, 5 and n=1, 2.

• Super-harmonic resonance: kΩωn, 3Ω2ωn, 2Ω3ωn, k= 2, 3, 4, 5 and n=1, 2.

Internal or secondary resonance: ω1sω2, ω2sω1, 2ω13ω2, 3ω12ω2s= 2, 3, 4, 5.

• Combined resonance: ω1±mω2±Ω,ω1±ω2±rΩ,ω1±2ω2±2Ω,ω1±2ω2±3Ω, ω1±3ω2±2Ω,2ω1±qω2±Ω,2ω1±ω2±pΩ,3ω1±nω2±Ω,3ω1±  ω2±2Ω,4ω1±ω2±Ω,m= 1, 2, 3, 4, r=2, 3, 4, q= 1, 2, 3, p= 2, 3, n=1, 2.

• Simultaneous or incident resonance.

Any combination of the above resonance cases is considered as simultaneous resonance.

4. Stability analysis

To describe quantitatively the closeness of the resonances, we introduce the detuning parameters σ1 and σ2 according to:

(21)
Ω   = ω 1 + ε σ 1 ,       ω 2 = 3 ω 1 + ε σ 2 .

Substituting Eq. (21) into Eqs. (11)-(12) and eliminating the secular and small divisor terms from x1 and y1, we get the following:

(22)
2 i ω 1 D 1 A = α ^ 1 + i α ^ 2 ω 1 A + 2 i α ^ 18 ω 1 + 2 α ^ 17 B B - A
              + α ^ 15 + i α ^ 16 ω 2 - i α ^ 20 ω 1 A - 2 B e x p i σ 2 T 1 + 3 α ^ 12 + i α ^ 14 ω 1 A 2 A - + ρ ^ 2 e x p i σ 1 T 1 ,
(23)
2 i ω 2 D 1 B = β ^ 3 + i β ^ 4 ω 2 B + 3 β ^ 13 + i β ^ 19 ω 2 B 2 B -
              + 2 β ^ 15 + 2 i β ^ 16 ω 2 A A - B + β ^ 12 + i β ^ 14 ω 1 A 3 e x p - i σ 2 T 1 .

To analyze the solution of Eqs. (22)-(23), it is convenient to express the A(T1) and B(T1) in the polar form as:

(24)
A T 1 = 1 2 a T 1 e i γ 1 T 1 ,         B T 1 = 1 2 b T 1 e i γ 2 T 1 ,

where a, b, γ1 and γ2 are unknown real-valued functions. Inserting Eq. (24) into Eqs. (22)-(23) and separating real and imaginary parts, we have:

(25)
a ' = α 2 2 a + α 14 8 a 3 + α 18 4 a b 2 + α 15 8 ω 1 a 2 b s i n θ 2 + α 16 ω 2 - α 20 ω 1 8 ω 1 a 2 b c o s θ 2
            + ρ 2 ω 1 s i n θ 1 ,
(26)
a γ 1 ' = - α 1 2 ω 1 a - 3 α 12 8 ω 1 a 3 - α 17 4 ω 1 a b 2 - α 15 8 ω 1 a 2 b c o s θ 2 + α 16 ω 2 - α 20 ω 1 8 ω 1 a 2 b s i n θ 2
            - ρ 2 ω 1 c o s θ 1 ,
(27)
b ' = β 4 2 b + β 19 8 b 3 + β 16 4 a 2 b - β 12 8 ω 2 a 3 s i n θ 2 + β 14 ω 1 8 ω 2 a 3 c o s θ 2 ,
(28)
b γ 2 ' = - β 3 2 ω 2 b - 3 β 13 8 ω 2 b 3 - β 15 4 ω 2 a 2 b - β 12 8 ω 2 a 3 c o s θ 2 - β 14 ω 1 8 ω 2 a 3 s i n θ 2 ,

where:

(29)
θ 1 = σ 1 T 1 - γ 1 ,         θ 2 = σ 2 T 1 + γ 2 - 3 γ 1 .

Then, it follows from Eq. (29) that:

(30)
γ 1 ' = σ 1 - θ 1 ' ,
(31)
γ 2 ' = 3 σ 1 - σ 2 + θ 2 ' - 3 θ 1 ' .

The control strategy is evaluated by analyzing the steady state response of the closed-loop system. Setting a'=b'=0 and θm'=0 (m= 1, 2) into Eqs. (25)-(28) we obtained:

(32)
α 2 2 a + α 14 8 a 3 + α 18 4 a b 2 + α 15 8 ω 1 a 2 b s i n θ 2
+ α 16 ω 2 - α 20 ω 1 8 ω 1 a 2 b c o s θ 2 + ρ 2 ω 1 s i n θ 1 = 0 ,
(33)
a σ 1 + α 1 2 ω 1 a + 3 α 12 8 ω 1 a 3 + α 17 4 ω 1 a b 2 + α 15 8 ω 1 a 2 b c o s θ 2 - α 16 ω 2 - α 20 ω 1 8 ω 1 a 2 b s i n θ 2
            + ρ 2 ω 1 c o s θ 1 = 0 ,
(34)
β 4 2 b + β 19 8 b 3 + β 16 4 a 2 b - β 12 8 ω 2 a 3 s i n θ 2 + β 14 ω 1 8 ω 2 a 3 c o s θ 2 = 0 ,
(35)
b ( 3 σ 1 - σ 2 ) + β 3 2 ω 2 b + 3 β 13 8 ω 2 b 3 + β 15 4 ω 2 a 2 b + β 12 8 ω 2 a 3 c o s θ 2 + β 14 ω 1 8 ω 2 a 3 s i n θ 2 = 0 .

Solving the resulting algebraic equations yields two possibilities for the fixed points for each case.

Case 1: In this case a 0,b= 0 and from Eqs. (32)-(35), the frequency response equation can be obtained in the form:

(36)
α 2 2 a + α 14 8 a 3 2 + a σ 1 + α 1 2 ω 1 a + 3 α 12 8 ω 1 a 3 2 - ρ 2 4 ω 1 2 = 0 .

Case 2: In this case a,b 0 and from Eqs. (32)-(35), the resulting two equations are:

(37)
α 2 2 a + α 14 8 a 3 + α 18 4 a b 2 2 + a σ 1 + α 1 2 ω 1 a + 3 α 12 8 ω 1 a 3 + α 17 4 ω 1 a b 2 2 - α 15 8 ω 1 a 2 b 2
            - ρ 2 4 ω 1 2 - α 16 ω 2 - α 20 ω 1 8 ω 1 2 a 4 b 2 - α 15 ρ 8 ω 1 2 a 2 b c o s ( θ 2 - θ 1 ) = 0 ,
(38)
β 4 2 b + β 19 8 b 3 + β 16 4 a 2 b 2 + b 3 σ 1 - σ 2 + β 3 2 ω 2 b + 3 β 13 8 ω 2 b 3 + β 15 4 ω 2 a 2 b 2
            - β 12 8 ω 2 2 a 6 - β 14 ω 1 8 ω 2 2 a 6 = 0 .

4.1. Non-linear solution

To determine the stability of the fixed points, one lets:

(39)
a = a 0 + a 1 ,       b = b 0 + b 1 ,         θ m = θ m 0 + θ m 1 ,       m = 1 ,   2 ,

where a0, b0 and θm0 are the solutions of Eqs. (32)-(35) and a1, b1, θm1 are perturbations which are assumed to be small compared to a0, b0 and θm0. Substituting Eq. (39) into Eqs. (25)-(28), using Eqs. (32)-(35) and keeping only the linear terms in a1, b1, θm1 we obtain:

1) For the solution (a 0,b= 0), we have:

(40)
a 1 ' = α 2 2 + 3 α 14 8 a 0 2 a 1 + ρ 2 ω 1 c o s θ 10 θ 11 ,
(41)
θ 11 ' = σ 1 a 0 + α 1 2 ω 1 a 0 + 9 α 12 8 ω 1 a 0 a 1 - ρ 2 a 0 ω 1 s i n θ 10 θ 11 .

The stability of a given fixed point to a disturbance proportional to exp(λt) is determined by the roots of:

(42)
A = α 2 2 + 3 α 14 8 a 0 2          ρ 2 ω 1 c o s θ 10 σ 1 a 0 + α 1 2 ω 1 a 0 + 9 α 12 a 0 8 ω 1         - ρ 2 a 0 ω 1 s i n θ 10 .

Consequently, a non-trivial solution is stable if and only if the real parts of both eigenvalues of the coefficient matrix Eq. (42) are less than zero. Solid / dotted lines denote stable / unstable solution on the response curves, respectively.

2) For the practical solution (a, b 0), the following are obtained:

(43)
a 1 ' = α 2 2 + 3 α 14 8 a 0 2 + α 18 4 b 0 2 + α 15 4 ω 1 a 0 b 0 s i n θ 20 + α 16 ω 2 - α 20 ω 1 4 ω 1 a 0 b 0 c o s θ 20 a 1
            + ρ 2 ω 1 c o s θ 10 θ 11 + α 18 2 a 0 b 0 + α 15 8 ω 1 a 0 2 s i n θ 20 + α 16 ω 2 - α 20 ω 1 8 ω 1 a 0 2 c o s θ 20 b 1
            + α 15 8 ω 1 a 0 2 b 0 c o s θ 20 - α 16 ω 2 - α 20 ω 1 8 ω 1 a 0 2 b 0 s i n θ 20 θ 21 ,
(44)
θ 11 ' = σ 1 a 0 + α 1 2 ω 1 a 0 + 9 α 12 8 ω 1 a 0 + α 17 4 ω 1 a 0 b 0 2 + α 15 4 ω 1 b 0 c o s θ 20 - α 16 ω 2 - α 20 ω 1 4 ω 1 b 0 s i n θ 20 a 1
- ρ 2 a 0 ω 1 s i n θ 10 θ 11 + α 17 2 ω 1 b 0 + α 15 8 ω 1 a 0 c o s θ 20 - α 16 ω 2 - α 20 ω 1 8 ω 1 a 0 s i n θ 20 b 1
- α 15 8 ω 1 a 0 b 0 s i n θ 20 - α 16 ω 2 - α 20 ω 1 8 ω 1 a 0 b 0 c o s θ 20 θ 21 ,
(45)
b 1 ' = β 16 2 a 0 b 0 - 3 β 12 8 ω 2 a 0 2 s i n θ 20 + 3 β 14 ω 1 8 ω 2 a 0 2 c o s θ 20 a 1 + β 4 2 + 3 β 19 8 b 0 2 + β 16 4 a 0 2 b 1
            - β 12 8 ω 2 a 0 3 c o s θ 20 + β 14 ω 1 8 ω 2 a 0 3 s i n θ 20 θ 21 ,
(46)
θ 20 ' = 3 σ 1 a 0 + 3 α 1 2 ω 1 a 0 + 27 α 12 8 ω 1 a 0 + 3 α 17 4 ω 1 a 0 b 0 2 + 3 α 15 4 ω 1 b 0 c o s θ 20 - 3 α 16 ω 2 - α 20 ω 1 4 ω 1 b 0 s i n θ 20 - β 15 2 ω 2 a 0 - 3 β 12 8 ω 2 b 0 a 0 2 c o s θ 20 - 3 β 14 ω 1 8 ω 2 b 0 a 0 2 s i n θ 20 a 1
- 3 ρ 2 a 0 ω 1 s i n θ 10 θ 11 + 3 α 17 2 ω 1 b 0 + 3 α 15 8 ω 1 a 0 c o s θ 20 - 3 α 16 ω 2 - α 20 ω 1 8 ω 1 a 0 s i n θ 20 + σ 2 - 3 σ 1 b 0 - β 3 2 ω 2 b 0 - 9 β 13 8 ω 2 b 0 - β 15 4 ω 2 b 0 a 0 2 b 1
+ - 3 α 15 8 ω 1 a 0 b 0 s i n θ 20 + 3 α 16 ω 2 - α 20 ω 1 8 ω 1 a 0 b 0 c o s θ 20 + β 12 8 ω 2 b 0 a 0 3 s i n θ 20 - β 14 ω 1 8 ω 2 b 0 a 0 3 c o s θ 20 θ 21 .

The stability of a particular fixed point with respect to perturbations proportional to exp(λt) depends on the real parts of the roots of the matrix. Thus, a fixed point given by Eqs. (43)-(46) is asymptotically stable if and only if the real parts of all roots of the matrix are negative.

5. Results and discussions

The differential equation of the rotor-seal system is solved numerically using Rung-Kutta fourth order method and MATLAB, MAPLE packages as shown in Figs. 2-4. Fig. 2 illustrates the response and the phase-plane for the non-resonant rotor-seal system (basic case) where Ωω1ω2 at some practical values of the equation parameters α1=0.03,α2=0.0001,α3= –0.05, α4= –0.2, α5=0.01, α6=0.02, α7=0.25, α8=0.004, α9=0.7, α10= –0.8, α11= –0.2, α12=0.004, α13=0.01, α14=0.02, α15=0.003, α16= –0.04, α17=0.001,α18= –0.5, α19= –0.03, α20= –0.04, α21=0.5, β1=0.05, β2=0.2, β3= –0.03, β4=0.0001, β5=0.002, β6=0.01, β7=0.2, β8=0.004, β9=0.8, β10=0.7, β11= –0.25, β12= –0.01, β13= –0.004, β14=0.03, β15= –0.001, β16=0.5, β17= –0.003, β18=0.04, β19= –0.02, β20= –0.5, β21=0.04, Ω=3, ω1= 2.8, ω2=3.25, ρ=1. It is observed from this Fig. 2 that the response of the first and second modes of the rotor-seal system start with increasing amplitude and becomes stable respectively and the steady state amplitudes x and y are about 0.3 and 0.7 respectively and the phase plane shows limit cycle, denoting that the system is free from chaos. Table 1 shows the results of the worst resonance conditions.

Fig. 2. Non-resonance rotor-seal system response (basic case)

 Non-resonance rotor-seal system response (basic case)
 Non-resonance rotor-seal system response (basic case)
 Non-resonance rotor-seal system response (basic case)
 Non-resonance rotor-seal system response (basic case)

For the simultaneous primary and internal resonance (Ωω1 and ω23ω1) which is one of the worst resonance cases as shown in Fig. 3, it can be shown that for the first and second modes of the rotor-seal system, the steady state amplitudes are increased to about 670 % and 1140 % respectively of that values shown in Fig. 2. The vibration of the first mode is increased and become stable with tuned oscillation. Also the response of the second mode is increased and becomes stale with some chaotic. The phase-plane shows a multi-limit cycle for the first and second modes.

Fig. 4 shows that the time response and phase-plane of the simultaneous primary and internal resonance (Ωω1 and ω2ω1) which another one of the worst resonance cases. From this figure we have that the amplitude of the first mode of the rotor-seal system is increased to about 670 % while the amplitude of the second mode is increased to about 570 % of the of that values shown in Fig. 2. The oscillations of the first and second modes become stable, chaotic and tuned and the phase plane shows multi-limit cycle.

5.1. Effects of different parameters

In this section, the effects of different parameters near the simultaneous primary and internal resonance case is investigated and studied. The frequency response equations are given by Eqs. (37) and (38) are solved numerically at the same values of the parameters shown in Fig. 2. In all figures, the solid lines stand for the stable solution and the dashed lines for the unstable solution.

Table 1. Summary of the worst resonance cases

Resonance cases
x (%)
y (%)
Remarks
Non-resonance case
100 %
100 %
Limit cycle
Ω ω 1
700 %
620 %
Multi limit cycle
Ω ω 1 2
560 %
570 %
Multi limit cycle
Ω ω 1 3
560 %
580 %
Multi limit cycle
Ω 2 ω 1 3
600 %
570 %
Multi limit cycle
Ω 3 ω 1 2
500 %
530 %
Multi limit cycle
Ω ω 2 2
560 %
570 %
Multi limit cycle
Ω ω 2 3
560 %
580 %
Multi limit cycle
Ω 2 ω 2 3
600 %
570 %
Multi limit cycle
Ω 3 ω 2 2
500 %
500 %
Multi limit cycle
Ω ω 1 ,   ω 2 ω 1
670 %
570 %
Multi limit cycle
Ω ω 1 ,   ω 2 2 ω 1
730 %
1000 %
Multi limit cycle
Ω ω 1 ,   ω 2 3 ω 1
670 %
1140 %
Multi limit cycle
Ω ω 1 ,   ω 2 4 ω 1
370 %
200 %
Multi limit cycle
Ω ω 1 ,   ω 2 5 ω 1
800 %
240 %
Multi limit cycle
Ω ω 1 ,   ω 2 ω 1 2
85 %
320 %
Multi limit cycle
Ω ω 1 ,   ω 2 ω 1 3
400 %
530 %
Multi limit cycle
Ω ω 1 ,   ω 2 ω 1 4
200 %
430 %
Multi limit cycle
Ω ω 1 ,   ω 2 2 ω 1 3
220 %
320 %
Multi limit cycle
Ω ω 2 ,   ω 1 ω 2
830 %
430 %
Multi limit cycle
Ω ω 2 ,  ω 1 2 ω 2
85 %
430 %
Multi limit cycle
Ω ω 2 ,  ω 1 3 ω 2
50 %
500 %
Multi limit cycle
Ω ω 2 ,  ω 1 4 ω 2
35 %
570 %
Multi limit cycle
Ω ω 2 ,  ω 1 5 ω 2
20 %
570 %
Multi limit cycle
Ω ω 2 ,   ω 1 ω 2 2
720 %
930 %
Multi limit cycle
Ω ω 2 ,   ω 1 2 ω 2 3
670 %
860 %
Multi limit cycle
Ω ω 2 ,   ω 1 3 ω 2 2
100 %
530 %
Multi limit cycle
Ω 2 ω 1 ,   ω 2 ω 1
670 %
290 %
Multi limit cycle
Ω 2 ω 1 ,   ω 2 4 ω 1
530 %
170 %
Multi limit cycle
Ω 3 ω 1 ,   ω 2 ω 1
430 %
290 %
Multi limit cycle
Ω 3 ω 1 ,   ω 2 3 ω 1
730 %
1190 %
Multi limit cycle

Fig. 3. Simultaneous primary and internal resonance case (Ωω1 and ω23ω1)

Simultaneous primary and internal resonance case (Ω≅ω1 and ω2≅3ω1)
Simultaneous primary and internal resonance case (Ω≅ω1 and ω2≅3ω1)
Simultaneous primary and internal resonance case (Ω≅ω1 and ω2≅3ω1)
Simultaneous primary and internal resonance case (Ω≅ω1 and ω2≅3ω1)

Fig. 4. Simultaneous primary combined and internal resonance case (Ωω1 and ω2ω1)

Simultaneous primary combined and internal resonance case (Ω≅ω1 and ω2≅ω1)
Simultaneous primary combined and internal resonance case (Ω≅ω1 and ω2≅ω1)
Simultaneous primary combined and internal resonance case (Ω≅ω1 and ω2≅ω1)
Simultaneous primary combined and internal resonance case (Ω≅ω1 and ω2≅ω1)

5.1.1. Response curves of the first mode

Fig. 5(a) shows the steady state amplitudes of the first mode against the detuning parameter σ1 at the practical case, where a0, b0. For increasing values of the nonlinear parameter α1 and α17, the curves are shifted to the left and the regions of stability are increased as shown in Figs. 5(b) and 5(j). Fig. 5(c) shows that the steady state amplitude of the first mode is a monotonic decreasing function in the linear damping coefficients α2. Fig. 5(d) shows that the steady state amplitude of the first mode is a monotonic increasing function in the nonlinear parameter α12 and the response curves are bent to the left leading to the occurrence of the jump phenomena, multi-valued amplitudes and the regions of stability are decreased. The steady state amplitude of the first mode is a monotonic increasing function in the nonlinear parameter α14 and the excitation amplitudes ρ and the regions of instability are increased as shown in Figs. 5 (e) and 5(o). For increasing values of the nonlinear parameter α15, α16 and α20 the steady state amplitude of the first mode is increased and the curves are diverged as shown in Figs. 5(f), 5(h), 5(m). For negative values of the nonlinear parameter α15 and α18, the steady state amplitude of the first mode is increased and the stability regions are increased as shown in Figs. 5(g) and 5(k). Figs. 5(i), 5(l) show that, for positive values of the nonlinear parameter α16 and α18, the system is unstable. Fig. 5(n) shows that the steady state amplitude of the first mode is a monotonic decreasing function in the natural frequencies ω1 and ω2.

Fig. 5. Effects of different parameters on the first mode

 Effects of different parameters on the first mode

a) Effects of the detuning parameter σ1

 Effects of different parameters on the first mode

b) Effects of the non-linear parameter α1

 Effects of different parameters on the first mode

c) Effects of the damping coefficient α2

 Effects of different parameters on the first mode

d) Effects of the non-linear parameter α12

 Effects of different parameters on the first mode

e) Effects of the non-linear parameter α14

 Effects of different parameters on the first mode

f) Effects of the non-linear parameter α15

 Effects of different parameters on the first mode

g) Effects of the non-linear parameter α15

 Effects of different parameters on the first mode

h) Effects of the non-linear parameter α16

 Effects of different parameters on the first mode

i) Effects of the non-linear parameter α16

 Effects of different parameters on the first mode

j) Effects of the non-linear parameter α17

 Effects of different parameters on the first mode

k) Effects of the non-linear parameter α18

 Effects of different parameters on the first mode

l) Effects of the non-linear parameter α18

 Effects of different parameters on the first mode

m) Effects of the non-linear parameter α20

 Effects of different parameters on the first mode

n) Effects of the natural frequencies ω1, ω2

 Effects of different parameters on the first mode

o) Effects of the excitation amplitude ρ

5.1.2. Response curves of the second mode

Fig. 6(a) shows the steady state amplitudes of the second mode against the detuning parameter σ2 at the practical case, where a0, b0. For increasing values of the nonlinear parameter β3 and β15, the curves are shifted to the right leading to the occurrence of the jump phenomena, multi-valued amplitudes as shown in Figs. 6(b) and 6(h). Figs. 6(c), 6(f), 6(k) show that the steady state amplitudes of the second mode are monotonic decreasing functions in the linear damping coefficients β4, the nonlinear parameter β14 and the natural frequencies ω1 and ω2. Figs. 6(d), 6(i), 6(j) shows that the steady state amplitudes of the second mode are monotonic increasing functions in the nonlinear parameters β12, β16 and β19. For negative and positive values of the nonlinear parameter β13, the steady state amplitude of the second mode is trivial leading to the occurrence of the saturation phenomena as shown in Fig. 6(e). For increasing positive values of the nonlinear parameter β14, the steady state amplitude of the second mode is increased and the system is unstable as shown in Fig. 6(g). The regions of stability are decreased for increasing values of the nonlinear parameters β3 and β12. The regions of stability are increased for increasing values of the nonlinear parameter β15 and the natural frequencies ω1 and ω2.

Fig. 6. Effects of different parameters on the second mode

 Effects of different parameters on the second mode

a) Effects of the detuning parameter σ2

 Effects of different parameters on the second mode

b) Effects of the non-linear parameter β3

 Effects of different parameters on the second mode

c) Effects of the damping coefficient β4

 Effects of different parameters on the second mode

d) Effects of the non-linear parameter β12

 Effects of different parameters on the second mode

e) Effects of the non-linear parameter β13

 Effects of different parameters on the second mode

f) Effects of the non-linear parameter β14

 Effects of different parameters on the second mode

g) Effects of the non-linear parameter β14

 Effects of different parameters on the second mode

h) Effects of the non-linear parameter β15

 Effects of different parameters on the second mode

i) Effects of the non-linear parameter β16

 Effects of different parameters on the second mode

j) Effects of the non-linear parameter β19

 Effects of different parameters on the second mode

k) Effects of the natural frequencies ω1, ω2

Fig. 7. a), b) Time responses and phase plane of the coupled system respectively; c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=1)

 a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=1)

a)

 a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=1)

b)

 a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=1)

c)

Fig. 8. a), b) Time responses and phase plane of the coupled system respectively; c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=0.4)

a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=0.4)

a)

a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=0.4)

b)

a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=3, ω1=2.8, ω2=3.25 and ρ=0.4)

c)

Fig. 9. a), b) Time responses and phase plane of the coupled system respectively; c) Poincare maps of the coupled system respectively; for (Ω=8.4, ω1=2.8, ω2=8.4 and ρ=1)

a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=8.4, ω1=2.8, ω2=8.4 and ρ=1)

a)

a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=8.4, ω1=2.8, ω2=8.4 and ρ=1)

b)

a), b) Time responses and phase plane of the coupled system respectively;  c) Poincare maps of the coupled system respectively; for (Ω=8.4, ω1=2.8, ω2=8.4 and ρ=1)

c)

6. Poincare map

Poincare maps are introduced via example using two-dimensional autonomous systems of differential equations. They are used extensively to transform complicated behavior in the phase space to discrete maps in a lower-dimensional space. Poincaré maps are used to determine stability and plot bifurcation diagrams. By simulating the movement of rotor-seal system (2) and Poincare mapping, the Fig. 7-9 are obtained. At ρ=1 with the non-resonant rotor-seal system where Ωω1ω2, there is a period-one harmonic solution, which is depicted as a closed curve in the phase plane and as and there exist one isolated points on Poincare map as shown in Fig. 7. When ρ=0.4 with the non-resonant rotor-seal system where Ωω1ω2, the system becomes chaotic and quasi-periodic motion appears and a close curve is observed on Poincare map as shown in Fig. 8. Fig. 9 shows that the time response of the rotor-seal system become stable and quasi-periodic motion appears at ρ=1 with the simultaneous sub-harmonic and internal resonance case (Ω3ω1, ω23ω1).

7. Conclusions

Multiple time scale perturbation method is useful to determine approximate solutions for the coupled nonlinear differential equations describing the rotor-seal system up to and including the second order approximation. It is quite clear that some of the simultaneous resonance cases are undesirable in the design of such system as they represent some of the worst behavior of the system. Both the frequency response equations and the phase-plane technique are applied to study the stability of the system. The effect of the different parameters of the system is studied numerically. From the above study the following may be concluded:

1) The amplitude of the first mode is increased to about 30 % of the maximum excitation forces amplitude ρ, while the amplitude of the second mode is increased to about 70 % at the non-resonant case (Ωω1ω2).

2) The amplitude of the first mode is increased to about 670 % while the amplitude of the second mode is increased to about 1140 % at resonance case (Ωω1 and ω23ω1).

3) The amplitude of the first mode is increased to about 670 % while the amplitude of the second mode is increased to about 570 % at resonance case (Ωω1 and ω2ω1).

4) The regions of stability are increased for increasing values of α1, α14 and α17.

5) For positive values of the nonlinear parameter α16 and α18, the system become unstable.

6) For increasing values of the nonlinear parameter β3 and β15, the curves are shifted to the right leading to the occurrence of the jump phenomena, multi-valued.

7) For negative and positive values of the nonlinear parameter β13, the steady state amplitude of the second mode is trivial leading to the occurrence of the saturation phenomena.

Acknowledgements

The authors would like to express their gratitude to the Editor and Referees for their encouragement and constructive comments in revising the paper.

References

  1. Dimarogonas A. D., Gomez-Mancilla J. C. Flow-excited turbine rotor instability. International Journal of Rotating Machinery, Vol. 1, Issue 1, 1994, p. 3-51. [Search CrossRef]
  2. Marquette O. R., Childs D. W., Andres L. S. Eccentricity effects on the rotor dynamic coefficients of plain annular seals theory versus experiment. ASME Journal of Tribology, Vol. 119, Issue 3, 1997, p. 443-447. [Search CrossRef]
  3. Klaus K. Dynamic coefficients of stepped labyrinth gas seals. Journal of Engineering for Gas Turbines and Power, Vol. 122, Issue 3, 2000, p. 473-477. [Search CrossRef]
  4. Dietzed F. J., Nordmann R. Calculating coefficients of seals by finite-difference techniques. SME Journal of Tribology, Vol. 109, Issue 3, 1987, p. 388-394. [Search CrossRef]
  5. Toshio H., Guo Z. L., Gordon K. R. Application of fluid dynamics analysis for rotating machinery – part II: labyrinth seal analysis. Journal of Engineering for Gas Turbines and Power, Vol. 127, Issue 4, 2005, p. 820-826. [Search CrossRef]
  6. Alford J. S. Protecting turbomachinery from self-excited rotor whirl. ASME Journal of Engineering for Power, Vol. 87, Issue 4, 1965, p. 333-344. [Search CrossRef]
  7. Ding Q., Zhang K. Order reduction and nonlinear behaviors of a continuous rotor system. Nonlinear Dynamics, Vol. 67, 2012, p. 251-262. [Search CrossRef]
  8. Chen Y. M., Meng G., Liu J. K. A new method for Fourier series expansions: applications in rotor-seal systems. Mechanics Research Communications, Vol. 38, 2011, p. 399-403. [Search CrossRef]
  9. Li W., Yang Y., Sheng D., Chen J. A novel nonlinear model of rotor/bearing/seal system and numerical analysis. Mechanism and Machine Theory, Vol. 46, 2011, p. 618-631. [Search CrossRef]
  10. Zhang H., Chen Y. Bifurcation analysis on full annular rub of a nonlinear rotor system. Science China, Technological Sciences, Vol. 54, Issue 8, 2011, p. 1977-1985. [Search CrossRef]
  11. Cheng M., Meng G., Jing J. Non-linear dynamics of a rotor-bearing-seal system. Archive of Applied Mechanics, Vol. 76, 2006, p. 215-227. [Search CrossRef]
  12. Ding Q., Cooper J. E., Leung A. Y. T. Hopf bifurcation analysis of a rotor-seal system. Journal of Sound and Vibration, Vol. 252, Issue 5, 2002, p. 817-833. [Search CrossRef]
  13. Li Z., Chen Y. Research on 1:2 subharmonic resonance and bifurcation of nonlinear rotor-seal system. Applied Mathematics and Mechanics, Vol. 33, Issue 4, 2012, p. 499-510. [Search CrossRef]
  14. Sayed M., Mousa A. A. Second-order approximation of angle-ply composite laminated thin plate under combined excitations. Communication in Nonlinear Science and Numerical Simulation, Vol. 17, 2012, p. 5201-5216. [Search CrossRef]
  15. Eissa M., Sayed M. A comparison between passive and active control of non-linear simple pendulum, part-I. Mathematical and Computational Applications, Vol. 11, 2006, p. 137-149. [Search CrossRef]
  16. Eissa M., Sayed M. A comparison between passive and active control of non-linear simple pendulum, part-II. Mathematical and Computational Applications, Vol. 11, 2006, p. 151-162. [Search CrossRef]
  17. Eissa M., Sayed M. Vibration reduction of a three DOF non-linear spring pendulum. Communication in Nonlinear Science and Numerical Simulation, Vol. 13, 2008, p. 465-488. [Search CrossRef]
  18. Sayed M. Improving the mathematical solutions of nonlinear differential equations using different control methods. Ph. D. Thesis, Menofia University, Egypt, 2006. [Search CrossRef]
  19. Amer Y. A., Bauomy H. S., Sayed M. Vibration suppression in a twin-tail system to parametric and external excitations. Computers and Mathematics with Applications, Vol. 58, 2009, p. 1947-1964. [Search CrossRef]
  20. Sayed M., Hamed Y. S. Stability and response of a nonlinear coupled pitch-roll ship model under parametric and harmonic excitations. Nonlinear Dynamics, Vol. 64, 2011, p. 207-220. [Search CrossRef]
  21. Hamed Y. S., El-Ganaini W. A., Kamel M. M. Vibration suppression in ultrasonic machining described by non-linear differential equations. Journal of Mechanical Science and Technology, Vol. 23, Issue 8, 2009, p. 2038-2050. [Search CrossRef]
  22. Hamed Y. S., El-Ganaini W. A., Kamel M. M. Vibration suppression in multi-tool ultrasonic machining to multi-external and parametric excitations. Acta Mechanica Sinica, Vol. 25, 2009, p. 403-415. [Search CrossRef]
  23. Hamed Y. S., El-Ganaini W. A., Kamel M. M. Vibration reduction in ultrasonic machine to external and tuned excitation forces. Applied Mathematical Modeling, Vol. 33, 2009, p. 2853-2863. [Search CrossRef]
  24. Hamed Y. S., Sayed M., Cao D.-X., Zhang W. Nonlinear study of the dynamic behavior of a string-beam coupled system under combined excitations. Acta Mechanica Sinica, Vol. 27, Issue 6, 2011, p. 1034-1051. [Search CrossRef]
  25. Sayed M., Kamel M. Stability study and control of helicopter blade flapping vibrations. Applied Mathematical Modelling, Vol. 35, 2011, p. 2820-2837. [Search CrossRef]
  26. Sayed M., Kamel M. 1:2 and 1:3 internal resonance active absorber for non-linear vibrating system. Applied Mathematical Modelling, Vol. 36, 2012, p. 310-332. [Search CrossRef]
  27. Sayed M., Hamed Y. S., Amer Y. A. Vibration reduction and stability of non-linear system subjected to external and parametric excitation forces under a non-linear absorber. International Journal of Contemporary Mathematical Sciences, Vol. 6, Issue 22, 2011, p. 1051-1070. [Search CrossRef]
  28. Sayed M., Mousa A. A. Vibration, stability, and resonance of angle-ply composite laminated rectangular thin plate under multi-excitations. Mathematical Problems in Engineering, 2013. [Search CrossRef]
  29. Poincaré H. Memory on the curves defined by a differential equation. Journal de Mathématiques Pures et Appliquées, Vol. 7, 1881, p. 375-422. [Search CrossRef]
  30. Nayfeh A. H. Introduction to Perturbation Techniques. John Wiley & Sons, Inc., New York, 1981. [Search CrossRef]
  31. Nayfeh A. H. Non-linear Interactions. Wiley-Inter-Science, New York, 2000. [Search CrossRef]