# Determination of the global responses characteristics of a piecewise smooth dynamical system with contact

Jun Jiang

► To cite this version:

Jun Jiang. Determination of the global responses characteristics of a piecewise smooth dynamical system with contact. Nonlinear Dynamics, 2009, 57 (3), pp.351-361. 10.1007/s11071-008-9446-z. hal-01433788

HAL Id: hal-01433788

https://hal.science/hal-01433788

Submitted on 13 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Determination of the global responses characteristics of a piecewise smooth dynamical system with contact

Jun Jiang

MOE Key Laboratory of Strength and Vibration, Xi'an
Jiaotong University, Xi'an, 710049, China
e-mail: jun.jiang@mail.xjtu.edu.cn

**Abstract** In this paper the global response characteristics of a piecewise smooth dynamical system with contact, which is specifically used to describe the rotor/stator rubbing systems, is studied analytically. A method to derive the global response characteristics of the model is proposed by studying each piece of the equations corresponding to different phases of the rotor motion, i.e., the phase without rubbing, the phase with rubbing and the phase of self-excited backward whirl. After solving the typical responses in each phase and deriving the corresponding existence boundaries in the parameter space, an overall picture of the global response characteristics of the model is obtained. As is shown, five types of the coexistences of the different rotor responses and deep insights into the interactive effect of parameters on the dynamic behavior of the model are gained.

**Keywords** Piecewise smooth systems · Global response characteristics · Stability analysis

## Nomenclature

A Amplitude of the synchronous no-rub motion

B Amplitude of synchronous full annular rub solution

C Real constant in the solution part of backward whirl motion

H Complex constant in the supposed solution

$R_0$ Non-dimensional clearance

R Non-dimensional deflection of the shaft center

$R_{\text{disk}}$ Non-dimensional radius of the disk

W Complex deflection of the rotor, $X + jY$

$X, Y$ Non-dimensional displacements of the shaft center

c Damping of the rotor

e Rotor mass eccentricity

$k_s, k_b$ Stiffness of the rotor shaft and the stator

m Mass of the rotor

$r_0$ Clearance between rotor and stator

r Displacement of the shaft geometric center

$r_{\text{disk}}$ Radius of the disk at the contact point

t Time

$x, y$ Horizontal and vertical displacements of the shaft center

$\alpha$ Real part of exponent of the solution part of backward whirl motion

$\beta$ Stiffness ratio of rotor-to-stator, $k_s/k_b$

$\varphi$ Phase angle

$\mu$ Coefficient of friction

$\tau$ Non-dimensional time

$\omega$ Rotating speed of the rotor

$\omega_0$ Backward whirl “natural” frequency of nonlinear coupled rotor/stator system

$\omega_2$ Natural frequency of the rotor system with zero clearance
$\omega_b$ Normalized whirl frequency of the rotor,
$\omega_b = \omega_w / \omega_2$

$\omega_w$ Whirl frequency of the rotor

$\Omega$ Normalized rotating speed of the rotor, $\omega/\omega_2$

$\Omega_e$ Critical rotating speed for the start of dry friction backward whirl

$\zeta$ Damping ratio of the rotor system

# 1 Introduction

Rotor-to-stator rubbing is a malfunction in rotating machinery that degrades the machine performance and may lead to a catastrophic failure of a whole machine through dry whip in the worst case. A large number of research papers deal with the dynamics of rub-related phenomena in rotor/stator contact systems. Some papers studied the possible rubbing responses induced through the rotor/stator rubbing and tried to derive the relationship between the rubbing responses and the system parameters, as summarized in [1]. Through the numerical and experimental studies, many rubbing behavior in rotor/stator contact systems have now become well known, for instance, the jump phenomenon and the synchronous full annular rubs [2–4], the partial rubs in sub- and super-synchronous whirl [5–7], the partial rubs in quasi-periodic whirl [8, 9], the chaotic motion [10, 11] as well as the dry friction backward whirl which includes two regimes: of dry whirl and dry whip [15–18].

The numerical and experimental studies are useful in detecting the possible rubbing responses in rotor/stator contact systems and in capturing some trends of dependence of the responses on system parameters. But a huge amount of work is inevitable in these ways in order to determine the boundaries of different rubbing responses of the rotor/stator contact systems. Furthermore, it was usually found that different researchers might present quite different response cascades and the inherent reason was not well explained. For instance, it was observed in tests of [12, 13] that the rotor under the unbalanced excitation would follow a normal behavior with the increase or decrease of the rotating speed, that is, the rotor follows the response cascade: no-rub motion → a synchronous full annular rub → no-rub motion, and vise versa. Jump phenomenon was detected in both run-up and run-down process. However, it was reported that when an additional outside disturbance (hitting on the shaft)

was applied on the rotor during the process, even at a very low rotor speed, the dry friction backward whirl with super-synchronous whirl frequency would be triggered from a synchronous response. A quite different response scenario with the increase of the rotor speed was observed in the tests of [3, 16] as: no-rub motion → a synchronous full annular rub → a partial rub → dry whip (heavy rub). So dry whip was triggered solely by the mass unbalance for the rotor/stator systems in this case, which is different from the case mentioned above where an outside disturbance was needed. This indicates the necessity to further explore the existence boundaries of different responses of rotor/stator contact systems so as to get an overall picture on the interactive effect of system parameters on the appearance of different response cascades. The recent analytical works [4, 17, 18] have provided some satisfactory explanation on the response characteristics of rotor/stator contact systems mentioned in above tests.

This paper aims to depict the global response characteristics of a specific rotor/stator contact system by determining the existence conditions of the most typical rotor responses of the piecewise smooth nonlinear system. As known, the piecewise smooth nonlinear model consists of three pieces of equations of motion governing three phases of rotor motion, namely, the phase without rub, the phase with rubs and the phase in self-excited dry friction backward whirl. As will be shown in this paper, the existence condition for the most often observed responses, i.e., synchronous no-rub motion, synchronous full annular rub, partial rub and dry friction backward whirl, can be derived solely from a single piece of the equations even though some responses, i.e., the partial rubs, run across two phases of the rotor motion. After incorporating the existence boundaries of all of the responses together into the parameter space, an overall picture of the response characteristics of the rotor/stator contact system is achieved. Deep insights into the interactive effect of different parameters on the system responses as well as the coexistence behavior are gained. Moreover, the analytical methods used in this paper can be also applied to the more complex rotor/stator systems, i.e., the model considering the stator dynamics and the contact stiffness at contact surfaces.

The paper is organized in the following ways. In Sect. 2, the equations of motion for each phase of motion of the rotor/stator contact system are first introduced. Then the equations for different phases of
motion are studied individually in Sect. 3 in order to derive the solutions and the existence boundaries of the corresponding responses. Section 4 will illustrate the global response characteristics of the present rotor/stator contact system in the parameter space through examples. Finally in Sect. 5, the conclusions from this work are drawn.

## 2 Mathematical model

The rotor/stator contact system consists of a Jeffcott rotor and a stator that is modeled as an added stiffness (see Fig. 1). This model physically describes the rotor in contact with a non-rotating susceptible circumferential stator or a mechanical seal. The rotor has a weightless shaft carrying a disk with mass $m$ at the middle of the span. The mass center of the rotor is located at a distance $e$ from its geometrical center. The stiffness of the shaft is $k_s$. The clearance between the rotor and the stator is denoted by $r_0$. The stator, which is rigidly fixed, has an elastic contact surface modeled as a symmetrical set of radial springs with isotropic

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 1 a Schematic diagram of Jeffcott rotor with the stator clearance. b The section of the rotor at the position of disk to show the forces applied on the rotor during the backward whirl of the rotor</figcaption></figure>

stiffness, $k_b$, and the (sliding) friction coefficient $\mu$. The whirl frequency of the shaft and the rotor radius at the axial position of contact are denoted by $\omega_w$ and $r_{disk}$ respectively. In the present analysis the gravity is neglected.

The governing equations of the rotor/stator contact system can be written as

$$ 
\begin{aligned}
m\ddot{x} + c\dot{x} + k_s x + \Theta k_b \left(1 - \frac{r_0}{r}\right) [x - \mu \text{sign}(v_{\text{rel}})y] &= m\omega^2 \cos \omega t \\
m\ddot{y} + c\dot{y} + k_s y + \Theta k_b \left(1 - \frac{r_0}{r}\right) [\mu \text{sign}(v_{\text{rel}})x + y] &= m\omega^2 \sin \omega t \\
v_{\text{rel}} &= \omega r_{\text{disk}} + \omega_w r
\end{aligned} 
\quad (1) $$

where $r = \sqrt{x^2 + y^2}$ and $v_{rel}$ represents the relative velocity between the rotor and the stator at the contact point. $\Theta$ is Heaviside function with $\Theta = 0$ when $r < r_0$ and $\Theta = 1$ when $r \ge r_0$. Equation (1) implies that only the sliding contact is assumed to occur between the rotor and the stator.

As can be seen from (1) the equations of motion, in fact, consist of three pieces that govern three phases of rotor motion: (i) the phase without no-rub when $\Theta = 0$; (ii) the phase with rubs $\Theta = 1$ (but $v_{rel} > 0$ all the time); (iii) the phase of self-excited backward whirl ($\Theta = 1$ and $v_{rel}$ fluctuates around zero where the dry friction effect must be taken into account). The first two phases of motion are quite familiar and widely used in the previous works to investigate the dynamics of the rotor/stator rubbing. The last phase of motion exhibits the well-known rubbing response: dry friction backward whirl, which contains two regimes—of dry whirl and dry whip—as defined in [16] and which is not really well explored.

Equation (1) will be first rewritten in the non-dimensional form:

$$ 
\begin{aligned}
& X'' + 2\zeta X' + \beta X + \Theta \left(1 - \frac{R_0}{R}\right) (X - \mu \text{sign}(V_{\text{rel}})Y) \\
& = \Omega^2 \cos \Omega \tau \\
& Y'' + 2\zeta Y' + \beta Y + \Theta \left(1 - \frac{R_0}{R}\right) (\mu \text{sign}(V_{\text{rel}})X + Y) \\
& = \Omega^2 \sin \Omega \tau \\
& V_{\text{rel}} = R_{\text{disk}}\Omega + R\omega_b
\end{aligned} 
\quad (2) $$
where the non-dimensional variables are defined as:

$$
\begin{align*}
X &= \frac{x}{e}, & Y &= \frac{y}{e}, & R_0 &= \frac{r_0}{e} \\
R_{\text{disk}} &= \frac{r_{\text{disk}}}{e}, & V_{\text{rel}} &= \frac{v_{\text{rel}}}{e}, & 2\zeta &= \frac{c}{\sqrt{k_b m}} \\
\beta &= \frac{k_s}{k_b}, & \omega_2 &= \sqrt{\frac{k_b}{m}}, & \Omega &= \frac{\omega}{\omega_2}, & \omega_b &= \frac{\omega_w}{\omega_2}
\end{align*}
$$

with $R = \sqrt{X^2 + Y^2}$ and $\Theta = 0$ when $R < R_0$ and $\Theta = 1$ when $R \ge R_0$. The prime indicates the differentiation with respect to the dimensionless time variable $\tau = \omega_2 t$.

By defining $W = X + jY$ and $|W| = \sqrt{X^2 + Y^2}$, (2) can be written in a complex form. Rather than writing a unified complex equation, we will write the equations into three pieces of equations according to the phases of motion.

The piece of equation governing the phase without rubbing ($|W| < R_0$) is given by

$$W'' + 2\zeta W' + \beta W = \Omega^2 e^{j\Omega\tau} \quad (3)$$

The piece of equation describing the phase with rubbing ($|W| \ge R_0$) can be written as

$$
\begin{align}
W'' + 2\zeta W' + \beta W + (1 + j\mu)\left(1 - \frac{R_0}{|W|}\right)W \nonumber \\
&= \Omega^2 e^{j\Omega\tau} \tag{4}
\end{align}
$$

The piece of equations describing the phase of self-excited backward whirl ($|W| \ge R_c$ with the sign of $V_{rel}$ alternating) can be written as

$$
\begin{align}
W'' + 2\zeta W' + \beta W + (1 + j\mu \operatorname{sign}(V_{\text{rel}})) \nonumber \\
\times \left(1 - \frac{R_0}{|W|}\right) W &= \Omega^2 e^{j\Omega\tau} \tag{5} \\
V_{\text{rel}} &= R_{\text{disk}}\Omega + |W|\omega_b \nonumber
\end{align}
$$

As is seen from the above discussion, the dynamics of the rotor/stator system in different phases of motion are governed by different equations of motion depending upon the amplitudes of the rotor deflection. The equation in each phase of motion is smooth. So we can solve each piece of equation(s) to get the corresponding solution and then carry out the stability analysis to determine its existence boundary. As will be shown, the most often observed responses in the rotor/stator systems, such as the synchronous no-rub motion, the synchronous full annular rub motion and the

dry friction backward whirl, can all be determined in this way. Although the quasi-periodic partial rub motion alternates between the phase without rubbing and the phase with rubbing, its existence boundaries can still be determined by investigating the piece of equation governing the phase with rubbing. After the existence boundaries of the solutions from all the pieces of equations are determined, the global response behavior of the system in the parameter space can be depicted.

### 3 The solutions and stability analysis

#### 3.1 The solution in the phase without rubbing

In this phase, the rotor is free from the contact with the stator, the synchronous no-rub motion is the only response and can be derived from (3). The periodic solution can be written as

$$W_0 = A e^{j(\Omega\tau + \varphi)} \\ \text{with } A = \frac{\Omega^2}{\sqrt{(\beta - \Omega^2)^2 + (2\zeta\Omega)^2}}, \quad \operatorname{tg}\varphi = \frac{-2\zeta\Omega}{\beta - \Omega^2} \quad (6)$$

The periodic solution is unconditionally stable. However, the solution is physically valid only when its amplitude is smaller than the clearance, namely $A < R_0$, due to the existence of stator. To find the parameter regions where the synchronous no-rub motion exists, let $A = R_0$ in (6) and we get the parameter equation defining the boundary of no-rub motion:

$$(R_0^2 - 1)\Omega^4 + R_0^2(4\zeta^2 - 2\beta)\Omega^2 + R_0^2\beta^2 = 0 \quad (7)$$

Equation (7) yields two pairs of real roots, each with an equal magnitude and opposite sign. Since the rotating speed is supposed to be positive and the synchronous no-rub response always undertakes forward whirl, two positive roots, $\Omega_l$ and $\Omega_u$, are chosen to define the boundaries of no-rub motion. That is, the synchronous no-rub motion exists when the rotating speed $\Omega$ is less than $\Omega_l$ or greater than $\Omega_u$.

#### 3.2 The solutions in the phase with rubbing

Synchronous full annular rub motion is a response by which the rotor is in the full annular contact with the stator and whirls forward with constant amplitude at
a frequency equal to the rotating speed of the rotor. Equation (4) can be used to determine this response by supposing a proper form of the solution according to the characteristics of the synchronous full annular rub.

Let us denote the periodic solution as $W_1 = Be^{j(\Omega\tau+\varphi)}$, two nonlinear equations of B and $\varphi$ are yielded after substituting $W_1$ into (4) and separating the real and imaginary parts of the obtained equation,

$$ (\beta - \Omega^2)B + (B - R_0) = \Omega^2 \cos \varphi \quad (8) \\ 2\zeta \Omega B + \mu(B - R_0) = -\Omega^2 \sin \varphi $$

To eliminate $\varphi$ from the above equations, we square both equations of (8) and then sum them up to get a second order polynomial equation of B in the form

$$ [(2\zeta\Omega + \mu)^2 + (1 + \beta - \Omega^2)^2]B^2 \\ - 2R_0[(1 + \beta - \Omega^2) + \mu(2\zeta\Omega + \mu)]B \\ + R_0^2(1 + \mu^2) - \Omega^4 = 0 \quad (9) $$

One can solve for the positive real roots of B from (9) and then substitute them back into either of the two equations of (8) to get $\varphi$ for the given set of parameters. From the condition with rubbing, B ought to be greater than $R_0$.

The response of the synchronous full annular rub corresponds to the stable periodic solution of (4). So the stability of the solution must be analyzed. To carry out the stability analysis on the periodic solution of the nonlinear system, the physical understanding of the response is of good use. It is known that a periodic solution in the inertial coordinate system will correspond to a fixed point in the rotating coordinate system at the same rotating speed as the period of the solution. To transform (4) into the rotating coordinate system, we use the transformation $W_1(\tau) = \rho_r(\tau)e^{j\Omega\tau}$. Substituting it into (4) and making some manipulations, yields

$$ \rho_r'' + 2(\zeta + j\Omega)\rho_r' + [(1 + \beta - \Omega^2) \\ + j(2\Omega\zeta + \mu)]\rho_r - (1 + j\mu)R_0 \frac{\rho_r}{|\rho_r|} = \Omega^2 \quad (10) $$

Equation (10) has a fixed point in the form of $\rho_{r0} = Be^{j\varphi}$ with B and $\varphi$ being defined by (8).

Further, the stability analysis on the fixed point, and thus on the stability of the synchronous full annular rub solution in the inertial coordinate system, will be

carried out. To accomplish this task, the complex amplitude of the response is written in the component form $\rho_r = X_r + jY_r$, where the subscript “r” indicates the variable in the rotating coordinate system.

By defining the state vector in the rotating coordinate system as $\underline{U} = \{X_r Y_r X'_r Y'_r\}^T$, (10) can be written in a time-invariant differential equation

$$ \underline{U}' = \underline{A} \underline{U} + \underline{g}(\underline{U}) = \underline{G}(\underline{U}) \quad (11) $$

where $\underline{A}$ and $\underline{g}(\underline{U})$ are respectively the coefficient matrix and the nonlinear vector term:

$$ \underline{A} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -(1 + \beta - \Omega^2) & 2\Omega\zeta + \mu & -2\zeta & 2\Omega \\ -(2\Omega\zeta + \mu) & -(1 + \beta - \Omega^2) & -2\Omega & -2\zeta \end{bmatrix} $$

and $\underline{g}(\underline{U}) = \{0 \ 0 \ H_X + \Omega^2 \ H_Y\}^T$ with $H_X = \frac{R_0}{R}(X_r - \mu Y_r)$, $H_Y = \frac{R_0}{R}(\mu X_r + Y_r)$.

Linearizing (11) about the solution of the fixed point $\underline{U}^0 = \{B \cos \varphi \ B \sin \varphi \ 0 \ 0\}$, yields a linear, time-invariant differential equation

$$ \underline{\delta U}' = \underline{J}(\underline{U}^0)\underline{\delta U} \quad (12) $$

where $\underline{\delta U} = \underline{U} - \underline{U}^0$ is the perturbation solution to the fixed point. $\underline{J}(\underline{U}^0)$ is the so-called Jacobian matrix which is given in the following form:

$$ \underline{[J(\underline{U}^0)]} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -1 - \beta + \Omega^2 & (1 - a)\mu + 2\zeta\Omega & -2\zeta & 2\Omega \\ -\mu - 2\zeta\Omega & -1 - \beta + \Omega^2 + a & -2\Omega & -2\zeta \end{bmatrix} \\ \text{with } a = R_0/B \quad (13) $$

From Jacobian matrix given in (13), one can set up the characteristic equation of the solution,

$$ s^4 + b_3 s^3 + b_2 s^2 + b_1 s + b_0 = 0 \quad (14) $$

with $b_3 = 4\zeta$, $b_2 = 2\Omega^2 + 4\zeta^2 - a + 2(1 + \beta)$, $b_1 = 2\mu(2 - a)\Omega - 2\zeta a + 4\zeta(1 + \beta + \Omega^2)$, $b_0 = \Omega^4 + (4\zeta^2 + a - 2)\Omega^2 + 2\mu\zeta(2 - a)\Omega + (1 + \beta)^2 - a(1 + \beta) - \mu a + \mu^2$.

The real parts of the characteristic roots decide the stability of the fixed point, and thus the stability of the synchronous full annular rub solution. Since the stability boundaries are of interest, bifurcation theory is being used as is done in [4]. The saddle-node bifurcation condition, by which a zero characteristic root exists, defines the existence boundary of the solution. The Hopf bifurcation condition, whereby a pair of
pure imaginary characteristic roots appears, gives the boundary between the quasi-periodic partial rub and the synchronous full annular rub.

Quasi-periodic partial rubbing is a response of the rotor/stator system, by which the rotor rubs the stator intermittently. There are generally more than one frequency component in the system responses during partial rubs. The equations governing the partial rubs should be a combination of (3) for the phase without rubbing and (4) for the phase with rubbing. Generally, we have difficulty in obtaining the analytical solutions of the piecewise smooth equations. Without the solution of partial rub, it seems impossible to determine their existence boundaries.

Fortunately, we can still sketch the existence boundaries of quasi-periodic partial rub based on the study of other responses and on the understanding of the physical reason for the loss of stability of the partial rub. For instance, it is easy to determine the boundary between no-rub motion and partial rub. Since the synchronous full annular rub bifurcates into the partial rub through Hopf bifurcation, the Hopf bifurcation condition of synchronous full annular rub solution will serve as the boundary between the synchronous full annular rub and the partial rub. As observed in experiment [3, 14] and simulations [4], the partial rub ceases to exist after the dry friction backward whirl was triggered by the imbalance with the increase rotating speed. In [17] the physical reason for the onset of dry friction backward whirl was exploited. It is known that the rotor in resonance at a backward whirl “natural” frequency of the coupled nonlinear rotor/stator system is the reason for the onset of the dry friction backward whirl from the partial rub. From this understanding, a method can be devised to determine the approximate onset condition—another existence boundary of the partial rub.

First, the backward whirl “natural” frequency of the coupled nonlinear rotor/stator system is sought. This frequency corresponds to the mode that the rotor whirls backwards. After setting the right-hand side of (4) to zero and then substituting supposed solution, $W = H \exp(j\omega_0\tau)$, into it, we solve the frequency with a negative value and the amplitude ($H$ should be greater than $R_0$) as

$$ \omega_0 = - \left( \frac{\zeta}{\mu} + \sqrt{\beta + \frac{\zeta^2}{\mu^2}} \right) $$

$$ H = \frac{R_0}{1 - 2\frac{\zeta}{\mu}\left(\frac{\zeta}{\mu} + \sqrt{\beta + \frac{\zeta^2}{\mu^2}}\right)} \\ \text{when } \frac{\zeta}{\mu} < \frac{1}{2\sqrt{1 + \beta}} \qquad (15) $$

Secondly, when the rotor runs at a rotating speed near the magnitude of the backward whirl “natural” frequency, (4) is transformed into a form that describes the dynamics of the rotor/stator system near the backward whirl “natural” frequency,

$$ \begin{aligned} W'' + 2\zeta W' + \beta W + (1 + j\mu)\left(1 - \frac{R_0}{|H|}\right)W \\ = R_0(1 + j\mu)\left(\frac{1}{|W|} - \frac{1}{|H|}\right)W + \Omega^2 e^{j\Omega\tau} \end{aligned} \qquad (16) $$

where the left-hand side of the complex system has the natural frequency equal to the above backward whirl “natural” frequency.

Thirdly, an approximated linear equation to (16) is set up after defining

$$ \varepsilon = D_f = \left( \frac{1}{|W|} - \frac{1}{|H|} \right) = O(\varepsilon) < 0 \qquad (17) $$

as a small parameter. The reason is that $|H|$ is the amplitude of the rotor in the backward whirl rubbing without imbalance excitation and $|W|$ is the amplitude of the rotor in the backward whirl rubbing with imbalance excitation. Since $\Omega < 1$ in the vicinity of the critical rotating speed of the onset of dry whip, the magnitude of the excitation $\Omega^2$ in (4) is thus a small quantity. So $|W|$ should be only a small modification to $|H|$ and the difference between the corresponding inverse relations should be also small. So (16) can be written approximately as:

$$ \begin{aligned} W'' + 2\zeta W' + \beta W + (1 + j\mu)\left(1 - \frac{R_0}{|H|}\right)W \\ = \varepsilon R_0(1 + j\mu)W + \Omega^2 e^{j\Omega\tau} \end{aligned} \qquad (18) $$

The approximated equation can be solved by, e.g., the Multiple Scale Method [19], since $W$ at the rotating speed near the magnitude of the backward whirl “natural” frequency, namely $\Omega = |\omega_0| + \varepsilon\sigma$ where $\sigma = O(1)$, is a small detuning.

Assume that the problem depends on many time scales:

$$ T_0 = \tau, \quad T_1 = \varepsilon\tau, \quad T_2 = \varepsilon^2\tau, \dots $$
The solution of (18) can be represented by an expansion having the form

$$W(\tau; \varepsilon) = W_0(T_0, T_1, \dots) + \varepsilon W_1(T_0, T_1, \dots) + \dots \quad (19)$$

The expansion to $O(\varepsilon^2)$ is carried out in this paper and the corresponding equations of $\varepsilon^0$ and $\varepsilon^1$ are

$$\varepsilon^0: \quad D_0^2 W_0 + 2\zeta D_0 W_0 + \kappa W_0 = \Omega^2 \exp(j\Omega T_0) \quad (20)$$

$$\varepsilon^1: \quad D_0^2 W_1 + 2\zeta D_0 W_1 + \kappa W_1 \\ = -2D_0 D_1 W_0 - 2\zeta D_1 W_0 + R_0(1 + j\mu)W_0 \quad (21)$$

where $\kappa = \beta + (1 + j\mu)(1 - R_0/|H|)$. By this moment, the range of the excitation frequency (the rotating speed) $\Omega$ is not specified as done in the analysis procedure of the Multiple Scale Method at the Secondary Resonance [19]. The solution of (20) yields

$$W_0 = a(T_1) \exp[j(\omega_0 T_0) + \theta(T_1)] + A \exp(j\Omega T_0) \quad (22)$$

where $a(T_1)$ and $\theta(T_1)$ are the real functions of $T_1$ and $A = \Omega^2/(\kappa - \Omega^2 + j2\zeta\Omega)$.

Substituting (22) into (21) leads to

$$D_0^2 W_1 + 2\zeta D_0 W_1 + \kappa W_1 \\ = [-2(j\omega_0 a' - \omega_0 a\beta') - 2\zeta(a' + ja\beta')] \\ \times \exp[j(\omega_0 T_0 + \theta)] \\ + R_0(1 + j\mu)a \exp[j(\omega_0 T_0 + \theta)] \\ + A R_0(1 + j\mu) \exp(j\Omega T_0) \quad (23)$$

where the prime denotes the derivative with respect to $T_1$.

Now, substituting the above expression of $\Omega$ into (23) and setting the coefficient of terms with $\exp(j\omega_0 T_0)$ on the right-hand side of (23) to zero in order to eliminate the secular term, leads to a complex equation

$$[-2(j\omega_0 a' - \omega_0 a\theta') - 2\zeta(a' + ja\theta') + R_0(1 + j\mu)a] \times \exp(j\theta) + A R_0(1 + j\mu) \exp(j\sigma T_1) = 0 \quad (24)$$

Separating the real and the imaginary parts and then solving $a'$ and $a\theta''$, gives

$$\begin{aligned} a' &= C_1(\mu\omega_0 + \zeta)a + C_1C_2\sin(\sigma T_1 - \theta + \alpha) \\ a\theta' &= C_1(\mu\zeta - \omega_0)a - C_1C_2\cos(\sigma T_1 - \theta + \alpha) \end{aligned} \quad (25)$$

where $C_1 = \frac{R_0}{2(\omega_0^2 + \zeta^2)}$, $C_2 = R_0\sqrt{(1 + \mu^2)(\omega_0^2 + \zeta^2)} \times |A|$, $\alpha = \arctan \frac{\zeta}{\omega_0} + \arctan \frac{\mu A_r + A_i}{A_r - \mu A_i}$. $A_r$ and $A_i$ are the real and imaginary parts of $A$.

Next, (25) is transformed into an autonomous system with $T_1$ not explicitly appearing by defining

$$\gamma = \sigma T_1 - \theta + \alpha \quad (26)$$

and substituting into (25) to get

$$\begin{aligned} a' &= C_1(\mu\omega_0 + \zeta)a + C_1C_2\sin(\gamma) \\ a\gamma' &= \sigma a - C_1(\mu\zeta - \omega_0)a + C_1C_2\cos(\gamma) \end{aligned} \quad (27)$$

The steady-state motions occur when $a' = \gamma' = 0$, which corresponds to the singular points of (27). They correspond to the solution of

$$\begin{aligned} C_1(\mu\omega_0 + \zeta)a_0 &= -C_1C_2\sin(\gamma_0) \\ \sigma a_0 - C_1(\mu\zeta - \omega_0)a_0 &= -C_1C_2\cos(\gamma_0) \end{aligned} \quad (28)$$

$|W_0|$ with $a$ and $\theta$ by (28) and (26) can serve as an approximation of $|W|$ in (17). By using $|H| = H$ in (15), $D_f$ in (17) can be calculated for the given parameters and $\varepsilon$.

Finally, the critical rotating speed for the onset of the dry friction backward whirl or the loss of stability of the partial rub can be solved after the deviation $\varepsilon\sigma$ to $|\omega_0|$ can be determined. To do this, we can first set $\sigma$ a reasonable constant number, i.e., $\sigma = 1$ as used in this paper or a value that depends on the system parameters as given in [17], then use a trial and error process on $\varepsilon$ in order to find the most proper $\varepsilon$ that makes $D_f$ (as given above) a largest negative value.

### 3.3 The solution in the phase of self-excited backward whirl

Dry friction backward whirl is a motion state of the rotor/stator system, by which the rotor is in continuous contact with the stator, rolling or slipping continuously on the contact surface and whirling backward at a super-synchronous frequency. Dry friction backward whirl is a self-excited vibration of the rotor/stator contact system due to the presence of dry friction effect at the contact surface. Although the mechanism of dry friction backward whirl is well understood, most of the research works used the governing equations, which
do not include the dry friction effect to explore this phenomenon. A little progress has been made in analytically determining the characteristics and the parameter region of dry friction backward whirl until now. Below, we will only sketch the solution procedure.

To make the reader to better understand the meaning of (5) and the methodology we adopted to approximate the amplitude of dry friction backward whirl, the mechanism for a rotor/stator system that performs dry friction backward whirl is first illustrated.

When rotor-to-stator contact occurs, the friction force at the contact point tends to drive the rotor into backward whirl. In case that the rotor whirls backward, the friction force will apply in the same direction as the rotor whirl and continuously turns the input rotational energy into the energy for the lateral vibration of the rotor system. The vibration amplitude of the rotor is driven higher and higher. The relative velocity at the contact point will change its sign when the rotor deflection is sufficiently large to make the whirl speed exceed the circumferential speed of the rotor at the contact point, i.e., $r = -\omega r_{\text{disk}}/\omega_w$ (refer to Fig. 1(b)) or $\omega_b|W| \ge \Omega R_{\text{disk}}$ (see (5)). At this time, the friction force also inverts its direction and is opposite to the rotor whirl to dissipate the energy. The rotor deflection will be driven down. When the rotor deflection decreases to make the relative velocity at the contact point to change sign again, the rotor deflection will be driven up again. In this way, the amplitude of rotor deflection at dry friction backward whirl is kept oscillating around the value with a zero relative velocity, that is, the rotor undertakes a dry friction induced self-excited vibration—dry friction backward whirl.

Now we can use the knowledge of dry friction backward whirl to derive the solution from (3). Since the fluctuation of the response amplitude around the zero relative velocity is much smaller than the total magnitude of the rotor deflection during dry friction backward whirl, the amplitude of dry friction backward whirl can be well approximated by the one at zero relative velocity. From the second equation in (5), it yields

$$|W| = -\Omega R_{\text{disk}}/\omega_b \qquad (29)$$

By substituting (29) into first equation of (5), we obtain the equation governing the motion of dry friction backward whirl in the following form:

$$\begin{aligned} W'' + 2\zeta W' + \beta W + (1 + j\mu) \left(1 + \frac{\omega_b R_0}{\Omega R_{\text{disk}}}\right) W \\ = \Omega^2 e^{j\Omega\tau} \end{aligned} \qquad (30)$$

To derive solution for dry friction backward whirl from (30), two important features should be taken into account. (a) Since the contribution of the forced vibration to the total response at dry friction backward whirl is generally too small to be considered, the supposed solution needs only to consider the frequency component with super-synchronous frequency. (b) Most importantly, when the governing equation changes into the form of (30), the solution corresponding to dry friction backward whirl should have a time increasing amplitude, that is, a positive real part in the exponent, because the dry friction effect does not explicitly take effect (its effect is taken into account in (29) to give a constant amplitude). So the solution can be assumed in the form

$$W_2 = C e^{(\alpha + j\omega_b)\tau} \qquad (31)$$

After substituting (31) into (29), one can get two equations of $\alpha$ and $\omega_b$, where $\omega_b$ indicates the whirl frequency at dry friction backward whirl and $\alpha$ the capability of the system to undertake self-excited vibration. When $\alpha$ is positive, self-excited vibration—dry friction backward whirl—can take place. When $\alpha$ is negative, the responses other than self-excited vibration occur. Therefore, $\alpha = 0$ is the critical condition for the existence of dry friction backward whirl [18]. In this way the boundary for dry friction backward whirl in the parameter space can be determined.

## 4 Global response characteristics

In this section, two examples are presented to illustrate the overall picture of the global dynamic response characteristics of the rotor/stator contact system with $\zeta = 0.05$, $R_0 = 2.00$ and $R_{\text{disk}} = 20R_0$. As can be seen from the above analysis the existence boundaries of no-rub motion, the synchronous full annular rub and the dry friction backward whirl are determined rigorously with very high accuracy. The numerical simulations have proved this point. For the onset condition of dry friction backward whirl from partial rub motion, some approximations are adopted. The validity of
**Fig. 2** Rotor response characteristics on parameter plane of Ω-μ with β = 0.04, where ζ = 0.05, R₀ = 2.0 and R<sub>disk</sub> = 20R₀. The regions of the five typical types of coexistence are identified. Ω<sub>l</sub> and Ω<sub>u</sub> are the critical rotating speeds of no-rub motion. SN and HP are the saddle-node and Hopf bifurcation boundaries of the synchronous full annular rub. DF and DW are the existence and onset boundaries of dry friction backward whirl

<figure><img src="image_4.png" /></figure>

the results is also checked by the numerical simulation as given in [17]. So the global dynamical characteristics of the rotor-to-stator contact system will be well demonstrated in the following analysis.

In the following figures, curves Ω<sub>l</sub> and Ω<sub>u</sub> indicate the rotating speed whereby the deflection of the linear rotor covers the clearance. Thus, the synchronous no-rub motion is in the rotating speed range less than Ω<sub>l</sub> and greater than Ω<sub>u</sub>. Curve HP and line SN are respectively the boundaries of Hopf bifurcation and saddle-node bifurcation of the synchronous full annular rub solution. Curve DW is the boundary by which dry whip is triggered by imbalance and the partial rub ceases to exist. Curve DF gives the existence boundary of the dry friction backward whirl. Also, the parameter regions marked by numbers represent different regions of coexistence, i.e., “0”—no-rub motion + dry friction backward whirl; “1”—synchronous full annular rub + dry friction backward whirl; “2”—partial rub + dry friction backward whirl; “3”—synchronous full annular rub + no-rub motion; “4”—no-rub motion + synchronous full annular rub + dry friction backward whirl.

Figure 2 depicts the global responses characteristics on the parameter planes of Ω-μ with β = 0.04. It is interesting to note that the dry friction backward whirl exists already from a very low rotating speed when the friction coefficient is larger than a certain value which satisfies the condition given in (15) as μ > 0.102. When the friction coefficient is below this

value, the only rubbing response is the synchronous full annular rub solution. Due to the coexistence of the no-rub motion and the synchronous full annular rub (see region “3”), the “jump” phenomenon will occur in this region. When the friction coefficient is above this value, one can find the regions that the dry friction backward whirl coexists with other responses. So it is quite possible to induce the more destructive dry friction backward whirl from a relatively mild rubbing response or even from no-rub motion when the rotor is subjected to a sudden outside disturbance. It is even worse that there is a region (enclosed by DW and Ω<sub>u</sub>) where dry whip is the sole stable response. Thus, it is important to avoid a large contact friction, especially that which exceeds the one defined by the criterion in (15). The coexistence of three responses can be also detected in a very small region—“4”.

Figure 3 shows the global response characteristics on the parameter planes of Ω-β with μ = 0.15. It is easily seen from Fig. 3 that with the increase of the stiffness ratio (it can be taken as the reduction of the stator stiffness but keeping the rotor stiffness unchanged, namely, using a soft stator), both Ω<sub>l</sub> and Ω<sub>u</sub> increase but the latter increases more rapidly. This makes the rotating speed range with the synchronous full annular rub larger and the speed range with the no-rub motion smaller. This seems to be unprofitable to the rotor response. However, it is found from Fig. 3 that when the stiffness ratio is larger than a certain value, i.e., greater than about 0.3 in this case, the
**Fig. 3** Rotor response characteristics on parameter plane of Ω-β with μ = 0.15, where ζ = 0.05, R₀ = 2.0 and Rₘₐrk = 20R₀. The regions of only three of the five typical types of coexistence are found in this case. Ωᵢ and Ω are the critical rotating speeds of no-rub motion. SN and HP are the saddle-node and Hopf bifurcation boundaries of the synchronous full annular rub. DF and DW are the existence and onset boundaries of dry friction backward whirl

<figure><img src="image_5.png" /></figure>

synchronous full annular rub will no longer bifurcate
into the more severe partial rub. More importantly, the
induction of the destructive dry whip through the im-
balance will not occur. In fact, when the stiffness ratio
is greater than about 2.4, the imbalance can no longer
trigger the dry whip although the rotor undertakes the
partial rub. From this sense, a larger stiffness ratio ben-
efits the rotor response. Since the criterion given in
(15) is met in the range of the stiffness ratio examined
in Fig. 3, the dry friction backward whirl exists in the
whole region right to line DF. This means that the dry
friction backward whirl coexists with all other rotor
responses we have examined. From the point of view
of the global responses, there are only three types of
coexistence regions in this case, except for region “3”
and region “4” mentioned above.

**5 Conclusions**

In this paper, the global response characteristics of
a piecewise smooth rotor/stator contact system are
investigated. An effective method to get the global
response characteristics of this specific model has
demonstrated that the piecewise smooth equations of
motion can be analyzed piece by piece so that the
existence boundaries of the typical responses in each
phases of motion, i.e., in the phase without rub, with
rub and of self-excited backward whirl, have been de-
termined. Then, they are incorporated into the same

parameter space to get the global response characteris-
tics of the rotor/stator contact system. The results show
deep insights into the interactive effect of different pa-
rameters on the system responses as well as the coex-
istence behavior. The main conclusions can be sum-
marized as follows.

From the point view of global responses, five typ-
ical types of coexistence are found for the piecewise
smooth rotor/stator contact system in this paper. They
are the coexistence of no-rub motion and synchronous
full annular rub, of no-rub motion and dry friction
backward whirl, of synchronous full annular rub and
dry friction backward whirl, of partial rub and dry fric-
tion backward whirl and of no-rub motion, synchro-
nous full annular rub and dry friction backward whirl.
Except for the first type of coexistence, the other types
of coexistence are not quite familiar and occur, in fact,
when the criterion given in (15) is met. From this cri-
terion and also from diagrams of the global response
characteristics in the parameter planes, it is known that
small friction on the contact surfaces and large stiff-
ness ratio of the rotor-to-stator stiffness will benefit the
rotor rubbing behavior by avoiding the occurrence of
the severe partial rub or the destructive dry whip.

**Acknowledgements** This work is supported by the National Natural Science Foundation of China (NSFC) under Grant 10472086 and by the Project for Excellent Young Teachers of MOE, China.
References

1. Muszynska, A.: Rotor-to-stationary element rub-related vibration phenomena in rotating machinery—literature survey. Sound Vib. Dig. **21**, 3–11 (1989)

2. Black, H.F.: Interaction of a whirling rotor with a vibrating stator across a clearance annulus. Int. J. Mech. Eng. Sci. **10**, 1–12 (1968)

3. Yu, J.J., Goldman, P., Bently, D.E., Muszynska, A.: Rotor/seal experimental and analytical study on full annular rub. ASME J. Eng. Gas Turbine Power **124**, 340–350 (2002)

4. Jiang, J., Ulbrich, H.: Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients. Nonlinear Dyn. **24**, 269–283 (2001)

5. Bently, D.E.: Forced subrotative speed dynamic action of rotating machinery. ASME, Paper No. 74-PET-16 (1974)

6. Childs, D.W.: Rub-induced parametric excitation in rotors. J. Mech. Des. **101**, 640–644 (1979)

7. Ehrich, F.F.: High order subharmonic response of high speed rotors in bearing clearance. ASME J. Vib. Acoust. Stress Reliab. Des. **110**, 9–16 (1988)

8. Day, W.B.: Asymptotic expansions in nonlinear rotor dynamics. Q. Appl. Math. **44**, 779–792 (1987)

9. Kim, Y.B., Noah, S.T.: Quasi-periodic response and stability analysis for a nonlinear Jeffcott rotor. J. Sound Vib. **190**, 239–253 (1996)

10. Goldman, P., Muszynska, A.: Chaotic behavior of rotor/stator system with rubs. J. Eng. Gas Turbines Power **116**, 693–701 (1994)

11. Karpenko, E.V., Wiercigroch, M., Cartmell, M.P.: Regular and chaotic dynamics of a discontinuously nonlinear rotor system. Chaos Solitons Fractals **13**, 1231–1242 (2002)

12. Choi, Y.S.: Investigation on the whirling motion of full annular rotor rub. J. Sound Vib. **258**, 191–198 (2002)

13. Lingener, A., Experimental investigation of reverse whirl of a flexible rotor. Trans. of IFToMM 3rd International Conference on Rotor Dynamics, Lyon, France, pp. 13–18 (1990)

14. Wu, F., Flowers, G.T.: An experimental study of the influence of disk flexibility and rubbing on rotor dynamics, ASME, DE vol. 60, pp. 19–26 (1993)

15. Zhang, W.: Dynamic instability of multi-degree-of-freedom flexible rotor systems due to full annular rub. In: IMechE C252/88, pp. 305–308 (1988)

16. Crandall, S.: From whirl to whip in rotordynamics. In: Proceedings of IFToMM 3rd International Conference on Rotor Dynamics, Lyon, France, pp. 19–26 (1990)

17. Jiang, J., Ulbrich, H.: The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems. ASME J. Acoust. Vib. **127**, 594–603 (2005)

18. Jiang, J.: The analytical solution and the existence condition of dry friction backward whirl in rotor-to-stator contact systems. ASME J. Acoust. Vib. **129**, 260–264 (2007)

19. Nayfeh, A.H., Mook, D.T.: Nonlinear Oscillations. Wiley, New York (1979)