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.

where Eq. (2) has been used. Moreover, the kinetic energy of the housing is $T(H)=12\u2211n=1nh[q\xb7n(H)]2$, and the potential energy of the housing is $V(H)=12\u222bVh[W(H),W(H)]dV=12\u2211n=1nh[\omega 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\u02dcY\u02dcZ\u02dc$ with $Z\u02dc$ being the spin axis in the virgin state. (Note that the rotor might wobble in actual motion; therefore, $Z\u02dc$ is not the spin axis in actual motion.) Also, one can define a rotating coordinate system $xyz$ with constant angular velocity $\omega 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\u02dcY\u02dcZ\u02dc$) coordinate system. Then

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

where $I\xaf3$ is the centroidal mass moment of inertia about the spin axis.

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),

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

### Vector Notation.

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

### Generalized Bearing Forces.

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

where Eq. (4) has been used.

where Eq. (10) has been used.

### Equations of Motion.

where $\omega H\u2261diag[{\omega 1(H)}2,{\omega 2(H)}2,\u2026,{\omega nh(H)}2]$, $\omega R=diag[{\omega 1(R)}2,{\omega 2(R)}2,\u2026,{\omega nr(R)}2]$, and $K\lambda (n,m)=\lambda nm$ defined in Eq. (14). $KB(t)$ in Eq. (39) is symmetric and has periodic coefficients with the spin speed $\omega 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.

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.

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

To obtain natural frequencies $\omega 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.

*constant*coefficients, i.e.,

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\u2260G'$ and $K\u2260K'$.) 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].

Complex $\lambda $ in Eq. (44) indicates instability of the system. After we transform the equation of motion to the ground-based observer, we have also calculated $\lambda $ according to Eq. (44). We observe that $\lambda $ 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.

### 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 $\omega 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 $\omega 3=0$. Theoretically, the solution of Eq. (39) will converge to the exact solution when infinitely many modes are retained. When $\omega 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 $\omega 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.

Mode | ANSYS | Matlab | % Relative error |
---|---|---|---|

1 | 148.55 | 149.41 | 0.58 |

2 | 148.55 | 149.41 | 0.58 |

3 | 294.22 | 296.10 | 0.64 |

4 | 313.30 | 313.22 | 0.03 |

5 | 423.82 | 423.99 | 0.04 |

6 | 423.82 | 423.99 | 0.04 |

7 | 452.99 | 452.89 | 0.00 |

8 | 440.45 | 456.42 | 3.63 |

9 | 472.59 | 472.59 | 0.02 |

10 | 925.25 | 928.78 | 0.38 |

11 | 925.27 | 929.01 | 0.40 |

12 | 1104.70 | 1109.39 | 0.42 |

13 | 1104.70 | 1109.39 | 0.42 |

14 | 1138.70 | 1139.08 | 0.03 |

Mode | ANSYS | Matlab | % Relative error |
---|---|---|---|

1 | 148.55 | 149.41 | 0.58 |

2 | 148.55 | 149.41 | 0.58 |

3 | 294.22 | 296.10 | 0.64 |

4 | 313.30 | 313.22 | 0.03 |

5 | 423.82 | 423.99 | 0.04 |

6 | 423.82 | 423.99 | 0.04 |

7 | 452.99 | 452.89 | 0.00 |

8 | 440.45 | 456.42 | 3.63 |

9 | 472.59 | 472.59 | 0.02 |

10 | 925.25 | 928.78 | 0.38 |

11 | 925.27 | 929.01 | 0.40 |

12 | 1104.70 | 1109.39 | 0.42 |

13 | 1104.70 | 1109.39 | 0.42 |

14 | 1138.70 | 1139.08 | 0.03 |

When the spin speed $\omega 3\u22600$, 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 $\omega $ and their bandwidth $\sigma $, where the combination resonance occurs at $\omega -\sigma /2<\omega 3<\omega +\sigma /2$.

Instabilities | Speed $\omega $ (rpm) | Bandwidth $\sigma $ (rpm) |
---|---|---|

1 | 64,075.78 | 1193.66 |

2 | 81,426.85 | 311.12 |

3 | 122,565.22 | 3761.47 |

4 | 135,695.50 | 629.20 |

Instabilities | Speed $\omega $ (rpm) | Bandwidth $\sigma $ (rpm) |
---|---|---|

1 | 64,075.78 | 1193.66 |

2 | 81,426.85 | 311.12 |

3 | 122,565.22 | 3761.47 |

4 | 135,695.50 | 629.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.

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\u22600$ into the equation of motion (Eq. (39)) in two steps.

In the first step, the rotor damping is absent (i.e., $\zeta (R)=0$), whereas the housing damping $\zeta (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.

In the second step, the viscous damping of the housing $\zeta (H)$ remains the same, but the rotor damping is introduced such that $\zeta (R)=0.05\zeta (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.

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.

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)
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\pi /\omega 3$.

- (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)
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)
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.