# 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.
Views 47
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 .

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  first derived the formula of the gas exciting force when he researched on the stability of aeroengine. Ding and Zhang  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.  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.  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  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.  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.  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  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  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  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.  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  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.  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.  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  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é  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 :

(1)
$\left(\begin{array}{cc}M+m& 0\\ 0& M+m\end{array}\right)\left(\begin{array}{c}\stackrel{¨}{x}\\ \stackrel{¨}{y}\end{array}\right)+\left(\begin{array}{cc}{D}_{xe}+D& 2m\lambda \mathrm{\Omega }\\ -2m\lambda \Omega & {D}_{ye}+D\end{array}\right)\left(\begin{array}{c}\stackrel{˙}{x}\\ \stackrel{˙}{y}\end{array}\right)$

where $M$ is the rotor mass, $m$ is the gas inertia effects, ${K}_{xe}$, ${K}_{ye}$, ${D}_{xe}$, ${D}_{ye}$ are the stiffness and external damping on the $x$ and $y$ directions respectively and $G$ is the gravitation effects of the rotor, $\rho$ is the amplitude of the excitation force.

Fig. 1. Model of rotor-seal system a) b)

#### 3. Mathematical analysis

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

(2a)
$\stackrel{¨}{x}+{\omega }_{1}^{2}x={f}_{x}\left(x,y,\stackrel{˙}{x},\stackrel{˙}{y}\right)+\rho \mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t\right),$
(2b)
$\stackrel{¨}{y}+{\omega }_{2}^{2}y={f}_{y}\left(x,y,\stackrel{˙}{x},\stackrel{˙}{y}\right)+\rho \mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\Omega }t\right),$

with initial condition $x\left(0\right)=\stackrel{˙}{x}\left(0\right)=y\left(0\right)=\stackrel{˙}{y}\left(0\right)=0$. Where $x$ and $y$ are the displacements of the shaft center, ${\omega }_{1}$ and ${\omega }_{2}$ the linear natural frequencies in $x$, $y$ directions, and $\mathrm{\Omega }$ the angular speed of the rotor, $\rho$ is the amplitude of the excitation force, ${f}_{x}\left(x,y,\stackrel{˙}{x},\stackrel{˙}{y}\right)$ and ${f}_{y}\left(x,y,\stackrel{˙}{x},\stackrel{˙}{y}\right)$ are power series function of the nonlinear force on $x$ and $y$ directions given by:

${f}_{x}\left(x,y,\stackrel{˙}{x},\stackrel{˙}{y}\right)={\alpha }_{1}x+{\alpha }_{2}\stackrel{˙}{x}+{\alpha }_{3}y+{\alpha }_{4}\stackrel{˙}{y}+{\alpha }_{5}{x}^{2}+{\alpha }_{6}{y}^{2}+{\alpha }_{7}x\stackrel{˙}{x}+{\alpha }_{8}xy+{\alpha }_{9}x\stackrel{˙}{y}$

${f}_{y}\left(x,y,\stackrel{˙}{x},\stackrel{˙}{y}\right)={\beta }_{1}x+{\beta }_{2}\stackrel{˙}{x}+{\beta }_{3}y+{\beta }_{4}\stackrel{˙}{y}+{\beta }_{5}{x}^{2}+{\beta }_{6}{y}^{2}+{\beta }_{7}x\stackrel{˙}{x}+{\beta }_{8}xy+{\beta }_{9}x\stackrel{˙}{y}$

.

All coefficient of ${\alpha }_{i}$, ${\beta }_{i}$ are defined in . The linear, non-linear coefficients and excitation amplitude are assumed to be:

(3)

where $\epsilon$ is a small perturbation parameter and $\text{0}<\epsilon <<\text{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)

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

(5)
$x\left(t,\text{\hspace{0.17em}}\epsilon \right)={x}_{0}\left({T}_{0},\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)+\epsilon \text{\hspace{0.17em}}{x}_{1}\left({T}_{0},\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)\text{\hspace{0.17em}}+{\epsilon }^{2}\text{\hspace{0.17em}}{x}_{2}\left({T}_{0},\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)\text{\hspace{0.17em}}+O\left({\epsilon }^{3}\right),$
(6)
$y\left(t,\text{\hspace{0.17em}}\epsilon \right)=\text{\hspace{0.17em}}{y}_{0}\left({T}_{0},\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)+\epsilon \text{\hspace{0.17em}}{y}_{1}\left({T}_{0},\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)\text{\hspace{0.17em}}+{\epsilon }^{2}\text{\hspace{0.17em}}{y}_{2}\left({T}_{0},\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)\text{\hspace{0.17em}}+O\left({\epsilon }^{3}\right).$

Derivatives with respect to $t$ are then transformed into:

(7)
$\frac{d}{dt}=\frac{\partial }{\partial {T}_{0}}\frac{\partial {T}_{0}}{\partial t}+\frac{\partial }{\partial {T}_{1}}\frac{\partial {T}_{1}}{\partial t}+\frac{\partial }{\partial {T}_{2}}\frac{\partial {T}_{2}}{\partial t}={D}_{0}+\epsilon \text{\hspace{0.17em}}{D}_{1}+{\epsilon }^{2}\text{\hspace{0.17em}}{D}_{2}\text{\hspace{0.17em}},$
(8)
$\frac{{d}^{2}}{d{t}^{2}}={D}_{0}^{2}+2\epsilon {D}_{0}\text{\hspace{0.17em}}{D}_{1}+{\epsilon }^{2}\left({D}_{1}^{2}+2{D}_{0}{D}_{2}\right)\text{\hspace{0.17em}},$

where ${D}_{n}=\partial /\partial {T}_{n}$, $n=$ 0, 1, 2. Terms of $O\left({\epsilon }^{3}\right)$ and higher orders are neglected. Substituting Eqs. (5)-(6) and (7)-(8) into Eqs. (2) and equating the coefficients of similar powers of $\epsilon$, one obtains the following set of ordinary differential equations:

Order ${\epsilon }^{0}$:

(9)
$\left({D}_{0}^{2}+{\omega }_{1}^{2}\right)\text{\hspace{0.17em}}{x}_{0}=0,$
(10)
$\left({D}_{0}^{2}+{\omega }_{2}^{2}\right)\text{\hspace{0.17em}}{y}_{0}=0.$

Order ${\epsilon }^{1}$:

(11)
$\left({D}_{0}^{2}+{\omega }_{1}^{2}\right)\text{\hspace{0.17em}}{x}_{1}=-2{D}_{0}{D}_{1}{x}_{0}+{\stackrel{^}{\alpha }}_{1}{x}_{0}+{\stackrel{^}{\alpha }}_{2}{D}_{0}{x}_{0}+{\stackrel{^}{\alpha }}_{3}{y}_{0}+{\stackrel{^}{\alpha }}_{4}{D}_{0}{y}_{0}+{\stackrel{^}{\alpha }}_{5}{x}_{0}^{2}+{\stackrel{^}{\alpha }}_{6}{y}_{0}^{2}$

(12)
$\left({D}_{0}^{2}+{\omega }_{2}^{2}\right)\text{\hspace{0.17em}}{y}_{1}=-2{D}_{0}{D}_{1}{y}_{0}+{\stackrel{^}{\beta }}_{1}{x}_{0}+{\stackrel{^}{\beta }}_{2}{D}_{0}{x}_{0}+{\stackrel{^}{\beta }}_{3}{y}_{0}+{\stackrel{^}{\beta }}_{4}{D}_{0}{y}_{0}+{\stackrel{^}{\beta }}_{5}{x}_{0}^{2}+{\stackrel{^}{\beta }}_{6}{y}_{0}^{2}$

Order ${\epsilon }^{2}$:

(13)
$\left({D}_{0}^{2}+{\omega }_{1}^{2}\right)\text{\hspace{0.17em}}{x}_{2}=-{D}_{1}^{2}{x}_{1}-2{D}_{0}{D}_{2}{x}_{0}-2{D}_{0}{D}_{1}{x}_{1}+{\stackrel{^}{\alpha }}_{1}{x}_{1}+{\stackrel{^}{\alpha }}_{2}\left({D}_{0}{x}_{1}+{D}_{1}{x}_{0}\right)+{\stackrel{^}{\alpha }}_{3}{y}_{1}$

(14)
$\left({D}_{0}^{2}+{\omega }_{2}^{2}\right)\text{\hspace{0.17em}}{y}_{2}=-{D}_{1}^{2}{y}_{1}-2{D}_{0}{D}_{2}{y}_{0}-2{D}_{0}{D}_{1}{y}_{1}+{\stackrel{^}{\beta }}_{1}{x}_{1}+{\stackrel{^}{\beta }}_{2}\left({D}_{0}{x}_{1}+{D}_{1}{x}_{0}\right)+{\stackrel{^}{\beta }}_{3}{y}_{1}$

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

(15)
${x}_{0}={A}_{1}\left(\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)\text{\hspace{0.17em}}\mathrm{e}\mathrm{x}\mathrm{p}\left(i{\omega }_{1\text{\hspace{0.17em}}}{T}_{0}\right)+cc,$
(16)
${y}_{0}={A}_{2}\left(\text{\hspace{0.17em}}{T}_{1},\text{\hspace{0.17em}}{T}_{2}\right)\text{\hspace{0.17em}}\mathrm{e}\mathrm{x}\mathrm{p}\left(i{\omega }_{2\text{\hspace{0.17em}}}{T}_{0}\right)+cc,$

where ${A}_{1}$ and ${A}_{2}$ are a complex function in ${T}_{1}$, ${T}_{2}$ 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)
$\text{\hspace{0.17em}}{x}_{1}={E}_{1}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(2i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{E}_{2}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(3i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{E}_{3}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(i{\omega }_{2}\text{\hspace{0.17em}}{T}_{0}\right)+{E}_{4}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(2i{\omega }_{2}\text{\hspace{0.17em}}{T}_{0}\right)$

(18)
$\text{\hspace{0.17em}}{y}_{1}={E}_{11}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{E}_{12}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(2i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{E}_{13}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(3i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{E}_{14}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(2i{\omega }_{2}\text{\hspace{0.17em}}{T}_{0}\right)$

where ${E}_{i}$ () are complex functions in ${T}_{1}$ and ${T}_{2}$. Substituting Eqs. (15)-(18) into Eqs. (13)-(14) and after eliminating the secular terms, the non-homogeneous solutions of Eqs. (13)-(14) are:

(19)

(20)
$\text{\hspace{0.17em}}{y}_{2}\text{\hspace{0.17em}}={G}_{1}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{G}_{2}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(2i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{G}_{3}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(3i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)+{G}_{4}\mathrm{e}\mathrm{x}\mathrm{p}\text{\hspace{0.17em}}\left(4i{\omega }_{1}\text{\hspace{0.17em}}{T}_{0}\right)$

where ${H}_{i}$, ${G}_{i}$ () are complex functions in ${T}_{1}$ and ${T}_{2}$. From the above derived solutions, the reported resonance cases are:

• Primary resonance:, $n=\text{1,\hspace{0.17em}2}$.

• Sub-harmonic resonance: $\mathrm{\Omega }\cong {k\omega }_{n}$, 2, 3, 4, 5 and $n=\text{1,\hspace{0.17em}2}$.

• Super-harmonic resonance: $k\mathrm{\Omega }\cong {\omega }_{n}$, $3\mathrm{\Omega }\cong {2\omega }_{n}$, $2\mathrm{\Omega }\cong 3{\omega }_{n}$, 2, 3, 4, 5 and $n=\text{1,\hspace{0.17em}2}$.

Internal or secondary resonance: ${\omega }_{1}\cong s\text{\hspace{0.17em}}{\omega }_{2}$, ${\omega }_{2}\cong s\text{\hspace{0.17em}}{\omega }_{1}$, $2{\omega }_{1}\cong 3\text{\hspace{0.17em}}{\omega }_{2}$, 2, 3, 4, 5.

• Combined resonance: ${\omega }_{1}\cong ±m\text{\hspace{0.17em}}{\omega }_{2}±\mathrm{\Omega }\text{,}$${\omega }_{1}\cong ±\text{\hspace{0.17em}}{\omega }_{2}±r\text{\hspace{0.17em}}\mathrm{\Omega }\text{,}$${\omega }_{1}\cong ±2{\omega }_{2}±2\mathrm{\Omega }\text{,}$${\omega }_{1}\cong ±2{\omega }_{2}±3\mathrm{\Omega }\text{,}$$2{\omega }_{1}\cong ±q\text{\hspace{0.17em}}{\omega }_{2}±\mathrm{\Omega }\text{,}$$2{\omega }_{1}\cong ±\text{\hspace{0.17em}}{\omega }_{2}±p\text{\hspace{0.17em}}\mathrm{\Omega }\text{,}$$3{\omega }_{1}\cong ±\text{\hspace{0.17em}}n\text{\hspace{0.17em}}{\omega }_{2}±\text{\hspace{0.17em}}\mathrm{\Omega }\text{,}$$3{\omega }_{1}\cong ±\text{\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}±2\text{\hspace{0.17em}}\mathrm{\Omega }\text{,}$$4{\omega }_{1}\cong ±\text{\hspace{0.17em}}{\omega }_{2}±\text{\hspace{0.17em}}\mathrm{\Omega }\text{,}$1, 2, 3, 4, $r=$2, 3, 4, 1, 2, 3, 2, 3, $n=\text{1,\hspace{0.17em}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 ${\sigma }_{1}$ and ${\sigma }_{2}$ according to:

(21)

Substituting Eq. (21) into Eqs. (11)-(12) and eliminating the secular and small divisor terms from ${x}_{1}$ and ${y}_{1}$, we get the following:

(22)
$2i{\omega }_{1}{D}_{1}A=\left({\stackrel{^}{\alpha }}_{1}+i{\stackrel{^}{\alpha }}_{2}{\omega }_{1}\right)A+\left[2i{\stackrel{^}{\alpha }}_{18}{\omega }_{1}+2{\stackrel{^}{\alpha }}_{17}\right]B\stackrel{-}{B}A\text{\hspace{0.17em}}$
(23)
$2i{\omega }_{2}{D}_{1}B=\left({\stackrel{^}{\beta }}_{3}+i{\stackrel{^}{\beta }}_{4}{\omega }_{2}\right)B+\left[3{\stackrel{^}{\beta }}_{13}+i{\stackrel{^}{\beta }}_{19}{\omega }_{2}\right]{B}^{2}\stackrel{-}{B}$

To analyze the solution of Eqs. (22)-(23), it is convenient to express the $A\left({T}_{1}\right)$ and $B\left({T}_{1}\right)$ in the polar form as:

(24)

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

(25)
${a}^{\text{'}}=\frac{{\alpha }_{2}}{2}a+\frac{{\alpha }_{14}}{8}{a}^{3}+\frac{{\alpha }_{18}}{4}a\text{\hspace{0.17em}}{b}^{2}+\frac{{\alpha }_{15}}{8{\omega }_{1}}{a}^{2}b\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2}+\left(\frac{{\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}}{8{\omega }_{1}}\right){a}^{2}b\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2}$
(26)
$a\text{\hspace{0.17em}}{\gamma }_{1}^{\text{'}}=-\frac{{\alpha }_{1}}{2{\omega }_{1}}a-\frac{3{\alpha }_{12}}{8{\omega }_{1}}{a}^{3}-\frac{{\alpha }_{17}}{4{\omega }_{1}}a{b}^{2}-\frac{{\alpha }_{15}}{8{\omega }_{1}}{a}^{2}b\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2}+\left(\frac{{\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}}{8{\omega }_{1}}\right){a}^{2}b\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2}$
(27)
${b}^{\text{'}}=\frac{{\beta }_{4}}{2}b+\frac{{\beta }_{19}}{8}{b}^{3}+\frac{{\beta }_{16}}{4}{a}^{2}b-\frac{{\beta }_{12}}{8{\omega }_{2}}{a}^{3}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2}+\frac{{\beta }_{14}{\omega }_{1}}{8{\omega }_{2}}{a}^{3}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2},$
(28)
$b{\gamma }_{2}^{\text{'}}=-\frac{{\beta }_{3}}{2{\omega }_{2}}b-\frac{3{\beta }_{13}}{8{\omega }_{2}}{b}^{3}-\frac{{\beta }_{15}}{4{\omega }_{2}}{a}^{2}b-\frac{{\beta }_{12}}{8{\omega }_{2}}{a}^{3}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2}-\frac{{\beta }_{14}{\omega }_{1}}{8{\omega }_{2}}{a}^{3}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2},$

where:

(29)

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

(30)
${\gamma }_{1}^{\text{'}}={\sigma }_{1}-{\theta }_{1}^{\text{'}},$
(31)
${\gamma }_{2}^{\text{'}}=3{\sigma }_{1}-{\sigma }_{2}+{\theta }_{2}^{\text{'}}-3{\theta }_{1}^{\text{'}}.$

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

(32)
$\frac{{\alpha }_{2}}{2}a+\frac{{\alpha }_{14}}{8}{a}^{3}+\frac{{\alpha }_{18}}{4}a\text{\hspace{0.17em}}{b}^{2}+\frac{{\alpha }_{15}}{8{\omega }_{1}}{a}^{2}b\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2}$
$+\left(\frac{{\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}}{8{\omega }_{1}}\right){a}^{2}b\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2}+\frac{\rho }{2{\omega }_{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{1}=0,$
(33)
$a\text{\hspace{0.17em}}{\sigma }_{1}+\frac{{\alpha }_{1}}{2{\omega }_{1}}a+\frac{3{\alpha }_{12}}{8{\omega }_{1}}{a}^{3}+\frac{{\alpha }_{17}}{4{\omega }_{1}}a{b}^{2}+\frac{{\alpha }_{15}}{8{\omega }_{1}}{a}^{2}b\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2}-\left(\frac{{\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}}{8{\omega }_{1}}\right){a}^{2}b\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2}$
(34)
$\frac{{\beta }_{4}}{2}b+\frac{{\beta }_{19}}{8}{b}^{3}+\frac{{\beta }_{16}}{4}{a}^{2}b-\frac{{\beta }_{12}}{8{\omega }_{2}}{a}^{3}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2}+\frac{{\beta }_{14}{\omega }_{1}}{8{\omega }_{2}}{a}^{3}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2}=0,$
(35)
$b\text{\hspace{0.17em}}\left(3{\sigma }_{1}-{\sigma }_{2}\right)+\frac{{\beta }_{3}}{2{\omega }_{2}}b+\frac{3{\beta }_{13}}{8{\omega }_{2}}{b}^{3}+\frac{{\beta }_{15}}{4{\omega }_{2}}{a}^{2}b+\frac{{\beta }_{12}}{8{\omega }_{2}}{a}^{3}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{2}+\frac{{\beta }_{14}{\omega }_{1}}{8{\omega }_{2}}{a}^{3}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{2}=0.$

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

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

(36)
${\left(\frac{{\alpha }_{2}}{2}a+\frac{{\alpha }_{14}}{8}{a}^{3}\right)}^{2}+{\left(a\text{\hspace{0.17em}}{\sigma }_{1}+\frac{{\alpha }_{1}}{2{\omega }_{1}}a+\frac{3{\alpha }_{12}}{8{\omega }_{1}}{a}^{3}\right)}^{2}-\frac{{\rho }^{2}}{4{\omega }_{1}^{2}}=0.$

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

(37)
${\left(\frac{{\alpha }_{2}}{2}a+\frac{{\alpha }_{14}}{8}{a}^{3}+\frac{{\alpha }_{18}}{4}a\text{\hspace{0.17em}}{b}^{2}\right)}^{2}+{\left(a\text{\hspace{0.17em}}{\sigma }_{1}+\frac{{\alpha }_{1}}{2{\omega }_{1}}a+\frac{3{\alpha }_{12}}{8{\omega }_{1}}{a}^{3}+\frac{{\alpha }_{17}}{4{\omega }_{1}}a\text{\hspace{0.17em}}{b}^{2}\right)}^{2}-{\left(\frac{{\alpha }_{15}}{8{\omega }_{1}}{a}^{2}b\right)}^{2}$
(38)
${\left(\frac{{\beta }_{4}}{2}b+\frac{{\beta }_{19}}{8}{b}^{3}+\frac{{\beta }_{16}}{4}{a}^{2}b\right)}^{2}+{\left(b\text{\hspace{0.17em}}\left(3{\sigma }_{1}-{\sigma }_{2}\right)+\frac{{\beta }_{3}}{2{\omega }_{2}}b+\frac{3{\beta }_{13}}{8{\omega }_{2}}{b}^{3}+\frac{{\beta }_{15}}{4{\omega }_{2}}{a}^{2}b\right)}^{2}$

#### 4.1. Non-linear solution

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

(39)

where ${a}_{0}$, ${b}_{0}$ and ${\theta }_{m0}$ are the solutions of Eqs. (32)-(35) and ${a}_{1}$, ${b}_{1}$, ${\theta }_{m1}$ are perturbations which are assumed to be small compared to ${a}_{0}$, ${b}_{0}$ and ${\theta }_{m0}$. Substituting Eq. (39) into Eqs. (25)-(28), using Eqs. (32)-(35) and keeping only the linear terms in ${a}_{1}$, ${b}_{1}$, ${\theta }_{m1}$ we obtain:

1) For the solution (0,0), we have:

(40)
${a}_{1}^{\text{'}}=\left(\frac{{\alpha }_{2}}{2}+\frac{3{\alpha }_{14}}{8}{a}_{0}^{2}\right){a}_{1}+\left(\frac{\rho }{2{\omega }_{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{10}\right){\theta }_{11},$
(41)
${\theta }_{11}^{\text{'}}=\left(\frac{{\sigma }_{1}}{{a}_{0}}+\frac{{\alpha }_{1}}{2{\omega }_{1}{a}_{0}}+\frac{9{\alpha }_{12}}{8{\omega }_{1}}{a}_{0}\right)\text{\hspace{0.17em}}{a}_{1}-\text{\hspace{0.17em}}\left(\frac{\rho }{2{a}_{0}{\omega }_{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{10}\right){\theta }_{11}.$

The stability of a given fixed point to a disturbance proportional to $\mathrm{e}\mathrm{x}\mathrm{p}\left(\lambda t\right)$ is determined by the roots of:

(42)
$A=\left[\begin{array}{cc}\frac{{\alpha }_{2}}{2}+\frac{3{\alpha }_{14}}{8}{a}_{0}^{2}& \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{\rho }{2{\omega }_{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{10}\\ \frac{{\sigma }_{1}}{{a}_{0}}+\frac{{\alpha }_{1}}{2{\omega }_{1}{a}_{0}}+\frac{9{\alpha }_{12}{a}_{0}}{8{\omega }_{1}}\text{\hspace{0.17em}}& \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}-\text{\hspace{0.17em}}\frac{\rho }{2{a}_{0}{\omega }_{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{10}\text{\hspace{0.17em}}\end{array}\right].$

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$, 0), the following are obtained:

(43)
${a}_{1}^{\text{'}}=\left(\frac{{\alpha }_{2}}{2}+\frac{3{\alpha }_{14}}{8}{a}_{0}^{2}+\frac{{\alpha }_{18}}{4}{b}_{0}^{2}+\frac{{\alpha }_{15}}{4{\omega }_{1}}{a}_{0}{b}_{0}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{20}+\frac{\left({\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}\right)}{4{\omega }_{1}}{a}_{0}{b}_{0}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{20}\right){a}_{1}$

(44)
${\theta }_{11}^{\text{'}}=\left(\begin{array}{c}\frac{{\sigma }_{1}}{{a}_{0}}+\frac{{\alpha }_{1}}{2{\omega }_{1}{a}_{0}}+\frac{9{\alpha }_{12}}{8{\omega }_{1}}{a}_{0}+\frac{{\alpha }_{17}}{4{\omega }_{1}{a}_{0}}{b}_{0}^{2}+\frac{{\alpha }_{15}}{4{\omega }_{1}}{b}_{0}cos{\theta }_{20}\\ -\frac{\left({\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}\right)}{4{\omega }_{1}}{b}_{0}sin{\theta }_{20}\end{array}\right)\text{\hspace{0.17em}}{a}_{1}$
$-\text{\hspace{0.17em}}\left(\frac{\rho }{2{a}_{0}{\omega }_{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{10}\right){\theta }_{11}+\left(\frac{{\alpha }_{17}}{2{\omega }_{1}}{b}_{0}+\frac{{\alpha }_{15}}{8{\omega }_{1}}{a}_{0}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{20}-\frac{\left({\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}\right)}{8{\omega }_{1}}{a}_{0}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{20}\right){b}_{1}$
$-\text{\hspace{0.17em}}\left(\frac{{\alpha }_{15}}{8{\omega }_{1}}{a}_{0}{b}_{0}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{20}-\frac{\left({\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}\right)}{8{\omega }_{1}}{a}_{0}{b}_{0}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{20}\right){\theta }_{21},$
(45)
${b}_{1}^{\text{'}}=\left(\frac{{\beta }_{16}}{2}{a}_{0}{b}_{0}-\frac{3{\beta }_{12}}{8{\omega }_{2}}{a}_{0}^{2}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{20}+\frac{3{\beta }_{14}{\omega }_{1}}{8{\omega }_{2}}{a}_{0}^{2}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{20}\right){a}_{1}+\left(\frac{{\beta }_{4}}{2}+\frac{3{\beta }_{19}}{8}{b}_{0}^{2}+\frac{{\beta }_{16}}{4}{a}_{0}^{2}\right){b}_{1}$
(46)
${\theta }_{20}^{\text{'}}=\left(\begin{array}{c}\frac{3{\sigma }_{1}}{{a}_{0}}+\frac{3{\alpha }_{1}}{2{\omega }_{1}{a}_{0}}+\frac{27{\alpha }_{12}}{8{\omega }_{1}}{a}_{0}+\frac{3{\alpha }_{17}}{4{\omega }_{1}{a}_{0}}{b}_{0}^{2}+\frac{3{\alpha }_{15}}{4{\omega }_{1}}{b}_{0}cos{\theta }_{20}\\ -\frac{3\left({\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}\right)}{4{\omega }_{1}}{b}_{0}sin{\theta }_{20}-\frac{{\beta }_{15}}{2{\omega }_{2}}{a}_{0}-\frac{3{\beta }_{12}}{8{\omega }_{2}{b}_{0}}{a}_{0}^{2}cos{\theta }_{20}-\frac{3{\beta }_{14}{\omega }_{1}}{8{\omega }_{2}{b}_{0}}{a}_{0}^{2}sin{\theta }_{20}\end{array}\right){a}_{1}$
$-\text{\hspace{0.17em}}\left(\frac{3\rho }{2{a}_{0}{\omega }_{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{10}\right){\theta }_{11}+\left(\begin{array}{c}\frac{3{\alpha }_{17}}{2{\omega }_{1}}{b}_{0}+\frac{3{\alpha }_{15}}{8{\omega }_{1}}{a}_{0}cos{\theta }_{20}-\frac{3\left({\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}\right)}{8{\omega }_{1}}{a}_{0}sin{\theta }_{20}\\ +\frac{\left({\sigma }_{2}-3{\sigma }_{1}\right)}{{b}_{0}}-\frac{{\beta }_{3}}{2{\omega }_{2}{b}_{0}}-\frac{9{\beta }_{13}}{8{\omega }_{2}}{b}_{0}-\frac{{\beta }_{15}}{4{\omega }_{2}{b}_{0}}{a}_{0}^{2}\end{array}\right){b}_{1}$
$+\left(\begin{array}{c}-\frac{3{\alpha }_{15}}{8{\omega }_{1}}{a}_{0}{b}_{0}sin{\theta }_{20}+\frac{3\left({\alpha }_{16}{\omega }_{2}-{\alpha }_{20}{\omega }_{1}\right)}{8{\omega }_{1}}{a}_{0}{b}_{0}cos{\theta }_{20}\\ +\frac{{\beta }_{12}}{8{\omega }_{2}{b}_{0}}{a}_{0}^{3}sin{\theta }_{20}-\frac{{\beta }_{14}{\omega }_{1}}{8{\omega }_{2}{b}_{0}}{a}_{0}^{3}cos{\theta }_{20}\end{array}\right){\theta }_{21}.$

The stability of a particular fixed point with respect to perturbations proportional to $\mathrm{e}\mathrm{x}\mathrm{p}\left(\lambda t\right)$ 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 $\mathrm{\Omega }\ne {\omega }_{1}\ne {\omega }_{2}$ at some practical values of the equation parameters ${\alpha }_{1}=\text{0.03,}$${\alpha }_{2}=\text{0.0001,}$–0.05, –0.2, ${\alpha }_{5}=\text{0.01}$, ${\alpha }_{6}=\text{0.02}$, ${\alpha }_{7}=\text{0.25}$, ${\alpha }_{8}=\text{0.004}$, ${\alpha }_{9}=\text{0.7}$, –0.8, –0.2, ${\alpha }_{12}=\text{0.004}$, ${\alpha }_{13}=\text{0.01}$, ${\alpha }_{14}=\text{0.02}$, ${\alpha }_{15}=\text{0.003}$, –0.04, ${\alpha }_{17}=\text{0.001,}$–0.5, –0.03, –0.04, ${\alpha }_{21}=\text{0.5}$, ${\beta }_{1}=\text{0.05}$, ${\beta }_{2}=\text{0.2}$, –0.03, ${\beta }_{4}=\text{0.0001}$, ${\beta }_{5}=\text{0.002}$, ${\beta }_{6}=\text{0.01}$, ${\beta }_{7}=\text{0.2}$, ${\beta }_{8}=\text{0.004}$, ${\beta }_{9}=\text{0.8}$, ${\beta }_{10}=\text{0.7}$, –0.25, –0.01, –0.004, ${\beta }_{14}=\text{0.03}$, –0.001, ${\beta }_{16}=\text{0.5}$, –0.003, ${\beta }_{18}=\text{0.04}$, –0.02, –0.5, ${\beta }_{21}=\text{0.04}$, $\text{\hspace{0.17em}}\mathrm{\Omega }=\text{3}$, , ${\omega }_{2}=\text{3.25}$, $\rho =\text{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)    For the simultaneous primary and internal resonance ($\mathrm{\Omega }\cong {\omega }_{1}$ and ${\omega }_{2}\cong 3{\omega }_{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 ($\mathrm{\Omega }\cong {\omega }_{1}$ and ${\omega }_{2}\cong {\omega }_{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 $\mathrm{\Omega }\cong {\omega }_{1}$ 700 % 620 % Multi limit cycle $\mathrm{\Omega }\cong \frac{{\omega }_{1}}{2}$ 560 % 570 % Multi limit cycle $\mathrm{\Omega }\cong \frac{{\omega }_{1}}{3}$ 560 % 580 % Multi limit cycle $\mathrm{\Omega }\cong \frac{2{\omega }_{1}}{3}$ 600 % 570 % Multi limit cycle $\mathrm{\Omega }\cong \frac{3{\omega }_{1}}{2}$ 500 % 530 % Multi limit cycle $\mathrm{\Omega }\cong \frac{{\omega }_{2}}{2}$ 560 % 570 % Multi limit cycle $\mathrm{\Omega }\cong \frac{{\omega }_{2}}{3}$ 560 % 580 % Multi limit cycle $\mathrm{\Omega }\cong \frac{2{\omega }_{2}}{3}$ 600 % 570 % Multi limit cycle $\mathrm{\Omega }\cong \frac{3{\omega }_{2}}{2}$ 500 % 500 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong {\omega }_{1}$ 670 % 570 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong 2{\omega }_{1}$ 730 % 1000 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong 3{\omega }_{1}$ 670 % 1140 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong 4{\omega }_{1}$ 370 % 200 % Multi limit cycle $\text{Ω}\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong 5{\omega }_{1}$ 800 % 240 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong \frac{{\omega }_{1}}{2}$ 85 % 320 % Multi limit cycle $\text{Ω}\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong \frac{{\omega }_{1}}{3}$ 400 % 530 % Multi limit cycle $\text{Ω}\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong \frac{{\omega }_{1}}{4}$ 200 % 430 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong \frac{2{\omega }_{1}}{3}$ 220 % 320 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{2}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{1}\cong {\omega }_{2}$ 830 % 430 % Multi limit cycle $\mathrm{\Omega }\cong \text{\hspace{0.17em}}{\omega }_{2}\text{,\hspace{0.17em}}{\omega }_{1}\cong 2{\omega }_{2}$ 85 % 430 % Multi limit cycle $\mathrm{\Omega }\cong \text{\hspace{0.17em}}{\omega }_{2}\text{,\hspace{0.17em}}{\omega }_{1}\cong 3{\omega }_{2}$ 50 % 500 % Multi limit cycle $\mathrm{\Omega }\cong \text{\hspace{0.17em}}{\omega }_{2}\text{,\hspace{0.17em}}{\omega }_{1}\cong 4{\omega }_{2}$ 35 % 570 % Multi limit cycle $\mathrm{\Omega }\cong \text{\hspace{0.17em}}{\omega }_{2}\text{,\hspace{0.17em}}{\omega }_{1}\cong 5{\omega }_{2}$ 20 % 570 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{2}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{1}\cong \frac{{\omega }_{2}}{2}$ 720 % 930 % Multi limit cycle $\mathrm{\Omega }\cong {\omega }_{2}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{1}\cong \frac{2{\omega }_{2}}{3}$ 670 % 860 % Multi limit cycle $\text{Ω}\cong {\omega }_{2}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{1}\cong \frac{3{\omega }_{2}}{2}$ 100 % 530 % Multi limit cycle $\text{Ω}\cong 2{\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong {\omega }_{1}$ 670 % 290 % Multi limit cycle $\mathrm{\Omega }\cong 2{\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong 4{\omega }_{1}$ 530 % 170 % Multi limit cycle $\mathrm{\Omega }\cong 3{\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong {\omega }_{1}$ 430 % 290 % Multi limit cycle $\mathrm{\Omega }\cong 3{\omega }_{1}\text{,\hspace{0.17em}\hspace{0.17em}}{\omega }_{2}\cong 3{\omega }_{1}$ 730 % 1190 % Multi limit cycle

Fig. 3. Simultaneous primary and internal resonance case ($\mathrm{\Omega }\cong {\omega }_{1}$ and ${\omega }_{2}\cong 3{\omega }_{1}$)    Fig. 4. Simultaneous primary combined and internal resonance case ($\mathrm{\Omega }\cong {\omega }_{1}$ and ${\omega }_{2}\cong {\omega }_{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 ${\sigma }_{1}$ at the practical case, where $a\ne \text{0}$, $b\ne \text{0}$. For increasing values of the nonlinear parameter ${\alpha }_{1}$ and ${\alpha }_{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 ${\alpha }_{2}$. Fig. 5(d) shows that the steady state amplitude of the first mode is a monotonic increasing function in the nonlinear parameter ${\alpha }_{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 ${\alpha }_{14}$ and the excitation amplitudes $\rho$ and the regions of instability are increased as shown in Figs. 5 (e) and 5(o). For increasing values of the nonlinear parameter ${\alpha }_{15}$, ${\alpha }_{16}$ and ${\alpha }_{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 ${\alpha }_{15}$ and ${\alpha }_{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 ${\alpha }_{16}$ and ${\alpha }_{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 ${\omega }_{1}$ and ${\omega }_{2}$.

Fig. 5. Effects of different parameters on the first mode a) Effects of the detuning parameter ${\sigma }_{1}$ b) Effects of the non-linear parameter ${\alpha }_{1}$ c) Effects of the damping coefficient ${\alpha }_{2}$ d) Effects of the non-linear parameter ${\alpha }_{12}$ e) Effects of the non-linear parameter ${\alpha }_{14}$ f) Effects of the non-linear parameter ${\alpha }_{15}$ g) Effects of the non-linear parameter ${\alpha }_{15}$ h) Effects of the non-linear parameter ${\alpha }_{16}$ i) Effects of the non-linear parameter ${\alpha }_{16}$ j) Effects of the non-linear parameter ${\alpha }_{17}$ k) Effects of the non-linear parameter ${\alpha }_{18}$ l) Effects of the non-linear parameter ${\alpha }_{18}$ m) Effects of the non-linear parameter ${\alpha }_{20}$ n) Effects of the natural frequencies ${\omega }_{1}$, ${\omega }_{2}$ o) Effects of the excitation amplitude $\rho$

#### 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 ${\sigma }_{2}$ at the practical case, where $\text{\hspace{0.17em}}a\ne 0$, $b\ne 0$. For increasing values of the nonlinear parameter ${\beta }_{\text{\hspace{0.17em}}3}$ and ${\beta }_{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 ${\beta }_{4}$, the nonlinear parameter ${\beta }_{\text{\hspace{0.17em}}14}$ and the natural frequencies ${\omega }_{1}$ and ${\omega }_{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 ${\beta }_{\text{\hspace{0.17em}}12}$, ${\beta }_{\text{\hspace{0.17em}}16}$ and ${\beta }_{\text{\hspace{0.17em}}19}$. For negative and positive values of the nonlinear parameter ${\beta }_{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 ${\beta }_{\text{\hspace{0.17em}}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 ${\beta }_{\text{\hspace{0.17em}}3}$ and ${\beta }_{\text{\hspace{0.17em}}12}$. The regions of stability are increased for increasing values of the nonlinear parameter ${\beta }_{\text{\hspace{0.17em}}15}$ and the natural frequencies ${\omega }_{1}$ and ${\omega }_{2}$.

Fig. 6. Effects of different parameters on the second mode a) Effects of the detuning parameter ${\sigma }_{2}$ b) Effects of the non-linear parameter ${\beta }_{3}$ c) Effects of the damping coefficient ${\beta }_{4}$ d) Effects of the non-linear parameter ${\beta }_{12}$ e) Effects of the non-linear parameter ${\beta }_{13}$ f) Effects of the non-linear parameter ${\beta }_{14}$ g) Effects of the non-linear parameter ${\beta }_{14}$ h) Effects of the non-linear parameter ${\beta }_{15}$ i) Effects of the non-linear parameter ${\beta }_{16}$ j) Effects of the non-linear parameter ${\beta }_{19}$ k) Effects of the natural frequencies ${\omega }_{1}$, ${\omega }_{2}$

Fig. 7. a), b) Time responses and phase plane of the coupled system respectively; c) Poincare maps of the coupled system respectively; for ($\mathrm{\Omega }=\text{3}$, ${\omega }_{1}=\text{2.8}$, ${\omega }_{2}=\text{3.25}$ and $\rho =\text{1}$) a) b) c)

Fig. 8. a), b) Time responses and phase plane of the coupled system respectively; c) Poincare maps of the coupled system respectively; for ($\mathrm{\Omega }=\text{3}$, ${\omega }_{1}=\text{2.8}$, ${\omega }_{2}=\text{3.25}$ and $\rho =\text{0.4}$) a) b) c)

Fig. 9. a), b) Time responses and phase plane of the coupled system respectively; c) Poincare maps of the coupled system respectively; for ($\mathrm{\Omega }=\text{8.4}$, ${\omega }_{1}=\text{2.8}$, ${\omega }_{2}=\text{8.4}$ and $\rho =\text{1}$) a) b) 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 $\rho =\text{1}$ with the non-resonant rotor-seal system where $\mathrm{\Omega }\ne {\omega }_{1}\ne {\omega }_{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 $\rho =\text{0.4}$ with the non-resonant rotor-seal system where $\mathrm{\Omega }\ne {\omega }_{1}\ne {\omega }_{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 $\rho =\text{1}$ with the simultaneous sub-harmonic and internal resonance case ($\mathrm{\Omega }\cong 3{\omega }_{1}$, ${\omega }_{2}\cong 3{\omega }_{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 $\rho$, while the amplitude of the second mode is increased to about 70 % at the non-resonant case ($\text{Ω}\ne {\omega }_{1}\ne {\omega }_{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 ($\mathrm{\Omega }\cong {\omega }_{1}$ and ${\omega }_{2}\cong 3{\omega }_{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 ($\mathrm{\Omega }\cong {\omega }_{1}$ and ${\omega }_{2}\cong {\omega }_{1}$).

4) The regions of stability are increased for increasing values of ${\alpha }_{1}$, ${\alpha }_{14}$ and ${\alpha }_{17}$.

5) For positive values of the nonlinear parameter ${\alpha }_{16}$ and ${\alpha }_{18}$, the system become unstable.

6) For increasing values of the nonlinear parameter ${\beta }_{\text{\hspace{0.17em}}3}$ and ${\beta }_{\text{\hspace{0.17em}}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 ${\beta }_{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.

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]