## Abstract

In machine systems where a rotor spins within a finite clearance space supported by bearings, contact between the rotor and its surround can result in persistent coupled vibration of the rotor and stator. When the vibration response is driven predominantly by friction forces, rotordynamic stability becomes a serious issue. This paper introduces a theory for model-based verification of dynamic stability in rotor systems with stator contact and rub. Generalized multi-degree-of-freedom linear models of rotor and stator lateral vibration are considered, combined with contact models that account for finite clearance and Coulomb friction. State-space conditions for global stability as well as stability of contact-free synchronous whirl responses are derived using Lyapunov's direct method. This leads to feasibility problems involving matrix inequalities that can be quickly verified using numerical routines for convex optimization. Parametric studies involving flexible rotor models indicate that tight bounds on regions of stability can be obtained. A case study involving a realistic machine model illustrates how design optimization based on the theory might be used to overcome instability problems in real machines.

## 1. Introduction

In turbomachinery and other rotordynamic systems where a rotor spins within a finite clearance space supported by bearings, contact between the rotor and its surround is ideally avoided. However, during machine operation, structural vibration can arise from various sources, including imbalance, misalignments, thermal bending, external disturbances and destabilizing fluid dynamic effects. In more extreme cases, this may lead to unwanted interaction between the rotor and the stator, typically involving sliding contact between the rotor and a machine part such as a housing, seal, bushing or, in magnetic bearing systems, a retainer bearing. The resulting vibration may be partly friction driven, with an associated potential for destructive instabilities. Accurate prediction of this behaviour using dynamic models can be hindered by uncertainties in the contact mechanics, as well as the complexity and lack of determinism that often characterize the behaviour of nonlinear dynamical systems.

In radially isotropic systems, solutions to the equations of motions can be obtained by assuming simple forms of steady-state vibration. However, the behaviour that occurs in actuality is sensitive to system disturbances and initial conditions. Numerical integration of the equations of motion can allow a fuller exploration of response and stability, but simulation of many test cases is both time consuming and does not allow for easy identification of contributing factors. Model-based methods that can be used to assess stability and thereby validate or refine a machine design more directly are therefore needed.

In this paper, a theory is developed for model-based verification of dynamic stability in rotor systems undergoing stator contact and rub. In the following sections, rotor–stator interaction models and response behaviour are first reviewed and some associated frequency-response-based results presented. Section 3 contains the main theoretical contribution of the present work, describing the derivation of a number of stability conditions based on generalized state-space models of the system dynamics. Section 4 looks at a selection of test cases for which parametric boundaries for stability are evaluated and compared. A design example is then given in §5 to show how parameter optimization based on the state-space stability conditions could be used to prevent unwantedvibration behaviour. Section 6 provides conclusions and an outlook on further work.

## 2. Rotor–stator interaction models and response phenomena

Basic predictions concerning rotor–stator interaction across a clearance annulus have been widely based on the polar receptance method, as introduced by Johnson (1962) and later extended by Black (1968) to account for friction effects. Black showed that, when rotor mass imbalance is comparable to the radial clearance, degenerate contacting and non-contacting whirl responses are both possible at certain rotational speeds. This leads to amplitude-jump phenomena and the possibility that vibration behaviour may differ depending on previous operating speeds and system disturbances. The potential for friction-driven backward whirl can also be predicted through the same approach of seeking periodic solutions to the equations of motion. More recently, experimental studies have addressed the issues of whether whirl responses predicted in this way are accurate or even relevant to the behaviour of real systems, including Lawen & Flowers (1999), Bartha (2000), Yu *et al*. (2000) and Williams (2004), among others. It is clear that multi-mode rotor and stator models can give rise to a multiplicity of periodic whirl solutions and that these solutions can be sensitive to model parameters. But provided local stability characteristics of the solutions are also determined, useful indications of possible system behaviour can be gleaned (Childs & Bhattacharya 2007; Cole 2008). However, owing to the potential for degenerate aperiodic vibration responses in some systems, it can still be argued that only by establishing global stability of the contact-free vibration response can firm guarantees on system behaviour be made.

### (a) Modelling approach

The approach to modelling will be introduced through reference to an illustrative case of a lumped-mass flexible rotor and stator system (figure 1). Throughout the paper, it will be assumed that both the rotor and the surround are circular in cross section and that the rotational frequency *ω* is constant. In this section, it will also be assumed that the system is radially isotropic, at least in terms of dynamic properties, so that complex representation can be usefully adopted. With , the equation describing transverse vibration of the rotor in the contact plane is(2.1)Typically, the synchronously rotating disturbance force *D* would represent rotor imbalance, in which case we may set *D*=*meω*^{2}, where *e* is the rotor mass-eccentricity. The equation of motion for the stator is(2.2)In this example, (2.1) and (2.2) are coupled only through the contact interaction force ** p**. In general, the rotor and stator dynamic equations may be further coupled through other machine components, such as bearings, seals and electromagnets. The only condition required in this analysis is that, with the contact force

**considered as an input, the stator and rotor dynamics form a linear model that, in the case of radially isotropic systems, can be conveniently expressed using the polar receptance.**

*p*Consider now the geometry associated with contact as shown in figure 2. Displacements are defined relative to the point *O* that is the equilibrium position for the rotor centre *C*. Assume for the moment that, with the rotor and stator centres positioned at *O*, there is a uniform radial clearance *c*. The contact force vector will be oriented to the contact normal at the friction angle . The normal component of the contact force ** f** can be related to the contact deflection

**through a contact stiffness coefficient**

*r**κ*according to(2.3)The free parameter

*θ*was introduced by Cole (2008) in order to ascertain how the two common assumptions either that

*θ*=

*ϕ*, consistent with gross translation of the surround in the direction of the contact force, or that

*θ*=0, which is representative of a local deformation in the direction of the contact normal, can lead to different conclusions about the existence or non-existence of limit-cycle solutions to the equations of motion. The deformation associated with each type of interaction is illustrated in figure 2. By transforming Black's conditions for the existence of circular rub solutions to conditions involving the Nyquist plot of the polar receptance, it was shown that negating existence was equivalent to the Nyquist plot avoiding certain critical regions of the complex plane determined by

*ϕ*,

*κ*and

*θ*. This method will be considered as a useful precursor to the more generalized stability theory developed in the present work and will first be outlined.

During circular whirl, the geometry associated with contact is unchanging in a frame that rotates in an anticlockwise sense about the point *O* at the whirl speed. Distinguishing stationary frame vectors ** p** and their rotating-frame counterparts

**=**

*P***e**

*p*^{−iωt}and then denoting the magnitude of both vectors from (2.3),(2.4)where the scalar

*F*is the normal component of the contact force. Closure of the vector triangle

*ABC*implies(2.5)and so, with , (2.1) and (2.2) give(2.6)Here,

**(**

*Q**ω*) is the response of the rotor without contact and

*g*(

*ω*) is the frequency-dependent polar receptance for the combined rotor–stator system. For the system shown in figure 1,(2.7)Combining (2.4), (2.5) and (2.6) gives(2.8)

### (b) Forward synchronous whirl with rub

Equation (2.8) defines a vector triangle, which from the cosine rule leads to two possible solutions for the contact force (Black 1968),(2.9)When *Q*>*c*, stator contact is inevitable. Only one solution for *P* will then be positive and therefore physically plausible. The case *Q*<*c* is perhaps more relevant to real situations, as typically a rotor will be sufficiently well balanced so that excursions beyond the clearance space are avoided during normal operation. Two real positive solutions for *P* may then exist, but only if and the non-contact whirl amplitude *Q* exceeds a critical value given by . Consequently, for any given system, it is possible to determine rotational frequencies at which there can be degenerate whirl responses with and without contact. If the phase of *g* is in the range (−*π*,0), as would be expected for a passive system, then requires(2.10)Accordingly, the frequency zones for sustained stator interaction may be established from a Nyquist plot of *g*(*ω*) as illustrated in figure 3. The frequency limits for the interaction zone follow from the intersection of the curve *g*(*ω*) with a line originating from the point −*λ*/*κ* with orientation −*π*/2−*ϕ*. Significantly, this line intersects the real axis at the point −1/*κ* and thus the points of intersection with *g*(*ω*) do not depend on *θ*. For the example shown, interaction can be expected for *ω*∈(*ω*_{1},*ω*_{2}). For rotational frequencies within this range, a contact-free orbit is considered unstable in the global sense as transgression to sustained interaction is possible, albeit for sufficiently small clearance. Using this graphical method, it can be clearly seen that, for a given rotational frequency, there is a corresponding value of *κ* which must be exceeded for sustained interaction to be possible. The minimum value of *κ* over all rotational speeds is denoted by *κ*_{1}, also indicated in figure 3.

### (c) Dry friction backward whirl

Backward whirl, or whip, occurs when friction forces are sufficient to drive the whirl in a reversed direction following stator contact. A backward whirl solution to the equations of motion can be found by assuming that external disturbances are negligible *Q*=0 and that a circular orbit involving continuous contact with slip occurs. The rotor and stator equations of motion then give(2.11)which, from (2.4) and (2.5), leads to(2.12)The condition for existence of a real positive solution for *P* follows as(2.13)where * indicates complex conjugation. In systems where the linear model has no *x*–*y* cross-coupling, . Existence of a backward whirl solution can then be determined from the Nyquist plot of *g*(*ω*) (or more generally *g*(−*ω*)^{*}), requiring that the curve intersects the line originating from the point −*λ*^{*}/*κ* with orientation −*π*+*ϕ*. The value *ω*=*Ω* corresponding to the intersection point gives the onset frequency for backward whirl. For given values of *κ* and *ϕ* and for *θ*∈(0,*ϕ*), the region of the complex plane that the curve *g*(*ω*) must avoid for stable operation can be constructed as shown in figure 4. For a given *g*(*ω*), it is possible to establish a maximum value of *κ*, denoted by *κ*_{2}, for which an unstable whirl solution does not exist. This value is essentially the upper limit on the contact stiffness for avoidance of backward whirl and it provides a useful indication of the robustness of the system design. In the state-space theory that follows, it is also a critical parameter in the determination of stability boundaries for any given system.

### (d) Whirl with asynchronous periodic bouncing

Proving that a friction-driven whirl response is theoretically plausible does not give any clear indication of its practical likelihood as it must be initiated by transient conditions. The transition from forward to backward whirl usually involves a series of bouncing motions. In some cases, these bouncing motions can persist indefinitely without the rotor progressing to continuous rub (Tamura *et al*. 2002; Keogh & Cole 2003). A method to predict steady bouncing whirl motions in radially isotropic rotor–stator systems has been proposed by Cole & Keogh (2003). By assuming that stator contacts involve regular short-duration impulses, a periodic Green's function can be established for the rotating-frame dynamics. The fractional energy dissipation resulting from each impact then determines the frequency and magnitude of sustained bouncing. Notably, the amplitude of bouncing increases with imbalance and so if imbalance is too large then the bouncing motion cannot be contained within the clearance space and transgression to full backward whirl is inevitable. This behaviour suggests that although imbalance may not have a significant effect when backward whirl is well established, it can play an important role in its instigation.

### (e) Aperiodic and subharmonic responses

Forced rotor vibration responses with subharmonic components generally involve repeated intermittent contact at certain locations on the surround. Such behaviour has been investigated by experimentation and model-based simulation for various cases of rotors with looseness, stator clearance, both with and without offsets, and bearing clearances (Childs 1979; Muszynska 1984; Ehrich 1988; Muszynska & Goldman 1993). With light damping, a rich pattern of dynamic behaviour is possible, with regimes of periodic and chaotic motions. The thrust of the current work is to derive mathematical conditions under which such response behaviours are eliminated, and so although some researchers have used numerical methods to directly determine periodic response solutions (Kim & Noah 1990; Flowers & Wu 1993; Von Groll & Ewins 2001), here we acknowledge that they can exist but seek only to determine when they do not.

### (f) Degenerate response examples

The previous survey has indicated how a number of different types of vibration response can be predicted using model-based methods. However, the determination of such solutions does not provide concrete indication of system behaviour, if degenerate responses can potentially occur under the same conditions. This can be the case for forward whirl, both with and without contact, backward whirl and aperiodic responses. To illustrate this point, simulation results for the system in figure 1 will now be presented for the parameter set 1 given in table 1. The contact stiffness is *κ*=20 N m^{−1} and the friction coefficient is *μ*=0.1. The response solutions shown in figure 5 are all produced under the same constant synchronous disturbance having angular frequency *ω*=2 rad s^{−1} (which is twice the rotor natural frequency). Only the initial conditions for the simulation runs differ.

## 3. State-space stability criteria: a Lyapunov approach

In the formulation of a generalized theory for stability, we seek to establish whether, for any given system, the following statements are true.

For the unforced system, the non-contact equilibrium state is a globally stable equilibrium state.

For a given forced synchronous whirl orbit satisfying the equations of motion, the response is both unique and stable.

Statement A is significant as it implies that the system cannot undergo an unstable (or limit-cycle) friction-driven response of the types described in §2 and therefore, following contact interaction, the unforced rotor system will always return to the non-contact equilibrium state. Statement B implies that, when a synchronous whirl orbit does not involve stator contact, a degenerate stable orbit involving contact interaction of any form is implausible. Therefore, following any transient disturbances leading to contact, the rotor will always return to the non-contacting orbit. Statement B also implies that sustained aperiodic responses and backward whirl are infeasible. Statement B is stronger than statement A as it implies not only stability but also the absence of any degeneracy in the forced response. Methods to determine the validity of each of these statements for a given system model will be presented in this section.

### (a) State-space system model

Multi-mode rotor and stator dynamic models are used to describe lateral vibration of complex machine structures. Following Nelson & McVaugh (1976), such models are usually based on consistent mass, stiffness and damping matrices derived by finite element (FE) discretization, which leads to a set of coupled second-order ODEs expressible in the following form.

Rotor dynamics:(3.1)

Stator dynamics:(3.2)The vectors *x*_{r}, *y*_{r} and *x*_{s}, *y*_{s} constitute nodal displacements of the rotor and stator elements, respectively. Disturbance force components acting on the rotor are embodied in *d*_{x} and *d*_{y}, while the contact force normal components (numbering *J*) are represented by (*f*_{x}, *f*_{y}). A rotor will also have a number of axial planes where interaction with the stator occurs in a normal fashion through bearings and seals, here indexed by the subscript *i*, and having associated interaction forces *g*_{x} and *g*_{y}. Linearized models of these interactions take the standard form(3.3)where the values of the stiffness and damping coefficients must be appropriate to the rotational speed under consideration. The displacement of the rotor relative to the stator (*z*_{x}, *z*_{y}) can be expressed as(3.4)Defining the state vectors , , a first-order state-space model follows in the form(3.5)The matrices within this equation can be calculated from those in (3.1)–(3.4) according to the formulae given in appendix A. The dynamics of the complete system can now be stated more concisely as(3.6)The nonlinear functions *β*_{j} relate the normal contact force components to the relative displacements in each of the *J* contact planes. Note that a dynamic model of this form will correspond to a certain fixed rotational speed only owing to speed-dependent terms in the FE model.

### (b) Stability analysis

Stability of system (3.6), in a global asymptotic sense, may be established through application of the Barbashin–Krasovskii theorem.

*Let* ** w**=0

*be an equilibrium point for*(

*3.6*).

*Let*

*be a continuously differentiable function such that*(3.7)(3.8)(3.9)

*then*

**=0**

*w**is a globally asymptotically stable equilibrium*.

A Lyapunov candidate *V*(** w**) that combines a quadratic term involving the system states with an energy-based storage function for the nonlinear rotor–stator interaction will be considered. To this end, consider a small change in rotor position relative to the stator in the plane of contact which produces a change in contact deflection as shown in figure 6

*a*. The change in stored elastic energy may be associated with the action of the contact force

**p**moving along a path in the tangential direction during which there is no change in the stored energy, followed by a displacement in the normal direction during which the surround is deflected. The work done on the surround is , but, as can be seen from figure 6

*a*, the normal components of d

**and d**

*r***are equal and so(3.10)A Lyapunov function candidate involving free parameters**

*z***and**

*P**ν*

_{j}follows as(3.11)If the matrix is positive definite and

*ν*

_{j}are positive scalars, then conditions (3.7) and (3.8) of theorem 3.1 are immediately satisfied. Noting from (3.10) that ,(3.12)If through particular choices of

**and**

*P**ν*

_{k}condition (3.9) of theorem 3.1 can also be satisfied, then

**=0 is a globally stable equilibrium. In order to establish whether this is the case, the geometry associated with the functions**

*w**β*

_{j}must be considered. It will be assumed that the deflection of the surround is aligned with the contact normal (corresponding to

*θ*=0 in (2.3)) as this is the more conservative assumption when we wish to establish stability. Two different situations will be dealt with in the following. First, the case where the rotor equilibrium position aligns with the centre of the clearance circle; and second, where there is an offset, as appropriate to situations of rotor misalignment or static deflection. The analysis of the second case will also be extended to establish global stability of circular whirl orbits.

#### (i) Case 1. Stability condition without rotor–stator offset

From (3.6), . However, from the form of (3.1)–(3.5), it follows that *C*_{j}*B*_{j}=0 and so (3.12) gives(3.13)From (2.3), with *γ*∈[0, 1) varying according to the amount of deflection. It then follows that(3.14)For to hold globally, the matrix within this quadratic form must be negative definite for all possible values of *γ*_{j}. The parameter space for *γ*_{j} is a *J*-dimensional unit box or polytope. As the matrix in (3.14) is linear in *γ*_{j}, this condition may be stated in terms of a set of 2^{J} matrix inequalities corresponding to each vertex of the box, which must all be negative definite. This provides a necessary and sufficient condition for stability.

**Stability condition 1**. *If there exist* *and* *ν*_{j}≥0 *such that*(3.15)*then the system defined by* (*3.6*) *is globally asymptotically stable at the equilibrium point* ** w**=0.

#### (ii) Case 2. Stability condition with rotor–stator offset

When the rotor equilibrium position is offset from the centre of the clearance circle but does not involve contact with the surround, the dynamics of the unforced system can be stated as(3.16)where *w*_{0} is the static equilibrium state. Defining and , the Lyaponov function candidate can be taken as(3.17)With , the stability condition follows as(3.18)With reference to figure 6*b*, depending on the rotor position, the deflection of the surround ** r** will be oriented at an angle

*α*to the rotor displacement from the equilibrium position

*A*. However, provided point

*A*is within the clearance circle, and the following constraint, which effects a limit on the degree of orientation, will always hold as follows:(3.19)From

**=−**

*f**κ*

**and then, reintroducing the contact plane index**

*r**j*,(3.20)It is now required to show that (3.18) is satisfied whenever constraint (3.20) holds. To establish whether this is the case, inequalities (3.18) and (3.20) can be combined into a single inequality via the ‘-procedure’ (Boyd

*et al*. 1994): the left-hand side of each inequality in (3.20) is multiplied by a free parameter

*σ*

_{j}≥0 and subtracted from to give the following sufficient condition:(3.21)Defining , and diagonal matrices and , (3.21) can be written more concisely as(3.22)That this must hold for all and

**leads to the following stability condition.**

*f***Stability condition 2**. *If there exist* , ** Y**≥0

*and Σ*≥0

*such that*(3.23)

*then the system defined by*(

*3.16*)

*is globally asymptotically stable at the equilibrium point*

**=**

*w*

*w*_{0}.

Note that stability condition 2 does not depend on the actual equilibrium state or clearance but owing to the use of (3.20) does require that this state is contact free.

### (c) Global stability of circular whirl motions

A modified version of stability condition 2 can be used to establish global asymptotic stability of synchronous whirl orbits in radially isotropic systems. Note that the system states and their rotating-frame counterparts are related by(3.24)A circular whirl, arising from imbalance or some other synchronous excitation, can be written as , which provides a periodic solution to (3.5). The perturbation from this critical point in the rotating frame may be expressed as(3.25)Substitution of (3.24) and (3.25) in the state equation (3.5) then gives(3.26)Here, it has been assumed that *x-* and *y*-axes are dynamically equivalent so that *A*_{xx}=*A*_{yy} and *A*_{yx}=*A*_{xy}. Consequently, the rotating-frame dynamics (3.26) are time invariant. The relationship between the contact force vector and the system states is unchanged by the transformation to a rotating frame and so these vectors also satisfy the bound (3.20),(3.27)The resulting stability condition takes the same form as stability condition 2 but with the matrix ** A** replaced by the whirl frequency-dependent matrix,When the stability condition is satisfied for a particular value of

*ω*, not only is the whirl motion stable but also there can be no degenerate response solutions involving contact. From an engineering perspective, it is clearly advantageous if a system can be designed such that this condition is satisfied for the entire range of operating speeds as this would imply that if rotor–stator contact does not occur during normal operation then, following any transient conditions leading to contact, the rotor will always return to a contact-free state. However, as stated at the beginning of the section, this is a stronger requirement than ensuring against backward whirl instability and thus may be harder to fulfil in a practical design.

## 4. Parametric boundaries for rotordynamic stability

Conditions for stability in each of the cases presented in §3 take the form of feasibility problems involving one or more linear matrix inequalities (LMIs). If a solution can be obtained such that all the eigenvalues of the parameter-dependent matrices are negative, then the stability condition is satisfied. As the dependence of these matrices on the free parameters ** P**,

*ν*

_{j}and

*σ*

_{j}is affine solutions can be found numerically using standard routines for convex optimization, treating the maximum eigenvalue as the objective function to be minimized (Boyd

*et al*. 1994).

### (a) Avoiding backward whirl

Although now equipped with a method to evaluate the stability of complex rotor–stator systems, it is useful to reconsider the basic model of figure 1 for which the results obtained from Black's original theory and the state-space stability conditions can now be compared. For a fixed stator system with equations of motion (2.1) and *z*_{s}=0, solutions to the backward whirl equation (2.12) can be obtained algebraically. With *θ*=0 in the contact model, the frequency of the backward whirl is independent of *κ* and given by (Zhang 1988; Muszynska 2002)(4.1)where and . As noted previously, there will be a maximum value of *κ* for which this unforced response solution exists, in this case given by . The condition for non-existence of such a solution also provides a necessary condition for global stability of the fixed stator system,(4.2)As might be expected, the value of *κ*_{2} for potential instability decreases as friction levels increase or rotor damping is reduced. Notably, the frequency of the whirl solution also drops and so it has been argued that backward whirl is then more easily instigated. It is of interest to compare this result with the stability boundary obtained through stability condition 1 for the same fixed stator system. The problem of maximizing *κ* subject to stability condition 1 is a generalized eigenvalue problem which can be directly solved using standard LMI solvers. For this example, the stability boundary obtained numerically is identical to (4.2), confirming that the state-space condition gives a tight bound on the stability limits. Moreover, if the circular whirl solution does not exist, then no other unstable response is possible. This optimization routine may have direct application to machine design, if the maximum value of *κ* obtained can be used to ensure stator structures are sufficiently compliant at potential contact interfaces.

In a more realistic model, stator motion may also be accounted for with a lumped parameter model as defined by (2.1) and (2.2). Two different systems will be considered, one having a ‘heavy’ and the other a ‘light’ stator with complete parameter sets given in table 1. Note that, within the overall system model, the contributions of the stator and rotor dynamics are equivalent and can be interchanged without affecting the stability results. Stability limits, in terms of the contact stiffness, have been calculated for a range of friction coefficient values and lead to the boundaries shown in figure 7.

For the case of the heavy stator, the state-space conditions exactly reproduce the boundary for existence of circular whirl solutions obtained using the Nyquist plot criterion on *g*(*ω*). For moderate levels of stator damping, this stability boundary is similar to the fixed stator case given by (4.2), which is also shown in figure 7*a*. Noting that the contact model parameters used to obtain the simulation results of figure 5 correspond to point *A*, by increasing the stator damping to the level given in parameter set 2 of table 1, the stability region can be expanded to encompass this point. As a result, both the unstable backward whirl and bouncing whirl response types can be eliminated. For this modified system, repeating the simulation conditions of case 3 in figure 5, the rotor instead settles into a steady rub condition as shown in figure 8.

For the light stator model, specified by parameter set 3 in table 1, the stability region is substantially diminished (figure 7). According to the frequency response analysis, this can be explained by the fact that, at high frequencies, the magnitude of *g*(*ω*) scales with 1/*m*_{s} and thus *g*(*ω*) extends further from the origin of the complex plane. The stability boundary obtained from the state-space condition closely follows the boundary for existence of circular whirl solutions obtained by the frequency response criterion (2.13) and deviates only in the vicinity where *g*(*ω*) forms a cusp on the complex plane.

### (b) Avoiding persistent rotor–stator interaction

For flexible rotors operating at supercritical speeds, coupled rotor–stator interaction is an important issue even if backward whirl is infeasible, as once contact has been instigated then rundown through critical speeds may be necessary to recover contact-free conditions and this could lead to machine damage. Considering again the single-mode rotor and stator system of figure 1, speeds at which contact-free whirl orbits are asymptotically stable can be determined according to the rotating-frame version of stability condition 2. Results for the parameter sets given in table 1 are shown in figure 9. Superimposed on the plots are the boundaries for the existence of degenerate circular rub solutions obtained from the frequency response criterion (2.10). Considering a finite speed range, the state-space conditions give a reasonably tight bound on the minimum value of *κ* for the existence of a degenerate rub response. These results indicate that the state-space condition may be usefully employed to check stability over a finite speed range, simply by checking stability at the maximum rotational speed. It is also found that, for low frequencies, where the forward rub solution does not exist, the stability boundary extends to the boundary for the existence of the backward whirl solution.

As an aside, a related behaviour revealed through simulation is that, when the system is also susceptible to backward whirl, a predicted forward circular whirl solution involving continuous rub may be unstable such that an initial forward whirl will progress to a bouncing whirl as seen in simulation case 2 of figure 5. If imbalance is large, then the bouncing will eventually result in the onset of full backward whirl; that is, a forward whirl with rub will spontaneously progress to full backward whirl.

## 5. Design example: vibration absorber optimization

There are a number of aspects to rotating machine design where the developed theory could be used to assist in the creation of systems that are less prone to contact-induced instabilities. As a final example, this section introduces a realistic high-speed flexible rotor model created using FE methods that account for shear, bending and gyroscopic effects. The rotor is mounted through bearings on a fixed foundation and contact can occur with a compliant housing that is also fixed to the foundation as shown schematically in figure 10. The rotor is symmetric and it is assumed that contact occurs mid-span, as would be the likely case when excitation of the first rotor flexural mode is predominant.

Although the design of a machine rotor is to a large extent predetermined by the requirements for machine function and efficiency, it may be quite feasible to purposefully tune the dynamic characteristics of the machine housing or foundation to avoid any propensity for contact-induced instability. It has already been seen that a rigid fixed housing gives poor protection against backward whirl. In this example, single-degree-of-freedom spring-mass-damper systems are attached to a compliant housing in orthogonal planes, thereby deliberately introducing low-frequency stator modes, the natural frequency and damping of which can be tuned to give improved response characteristics of the coupled rotor–stator system. Although the case considered here could be viewed as representative of a multi-stage compressor operating above the first critical speed, the results presented should be regarded as a proof of concept rather than a true design exercise.

The rotor is 1 m in length and has a mass of 63 kg. For lateral vibration of the non-rotating rotor, the natural frequency corresponding to the first flexural mode is 1060 rad s^{−1}, while the second and third modes have natural frequencies of 2730 and 4460 rad s^{−1}, respectively. The stiffness of the rotor mid-span is 42 MN m^{−1}, while the housing stiffness is 7.2 MN m^{−1} in the same plane. The housing dynamics are modelled by a lumped mass of 20 kg at the mid-span location, resulting in a natural frequency of 600 rad s^{−1}. This mass is additionally coupled to reaction masses of 10 kg through springs and dampers connected in orthogonal directions. The full-order state-space model has 128 states. The nominal rotational speed is 1500 rad s^{−1}. A friction coefficient of 0.2 is assumed for the purpose of optimization.

The design objective is to select optimal stiffness and damping characteristics of the vibration absorbers, *k*_{a} and *b*_{a}, respectively, to minimize the possibility of backwards whirl instability. According to the preceding theory, we may interpret this as maximizing the value of the contact stiffness *κ* over the design parameters subject to stability condition 1. For a single contact plane, stability condition 1 is homogeneous in the free variables and so *ν* can be set to 1. The dependency of the state-space model on the design variables may be expressed as , where . The resulting optimization problem is then as follows.

*Design optimization problem*. Maximize *κ* over , subject to

To tackle an optimization problem of this form, which is bilinear in the free variables ** P** and and therefore non-convex, algorithms developed for similar optimal control problems may be employed. The branch-and-bound technique, as described in the review paper by Van Antwerp & Braatz (2000), is one such approach in which LMI relaxations are used iteratively to obtain lower bounds on the objective function and thereby close in on the global optimum. To expedite the calculations in the current problem, a reduced-order rotor model retaining only the first rotor bending modes (both forward and backward) was used to obtain a first-cut solution: for the unmodified system, calculation of the polar receptance indicated that the first rotor mode dominates in the critical frequency range, where crosses the −

*π*+

*ϕ*phase boundary. This solution was then refined using a higher order model retaining the first three rotor modes. Figure 11 shows the critical part of the polar receptance plot for the un-optimized and optimized systems. For the optimized system, the phase of

*g*(

*ω*) is shifted in a positive sense over the critical frequency range. This means that the −

*π*+

*ϕ*phase line is crossed at a much higher frequency and with a higher corresponding value of contact stiffness . The optimized system is not only stable for larger values of

*κ*, but also the critical frequency

*Ω*for the onset of backward whirl is increased. The results are summarized in table 2.

## 6. Conclusions

This paper has introduced a theory for evaluating stability of rotordynamic systems undergoing rotor–stator contact interaction. It has been shown that conditions for the global stability of the unforced system, as well as certain synchronous whirl solutions, lead to feasibility problems involving LMIs that can be solved numerically using standard routines for convex optimization. The principal advantages of the new theory that may be further exploited in order to overcome instability problems in rotor–stator systems are as follows.

Numerical routines used to establish stability can involve system models that are dependent on a number of design parameters. These additional free variables may then be optimized within the routine to achieve improved robustness of the machine design.

The theory can be applied to systems with active vibration control elements and thus used as a basis for control optimization.

Guaranteed global stability can be established in situations where the system dynamics are anisotropic and when there are multiple potential contact planes.

Possibility for further development of the theory to solve machine design problems of these types is currently being explored. The research also opens up the possibility that similar techniques may be used to assess the stability of rotor systems subject to other forms of nonlinear and non-conservative force characteristics. These may arise, for example, from machine elements, such as fluid film bearings and seals, or from higher order rotordynamic effects, such as coupling of torsional and lateral shaft vibration. The issue of model uncertainty and the associated implications for any model-based indications of stability are also worthy of further consideration. It is foreseeable that modified versions of the stability conditions described here may be derived which can account directly for uncertainty in model parameters. Otherwise, a thorough sensitivity analysis for uncertain parameters would be recommended in an evaluation of any real machine design.

## Acknowledgments

This work was partly supported by a Visiting Research Fellowship at the University of Bath, UK, and with finance under EPSRC Platform Grant GR/S64448/01.

## Footnotes

- Received June 7, 2008.
- Accepted July 31, 2008.

- © 2008 The Royal Society