This paper is meant to model free vibration of a coupled rotor-bearing-housing system. In particular, the rotor is cyclic symmetric and spins at constant speed while the housing is stationary and flexible. The rotor and housing are assembled via multiple, linear, elastic bearings. A set of equations of motion is derived using component mode synthesis, in which the rotor and the housing each are treated as a component. The equations of motion take the form of ordinary differential equations with periodic coefficients. Analyses of the equations of motion indicate that instabilities could appear at certain spin speed in the form of combination resonances of the sum type. To demonstrate the validity of the formulation, two numerical examples are studied. For the first example, the spinning rotor is an axisymmetric disk, and the housing is a square plate with a central shaft. The rotor and the housing are connected via two linear elastic bearings. For the second example, the rotor is cyclic symmetric in the form of a disk with four evenly spaced radial slots. The housing and bearings remain the same. In both examples, instability appears as a combination resonance of the sum type between a rotor mode and an elastic housing mode. The cyclic symmetric rotor, however, has more instability zones. Finally, effects of damping are studied. Damping of the housing widens the instability zones, whereas the damping of the rotor does the opposite.

## Introduction

A spinning cyclic symmetric rotor mounted on a stationary housing via multiple bearings is a very common platform used in modern rotary machinery. Representative examples include propellers, wind turbines, bladed turbine disks, and compressors. Traditionally, the housing is massive and rigid, so vibration tends to be more pronounced at the rotor. Therefore, most literature on rotor dynamics focuses on vibration analyses of spinning rotors.

In the last decade, emergent applications started to call for designs that employ a lighter housing or a larger rotor. A natural consequence is that vibration of the rotor and housing becomes coupled through the bearings. As an example, wind turbines now adopt long, huge blades, whereas the housing is reduced in weight to lower material costs. Research has shown that wind turbines experience considerable coupled vibration between the blades and the tower [1–5]. Similarly, in developing greener engines with higher fuel efficiency and less emission, the turbine engine industry seeks to reduce the weight (and hence rigidity) and to increase the bypass ratio of turbine engines. This design trend tends to cause considerable coupled vibration between the rotating parts and the stationary parts and motivates vibration analyses of entire rotor-bearing-housing systems [6,7].

Although these early publications have attempted to model a rotor-bearing-housing system, the rotor is assumed to take a specific geometry (e.g., a beam). A universal model for cyclic symmetric rotors and housing with arbitrary geometry remains open. Moreover, characteristic responses of cyclic symmetric rotor-bearing-housing systems have not been thoroughly studied and understood.

Motivated by industrial needs and academic challenges, this paper is meant to develop a formulation that accommodates arbitrary geometry for the spinning rotor and the stationary housing. Such a formulation is universal and will be valid for various cyclic symmetric rotors, ranging from wind turbines to bladed turbine disks. Characteristics drawn from this model will then be applied to any cyclic symmetric rotor mounted on a flexible housing.

Developing such a universal formulation, however, is a nontrivial task. There are two major challenges to overcome. The first challenge is that traditional vibration analyses of cyclic symmetric rotors assume no bearings or housing [8–19]. A systematic way to incorporate a stationary housing and multiple bearings to a spinning cyclic symmetric rotor is yet to be developed. The only knowledge thus far in this context is that the presence of the bearings may cause some rotor modes to be coupled with the bearings while other modes may not [20]. More specifically, if a rotor vibration mode presents unbalanced inertia forces or moments as the rotor vibrates, the rotor vibration mode will couple with the housing. These vibration modes are known as unbalanced modes [20].

The second challenge is that equations of motion could possess extremely high order if a rotor-bearing-housing system is not properly formulated. Theoretically, such a system is governed by partial differential equations with time-varying coefficients because the rotor's mass and stiffness distributions change periodically with respect to the housing [3–5]. If a finite element method is used to discretize the system, the order will be so high that numerical solutions become impractical. Therefore, novel ways must be sought to significantly reduce the order to make numerical predictions feasible.

To overcome these challenges, we employ two coordinate systems in the formulation. A rotor-based coordinate system is used for the spinning rotor to accommodate its cyclic symmetry [16–19]. A ground-based coordinate is used to describe the response of the stationary housing [21]. The rotor and the housing are then discretized in respective coordinate systems, for example via mode shapes, to reduce the order. Finally, the rotor, housing, and bearings are assembled together to derive the equation of motion with a reduced order.

Specifically, derivation of the equation of motion with a reduced order consists of the following steps. The first step is to extract mode shapes of the rotors and housings. Since the two coordinates are fixed to the rotor and the housing, respectively, the rotors and housings could be considered stationary in the respective coordinates. Then a finite element analysis is conducted on the rotor and housing to obtain their natural frequencies and mode shapes. Since the mode shapes are orthogonal in respective coordinate systems, they form two complete sets of basis.

The second step is discretization through use of the two coordinate systems. Since the mode shapes are complete, the rotor-based response of the spinning rotor can be represented in terms of its mode shapes and their modal response (also known as generalized coordinates). Likewise, ground-based response of the stationary housing can be represented in terms of its mode shapes and their modal response. Therefore, kinetic and potential energies of the rotor and housing can be discretized in terms of these two sets of generalized coordinates. Similarly, bearing deformation resulting from relative displacements between the rotor and housing can also be discretized. (Since two coordinate systems are present, the relative displacement is only meaningful after a coordinate transformation is performed.) Through the use of the Lagrange equation, a set of equations of motion with periodic coefficients is obtained, where the periodic coefficients come from the bearing deformation and the coordinate transformation. Solution of the equation of motion can be obtained via the method of multiple scales or numerical integrations.

For the rest of the paper, we will first derive the equation of motion with a reduced order following the outline described above. Next, we will have a brief discussion on the convergence rate of this formulation. Then, we will use the equations of motion to conduct two numerical examples. The first example has an axisymmetric disk as the spinning rotor, while the second example has a cyclic symmetric rotor in the form of a circular disk with evenly spaced radial slots. In both examples, the rotor is supported on a stationary elastic plate (i.e., the housing) via two elastic bearings. We then analyze stability of the rotor-bearing-housing system as a function of the spin speed to study how cyclic symmetry affects the stability. Finally, we demonstrate how damping affects the stability of a coupled rotor-bearing-housing system.

## Mathematical Formulations

This section includes formulations of the stationary housing, spinning rotor, bearing deformations, generalized bearing forces, and the equations of motion. In addition, methods of solving the equations of motion are briefly discussed.

### Formulation of the Stationary Housing.

The stationary housing is modeled as an elastic structure with arbitrary shapes. Let $O$ be a convenient reference point on the housing. With $O$ being the origin, one can define an inertia frame $XYZ$ with unit vectors $I$, $J$, and $K$. Now consider an arbitrary point $P$ in the housing with position vector
$r∧≡r∧xI+r∧yJ+r∧zK$
(1)
When the housing undergoes an elastic deformation $W(H)(r∧,t)$, where the superscript $(H)$ refers to the stationary housing, it can be approximated in terms of $nh$ vibration modes via
$W(H)(r∧,t)≈∑n=1nhWn(H)(r∧)qn(H)(t)$
(2)
In Eq. (2),
$Wn(H)(r∧)≡Wxn(H)(r∧)I+Wyn(H)(r∧)J+Wzn(H)(r∧)K$
(3)
is the $n$-th vibration mode shape, and $qn(H)(t)$ is the corresponding generalized coordinate. In addition, the mode shapes satisfy the orthonormality conditions $∫Wm(H)(r∧)·Wn(H)(r∧)dm=δmn$ and $∫VH[Wm(H)(r∧),Wn(H)(r∧)]dV=[ωn(H)]2δmn$, where $δmn$ is the Kronecker δ, $ωn(H)$ is the natural frequency of the $n$-th vibration mode, and $VH[•,•]$ is a potential energy inner-product operator. Since the housing is stationary, the motion of point $P$ results entirely from the elastic deformation $W(H)(r∧,t)$. Hence, the displacement of $P$ is
$RP≡W(H)(r∧,t)$
(4)
and the velocity of point $P$ is
$R·P=∑n=1nhWn(H)(r∧)q·n(H)(t)$
(5)

where Eq. (2) has been used. Moreover, the kinetic energy of the housing is $T(H)=12∑n=1nh[q·n(H)]2$, and the potential energy of the housing is $V(H)=12∫Vh[W(H),W(H)]dV=12∑n=1nh[ωn(H)qn(H)(t)]2$.

### Formulation of the Rotor.

The rotor is modeled as an elastic solid of arbitrary geometry with its centroid $G$ located on the spin axis (i.e., rotor is perfectly balanced). Moreover, the rotor is subject to free-body boundary conditions and can undergo rigid-body motion as well as elastic deformation. To formulate the motion of the rotor, let us first define a virgin state for reference. In the virgin state, the rotor-bearing-housing system experiences no elastic deformation and no rigid-body motion except the spin of the rotor. Moreover, let $G'$ denote the location of the centroid in the virgin state; see Fig. 1. (Note that $G$ and $G'$ coincide in the virgin state but will correspond to two different points in actual motion.) With $G'$ as the origin, one can define an inertia Cartesian coordinate system $X˜Y˜Z˜$ with $Z˜$ being the spin axis in the virgin state. (Note that the rotor might wobble in actual motion; therefore, $Z˜$ is not the spin axis in actual motion.) Also, one can define a rotating coordinate system $xyz$ with constant angular velocity $ω3$. It is also convenient to assume that $xyz$ coincides with the principal axes of the rotor in the virgin state. Note that $xyz$ axes do not attach to the rotor because the rotor can rock and translate, but $xyz$ axes cannot. Nonetheless, $xyz$ axes are the best coordinate system to describe the motion of the rotor. Let $i$, $j$, and $k$ be the unit vectors of the $xyz$ coordinate system and $I$, $J$, and $K$ be the unit vectors of the inertia frame $XYZ$ (or $X˜Y˜Z˜$) coordinate system. Then

Fig. 1
Fig. 1
Close modal
$(ijk)=[cosω3tsinω3t0-sinω3tcosω3t0001](IJK)$
(6)
Now let us consider the motion of a point $P$ on the rotor. In the virgin state, position of $P$ relative to $G'$ is defined by vector $r˜$ (Fig. 1), where $r˜≡r˜xi+r˜yj+r˜zk$. In the deformed state, the position of $P$ with respect to the inertia frame $XYZ$ becomes
$RP(R)(r˜,t)=OG'→+r˜+W(R)(r˜,t)$
(7)

where $W(R)(r˜,t)$ is the displacement of $P$ that includes rigid-body motion (excluding the spin) and elastic deformation. The superscript $(R)$ refers to the rotor.

Under the coordinate system $xyz$, several sets of complete mode shapes are available. For example, when the rotor is stationary, undamped, and subjected to free boundary conditions, the rotor will have a set of complete mode shapes
$Wn(R)(r˜)≡Wxn(R)(r˜)i+Wyn(R)(r˜)j+Wzn(R)(r˜)k$
(8)
with natural frequencies $ωn(R)$, $n=0,1,2,…$. (For rotors with complicated geometry, $Wn(R)(r˜)$ and $ωn(R)$ can be calculated through finite element analyses; for example.) Moreover, these mode shapes are orthogonal satisfying $∫Wm(R)(r˜)·Wn(R)(r˜)dm=δmn$ and $∫VR[Wm(R)(r˜),Wn(R)(r˜)]dV=[ωn(R)]2δmn$, where $VR[•,•]$ is a potential energy inner-product operator of the rotor. Also note that the mode shapes $Wn(R)(r˜)$ include six infinitesimal rigid-body modes with zero natural frequencies. Let us label the zero-th mode as the infinitesimal rigid-body spin of the rotor, i.e.,
$W0(R)(r˜)≡1I¯3(-y˜i+x˜j)$
(9)

where $I¯3$ is the centroidal mass moment of inertia about the spin axis.

Since the vibration mode shapes are complete, the position vector $RP(R)(r˜,t)$ of the rotor in Eq. (7) can be approximated as
$RP(R)(r˜,t)≈OG'→+r˜+∑n=1nrWn(R)(r˜)qn(R)(t)$
(10)
where $qn(R)(t)$ is the generalized coordinates whose response is to be determined, and $nr$ is the number of modes retained in the series for approximation. Note that the rigid-body spin of the rotor $W0(R)(r˜)$ does not appear in $RP(R)(r˜,t)$ because the spin of the rotor is redundant and has been described by the rotation of the coordinate system $xyz$. With the displacement in Eq. (10), the potential energy of the rotor is
$V(R)=12∫VR[W(R),W(R)]dV=12∑n=1nr[ωn(R)qn(R)(t)]2$
(11)
According to Eq. (10), the velocity of point $P$ is
$R·P(R)(r˜,t)=ω×r˜+∑n=1nrWn(R)(r˜)q·n(R)(t)+∑n=1nr[ω×Wn(R)(r˜)]qn(R)(t)$
(12)
where $ω=ω3k$. With the velocity field in Eq. (12), the kinetic energy of the rotor is
$T(R)=12I¯3ω32+12∑n=1nr[q·n(R)]2+12ω32∑n=1nr∑m=1nrλmnqm(R)(t)qn(R)(t)+ω3∑n=1nrJbn(R)q·n(R)(t)+ω32∑n=1nrJcn(R)qn(R)(t)+ω3∑n=1nr∑m=1nrgnmqm(R)(t)q·n(R)(t)(13)$
(13)
where
$λmn=λnm≡δmn-∫Wmz(R)(r˜)Wnz(R)(r˜)dm$
(14)
$Jbn(R)≡k·∫r˜×Wn(R)(r˜)dm$
(15)
$Jcn(R)≡∫[r˜·Wn(R)(r˜)-r˜zWnz(R)(r˜)]dm$
(16)
$gnm=-gmn≡k·∫[Wm(R)(r˜)×Wn(R)(r˜)]dm$
(17)

For Eqs. (14)–(17), the integration is carried out with respect to the rotor.

### Bearing Deformations.

The purpose of this section is to derive bearing deformations in terms of $qn(H)(t)$ and $qn(R)(t)$. Let $A$ and $A'$ be two mating surfaces of a bearing (e.g., inner race and outer race), where $A$ is on the rotor and $A'$ is on the stationary housing. When the spindle vibrates, the two bearing surfaces move relatively to each other. In addition, the relative motion could be linear or angular.

To determine the linear bearing deformation, let $rA$ and $rA'$ be the position vector of $A$ and $A'$, respectively. Then $rA-rA'$ defines the linear bearing deformation. According to Fig. 2 and Eq. (10),

Fig. 2
Fig. 2
Close modal
$rA=OG'→+G'A→+∑n=1nrWn(R)(r˜A)qn(R)(t)$
(18)
According to Eq. (4),
$rA'=∑n=1nhWn(H)(r∧A)qn(H)(t)$
(19)
Based on Eqs. (18) and (19), the linear bearing deformation is $rA-rA'≡ΔxI+ΔyJ+ΔzK$, where
$Δx≡∑n=1nr[Wxn(R)(r˜A)cosω3t-Wyn(R)(r˜A)sinω3t]qn(R)(t)-∑n=1nhWxn(H)(r∧A)qn(H)(t)$
(20)
$Δy≡∑n=1nr[Wxn(R)(r˜A)sinω3t+Wyn(R)(r˜A)cosω3t]qn(R)(t)-∑n=1nhWyn(H)(r∧A)qn(H)(t)$
(21)
and
$Δz≡∑n=1nrWzn(R)(r˜A)qn(R)(t)-∑n=1nhWzn(H)(r∧A)qn(H)(t)$
(22)

Note that the coordinate transformation in Eq. (6) is used because $Wn(R)(r˜A)$ in Eq. (18) is represented in coordinates $xyz$ and must be transformed into the inertia frame $XYZ$.

For the angular deformation of the bearing, the elastic deformation of the stationary housing will induce an infinitesimal rigid-body rotation
$α≡αxI+αyJ+αzK=12∇×∑n=1nhWn(H)(r∧A)qn(H)(t)$
(23)
where $r∧A$ is the position vector of $A$ in the $XYZ$ frame. Similarly, the rotor will experience an infinitesimal rigid-body rotation
$β≡βxi+βyj+βzk=12∇×∑n=1nrWn(R)(r˜A)qn(R)(t)$
(24)
Note that $β$ in Eq. (24) is written in the rotating frames $xyz$. Transforming $β$ back to the inertia coordinates $XYZ$ via Eq. (6) results in
$β=∑n=1nr(βxncosω3t-βynsinω3t)qn(R)(t)I+∑n=1nr(βxnsinω3t+βyncosω3t)qn(R)(t)J+∑n=1nrβznqn(R)(t)K$
(25)
The angular displacements of $A$ relative $A'$ are $ξ≡ξxI+ ξyJ+ξzK$, where
$ξx=∑n=1nr[βxncosω3t-βynsinω3t]qn(R)(t)-∑n=1nhαxnqn(H)(t)$
(26)
$ξy=∑n=1nr[βxnsinω3t+βyncosω3t]qn(R)(t)-∑n=1nhαynqn(H)(t)$
(27)
and
$ξz=∑n=1nrβznqn(R)(t)-∑n=1nhαznqn(H)(t)$
(28)

### Vector Notation.

To keep track of the lengthy derivation, let us define the vector of generalized coordinates as
$q(t)=[qHT(t),qRT(t)]T$
(29)
where $qH(t)=[q1(H),q2(H),…,qnh(H)]T$ and $qR(t)=[q1(R),q2(R),…,qnr(R)]T$. With the vector notation, the linear and angular bearing deformations can be rewritten as
$Δb≡(Δx,Δy,Δz,ξx,ξy)T=Bb(t)q(t)$
(30)

where $Bb(t)≡Bb(0)+Bb(1)cosω3t+Bb(2)sinω3t$ and the subscript $b$ refers to the bearings.

### Generalized Bearing Forces.

For each bearing, the bearing forces consist of three force components and two moment components. Moreover, the bearing forces are linear combination of spring forces described by
$Fb≡(Fxb,Fyb,Fzb,Mxb,Myb)T=-KbΔb$
(31)
where $Kb$ is a $5×5$ stiffness matrix. The virtual work done by the bearing forces is
$δWb=∑FbTδΔb=-δqTQb$
(32)

where the summation sums over all the bearings and $Qb$ is the generalized bearing force vector. Substitution of Eqs. (30) and (31) into Eq. (32) yields $Qb=KB(t)q$, where $KB(t)≡∑BbT(t)KbBb(t)$.

### Generalized Damping Forces.

For rotating machines, damping is very difficult to model accurately and effectively. If the damping model is too simple, predictions on forced response amplitude will be inaccurate. If the damping model is too complicated, it could become impractical to apply. As a compromise, the damping model used in this section assumes proportional damping and modal viscous damping factors.

For the stationary housing, let us assume a damping force density (damping force per unit volume) as
$FH(r∧,t)=-2ρH∑n=1nHζn(H)ωn(H)Wn(H)(r∧)q·n(H)(t)$
(33)
where $ρH$ is the density of the stationary housing, and $ζn(H)$ is the viscous damping factor of each housing mode. In this case, the model allows some flexibility through $ζn(H)$ while keeping the complexity of the model under control. The virtual work done by the damping force is
$∫VH(FH·δRP)dVH=-2∑n=1nHζn(H)ωn(H)q·n(H)(t)δqn(H)(t)$
(34)

where Eq. (4) has been used.

For the rotating part, there are two types of damping: internal damping and external damping. For asymmetric rotors, external damping is very difficult to model accurately and remains largely open. Therefore, we will only focus on internal damping only. For internal damping, let us assume
$FR(r˜,t)=-2ρR∑n=1nRζn(R)ωn(R)Wn(R)(r˜)q·n(R)(t)$
(35)
where $ρR$ is the density of the rotor, and $ζn(R)$ is the internal viscous damping factor of each rotor mode. The virtual work done by the internal damping force is
$∫VR(F·RδR(R))dVR=-2∑n=1nRζn(R)ωn(R)q·n(R)(t)δqn(R)(t)$
(36)

where Eq. (10) has been used.

Finally, the virtual work done by all the damping forces is arranged in a matrix form
$δWd=-δqTDq·$
(37)
where
$D=2×diag[ζ1(H)ω1(H),…,ζnH(H)ωnH(H),ζ1(R)ω1(R),…,ζnR(R)ωnR(R)]$
(38)

### Equations of Motion.

With the kinetic energy, potential energy, and generalized bearing forces of the rotor-bearing-housing system, one can derive the equations of motion through use of Lagrangian equations. The resulting equation of motion is
$q··+[G+D]q·+[K+KB(t)]q=0$
(39)
In Eq. (39),
$G=[000G22(n,m)]$
(40)
with $G22(n,m)=2ω3gnm$ given in Eq. (17) and $D$ is defined in Eq. (38). Also,
$K=[ωH00ωR-ω32Kλ(n,m)]$
(41)

where $ωH≡diag[{ω1(H)}2,{ω2(H)}2,…,{ωnh(H)}2]$, $ωR=diag[{ω1(R)}2,{ω2(R)}2,…,{ωnr(R)}2]$, and $Kλ(n,m)=λnm$ defined in Eq. (14). $KB(t)$ in Eq. (39) is symmetric and has periodic coefficients with the spin speed $ω3$.

There are two things worth noting regarding the equation of motion (Eq. (39)). First, not all rotor modes need to be retained in Eq. (39). According to Kim et al. [20], the presence of bearings may cause some rotor modes to couple with the bearings, while other modes may not [20]. If a rotor mode presents unbalanced inertia forces or moments as the rotor vibrates, the rotor mode is called an unbalanced mode. The unbalanced inertia forces or moments will result in deforming the bearings. Therefore, only unbalanced modes will be coupled to the housing via the bearings, and only unbalanced modes need to be retained in Eq. (39).

Second, instability of the rotor-bearing-housing system can be found analytically or numerically. For analytical solutions, one can use any perturbation method, such as the method of multiple scales, to Eq. (39). Since $KB(t)$ is symmetric in Eq. (39), the rotor-bearing-housing system will present instability in the form of combination resonances of the sum type, if the gyroscopic matrix $G$ and the damping matrix $D$ are not present. Moreover, the range of spin speed in which the combination resonance occurs can be predicted analytically in terms of the coefficients in $K$ and $KB(t)$.

The instability can also be found by numerically integrating the equation of motion (Eq. (39)). For example, Eq. (39) can be written in a state-space form and integrated over one period of time to find the fundamental matrix. If any eigenvalue of the fundamental matrix has a magnitude greater than unity, the rotor-bearing-housing system is unstable.

### Convergence Rate.

The essence of this mathematical model is component-mode synthesis, in which the rotor and housing are each treated as an individual component. Therefore, the shape functions used to discretize each component will affect the convergence rate of the model.

In this model, natural boundary conditions are imposed on each component at the bearings to obtain mode shapes as the shape functions. Such an approach provides great convenience because mode shapes form a complete set of admissible functions. Moreover, this approach allows a simpler formulation of the centrifugal softening, gyroscopic effects, and bearing deformation in general. Nevertheless, a trade-off is the convergence rate because the natural boundary conditions may not be satisfied after the components are assembled.

When modeling flexible multibody systems, Meirovitch and Kwak [22] point out that convergence can suffer if the admissible functions (mode shapes in our case) of substructures are not chosen properly. The poor convergence results from the fact that a finite linear combination of the admissible functions is not able to satisfy the natural boundary conditions. As a result, a relatively large number of degrees of freedom are needed for a reasonable convergence.

To prevent poor convergence rate, we use the concept of balanced and unbalanced modes [20] to minimize the number of mode shapes that are capable of satisfying all natural boundary conditions. (Since unbalanced modes include both unbalanced inertia forces and moments, their linear combination comprises the minimum number of mode shapes to satisfy all possible natural boundary conditions.) We also conducted a convergence test to ensure that the equations of motion have a decent convergence rate.

## Example 1: Axisymmetric Rotor

In this example, the rotor is a circular disk with a rigid hub; see Fig. 3. The disk has an outer diameter of 152.8 mm, an inner diameter of 37.28 mm, and a thickness of 0.79 mm. The disk is elastic with a Young's modulus of 280 GPa and Poisson's ratio of 0.33 and has a density of 2700 kg/$m3$. The disk has a rigid hub at the center with an outer diameter of 37.28 mm, an inner diameter of 25 mm, and a thickness of 4.74 mm. Also, the hub is symmetric with respect to the disk plane. Based on the rotor geometry, gyroscopic term will be negligible, i.e., $G=0$. Moreover, damping is assumed to be zero, i.e., $D=0$, in this example to serve as a reference.

Fig. 3
Fig. 3
Close modal

The housing consists of a square plate and a shaft; see Fig. 4. The length and thickness of the square plate are 180 mm and 3.95 mm, respectively. The plate is elastic with a Young's modulus of 280 GPa and a Poisson's ratio of 0.33 and has a density of 2700 kg/$m3$. Moreover, the square plate is simply supported at the four corners. The shaft is rigid and is located at the center of the plate. The shaft has an outer diameter of 25 mm and length of 5.925 mm.

Fig. 4
Fig. 4
Close modal
The rotor and housing are connected via two ball bearings. The outer race of the upper bearing is located at 0.79 mm above the disk plane on the rigid hub. The outer race of the lower bearing is located at 0.79 mm below the disk plane on the rigid hub. The inner races of the bearings are located on the shaft of the housing. Moreover, the bearing stiffness coefficients are given as
$Kb=diag[2.16×107,2.16×107,5.76×106,58.5,58.5]$
(42)

where the coefficients in Eq. (42) are in mks units.

To obtain natural frequencies $ωn(R)$ and mode shapes $Wn(R)(r)$ of the rotor, we first note that the rotor is subjected to free boundary conditions. Therefore, the rotor will present both rigid-body modes and elastic modes. The rigid-body modes are obtained analytically (cf. Eq. (9)), while the elastic modes are obtained via a finite element analysis in ANSYS. The elastic modes basically appear as disk vibration with nodal diameters and nodal circles. (They are not shown here due to space limitation.) A quick calculation via Ref. [20] indicates that only zero- and one-nodal-diameter modes are unbalanced modes.

For this example, the following modes are retained in Eq. (39): (a) three rigid-body translational modes of the rotor along the $x$, $y$, and $z$ axes, (b) two rigid-body rocking modes of the rotor around the $x$ and $y$ axes, (c) unbalanced modes of the rotors including zero- and one-nodal-diameter modes whose natural frequencies are lower than 4000 Hz, and (d) 14 elastic modes from the housing.

To predict the stability of the rotor-bearing-housing system, the equation of motion (Eq. (39)) is rewritten in a state-space form and numerically integrated over a period of time to obtain its fundamental matrix. If any eigenvalue has a magnitude greater than 1, the system is unstable. Figure 5 shows the largest magnitude of all eigenvalues from the fundamental matrix. As shown in Fig. 5, there are three instability zones and each instability zone corresponds to a combination resonance of the sum type. Each combination resonance involves a vibration mode from the housing and a vibration mode from the rotor. The vibration mode from the housing is the same for all instability zones. As far as the rotor modes, the first and the last instability zones are dominated by rigid-body rocking modes and rigid-body translational modes, respectively. The second instability zone is dominated by a one-nodal-diameter unbalanced mode with natural frequency of 1325.90 Hz. For the speed range shown in Fig. 5, the stiffness matrix $K$ remains positive definite. If the spin speed increases further, $K$ will lose its positive definiteness (cf. Eq. (41)) resulting in a divergence instability.

For this example and its instability, there are several issues worth further exploration. First, the equation of motion of an axisymmetric rotor-bearing-housing system can be written in a different form than Eq. (39). Since the spinning rotor is axisymmetric, one can use a ground-based observer to describe motion of the spinning axisymmetric rotor. By doing so, Tseng et al. [21] shows that the equation of motion of an axisymmetric rotor-bearing-housing system can take the form of ordinary differential equations with constant coefficients, i.e.,
$q'··+G'q'·+K'q'=0$
(43)

where $q'$ is a vector of generalized coordinates from the ground-based observer only, and $G'$ and $K'$ are corresponding constant gyroscopic and stiffness matrices observed in the ground-based coordinates. (Note that $G≠G'$ and $K≠K'$.) For systems with an axisymmetric rotor, equations of motion (Eqs. (39) and (43)) are equivalent. One equation of motion can be derived from the other if a coordinate transformation is introduced between the rotor-based coordinates and the ground-based coordinates. Indeed, we have verified it but cannot afford to present it here due to the limited space. Such transformation has also been discussed in Ref. [23].

Second, the stability predicted from Eqs. (39) and (43) should be identical. According to Huseyin [24], the stability of the rotor-bearing-housing system can be predicted through the use of the characteristic equation of Eq. (43). Let $λ$ be eigenvalues of the system, the characteristic equation of Eq. (43) can be written as a function of $λ2$, i.e.,
$f(λ2)≡|-λ2-iλG'+K'|=0$
(44)

Complex $λ$ in Eq. (44) indicates instability of the system. After we transform the equation of motion to the ground-based observer, we have also calculated $λ$ according to Eq. (44). We observe that $λ$ becomes complex at the same spin speed where the rotor-bearing-housing system experiences instability as shown in Fig. 5. Therefore, the stability predicted from Eqs. (39) and (43) are indeed the same.

Third, instability does appear in gyroscopic systems under certain conditions. Meirovitch [25] shows that eigenvalues of undamped gyroscopic systems are purely imaginary implying no instability. That conclusion, however, is achieved under the assumption that the stiffness matrix $K'$ is positive definite. Huseyin [24] indicates that the instability can occur when a gyroscopic system loses its positive definiteness. According to our calculation, all instabilities in Fig. 5 occur after the stiffness matrix $K'$ in Eq. (43) has lost its positive definiteness. This finding is consistent with the stability predictions by Huseyin [24] and Meirovitch [25].

## Example 2: Cyclic Symmetric Rotor

In this example, the rotor is a circular disk with four radial slots evenly spaced in the circumferential direction; see Fig. 6. The slotted disk also has a rigid hub. The disk has the same dimensions and material as in example 1, except that the slots are 4.8 mm wide and 52 mm long. The housing remains the same, and the rotor is mounted on the housing through the same two bearings. Again, there is no gyroscopic term, i.e., $G=0$. Both undamped and damped systems are considered.

Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal

### Undamped System.

In this part of the example, the damping term is assumed to be zero, i.e., $D=0$. For the rotor, its rigid-body modes are calculated analytically. For the rotor's elastic modes, their natural frequencies $ωn(R)$ and mode shapes $Wn(R)(r)$ are obtained via ANSYS. For this example, the following modes are retained in Eq. (39): (a) three rigid-body translational modes of the rotor along the $x$, $y$, and $z$ axes, (b) two rigid-body rocking modes of the rotor around the $x$ and $y$ axes, (c) unbalanced modes of the rotors with natural frequencies lower than 4000 Hz, and (d) 14 elastic modes from the housing. Please note that the unbalanced modes for the cyclic symmetric rotor are significantly different from those of the axisymmetric rotor in example 1. Altogether, 34 modes are retained in Eq. (39) for the cyclic symmetric rotor-bearing-housing system due to the large number of unbalanced modes under 4000 Hz resulting from the slots.

To ensure that enough number of modes is retained, we check the convergence by considering a stationary rotor with $ω3=0$. Theoretically, the solution of Eq. (39) will converge to the exact solution when infinitely many modes are retained. When $ω3=0$, the exact solution can be approximated fairly well by using a finite element analysis of the entire rotor-bearing-housing system. Table 1 compares the first 14 natural frequencies between these two sets of solutions. The frequencies in the column labeled “ANSYS” are natural frequencies predicted from the finite element analysis of the entire rotor-bearing-housing system with $ω3=0$. The frequencies in the column labeled “matlab” are natural frequencies predicted from Eq. (39) with 34 modes retained. According to Table 1, the largest relative error between these two solutions is 3.63%. Therefore, 34 modes are enough to ensure convergence in predicting an accurate response of the cyclic symmetric rotor-bearing-housing system.

Table 1

Natural frequency (Hz): ANSYS versus Matlab

ModeANSYSMatlab% Relative error
1148.55149.410.58
2148.55149.410.58
3294.22296.100.64
4313.30313.220.03
5423.82423.990.04
6423.82423.990.04
7452.99452.890.00
8440.45456.423.63
9472.59472.590.02
10925.25928.780.38
11925.27929.010.40
121104.701109.390.42
131104.701109.390.42
141138.701139.080.03
ModeANSYSMatlab% Relative error
1148.55149.410.58
2148.55149.410.58
3294.22296.100.64
4313.30313.220.03
5423.82423.990.04
6423.82423.990.04
7452.99452.890.00
8440.45456.423.63
9472.59472.590.02
10925.25928.780.38
11925.27929.010.40
121104.701109.390.42
131104.701109.390.42
141138.701139.080.03

When the spin speed $ω3≠0$, the equation of motion (Eq. (39)) can be analyzed via the method of multiple scales to predict the speed range in which combination resonance may occur. Table 2 shows the first four predicted speeds $ω$ and their bandwidth $σ$, where the combination resonance occurs at $ω-σ/2<ω3<ω+σ/2$.

Table 2

Instabilities: the method of multiple scales

InstabilitiesSpeed $ω$ (rpm)Bandwidth $σ$ (rpm)
164,075.781193.66
281,426.85311.12
3122,565.223761.47
4135,695.50629.20
InstabilitiesSpeed $ω$ (rpm)Bandwidth $σ$ (rpm)
164,075.781193.66
281,426.85311.12
3122,565.223761.47
4135,695.50629.20

At the same time, Eq. (39) can be numerically integrated as explained above. Figure 7 shows the largest magnitude of eigenvalues from the fundamental matrix. (The numerical integration was only performed within 64,000–135,000 rpm with a resolution of 100 rpm.) As shown in Fig. 7, there are five instability zones. The five instability zones basically have the same mechanism as those shown in example 1. Each instability zone results from a combination resonance between a housing mode and a rotor mode. The housing mode involved is the same from example 1. The first instability zone involves rigid-body rocking of the rotor. The second, third, and fourth instability zones each involves an unbalanced mode of the rotor. Finally, the last instability zone involves rigid-body translation in the radial ($x$ and $y$) directions.

Fig. 7
Fig. 7
Close modal

There are several things worth discussing. First, the predictions from the method of multiple scales (cf. Table 2) are quite accurate. According to Fig. 7, the first four instability zones center around 64,000 rpm, 81,000 rpm, 122,000 rpm, and 135,000 rpm. These predictions, in general, agree well with those shown in Table 2.

Second, combination resonances resulting from the rigid-body translation and rocking modes are almost identical for the axisymmetric and cyclic symmetric rotors. Combination resonances resulting from the unbalanced modes, however, are very different. For the axisymmetric rotor, only one unbalanced mode with a natural frequency of 1325.90 Hz leads to an instability zone. For the cyclic symmetric rotor, it has three unbalanced modes under 4000 Hz (i.e., 419.24 Hz, 1080.30 Hz, and 1272.40 Hz, respectively). These unbalanced modes lead to the second, third, and fourth instability zones shown in Fig. 7.

### Effects of Damping.

To examine how damping affects the instability, let us incorporate the damping terms $D≠0$ into the equation of motion (Eq. (39)) in two steps.

In the first step, the rotor damping is absent (i.e., $ζ(R)=0$), whereas the housing damping $ζ(H)$ is assumed constant. Figure 8 is the blowup of the fourth instability zone in Fig. 7 in the interval of 136,000–152,000 rpm after damping of the housing modes is included. As seen in Fig. 8, the width of the instability zone increases as the viscous damping factor of the housing increases (following the arrow). In other words, the housing damping has a destabilizing effect.

Fig. 8
Fig. 8
Close modal

In the second step, the viscous damping of the housing $ζ(H)$ remains the same, but the rotor damping is introduced such that $ζ(R)=0.05ζ(H)$. Figure 9 is the same blowup of the fourth instability zone with damping of both the housing and rotor included. As seen in Fig. 9, the width of the instability zone decreases as the rotor damping increases. In other words, the rotor damping has a stabilizing effect.

Fig. 9
Fig. 9
Close modal

The numerical results shown in Figs. 8 and 9 may seem somewhat counterintuitive. There are, however, theoretical supports that housing damping could destabilize the system and rotor damping could do the opposite. Let us investigate this in more detail as follows.

According to Eq. (38), the damping matrix $D$ has nonuniform diagonal values resulting from different natural frequencies and viscous damping factors of the housing and rotor modes. The effect of nonuniform damping on parametrically excited systems, whose equations of motion are similar to Eq. (39), is particularly intriguing [26–28]. According to Refs. [26,27], nonuniform damping could worsen combination resonances by increasing the width of instability zones under certain circumstances. The particular circumstances depend highly on the structure of the periodic matrices as well as the damping ratio of each mode. Therefore, effects of nonuniform damping may vary significantly from one problem to another.

On the other hand, this phenomenon can be viewed from a different perspective. If there exists a coordinate transformation such that Eq. (39) can be transformed into an equation with constant coefficients (as in example 1), it would take the form
$q'··+[G'+D'R+D'H]q'·+[K'+S'R]q'=0$
(45)

In Eq. (45), $DR'$ and $SR'$ are constant damping and circulatory matrices resulting from the damping of the rotor modes, whereas $DH'$ results from the damping of the stationary housing. Equation (45) is a standard form of damped gyroscopic systems with nonconservative forces [29–32].

For such systems, damping can destabilize the otherwise stable gyroscopic systems, and yet combination of damping and circulatory matrices can do the opposite [29,31,32]. In the case of asymmetric rotors, such a coordinate transformation may not even exist, and Eq. (45) is not necessarily achievable. Nonetheless, the criteria provided in Refs. [29,31,32] are still good guidelines to investigate the qualitative behavior of Eq. (39) when damping is present.

Based on the results from Refs. [29,31,32], the presence of housing damping will introduce only $DH'$, which tends to destabilize the coupled rotor-bearing-housing system. This is exactly what the numerical results show in Fig. 8. In contrast, the damping of the rotor modes contributes to both damping and circulatory matrices in Eq. (45), which can stabilize the coupled rotor-bearing-housing system. This is exactly what the numerical results show in Fig. 9. Therefore, the numerical results in Figs. 8 and 9 are in line with theoretical predictions from the literature [29,31,32].

## Conclusions

Through the derivation and simulations shown in this paper, we reach the following conclusions:

1. (1)

For a spinning, cyclic symmetric rotor supported with linear bearings and an elastic housing, its equation of motion takes the form of a set of ordinary differential equation with periodic coefficient as shown in Eq. (39). The period is $2π/ω3$.

2. (2)

When the rotor is axisymmetric, the rotor-bearing-housing system could become unstable due to combination resonance of the sum type. The combination resonance involves a vibration mode of the housing and a vibration mode from the rotor. The vibration mode from the rotor can be a rigid-body translation mode, a rigid-body rocking mode, or an elastic rotor mode that has unbalanced inertia forces or moments (i.e., an unbalanced rotor mode). This combination resonance only occurs when the stiffness matrix, viewed from the ground-base observer, loses its positive definiteness.

3. (3)

When the rotor is cyclic symmetric, the rotor-bearing-housing system could become unstable due to combination resonance of the sum type. The combination resonance has the same characteristics as in the case of an axisymmetric rotor. Since cyclic symmetric rotors have more unbalanced modes, more combination resonances appear in spinning cyclic symmetric rotor-bearing-housing systems.

4. (4)

When damping is present in a coupled rotor-bearing-housing system, qualitative behavior of the instability is similar to that of a damped gyroscopic system with nonconservative forces. Damping of the housing widens the instability zones further destabilizing the system. Rotor damping tends to stabilize the system by reducing the width of the instability zones.

## Acknowledgment

This material is based upon work supported by the National Science Foundation under Grant No. CMMI-0969024.

## References

1.
Murtagh
,
P. J.
,
Ghosh
,
A.
,
Basu
,
B.
, and
Broderick
,
B. M.
,
2008
, “
Passive-Control of Wind Turbine Vibrations Including Blade/Tower Interaction and Rotationally Sampled Turbulence
,”
Wind Energy
,
11
(
4
), pp.
305
317
.10.1002/we.249
2.
Murtagh
,
P. J.
,
Basu
,
B.
, and
Broderick
,
B. M.
,
2005
, “
,”
Eng. Struct.
,
27
, pp.
1209
1219
.10.1016/j.engstruct.2005.03.004
3.
Hansen
,
M. H.
,
2003
, “
Improved Modal Dynamics of Wind Turbines to Avoid Stall-Induced Vibrations
,”
Wind Energy
,
6
, pp.
179
195
.10.1002/we.79
4.
Lee
,
D.
,
Hodges
,
D. H.
, and
Patil
,
M. J.
,
2002
, “
Multi-Flexible-Body Dynamic Analysis of Horizontal Axis Wind Turbines
,”
Wind Energy
,
5
, pp.
281
300
.10.1002/we.66
5.
Malcolm
,
D. J.
,
2002
, “
Modal Responses of 3-Bladed Wind Turbines
,”
ASME J. Sol. Energy Eng.
,
124
, pp.
372
377
.10.1115/1.1509479
6.
Santiago
,
O. D.
, and
Abraham
,
E.
,
2008
, “
Rotordynamic Analysis of a Power Turbine Including Support Flexibility Effects
,”
Proceedings of the ASME Turbo Exposition
,
Berlin
, June 9–13,
ASME
Paper No. GT2008-50900.10.1115/GT2008-50900
7.
Bonello
,
P.
, and
Hai
,
P. M.
,
2009
, “
Computational Studies of the Unbalance Response of a Whole Aero-Engine Model With Squeeze-Film Bearings
,”
Proceedings of the ASME Turbo Expo
,
Orlando, FL
, June 8–12,
ASME
Paper No. GT2009-59687.10.1115/GT2009-59687
8.
Tseng
,
J. G.
, and
Wickert
,
J. A.
,
1994
, “
On the Vibration of Bolted Plate and Flange Assemblies
,”
ASME J. Vibr. Acoust.
,
116
, pp.
468
473
.10.1115/1.2930450
9.
Shen
,
I. Y.
,
1994
, “
Vibration of Rotationally Periodic Structures
,”
J. Sound Vib.
,
172
, pp.
459
470
.10.1006/jsvi.1994.1189
10.
Chang
,
J. Y.
, and
Wickert
,
J. A.
,
2001
, “
Response of Modulated Doublet Modes to Travelling Wave Excitation
,”
J. Sound Vib.
,
242
, pp.
69
83
.10.1006/jsvi.2000.3363
11.
Chang
,
J. Y.
, and
Wickert
,
J. A.
,
2002
, “
Measurement and Analysis of Modulated Doublet Mode Response in Mock Bladed Disks
,”
J. Sound Vib.
,
250
, pp.
379
400
.10.1006/jsvi.2001.3942
12.
Tang
,
J.
, and
Wang
,
K. W.
,
1999
, “
Vibration Control of Rotationally Periodic Structures Using Passive Piezoelectric Shunt Networks and Active Compensation
,”
ASME J. Vibr. Acoust.
,
121
, pp.
379
390
.10.1115/1.2893991
13.
Thomas
,
D. L.
,
1979
, “
Dynamics of Rotationally Periodic Structures
,”
Int. J. Numer. Methods Eng.
,
14
, pp.
81
102
.10.1002/nme.1620140107
14.
Wildheim
,
J.
,
1981
, “
Excitation of Rotating Circumferentially Periodic Structures
,”
J. Sound Vib.
,
75
, pp.
397
416
.10.1016/0022-460X(81)90386-2
15.
Wildheim
,
J.
,
1981
, “
Vibration of Rotating Circumferentially Periodic Structures
,”
Q. J. Mech. Appl. Math.
,
34
, pp.
213
229
.10.1093/qjmam/34.2.213
16.
Jacquet-Richardet
,
G.
,
Ferraris
,
G.
, and
Rieutord
,
P.
,
1996
, “
Frequencies and Modes of Rotating Flexible Bladed Disc-Shaft Assemblies: A Global Cyclic Symmetry Approach
,”
J. Sound Vib.
,
191
, pp.
901
915
.10.1006/jsvi.1996.0162
17.
Nikolic
,
M.
,
Petrov
,
E. P.
, and
Ewins
,
D. J.
,
2007
, “
Coriolis Forces in Forced Response Analysis of Mistuned Bladed Disks
,”
ASME J. Turbomach.
,
129
, pp.
730
739
.10.1115/1.2720866
18.
Kim
,
H.
, and
Shen
,
I. Y.
,
2009
, “
Ground-Based Vibration Response of a Spinning, Cyclic Symmetric Rotor With Gyroscopic and Centrifugal Softening Effects
,”
ASME J. Vibr. Acoust.
,
131
, p.
021007
.10.1115/1.3025847
19.
Shen
,
I. Y.
, and
Kim
,
H.
,
2006
, “
A Linearized Theory on Groung-Based Vibration Response of Rotating Asymmetric Flexible Structures
,”
ASME J. Vibr. Acoust.
,
128
(
2
), pp.
375
384
.10.1115/1.2172265
20.
Kim
,
H.
,
Colonnese
,
N. T. K.
, and
Shen
,
I. Y.
,
2009
, “
Mode Evolution of Cyclic Symmetric Rotors Assembled to Flexible Bearings and Housing
,”
ASME J. Vibr. Acoust.
,
131
, p.
051008
.10.1115/1.3147167
21.
Tseng
,
C. W.
,
Shen
,
J. Y.
,
Kim
,
H.
, and
Shen
,
I. Y.
,
2005
, “
A Unified Approach to Analyze Vibration of Axisymmetric Rotating Structures With Flexible Stationary Parts
,”
ASME J. Vibr. Acoust.
,
127
, pp.
125
138
.10.1115/1.1857917
22.
Meirovitch
,
L.
, and
Kwak
,
M. K.
,
1991
, “
Rayleigh-Ritz Based Substructure Synthesis for Flexible Multibody Systems
,”
AIAA J.
,
29
, pp.
1709
1719
.10.2514/3.10794
23.
Huseyin
,
K.
,
Hagedorn
,
P.
, and
Teschner
,
W.
,
1983
, “
On the Stability of Linear Conservative Gyroscopic Systems
,”
Z Angew. Math. Phys.
,
34
(
6
), pp.
807
815
.10.1007/BF00949057
24.
Huseyin
,
K.
,
1976
, “
Standard Forms of Eigenvalue Problems Associated With Gyroscopic Systems
,”
J. Sound Vib.
,
45
, pp.
29
37
.10.1016/0022-460X(76)90665-9
25.
Meirovitch
,
L.
,
1974
, “
New Method of Solution of Eigenvalue Problem for Gyroscopic Systems
,”
AIAA J.
,
12
, pp.
1337
1342
.10.2514/3.49486
26.
Iwatsubo
,
T.
,
Sugiyama
,
Y.
, and
Ogino
,
S.
,
1974
, “
Simple and Combination Resonances of Columns Under Periodic Axial Loads
,”
J. Sound Vib.
,
33
, pp.
211
221
.10.1016/S0022-460X(74)80107-0
27.
Takahashi
,
K.
,
1982
, “
Instability of Parametric Dynamic Systems With Non-Uniform Damping
,”
J. Sound Vib.
,
85
, pp.
257
262
.10.1016/0022-460X(82)90520-X
28.
Kirillov
,
O. N.
,
2010
, “
Paradoxes of Dissipation-Induced Destabilization or Who Opened Whitneys Umbrella?
,”
ZAMM
,
90
, pp.
462
488
.10.1002/zamm.200900315
29.
Huseyin
,
K.
,
1978
,
Vibrations and Stability of Multiple Parameter Systems
,
Springer
,
New York
, pp.
171
203
.
30.
Meirovitch
,
L.
, and
Ryland
,
G.
,
1985
, “
A Perturbation Technique for Gyroscopic Systems With Small Internal and External Damping
,”
J. Sound Vib.
,
100
, pp.
393
408
.10.1016/0022-460X(85)90295-0
31.
Yang
,
S. M.
, and
Mote
,
C. D.
,
1991
, “
Stability of Non-Conservative Linear Discrete Gyroscopic Systems
,”
J. Sound Vib.
,
147
(
3
), pp.
453
464
.10.1016/0022-460X(91)90493-4
32.
Kirillov
,
O. N.
,
2007
, “
Gyroscopic Stabilization in the Presence of Nonconservative Forces
,”
Dokl. Math.
,
76
, pp.
780
785
.10.1134/S1064562407050353