# A finite element-based algorithm for rubbing induced vibration prediction in rotors

Mehdi Behzad, Mehdi Alvandi, David Mba, Jalil Jamali

► To cite this version:

Mehdi Behzad, Mehdi Alvandi, David Mba, Jalil Jamali. A finite element-based algorithm for rubbing induced vibration prediction in rotors. Journal of Sound and Vibration, 2013, 332, pp.5523 - 5542. 10.1016/j.jsv.2013.05.016. hal-01396072

HAL Id: hal-01396072

https://hal.science/hal-01396072

Submitted on 13 Nov 2016

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
# A finite element-based algorithm for rubbing induced vibration prediction in rotors

Mehdi Behzad<sup>a</sup>, Mehdi Alvandi<sup>a</sup>, David Mba<sup>b</sup>, Jalil Jamali<sup>c</sup>

<sup>a</sup> Department of Mechanical Engineering, Sharif University of Technology, Tehran, Iran

<sup>b</sup> School of Engineering, Cranfield University, Cranfield, Bedfordshire MK43 0AL, UK

<sup>c</sup> Mechanical Department, Shoushtar Branch, Islamic Azad University, Shoushtar, Iran

In this paper, an algorithm is developed for more realistic investigation of rotor to stator rubbing vibration, based on finite element theory with unilateral contact and friction conditions. To model the rotor, cross sections are assumed to be radially rigid. A finite element discretization based on traditional beam theories which sufficiently accounts for axial and transversal flexibility of the rotor is used. A general finite element discretization model considering inertial and viscoelastic characteristics of the stator is used for modeling the stator. Therefore, for contact analysis, only the boundary of the stator is discretized. The contact problem is defined as the contact between the circular rigid cross section of the rotor and "nodes" of the stator only. Next, Gap function and contact conditions are described for the contact problem. Two finite element models of the rotor and the stator are coupled via the Lagrange multipliers method in order to obtain the constrained equation of motion. A case study of the partial rubbing is simulated using the algorithm. The synchronous and subsynchronous responses of the partial rubbing are obtained for different rotational speeds. In addition, a sensitivity analysis is carried out with respect to the initial clearance, the stator stiffness, the damping parameter, and the coefficient of friction. There is a good agreement between the result of this research and the experimental result in the literature.

## 1. Introduction

Current industrial needs are continuously demanding higher performance from current turbomachinery. To achieve high performance, losses and hence clearances should be minimized. However, tighter clearances increase the possibilities of contact between moving and stationary parts. This contact is called "rubbing" and can lead to different physical phenomena such as rubbing vibration, thermal bow, power dissipation due to friction forces, etc. Therefore, to achieve high efficiency, more engineering knowledge of the rubbing phenomena should be incorporated into the design process. Data from the early 1970s show the losses at clearances can be of approximately 2.6 percent for a 600 MW steam turbine; and the additional expenses due to worn out seals can be up to 1 percent of the total cost of fuel for a single aircraft turbine [1].

Rubbing contacts are divided into two classes: blade to shroud and rotor to stator. In the latter, which is the focus of this paper, the rotor geometry is considered circular. The stator can be a seal, a bearing, casing, fixed blades, etc.
The rotor to stator rubbing vibration is recognized by two motions types: full annular rubbing wherein the rotor maintains contact with the stator during every whirling motion; and partial rubbing wherein the rotor occasionally makes contact with the stator.

Over the last few decades, many studies have been conducted on the rubbing phenomena, its different dynamics regimes, and its consequences on the machine performance and safety. Muszynska [2] compiled different aspects of rubbing phenomena and some simple models for investigation of rubbing characteristics, from the works which had been previously presented in the literature. Muszynska stated that modeling of the rotor to stator rubbing, as well as definitions of practical machine design criteria to prevent occurrence of unwelcome rub related events, is still a considerable challenge. Choy and Padovan [3] developed a basic model of the transient motion of a rotor bearing casing system during rubbing contact. They employed a simplistic model where the rotor bearing assembly was assumed to be a Jeffcott rotor model; the casing deformation effects were approximated by a rigid casing supported by radial springs; only light rub conditions were considered such that the tangential force generated from friction varied linearly with the radial force exerted by the rotor. Despite the simplicity, this model has been used on numerous occasions to evaluate the effects of various parameters (e.g. stiffness, mass and damping of the rotor or the stator, clearances, rotational speed, etc.) on the machine behavior during rubbing conditions. Furthermore, steady state, transient, or chaotic responses of the rotor due to the rubbing have also studied by this simple model [1,4 18].

However, Roques et al. [19] initiatively employed finite element method to model the diaphragm of a turbogenerator as an actual stator in order to investigate speed transients with rotor to stator rubbing caused by an accidental blade off imbalance. They solved highly nonlinear equations due to contact conditions through a time marching procedure combined with the Lagrange multiplier approach dealing with a "node to line" contact strategy. They assumed that only one point of the cross section boundary of the rotor comes into contact with one curved element of the diaphragm. This approach seems to poorly predict the interaction between the rotor and the stator when the contact area cannot be approximated by a point due to local deformation of the stator.

Using finite element (FE) theory with unilateral contact and friction conditions, the objective of this work is to develop an algorithm in which the model of the stator is extended from an annulus rigid block to a deformable body with any arbitrary geometry. Assuming that two FE models of the flexible rotor and the deformable stator are independently available, a quasi "node to node" contact strategy is used for treating the interaction of bodies at the contact interface. This strategy efficiently works when the contact area widens due to local deformation. Then, two FE models are coupled via the Lagrange multiplier technique in order to obtain the constrained equation of motion.

## 2. Finite element modeling

Every turbomachine has a rotating part, the "rotor", which generally exhausts energy from a working fluid. The rotor is usually composed of a shaft, disks, spacers and blades. This rotor is nested in and enclosed by a casing through a set of bearings. The leakage of working fluids is suppressed via the seals at the ends of the rotor. The seals, bearings, and any stationary part located around the rotor are called the "stator". The schematic configuration of a typical rotor, associated with the bearings and the seals, is depicted in Fig. 1. The geometry of rotor cross sections in stations which are prone to rub with the stator is usually circular. For this reason, the contact algorithm is aimed to treat dynamics of the contact occurring at any circular cross section of the rotor around which the stator is located.

In order to have a simple model for the rotor, cross sections are assumed to be radially rigid. Thus, a finite element discretization based on a traditional beam theory would adequately accounts for axial and transversal flexibility of the rotor. For example, the shaft can be axially discretized using Timoshenko's beam theory, and the bearings, disks, blades and unbalances effects are then included in the model. A simple Jeffcott model can also be used as the rotor model. Thus, the

<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic configuration of a rotor: (a) axial view and (b) sectional view.</figcaption></figure>
rotor can be represented either by only one point, as in the Jeffcott model, or by a set of nodes, as in the FE model. Moreover, to model the stator, the body of the stator is also discretized and inertial and viscoelastic characteristics can be considered in the model. Depending on how the real stator is modeled, the stator geometry can be a continuous annulus, segmented arcs or a set of separated nodes or a single node.

## 2.1. Assumptions and simplifications

The following assumptions were employed for the finite element model:

* Displacements are infinitesimal either in the rotor or in the stator.

* Both the material of the rotor and stator behave linearly and elastically.

* Despite axial and transversal flexibility, the rotor is radially rigid.

* Angular rotation of the cross section is discarded.

* Interaction of the fluid between the rotor and the stator is ignored.

* To model contact friction forces, it is assumed that the rotor never sticks to the stator and always slips due to high rotational speed.

The rotor and stator modeling, and dynamics of the contact between them are presented in the next sections.

## 2.2. Modeling of the rotor

To model the rotor system all local matrices of shaft elements, disks, added masses and unbalances should be calculated and assembled into global matrices. Timoshenko's beam theory can be used to model the shaft. Considering gyroscopic effects, the FE equation for a typical rotor is as follows [20]:

$$ M_r u_r + (D + G(\Omega)) \dot{u}_r + K_r u_r = f_r(t), \quad (1) $$

where $u_r$ is the vector of degrees of freedom, $f_r$ is the vector of external forces, $K_r$ is the stiffness matrix, $M_r$ is the mass matrix, $D$ is the damping matrix due to bearings, and $G(\Omega)$ is the gyroscopic matrix of the rotor with the angular velocity $\Omega$. Letting total damping matrix of the rotor be $C_r$, Eq. (1) can be rewritten as:

$$ M_r u_r + C_r \dot{u}_r + K_r u_r = f_r(t). \quad (2) $$

## 2.3. Modeling of the stator

The stator can be an obstacle, a protrusion, circular casing, turbine diaphragm, segmented labyrinth seal, retainer bearings, sliding bearings, or fixed blades of the machine. To model these stators, various element types such as bar, beam, frame, curved beam, quadrilateral, or shell elements can be used in combination with each other. Any of these stators can be connected to the casing and foundation in order to study the effect of rubbing forces on the casing or the foundation. The FE equation of a linear stator is constructed as follows:

$$ M_s u_s + C_s \dot{u}_s + K_s u_s = f_s(t), \quad (3) $$

Eq. (3) is similar to Eq. (2) whereas the subscript "s" stands for the stator.

## 2.4. Modeling of contact dynamics

The aim of this section is to illustrate how two FE models are combined and how the constrained equation of motion is obtained on the basis of FE theory with unilateral contact and friction conditions. Many powerful approaches have been already developed for this purpose. The solution of a contact problem generally includes: (a) definition of appropriate contact conditions, (b) finding constrained equation of motion via a suitable formulation, (c) using a search algorithm to identify points on a boundary that interact with those on another boundary and (d) the insertion of the contact conditions in order to prevent the penetration, and correctly transmit the traction between the bodies [21]. Among different methods for formulating constrained equations of motion and contact conditions, the Lagrange multiplier method is best favored for rubbing contact problem [19].

The first step in implementing the Lagrange multiplier method is to define gap function for the contact problem. The contact conditions ensuring impenetrability of the boundaries can then mathematically be described. For the next step, the contact conditions are enforced to the equation of motion through a variational procedure. Finally, for integration one can use either an implicit or an explicit numerical method.
### 2.4.1. "Node to nodes" contact strategy

In practice, outer boundary of the rotor smoothly moves along the inner boundary of the stator. On the other hand, the boundaries are discretized for obtaining the numerical solution. According to Fig. 2, we propose three methods for discretization of the boundaries:

I. Both boundaries of the rotor and the stator are discretized, as shown in Fig. 2(a). The nodes of the rotor boundary are considered as slave nodes, and element edges of the stator boundary as master segments. The contact conditions should ensure that none of the slave nodes on the rotor boundary penetrates into the master segments of the stator boundary.

II. The rotor boundary is not discretized. The cross section is represented by its center and radius (Fig. 2(b)). The contact conditions need to ensure that the distance between the rotor center as slave node and the master segments of the stator boundary does not exceed the rotor radius.

III. The rotor boundary is not discretized. The cross section is represented by its center and radius. The rotor boundary is allowed to penetrate into the element edges of the stator boundary. The contact conditions need to guarantee that the distance between the rotor center and only the "nodes" of the stator does not exceed the rotor radius.

Method I is computationally time consuming and more complex than Method II. Method III is numerically easier than the other methods since the search algorithm is limited only to determination of which "nodes" of the stator penetrates into the rotor boundary. Neglecting the edges of the stator boundary may lead to a non smooth response if discretization size is large. An acceptable smoothness will be provided if discretization size is small enough. Method III is employed in this work as it treats the contact problem with a quasi "node to node" strategy. In contact mechanism, the rotor first comes into contact with only one node of the stator boundary. As a result, it deforms the stator and can interact with more nodes of the stator depending on the severity of the contact. Because the rotor may contact with several nodes of the stator at any instant, we call this strategy "node to nodes".

<figure><img src="image_3.png" /><figcaption>Fig. 2. Methods for discretization: (a) boundary of the rotor cross section is discretized and (b) rotor cross section represented by its central node.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Gap function for a rigid rotor and deformable stator.</figcaption></figure>
### 2.4.2. Gap function definition

Points C and A in Fig. 3 are the central node of the rotor cross section and one node of the stator, respectively. Let $\mathbf{x}_C$ and $\mathbf{x}_A$ be the current position vectors, and $\mathbf{u}_C$ and $\mathbf{u}_A$ the current displacement vectors of the points C and A, respectively. If C and A are considered to be one nodal contact pair, the corresponding normal gap function can be written as:

$$g_N = \mathbf{n}^T(\mathbf{x}_A - \mathbf{x}_C) \ r, \qquad (4)$$

where $g_N$ is the normal gap function, $\mathbf{n}$ is the unit normal vector pointing to A, and $r$ is the radius of the rotor cross section. Eq. (4) can be written in matrix form as:

$$g_N = \begin{bmatrix} n_x & n_y & n_x & n_y \end{bmatrix} \begin{bmatrix} \mathbf{x}^C \\ \mathbf{y}^C \\ \mathbf{x}^A \\ \mathbf{y}^A \end{bmatrix} = r \begin{bmatrix} \mathbf{G}_L \end{bmatrix}_{1 \times 4} \begin{bmatrix} \mathbf{x}^C \\ \mathbf{y}^C \\ \mathbf{x}^A \\ \mathbf{y}^A \end{bmatrix}, \qquad (5)$$

in which $\mathbf{G}_L$ is called local contact matrix, and $\mathbf{x}$ and $\mathbf{y}$ are the components of the position vector $\mathbf{x}$ in horizontal and vertical directions, respectively. Contribution of all local contact matrices can be assembled into a global contact matrix which contains information of all nodal contact pairs. Thus, if $N_r, N_s$ and $q$ denote number of the rotor axial nodes, number of the stator nodes, and number of activated nodal contact pairs at which contact have occurred, then the gap function can be generalized as:

$$\begin{bmatrix} g_N \\ \vdots \\ g_N^q \end{bmatrix}_{q \times 1} = [\mathbf{G}_N]_{q \times 2(N_r+N_s)} \begin{bmatrix} \mathbf{x}_1^r \\ \mathbf{y}_1^r \\ \vdots \\ \mathbf{x}_1^s \\ \mathbf{y}_1^s \\ \vdots \end{bmatrix}_{2(N_r+N_s) \times 1} \begin{bmatrix} r_1 \\ \vdots \\ r_q \end{bmatrix}_{q \times 1}, \qquad (6)$$

where $g_N^i$ is the normal gap function associated with the ith nodal contact pair, and $\mathbf{G}_N$ is the global normal contact matrix which is assembled as below:

$$[\mathbf{G}_N]_{q \times 2(N_r+N_s)} = [\mathbf{G}_L^1]^T [\mathbf{G}_L^2]^T \dots [\mathbf{G}_L^q]^T, \qquad (7)$$

in Eq. (7), $\mathbf{G}_L^i$ is the ith local contact matrix associated with the ith nodal contact pair. It is evident that $r_i$ is identical for the stator nodes located around an identical cross section of the rotor. Eq. (6) can be rewritten in compact form as:

$$\mathbf{g}_N = \mathbf{G}_N \mathbf{x} \ \mathbf{r} \text{ or } \mathbf{g}_N(\mathbf{u}) = \mathbf{G}_N(\mathbf{x}^0 + \mathbf{u}) \ \mathbf{r}, \qquad (8)$$

where $\mathbf{x}^0$ and $\mathbf{u}$ are the initial position and the current displacement vectors of nodes, respectively. The rotor boundary effect is accounted for by adding $\mathbf{r}$ in the equation of the gap function.

Similarly, one can write for the relative tangential motion between the points C and A along the rotor boundary as follows:

$$\dot{\mathbf{g}}_T = \mathbf{a}^T(\dot{\mathbf{u}}_A - \dot{\mathbf{u}}_C), \qquad (9)$$

where $\dot{\mathbf{g}}_T$ is the relative tangential velocity or time derivative of tangential gap function between the points C and A, $\mathbf{a}$ is the unit tangent vector being perpendicular to $\mathbf{n}$ and in the same direction as the circumferential velocity of the rotor boundary, and $\dot{\mathbf{u}}_A$ and $\dot{\mathbf{u}}_C$ are the velocity vectors of the points C and A, respectively. As shown in Fig. 3, the unit tangent vector is defined by:

$$\mathbf{a} = \Omega \times \mathbf{n}, \qquad (10)$$

where $\Omega$ is the vector of angular velocity of the rotor.

Analogously, Eq. (9) can be rewritten in matrix form as:

$$\dot{\mathbf{g}}_T = \begin{bmatrix} a_x & a_y & a_x & a_y \end{bmatrix} \begin{bmatrix} \dot{\mathbf{u}}_x^C \\ \dot{\mathbf{u}}_y^C \\ \dot{\mathbf{u}}_x^A \\ \dot{\mathbf{u}}_y^A \end{bmatrix}, \qquad (11)$$
the relative tangential velocity for all activated nodal contact pairs can be generalized as follows:

$$ \begin{bmatrix} \dot{g}_T \\ \vdots \\ \dot{g}_T^q \end{bmatrix}_{q \times 1} = [\mathbf{G}_T]_{q \times 2(N_r+N_s)} \begin{bmatrix} \dot{u}_1^{rx} \\ \dot{u}_1^{ry} \\ \vdots \\ \dot{u}_1^{sx} \\ \dot{u}_1^{sy} \\ \vdots \end{bmatrix}_{2(N_r+N_s) \times 1} \quad \text{or} \quad \dot{\mathbf{g}}_T = \mathbf{G}_T \dot{\mathbf{u}}, \qquad (12) $$

where $\mathbf{G}_T$ is called global tangential contact matrix, and $\dot{u}_1^{rx}$ is, for example, the velocity of jth node of the rotor in horizontal direction.

### 2.4.3. Contact conditions

The contact conditions should ensure that: the penetration of the stator nodes into the rotor boundary will not occur; the normal contact force at contact zone will be zero or compressive; the normal contact force will be zero when the gap is positive; and the gap will be zero when the normal force is negative. All these, which are also known as Hertz Signorini Moreau conditions [22], are described as:

$$ \left\{ \begin{array}{ll} g_N^i \ge 0 & \\ p_N^i \le 0 & \text{for } i = 1, 2, \dots, q, \\ p_N^i g_N^i = 0 & \end{array} \right. \qquad (13) $$

where $p_N^i$ is the normal contact pressure associated with the ith nodal contact pair.

In practice, friction always exists. Thus, a model considering the frictional effects must be introduced into the constrained equation of motion. An appropriate model is Coulomb's law which relates the tangential surface traction to the normal pressure as follows:

$$ \|\mathbf{t}_T^i\| \le \mu |p_N^i|, \qquad (14) $$

in which $\mathbf{t}_T^i$ is the tangential surface traction at contact zone, and $\mu$ is the coefficient of friction.

In this work, it is assumed that the rotor has a constant angular velocity, and the energy dissipation due to the friction is compensated by a driven torque. Thus, the rotor boundary steadily slides over the stator boundary and the tangential surface traction will be proportional to the normal pressure as given by:

$$ \mathbf{t}_T^i = \mu |p_N^i| \mathbf{a}, \qquad (15) $$

noting that the normal pressure is to be negative or zero, Eq. (15) is obtained in scalar form as:

$$ \mathbf{t}_T^i = \mu p_N^i, \qquad (16) $$

if the contact area corresponding to the ith nodal contact pairs is denoted by $A^i$, one can define:

$$ \left\{ \begin{array}{ll} \lambda_N^i = p_N^i A^i & \\ \lambda_T^i = t_T^i A^i & \text{for } i = 1, 2, \dots, q, \end{array} \right. \qquad (17) $$

where $\lambda_N^i$ and $\lambda_T^i$ are the normal contact force and the scalar value of the tangential contact force associated with the ith nodal contact pair, respectively. So, Eqs. (13) and (16) can be rewritten in terms of these contact forces as:

$$ \left\{ \begin{array}{ll} g_N^i \ge 0 & \\ \lambda_N^i \le 0 & \text{for } i = 1, 2, \dots, q, \\ \lambda_N^i g_N^i = 0 & \end{array} \right. \qquad (18) $$

and

$$ \lambda_T^i = \mu \lambda_N^i, \qquad (19) $$

respectively.

### 2.4.4. Constrained equation of motion

The Lagrange multiplier method is used to add the contact conditions as constraints to weak form of the equation of motion. As mentioned, a steady slip motion between the bodies is assumed. The variation of the energy related to the contact interface $\Gamma_c$ at which the slip motion takes place is given by:

$$ \delta\Pi_c^{LM} = \int_{\Gamma_c} (p_N \delta g_N + \mathbf{t}_T \delta \mathbf{g}_T) dA + \int_{\Gamma_c} \delta p_N g_N dA, \qquad (20) $$
where the first integral is the virtual work of the contact forces along the variation of gap functions, and the second one enforces the constraints [22]. In the second integral, the tangential gap function does not appear because it is not needed to be zero in a slip motion. It is stressed that the contact pressure $p_N$ solely acts as the Lagrange multiplier because the tangential surface traction is dependent on the normal pressure.

Since the tangential surface traction and the virtual tangential gap are parallel in two dimensional contact problem, the expression $t_T \delta g_T$ in Eq. (20) is simply reduced to the multiplication of the scalar values, i.e. $t_T \delta g_T$. Using node to node contact strategy in FE discretization, Eq. (20) is simplified for the q activated nodal contact pairs as follows [22]:

$$ \int_{t_c} (p_N \delta g_N + t_T \delta g_T) dA + \int_{t_c} \delta p_N g_N dA = \sum_{i}^{q} (p_N^i A^i \delta g_N^i + t_T^i A^i \delta g_T^i) + \sum_{i}^{q} \delta p_N^i A^i g_N^i, \quad (21) $$

using Eqs. (17) and (21), Eq. (20) can be written as:

$$ \delta\Pi_c^{LM} = \sum_{i}^{q} (\lambda_N^i \delta g_N^i + \lambda_T^i \delta g_T^i) + \sum_{i}^{q} \delta\lambda_N^i g_N^i. \quad (22) $$

Hereafter, we regard $\lambda_N^i$ as the Lagrange multiplier associated with the ith nodal contact pair. Eq. (22) can be rewritten in matrix form as:

$$ \delta\Pi_c^{LM} = \delta g_N^T \lambda_N + \delta g_T^T \lambda_T + (\delta\lambda_N)^T g_N, \quad (23) $$

where $\lambda_N = [\lambda_N^1, \lambda_N^2, \dots, \lambda_N^q]^T$, $\lambda_T = [\lambda_T^1, \lambda_T^2, \dots, \lambda_T^q]^T$, and $g_N = [g_N^1, g_N^2, \dots, g_N^q]^T$.

By applying FE discretization procedure to the linear elastic bodies of the rotor and the stator, the potential energy due to static deformation is given by [22]:

$$ \Pi = \frac{1}{2} \mathbf{u}^T \mathbf{K} \mathbf{u} \quad \mathbf{u}^T \mathbf{f}, \quad (24) $$

where $\mathbf{u}$ is the vector of nodal displacements, $\mathbf{K}$ is the stiffness matrix, and $\mathbf{f}$ is the vector of body forces and surface tractions. These variables contain information of both bodies of the rotor and the stator. Using Eqs. (2) and (3), the detailed form of Eq. (24) is:

$$ \Pi = \frac{1}{2} \begin{bmatrix} \mathbf{u}_r \\ \mathbf{u}_s \end{bmatrix}^T \begin{bmatrix} \mathbf{K}_r & 0 \\ 0 & \mathbf{K}_s \end{bmatrix} \begin{bmatrix} \mathbf{u}_r \\ \mathbf{u}_s \end{bmatrix} \begin{bmatrix} \mathbf{u}_r \\ \mathbf{u}_s \end{bmatrix}^T \begin{bmatrix} \mathbf{f}_r \\ \mathbf{f}_s \end{bmatrix}, \quad (25) $$

in which deformations of the rotor and the stator are not yet coupled. Taking variation of Eq. (24) with respect to the displacements gives the weak form of the equilibrium equation:

$$ \delta\Pi = \delta\mathbf{u}^T \mathbf{K} \mathbf{u} \quad \delta\mathbf{u}^T \mathbf{f}, \quad (26) $$

to obtain the constrained equilibrium equation, the variation of the energy related to the contact interface must be added to the variation of the potential energy related to the elastic deformation of the bodies. This means that:

$$ \delta\Pi^{LM} = \delta\Pi + \delta\Pi_c^{LM}, \quad (27) $$

where $\delta\Pi^{LM}$ is the variation of the total energy. Equating Eq. (27) to zero yields the constrained equilibrium equation. Prior to this, $\delta\Pi_c^{LM}$ is required to be expressed in terms of the virtual displacements. Using Eqs. (8) and (12), one can rewrite Eq. (23) as:

$$ \delta\Pi_c^{LM} = \delta\mathbf{u}^T \mathbf{G}_N^T \lambda_N + \delta\mathbf{u}^T \mathbf{G}_T^T \lambda_T + (\delta\lambda_N)^T \mathbf{g}_N = \delta\mathbf{u}^T \mathbf{G}^T \lambda + \delta\lambda_N^T \mathbf{g}_N, \quad (28) $$

where $\lambda^T = [\lambda_N^T \lambda_T^T]$ and $\mathbf{G}^T = [\mathbf{G}_N^T \mathbf{G}_T^T]$.

Inserting Eqs. (26) and (28) into Eq. (27) and then equating the result to zero yields:

$$ \delta\Pi^{LM} = \delta\mathbf{u}^T \mathbf{K} \mathbf{u} \quad \delta\mathbf{u}^T \mathbf{f} + \delta\mathbf{u}^T \mathbf{G}^T \lambda + \delta\lambda_N^T \mathbf{g}_N = 0, \quad (29) $$

considering the contact conditions stated in Eq. (18) and also noting that Eq. (29) must hold for any arbitrary $\delta\mathbf{u}$ and $\delta\lambda_N$, one can obtain the equation:

$$ \left\{ \begin{array}{l} \mathbf{K} \mathbf{u} + \mathbf{G}^T \lambda = \mathbf{f} \\ \mathbf{g}_N \ge 0; \lambda_N \le 0; \lambda_N^i g_N^i = 0 \text{ for } i = 1, \dots, q \end{array} \right. \quad (30) $$

which treats the static deformation of the elastic bodies being in contact.

Considering inertial forces to study the contact dynamics between the bodies of the rotor and the stator, Eq. (30) is rewritten as:

$$ \left\{ \begin{array}{l} \mathbf{M} \mathbf{u} + \mathbf{C} \dot{\mathbf{u}} + \mathbf{K} \mathbf{u} + \mathbf{G}^T \lambda = \mathbf{f} \\ \mathbf{G}_N (\mathbf{x}^0 + \mathbf{u}) \mathbf{r} \ge 0; \lambda_N \le 0; \lambda_N^i g_N^i = 0 \text{ for } i = 1, \dots, q \end{array} \right. \quad (31) $$

where $\mathbf{u}$ and $\dot{\mathbf{u}}$ are the vectors of nodal acceleration and velocity, respectively. Moreover, $\mathbf{M}$ is the mass matrix, $\mathbf{C}$ is the damping matrix, and $\mathbf{f}$ is external force exerted upon the bodies. Eq. (31) is the constrained equation of motion for the rotor and the stator which should be solved by numerical integration.
Recalling Eq. (19), it is evident that $\lambda_T = \mu\lambda_N$. So, the simpler form of Eq. (31) is given by:

$$ \left\{ \begin{array}{l} \mathbf{M}u + \mathbf{C}\dot{u} + \mathbf{K}u + \mathbf{G}_{NT}^T \lambda_N = \mathbf{f} \\ \mathbf{G}_N(x^0 + \mathbf{u}) \mathbf{r} \ge 0; \lambda_N \le 0; \quad \lambda_N^i g_N^i = 0 \text{ for } i = 1, \dots, q \end{array} \right. \qquad (32) $$

in which $G_{NT}$ is here called the normal tangential contact matrix and given by:

$$ G_{NT} = G_N + \mu G_T. \qquad (33) $$

## 2.5. Numerical solution

The numerical solution procedure is adopted from “forward increment Lagrange multipliers method” formulated in [23]. The constrained equation of motion is expressed by:

$$ \left\{ \begin{array}{l} \mathbf{M}u^n + \mathbf{C}\dot{u}^n + \mathbf{K}u^n + (\mathbf{G}_{NT}^{n+1})^T \lambda_N^n = \mathbf{f}^n \\ \mathbf{g}_N(\mathbf{u}^{n+1}) \ge 0. \end{array} \right. \qquad (34) $$

it is initially assumed that all information associated with the time $t^n$ are known, and no contact occurs in the current time step, i.e. $t^n$ to $t^{n+1}$. In other words, $\lambda_N^n$ is assumed to be zero. Thus, based on the finite difference method, the “usual explicit displacement update” at the end of the current time step is calculated as [23,24]:

$$ \mathbf{u}_*^{n+1} = \bar{\mathbf{M}}^{-1} \tilde{\mathbf{f}}^n, \qquad (35) $$

where

$$ \tilde{\mathbf{M}} = \frac{1}{\Delta t^2} \mathbf{M} + \frac{1}{2\Delta t} \mathbf{C}, \quad \tilde{\mathbf{f}}^n = \mathbf{f}^n \left( \mathbf{K} \frac{2}{\Delta t^2} \mathbf{M} \right) \mathbf{u}^n \left( \frac{1}{\Delta t^2} \mathbf{M} \frac{1}{2\Delta t} \mathbf{C} \right) \mathbf{u}^{n-1}, \qquad (36) $$

in Eq. (36), $\Delta t$ is the time step size. Now, occurrence of the penetration should be examined using $\mathbf{u}_*^{n+1}$ and initial positions of the boundaries. If the calculated gap function is positive, no contact has occurred and the total displacement $\mathbf{u}^{n+1}$ is the

<figure><img src="image_5.png" /><figcaption>Fig. 4. Case study of partial rubbing.</figcaption></figure>

<table><caption>**Table 1**

Data for numerical model.</caption><thead><tr><th>Characteristics</th><th>Quantity</th></tr></thead><tbody><tr><td>Shaft elastic modulus</td><td>200 GPa</td></tr><tr><td>Shaft density</td><td>7800 kg/m<sup>3</sup></td></tr><tr><td>Shaft length</td><td>40 cm</td></tr><tr><td>Shaft diameter</td><td>1 cm</td></tr><tr><td>Disk diameter</td><td>12 cm</td></tr><tr><td>Disk density</td><td>7800 kg/m<sup>3</sup></td></tr><tr><td>Disk thickness</td><td>3 cm</td></tr><tr><td>Rod elastic modulus</td><td>3 GPa</td></tr><tr><td>Rod length</td><td>4 cm</td></tr><tr><td>Rod area</td><td>25 mm<sup>2</sup></td></tr><tr><td>Unbalance mass</td><td>0.008 kg</td></tr><tr><td>Radius of unbalance mass</td><td>4.5 cm</td></tr></tbody></table>
same as the displacement update $\mathbf{u}_*^{n+1}$. Otherwise, the Lagrange multipliers vector $\lambda_N^n$ needs to be determined and included into the computations so that the contact condition is satisfied. Finally the total displacement should be corrected as follows:

$$ \mathbf{u}^{n+1} = \mathbf{u}_*^{n+1} + \mathbf{u}_\epsilon, \quad (37) $$

<table><caption>**Table 2**
Non-dimensional parameters.</caption><thead><tr><th>Parameter</th><th>Remarks</th></tr></thead><tbody><tr><td>δ<sub>nd</sub></td><td>Non-dimensional initial clearance; δ initial clearance; |Z<sub>0</sub>| magnitude of steady state response of the rotor without presence of the stator</td></tr><tr><td>ω<sub>n</sub>/ω<sub>0</sub></td><td>Non-dimensional rotational speed; ω<sub>n</sub> angular velocity corresponding to natural frequency of the rotor</td></tr><tr><td>ε</td><td>Damping parameter; modal damping is calculated by c<sub>r</sub> 2ε√m<sub>r</sub>k<sub>r</sub></td></tr><tr><td>ε</td><td>Non-dimensional stiffness of the stator; k<sub>s</sub> linear stiffness of the rod (stator)</td></tr><tr><td>μ</td><td>Coefficient of friction</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Orbits of the rotor center displacement for δ<sub>nd</sub> 0.6, ε 25, ε 0.2, μ 0.1, and (a) (ω/ω<sub>n</sub>) 0.5, (b) (ω/ω<sub>n</sub>) 0.7, (c) (ω/ω<sub>n</sub>) 0.9, (d) (ω/ω<sub>n</sub>) 1.1, (e) (ω/ω<sub>n</sub>) 1.3, (f) (ω/ω<sub>n</sub>) 1.5, (g) (ω/ω<sub>n</sub>) 1.7, (h) (ω/ω<sub>n</sub>) 1.9 and (i) (ω/ω<sub>n</sub>) 2.1: (-) initial clearance δ<sub>nd</sub>.</figcaption></figure>
where $\mathbf{u}_c$ is a corrective incremental displacements produced by the Lagrange multipliers vector $\lambda_N^n$. To find the proper $\lambda_N^n$, the impenetrability condition is used. So, if the initial position vector $\mathbf{x}^0$ is added to Eq. (37), it follows that:

$$ \mathbf{x}^0 + \mathbf{u}^{n+1} = \mathbf{x}^0 + \mathbf{u}_*^{n+1} + \mathbf{u}_c, \qquad (38) $$

multiplying Eq. (38) by the contact matrix $\mathbf{G}_N^{n+1}$, and then subtracting $\mathbf{r}$ from both sides, one has:

$$ \mathbf{G}_N^{n+1}(\mathbf{x}^0 + \mathbf{u}^{n+1}) \mathbf{r} = \mathbf{G}_N^{n+1}(\mathbf{x}^0 + \mathbf{u}_*^{n+1}) \mathbf{r} + \mathbf{G}_N^{n+1}\mathbf{u}_c, \qquad (39) $$

considering impenetrability condition, the gap function has to be zero at the time $t^{n+1}$:

$$ \mathbf{G}_N^{n+1}(\mathbf{x}^0 + \mathbf{u}^{n+1}) \mathbf{r} = \mathbf{g}_N(\mathbf{u}^{n+1}) = 0, \qquad (40) $$

combining Eq. (39) and Eq. (40) results in:

$$ \mathbf{G}_N^{n+1}\mathbf{u}_c = \mathbf{G}_N^{n+1}(\mathbf{x}^0 + \mathbf{u}_*^{n+1}) \mathbf{r}, \qquad (41) $$

similar to Eq. (35), one can write for the corrective displacements vector:

$$ \tilde{\mathbf{M}}\mathbf{u}_c = \tilde{\mathbf{f}}_c^n, \qquad (42) $$

where $\tilde{\mathbf{f}}_c$ is the vector of generalized contact forces. So, Eq. (41) can be described as:

$$ \mathbf{G}_N^{n+1} \tilde{\mathbf{M}}^{-1} \tilde{\mathbf{f}}_c^n = \mathbf{g}_N(\mathbf{u}_*^{n+1}), \qquad (43) $$

<figure><img src="image_7.png" /><figcaption>Fig. 6. Orbits of the rotor center displacement for $\delta_{nd}$: 0.6, ε 25, ξ 0.2, μ 0.1, and (a) $(\omega/\omega_n)$ 2.3, (b) $(\omega/\omega_n)$ 2.5, (c) $(\omega/\omega_n)$ 2.7, (d) $(\omega/\omega_n)$ 2.9, (e) $(\omega/\omega_n)$ 3.1, (f) $(\omega/\omega_n)$ 3.3, (g) $(\omega/\omega_n)$ 3.5, (h) $(\omega/\omega_n)$ 3.7 and (i) $(\omega/\omega_n)$ 4: (-) initial clearance $\delta_{nd}$.</figcaption></figure>
on the other hand, the vector of generalized contact forces is given by:

$$ \tilde{\mathbf{f}}_c^n = (\mathbf{G}_{NT}^{n+1})^T \lambda_N^n, \quad (44) $$

substituting Eq. (44) into Eq. (43) yields:

$$ \mathbf{G}_N^{n+1} \tilde{\mathbf{M}}^{-1} (\mathbf{G}_{NT}^{n+1})^T \lambda_N^n = \mathbf{g}_N (\mathbf{u}_*^{n+1}), \quad (45) $$

so, the Lagrange multipliers can be determined from Eq. (45) as:

$$ \lambda_N^n = (\mathbf{G}_N^{n+1} \tilde{\mathbf{M}}^{-1} (\mathbf{G}_{NT}^{n+1})^T)^{-1} \mathbf{g}_N (\mathbf{u}_*^{n+1}), \quad (46) $$

the total displacement is finally obtained using Eq. (37) and the Lagrange multipliers from Eq. (46):

$$ \mathbf{u}^{n+1} = \mathbf{u}_*^{n+1} \tilde{\mathbf{M}}^{-1} (\mathbf{G}_{NT}^{n+1})^T \lambda_N^n, \quad (47) $$

it should be noted that both $\mathbf{G}_N^{n+1}$ and $\mathbf{G}_{NT}^{n+1}$ which are computed using the explicit position update $\mathbf{x}^0 + \mathbf{u}_*^{n+1}$ are assumed to remain constant within the time step and are updated for the next time step.

### 2.5.1. Solution for elastic stators

If the stator has a massless body without any external forces, a simpler solution can be obtained. First, it is assumed that the elastic stator is “undeformed” at the beginning of the current time step. Then, the usual explicit displacement update for

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Orbits of the rotor center displacement for $\delta_{nd}$ 0.6, $\epsilon$ 25, $\xi$ 0.2, $\mu$ 0.1, and (a) $(\omega/\omega_n)$ 4.3, (b) $(\omega/\omega_n)$ 4.5, (c) $(\omega/\omega_n)$ 4.7, (d) $(\omega/\omega_n)$ 4.9, (e) $(\omega/\omega_n)$ 5.1, (f) $(\omega/\omega_n)$ 5.3, (g) $(\omega/\omega_n)$ 5.5, (h) $(\omega/\omega_n)$ 5.7 and (i) $(\omega/\omega_n)$ 5.9: (-) initial clearance $\delta_{nd}$.</figcaption></figure>
the rotor is calculated by:

$$ \mathbf{u}_{r*}^{n+1} = \tilde{\mathbf{M}}_r^{-1} \tilde{\mathbf{f}}_r^n, \qquad (48) $$

now, in case of penetration, the Lagrange multipliers vector $\lambda_N^n$ should be calculated such that the resulting contact force $\tilde{\mathbf{f}}_c^n$ deforms the stator and, therefore, moves back the rotor up to where the impenetrability condition is satisfied. Actually, the corrective displacements are redefined as follow:

$$ \begin{bmatrix} \tilde{\mathbf{M}}_r & 0 \\ 0 & \mathbf{K}_s \end{bmatrix} \begin{bmatrix} \mathbf{u}_c^r \\ \mathbf{u}_c^s \end{bmatrix} = \tilde{\mathbf{f}}_c^n, \qquad (49) $$

where $\mathbf{u}_c^r$ and $\mathbf{u}_c^s$ are the corrective displacements associated with the rotor and the stator, respectively. Consequently, Eqs. (46) and (47) will be simply applicable by replacing $\tilde{\mathbf{M}}$ by $\begin{bmatrix} \tilde{\mathbf{M}}_r & 0 \\ 0 & \mathbf{K}_s \end{bmatrix}$.

### 3. Simulation of the partial rubbing vibration

The case study for simulation consists of a simple Jeffcott rotor, with two degrees of freedom, and an elastic rod which disturbs the whirling motion of the rotor [2,25], as shown in Fig. 4. The rod is generally located close enough to the rotor in

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Orbits and frequency spectra of vertical component of the rotor center for $\delta_{nd}$ 0.6, $\epsilon$ 25, $\xi$ 0.2, $\mu$ 0.1, and (a) $(\omega/\omega_n)$ 1.7, (b) $(\omega/\omega_n)$ 3.1 and (c) $(\omega/\omega_n)$ 4.9: (-) initial clearance $\delta_{nd}$ (-) the synchronous frequency of the rotational speed.</figcaption></figure>
order that the rubbing is initiated, i.e., the clearance value is chosen less than the magnitude of the steady state vibration displacement.

Considering isotropy and one lateral mode, the equation of motion of the Jeffcott rotor can be written in matrix form as:

$$ \begin{bmatrix} m_r & 0 \\ 0 & m_r \end{bmatrix} \begin{bmatrix} u_{rx} \\ u_{ry} \end{bmatrix} + \begin{bmatrix} c_r & 0 \\ 0 & c_r \end{bmatrix} \begin{bmatrix} \dot{u}_{rx} \\ \dot{u}_{ry} \end{bmatrix} + \begin{bmatrix} k_r & 0 \\ 0 & k_r \end{bmatrix} \begin{bmatrix} u_{rx} \\ u_{ry} \end{bmatrix} = \begin{bmatrix} f_{rx}(t) \\ f_{ry}(t) \end{bmatrix}, \quad (50) $$

where $m_r$, $c_r$ and $k_r$ are the modal mass, modal damping and modal stiffness, respectively. In addition, $u_{rx}$ and $u_{ry}$ are the displacements of the rotor center in horizontal and vertical direction, respectively. The excitation force generated by an unbalance mass is defined as:

$$ \begin{bmatrix} f_{rx}(t) \\ f_{ry}(t) \end{bmatrix} = \begin{bmatrix} m_u e^{\omega^2 \cos(\omega t)} \\ m_u e^{\omega^2 \sin(\omega t)} \end{bmatrix}, \quad (51) $$

where $m_u$, $e$, and $\omega$ are the unbalance mass, the radius of eccentricity of the unbalance mass, and the rotational speed of the rotor, respectively.

The rod is discretized by two dimensional bar element in FE model of the stator. The mass and the damping of the rod are neglected. The geometrical characteristics and the material properties of the model are presented in Table 1.

For analysis, it is convenient to use non-dimensional parameters as introduced in Table 2. If the linear stiffness of the rod is computed and shown by $k_s$, the non-dimensional stiffness of the stator will be calculated $\epsilon \approx 25$ (see Table 2). The damping

<figure><img src="image_10.png" /><figcaption>Fig. 9. Vertical component of the contact force normalized to magnitude of the unbalance force for $\delta_{nd}$ 0.6, $\epsilon$ 25, $\xi$ 0.2, $\mu$ 0.1, and (a) $(\omega/\omega_n)$ 1.7, (b) $(\omega/\omega_n)$ 3.1 and (c) $(\omega/\omega_n)$ 4.9.</figcaption></figure>
parameter is set as high as $\xi=0.2$ in order to avoid chaos phenomena which is beyond the scope of this work. The non-dimensional initial clearance $\delta_{nd}$ is also set to be less than unity at all rotational speeds for the rubbing initiation and continuation. The coefficient of friction is also assumed to be $\mu = 0.1$.

In numerical computation, a proper time step size should be chosen for the stability. The time step size needs to be smaller than $min((T_{min}/\pi), (2/\omega))$ where $T_{min}$ is the smallest period of frequencies in the rotor stator system, and $\omega$ is the frequency of excitation force [24]. In practice, the superharmonics of exciting frequency are present in the response due to the impacts. Thus, we set the step size to $min((T_{min}/32n), (2\pi/32n\omega))$ in order to obtain the response with frequency content up to $max((2n\pi/T_{min}), (n\omega))$.

# 4. General results and discussions

In this section, the synchronous and subsynchronous vibration responses of the partial rubbing model are presented with respect to change in the rotational speed. Furthermore, a sensitivity analysis with respect to the initial clearance, the stator stiffness, the damping parameter, and the coefficient of friction is performed. The response is shown in form of orbit and frequency spectrum. Both horizontal and vertical components of the rotor center displacement are normalized to $|Z_0|$. All frequencies are also normalized to the natural frequency of the rotor. The transient part of the response is cut off before drawing the orbits and performing the frequency analysis. For obtaining a Fast Fourier Transform (FFT), a Hanning window was employed to prevent leakage problem.

## 4.1. Rotor speed effect on the response

We intend to investigate $1 \times$ synchronous, $1/2 \times$, and $1/3 \times$ subsynchronous vibrations regimes of the response. For this purpose, the rotational speed of the rotor is set to change from $\omega/\omega_n = 0.5$ to 6.0.

The orbits of the rotor center motion with respect to change in the rotational speed from $\omega/\omega_n = 0.5$ to 2.1 with increment of 0.2 are shown in Fig. 5. At the rotational speed $\omega/\omega_n = 0.5$, according to Fig. 5(a), the rubbing contact occurs twice per revolution of the rotor because the orbit exceeds the initial clearance twice. This rotational speed is low enough that the unbalance force pushes the rotor towards the stator twice within one rotor revolution. It is evident that several contacts may occur per revolution for rotational speeds below $\omega/\omega_n = 0.5$. For the speeds above $\omega/\omega_n = 0.5$, as can be seen

<figure><img src="image_11.png" /><figcaption>Fig. 10. Stress and strain at the stator for $\delta_{nd}$ 0.6, $\varepsilon$ 25, $\xi$ 0.2, $\mu$ 0.1, and (a) $(\omega/\omega_n)$ 1.7, (b) $(\omega/\omega_n)$ 3.1 and (c) $(\omega/\omega_n)$ 4.9.</figcaption></figure>
in Fig. 5, all orbits are single loop and exceeds the initial clearance once per revolution. It means that the rubbing contact
occurs once per revolution of the rotor. The shape of the orbit is first an oblique oval. By increasing the rotational speed, the
orbit becomes sharper and orients more vertically. This is due to increase in the phase delay of the response lagging
the unbalance excitation force. Actually, the vertical component of the contact force is intensified by downward effect of the
rotating unbalance force. It can also be observed that the size of orbit first grows and then decreases while the rotational
speed increases continuously. This is due to passing a resonance speed about ω/ωn = 1.7. In the range of ω/ωn = 0.5 to 2.1,
the vibration is said to be 1× synchronous because the motion repeats itself once per revolution of the rotor even though
two or more contact occur at low rotational speeds.

The orbits of the rotor center motion for the rotational speed range of ω/ωn = 2.3 to 4.0 in increment of 0.2 are depicted in Fig. 6. In this range, the orbits become double loops. It is obvious that the rubbing motion repeats itself per two revolutions of the rotor. It means that the vibration regime is here 1/2× subsynchronous, and the rubbing contact occurs once per two rotor revolutions. The orbit, as can be seen, exceeds the initial clearance only once. It is noted that irregularities such as chaos phenomena (see Fig. 6(b)) are observed in the beginning of the speed range during transition from 1× synchronous vibration regime to 1/2× subsynchronous one. The size of the orbit increases and then decreases while the rotational speed increases. This is due to passing through another resonance about ω/ωn = 3.1. It can be concluded that a "quasi resonance speed" exists within the rotational speed range associated with 1/2× subsynchronous vibration regime (i.e. ω/ωn = 2.3 to ω/ωn = 4.0).

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Spectrum cascade plot of vertical component of the rotor center displacement for ε = 25, ξ = 0.2, μ = 0.1, (ω/ωn) = 0.5 to 5.0 in increment of 0.1, and (a) δnd = 0.9, (b) δnd = 0.7 and (c) δnd = 0.5: (-) 1× synchronous, (---) 1/2× subsynchronous and (...) 1/3× subsynchronous vibration lines.</figcaption></figure>
The orbits of the rotor center motion for the rotational speed range $\omega/\omega_n = 4.3$ to 5.9 in increment of 0.2 are presented in Fig. 7. In this range, the orbits are triple loops. So, the rubbing motion repeats itself per three revolutions of the rotor. It means that the vibration regime is $1/3 \times$ subsynchronous, and the rubbing contact mainly occurs once per three revolution of the rotor. Similarly, a transition accompanied with some irregularities from $1/2 \times$ subsynchronous vibration regime to $1/3 \times$ one happens in the beginning of the speed range. Grow and decay in size of the orbits is also observed here. This is due to passing through another resonance about $\omega/\omega_n = 4.9$. Therefore, there is another "quasi resonance speed" within the rotational speed range associated with $1/3 \times$ subsynchronous vibration regime (i.e. from $\omega/\omega_n = 4.3$ to $\omega/\omega_n = 5.9$).

For further investigation, FFT analysis is performed on the vertical component of the rotor center displacement at each rotational speed. The increment of the rotational speed is set to be $\omega/\omega_n = 0.1$. All frequency spectra corresponding to all rotational speed are assembled into a plot called "spectrum cascade". The spectrum cascade plot for the response of the model is given in Fig. 12(b). The figure shows that the $1 \times$ synchronous component is present in all rotational speeds due to the unbalance synchronous excitation. However, it significantly appears at the rotational speed range associated with the $1 \times$ synchronous vibration regime due to the rubbing contact effects. The axis of the rotational speed in the plots can be divided to successive distinctive zones each of them representing a fractional subsynchronous vibration regime. The $1 \times$ synchronous, $1/2 \times$, and $1/3 \times$ subsynchronous vibration regimes, the transition between regimes, the irregularities during the transition, and the relative peaks denoting the quasi resonance speeds are clearly illustrated in the spectrum cascade plot.

The orbits and frequency spectra at the three mentioned resonance speeds are shown in Fig. 8. It is observed that both frequencies of $1/2 \times$ and $1/3 \times$ subsynchronous components at second and third resonance speeds, respectively, are close to that of the $1 \times$ synchronous component at the first resonance speed depicted in Fig. 8(a). Investigation on the spectrum cascade plot reveals that the lowest dominant frequency in the synchronous and subsynchronous regimes does not tend to go much beyond the frequency of $1 \times$ synchronous component at the first resonance speed. In fact, successive transitions in the vibration regimes and, consequently, sudden shift downs of the lowest dominant frequency happen when the rotational speed increases (see Fig. 12(b)).

The vertical component of the contact force which is normalized to the magnitude of the unbalance force for the three resonance speeds is presented in Fig. 9. The normalized contact force reaches its steady state value after several contacts. The normalized contact force reaches its maximum value at the first resonance speed, so it has smaller values at other resonance speeds corresponding to the $1/2 \times$ and $1/3 \times$ subsynchronous vibration regimes.

<figure><img src="image_13.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Spectrum cascade plot of vertical component of the rotor center displacement for $\delta_{nd}$ 0.6, $\xi$ 0.2, $\mu$ 0.1, $(\omega/\omega_n)$ 0.5 to 6.0 (5.0 for (a)) in increment of 0.1, and (a) $\epsilon$ 5, (b) $\epsilon$ 25, (c) $\epsilon$ 50 and (d) $\epsilon$ 100: (-) $1 \times$ synchronous, (- -) $1/2 \times$ subsynchronous and (...) $1/3 \times$ subsynchronous vibration lines.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Spectrum cascade plot of vertical component of the rotor center displacement for $\delta_{nd}$ 0.6, $\epsilon$ 25, $\mu$ 0.1, $(\omega/\omega_n)$ 0.5 to 6.0 in increment of 0.1, and (a) $\xi$ 0.1, (b) $\xi$ 0.2 and (c) $\xi$ 0.3: (-) 1× synchronous (---), 1/2× subsynchronous and (...) 1/3× subsynchronous vibration lines.</figcaption></figure>

The stress and strain in the stator are presented in Fig. 10. Unlike the normalized contact force, the stress and strain have greater value at next resonance speed. Thus, higher rotational speed intensifies the rubbing condition considering the contact forces and the subsequent stresses and strains. It is worth to note that the strain range verifies the assumption of linear behavior. The results are obtained by keeping the non dimensional initial clearance and other parameters constant.

## 4.2. Sensitivity analysis

The spectrum cascade plots of the response for the three initial clearances are shown in Fig. 11. By decreasing the non dimensional initial clearance, the curves in the spectra and relative peaks representing resonance speeds shift towards higher rotational speeds. In fact, reducing the clearance increases the potential engagement of the rotor and stator and thus effective stiffness of the system. The transitions also occur at higher rotational speeds while the size of synchronous and subsynchronous components remains unchanged.

The spectrum cascade plots of the response for different stator stiffnesses are presented in Fig. 12. By increasing the stator stiffness, the curves in the spectra, the transition speeds, and the relative peaks representing the resonance speeds will shift toward higher rotational speeds. This is evidently due to increase in total stiffness of the system. In addition, more irregularities appear during the transition. The size of fractional subsynchronous components grows without any considerable change in the synchronous component.

The spectrum cascade plots of the response for some damping parameters are depicted in Fig. 13. One can see that when the damping parameter increases, the curves in the spectra, the transition speeds, and the relative peaks representing resonance speeds shift towards the lower rotational speeds. The irregularities associated with the transitions also disappear
<figure><img src="image_18.png" /></figure>

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

<figure><img src="image_20.png" /><figcaption>**Fig. 14.** Spectrum cascade plot of horizontal component of the rotor center displacement for $\delta_{nd}$ 0.6, $\varepsilon$ 25, $\xi$ 0.2, $(\omega/\omega_n)$ 0.5 to 6.0 in increment of 0.1, and (a) $\mu$ 0.1, (b) $\mu$ 0.3 and (c) $\mu$ 0.5: (-) 1× synchronous, (---) 1/2× subsynchronous and (...) 1/3× subsynchronous vibration lines.</figcaption></figure>

progressively. Furthermore, the size of synchronous and fractional subsynchronous components significantly reduces because of the damping effect.

The spectrum cascade plots of the horizontal component of the rotor center displacement for several coefficients of friction are shown in Fig. 14. When the coefficient of friction increases, the transition speeds, and the rotational speeds of the relative peaks, and the irregularities do not change. However, the fractional subsynchronous components significantly grow because higher friction intensifies the horizontal contact forces. The orbits of the rotor center motion for different coefficient of friction are presented in Fig. 15. By increasing the coefficient of friction, the shape of the orbit does not change but it turns in the same wise the rotor rotates.

Considering the presented results it can be concluded that change in rub related parameters significantly affects the rubbing phenomena. It actually depends on the subsequent shift of the curves in the spectrum cascades towards higher rotational speeds or lower ones, and also on whether the rotor speed is under or above the nearest critical speed of the rotor. For example, if the change leads to a shift of the curves towards the higher rotational speeds, and the rotor speed is under the resonance speed, the rubbing condition will be alleviated.

### 4.3. Comparison of the results with the literature

The obtained numerical results is compared with experimental data published by Muszynska [2]. In the experimental study, a rubbing test was performed on a rotor rig composed of a steel rotor and a brass screw mounted vertically in a stationary holder. Fig. 16 shows the results of the experimental rubbing responses which were presented in the format of spectrum cascade plot together with samples of the rotor orbits. The orbits corresponds to 1× synchronous and 1/2× and 1/3× subsynchronous vibrations regimes.
<figure><img src="image_21.png" /><figcaption>**Fig. 15.** Orbits of the rotor center motion for $\delta_{nd}$ 0.6, $\xi$ 0.2, $\epsilon$ 25, $\mu$ {0.1, 0.3, 0.5} (increasingly left to right), and (a) $(\omega/\omega_n)$ 1.7, (b) $(\omega/\omega_n)$ 2.8 and (c) $(\omega/\omega_n)$ 4.9.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 16.** Rotor response spectrum cascade plot and associated orbits for an experimental setup designed for the partial rubbing investigation [2].</figcaption></figure>
By comparing these orbits with those presented in Fig. 8, and the experimental spectrum cascade plot with the numerical one in Fig. 12(a) or (b), one can conclude that there is an excellent agreement between these results.

## 5. Conclusions

In this paper, an algorithm is developed for investigating the rotor to stator rubbing vibration, based on FE theory with unilateral contact and friction conditions. This algorithm can describe the rubbing contact between a circular rigid cross section of a rotor and any deformable stator with any arbitrary geometry. Different methods are proposed for discretizing the boundaries of the rotor and the stator. In the best chosen method, only the boundary of the stator is discretized. The contact problem is defined as the contact between the circular cross section of the rotor and only "nodes" of the stator. Next the constrained equation of motion is obtained by coupling the finite element models of the rotor and stator via the Lagrange multiplier technique and solved by numerical computations.

By using this algorithm, partial rubbing is simulated for a sample rotor stator system. The synchronous and subsynchronous responses are obtained for different rotational speeds. Sensitivity analysis is carried out considering parameters such as initial clearance, stator stiffness, damping parameter, and the coefficient of friction. A qualitative verification is performed on the numerical results using experimental data in the literature, and excellent agreement is found.

This algorithm is able to more realistically model the stationary parts of turbomachinery such as different seals and bearings.

## Acknowledgment

The authors are thankful to Iranian National Gas Transmission Company for providing financial and technical support for this research.

## References

[1] A.R. Bartha, *Dry Friction Backward Whirl of Rotor*, ETH, Zurich, Switzerland, 2000 (Ph.D. thesis).

[2] A. Muszynska, *Rotordynamics*, CRC Taylor and Francis Group, New York, 2005.

[3] F.K. Choy, J. Padovan, Non-linear transient analysis of rotor-casing rub events, *Journal of Sound and Vibration* 113 (1987) 529-545.

[4] B.O. Al-bedoor, Transient torsional and lateral vibrations of unbalanced rotors with rotor-to-stator rubbing, *Journal of Sound and Vibration* 229 (2000) 627-645.

[5] D.W. Childs, Rub-induced parametric excitation in rotors, *Journal of Mechanical Design, Transactions of the ASME* 101 (1979) 640-644.

[6] F. Chu, W. Lu, Determination of the rubbing location in a multi-disk rotor system by means of dynamic stiffness identification, *Journal of Sound and Vibration* 248 (2001) 235-246.

[7] F. Chu, W. Lu, Stiffening effect of the rotor during the rotor-to-stator rub in a rotating machine, *Journal of Sound and Vibration* 308 (2007) 758-766.

[8] X. Dai, X. Zhang, X. Jin, The partial and full rubbing of a flywheel rotor-bearing-stop system, *International Journal of Mechanical Sciences* 43 (2001) 505-519.

[9] F. Ehrich, Observations of subcritical superharmonic and chaotic response in rotordynamics, *Journal of Vibration, Acoustics, Stress, and Reliability in Design* 114 (1992) 93-100.

[10] Z.C. Feng, X.Z. Zhang, Rubbing phenomena in rotor-stator contact, *Chaos, Solitons and Fractals* 14 (2002) 257-267.

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

[12] E.E. Pavlovskaia, E.V. Karpenko, M. Wiercigroch, Non-linear dynamic interactions of a Jeffcott rotor with preloaded snubber ring, *Journal of Sound and Vibration* 276 (2004) 361-379.

[13] Z.K. Peng, F.L. Chu, P.W. Tse, Detection of the rubbing-caused impacts for rotor-stator fault diagnosis using reassigned scalogram, *Mechanical Systems and Signal Processing* 19 (2005) 391-409.

[14] P. Pennacchi, N. Bachschmid, E. Tanzi, Light and short arc rubs in rotating machines: experimental tests and modelling, *Mechanical Systems and Signal Processing* 23 (2009) 2205-2227.

[15] W. Qin, G. Chen, G. Meng, Nonlinear responses of a rub-impact overhung rotor, *Chaos, Solitons and Fractals* 19 (2004) 1161-1172.

[16] G. Von Groll, D.J. Ewins, The harmonic balance method with arc-length continuation in rotor/stator contact problems, *Journal of Sound and Vibration* 241 (2001) 223-233.

[17] J.J. Yu, P. Goldman, D.E. Bently, A. Muszynska, Rotor/seal experimental and analytical study on full annular rub, *Journal of Engineering for Gas Turbines and Power* 124 (2002) 340-350.

[18] W.-M. Zhang, G. Meng, Stability, bifurcation and chaos of a high-speed rub-impact rotor system in MEMS, *Sensors and Actuators A* 127 (2006) 163-178.

[19] S. Roques, M. Legrand, P. Cartraud, C. Stoisser, C. Pierre, Modeling of a rotor speed transient response with radial rubbing, *Journal of Sound and Vibration* 329 (2010) 527-546.

[20] M. Lalanne, G. Ferraris, *Rotordynamics Prediction in Engineering*, Wiley, New York, 1998.

[21] O.C. Zienkiewicz, R.L. Taylor, *The Finite Element Method for Solid and Structural Mechanics*, sixth ed. Butterworth-Heinemann, Oxford, 2005.

[22] P. Wriggers, *Computational Contact Mechanics*, Wiley, New York, 2002.

[23] N.J. Carpenter, R.L. Taylor, M.G. Katona, Lagrange constraints for transient finite element surface contact, *International Journal for Numerical Methods in Engineering* 32 (1991) 103-128.

[24] K.J. Bathe, *Finite Element Procedure*, Prentice-Hall, New Jersey, 1996.

[25] Y.S. Choi, On the contact of partial rotor rub with experimental observations, *KSME International Journal* 15 (2001) 1630-1638.