## Abstract

We consider the synchronous whirl of arbitrary axisymmetric rotors supported on rigid bearings. Prior computational treatments of this problem were based on adding element-level gyroscopic terms to the governing equations. Here, we begin with a direct continuum formulation wherein gyroscopic terms need not be added on separately and explicitly: all gyroscopic effects are captured implicitly within the continuum elastodynamics. We present two new methods for obtaining the whirl speed, where we project the dynamic equilibrium equations of the rotor on to a few of its non-spinning vibration mode shapes. The first modal projection method is direct and more accurate, but requires numerical evaluation of more demanding integrals. The second method is iterative and involves a small approximation, but is simpler. Both the methods are based on one new insight: the gyroscopic terms used in other treatments are essentially the result of a prestress in the rotor caused by the non-zero spin rate, and may be incorporated as such in the continuum formulation. The accuracy of the results obtained, for several examples, is verified against detailed calculations with a commercial finite-element package, against our own nonlinear finite-element code or against analytical estimates. For further verification and illustration, a closed-form analytical solution for a simple problem, obtained using our method, matches the results obtained with explicit gyroscopic terms.

## 1. Introduction

Whirling of rotors is an important topic in the dynamics of machines. In many applications, rotors consist of shafts, possibly slender, and with or without discs or other masses attached to them. Practical instabilities or related problems in rotordynamics can involve details of the behaviour of bearings, material damping effects, unbalance and other complications (e.g. Den Hartog 1956; Dimentberg 1961; Rao 1983). In this paper, however, we concentrate on an idealized problem.

We study the synchronous whirl of an arbitrary axisymmetric rotor supported on ideal (rigid, massless and frictionless) bearings. The rotor is *not* assumed to be a shaft (i.e. a one-dimensional continuum) with or without objects attached to it. We consider genuine three-dimensional treatments for arbitrary axisymmetric cross sections, i.e. shapes such as hollow cones, bells and bottles. In particular, we consider a continuum formulation and develop two new modal projections for calculating the whirl speed. These modal projections can be used in conjunction with widely available commercial finite-element packages, with external numerical evaluation of some volume integrals.

We focus on the synchronous (forward) whirl, both for its simplicity and importance. We further restrict our attention to rotors where the spin-induced radial expansion effects produce a negligible change in the dimensions of the rotor. To help a general reader fix ideas about the kinematics, we note that the synchronous whirl is a special motion at a special speed. When viewed in a rotating coordinate system that spins about the undeformed axis at that special speed, the synchronous whirl gives a static, bent configuration. In the engineering practice, the synchronous whirl speed is usually the most important among the rotor's critical speeds.

We mention that our approach can in principle be extended to study non-rigid bearings, damping effects, asynchronous whirl, vibrations at non-critical rotation speeds and other important considerations. However, in the interest of clarity, brevity and focus, we have not discussed those issues here except for a single simple analytical example of asynchronous whirl in §8.

## 2. Contribution of this paper

Conventional formulations of rotordynamics use what are called gyroscopic terms (e.g. Genta 2005). A key aspect of our work is that we need not explicitly incorporate these gyroscopic terms, but still obtain correct results, as explained below.

The general governing equation for nonlinear elastodynamics of an arbitrarily moving body must remain true whether or not the body is a spinning rotor. Therefore, a correct three-dimensional continuum formulation for a spinning elastic rotor must implicitly capture any and all effects of the so-called gyroscopic terms commonly encountered in rotor-specific formulations. We use such a continuum formulation and do not convert our equations into typical rotor-specific forms (with gyroscopic terms). Neither do we confine our attention to problems where such gyroscopic terms play a small role. It will be clear that gyroscopic *effects* are duly and correctly accounted for.

We consider the continuum formulation and project the governing equations on to a few of the rotor's lateral, non-spinning, vibration mode shapes. The modal projections themselves are based on some physical insights, which we will discuss in detail. These modal projections constitute our main contribution.

We also use finite-element solutions of the continuum governing equations. Their role, however, will be solely to validate our modal projections.

Finally, we use Ansys (v. 8), which is not set up to directly and easily solve our rotor whirl problem. Our primary use of Ansys is restricted to mesh generation, modal analysis (non-spinning) and other non-whirling calculations required for our modal projections. We also use Ansys in a laborious, load-stepping calculation to further check the accuracy of our modal projection results (see electronic supplementary material).

## 3. Literature review

The literature on critical speeds of rotors, while extensive, focuses on shafts treated as beams (spatially one-dimensional). Several classical papers considered the problem of rigid discs on flexible shafts and bearings, with shafts modelled as the Euler–Bernoulli beams or simply as springs (e.g. Dunkerley 1894; Smith 1933). Many subsequent papers used more sophisticated theories and considered greater loading complexities but still modelled the rotor using beam or shaft elements. Eshleman & Eubanks (1969) included the effects of transverse shear, rotary inertia and gyroscopic moments and investigated the effects of an externally applied axial torque on the critical speeds. Similar considerations were made by Joshi & Dange (1976). The effects of viscous and hysteretic damping were also included in later work (e.g. Özgüven & Özkan 1984). Prohl (1945) used the transfer matrix method, based on the beam theory, for determining the critical speeds of rotors. The method has since been used widely for rotordynamics analysis (e.g. Flack & Rooke 1980; Sakata *et al*. 1983; Hsieh *et al*. 2006). Several papers have concentrated on finite beam elements of various kinds (e.g. conical or tapered), which include Rouch & Kao (1979), Nelson (1980), Genta (1985), Greenhill *et al*. (1985), Genta & Gugliotta (1988), Edney *et al*. (1990) and Gmür & Rodrigues (1991). Several other papers also studied Jeffcott rotors (a rotor with a point mass attached to a massless shaft) or simple variations thereof. Among them, Hassenpflug *et al*. (1980) experimentally studied the critical speed response of a Jeffcott rotor with acceleration. Genta & Delprete (1995) considered the effect of angular acceleration in the dynamic response of Jeffcott and other rotors passing through critical speeds. Ota *et al*. (1996) analysed the whirling and torsional vibration of a rotating massless flexible shaft, with a rigid disc, beyond its first critical speed. Muszynska (1996) studied a vertical, overhung, imbalanced rotor supported by flexible, anisotropic bearings. Genta & Brusa (2000) considered the role of non-synchronous damping on the free whirling and unbalance response of a Jeffcott rotor model. A noteworthy point in the above representative, if incomplete, review of the literature is that all these works involved modelling the shaft either simply as a massless spring or at most as being composed of elements of beams (spatially one-dimensional).

However, while such one-dimensional rotors are both common and important, there are rotor systems, such as specialized centrifuges, where a genuinely three-dimensional calculation is desirable. Such calculations can presently be done using commercial codes (with some effort; see electronic supplementary material) or specialized approaches (not widely available). Here, we present two simple modal projection-based methods for computing the whirl speeds of such rotors.

In prior work, there have been finite-element analyses of arbitrary rotor whirl with accounting of inertia terms arising out of the coupling between rotor deformation and spin rate (usually referred to as gyroscopic terms).

For example, Nandi & Neogy (2001) initially presented governing equations (their eqn (4)), which, in principle, match our own initially attempted (incorrect) modal projection for this problem, which we will discuss below. However, they subsequently added on an inertial term which we believe can give good results, but which is *ad hoc* in that it is not derived from underlying continuum equations; their added term also takes their approach away from ours, and a meaningful comparison between the two seems difficult without making specific kinematic assumptions which we have avoided in this paper.

Stephenson & Rouch (1989, 1993), who used axisymmetric harmonic elements, added a separately calculated matrix of gyroscopic terms (much in the spirit of the Timoshenko rotor analysis of Choi *et al*. (1992)). Our approach differs from that of adding on element-level gyroscopic terms (as in the work of Stephenson & Rouch) in several significant ways. The idea of adding on such gyroscopic (inertial) terms is shown to be correct using a verifying finite-strain calculation for Timoshenko rotors by Choi *et al*. (1992). However, as far as we know, such terms have not been derived starting from continuum-level formulations; rather, these terms are added based on engineering insights, and lead to specialized finite-element formulations that are not included with typical commercial finite-element packages. Here, in contrast to Stephenson & Rouch, we offer a continuum point-level, stress-based treatment that both keeps track of all terms dropped or retained and provides some fresh physical insights. Moreover, while they adopt a matrix reduction technique based on an explicit choice of master and slave degrees of freedom, we adopt a simpler strategy of projecting the governing equations directly onto a small number of lateral vibration modes. While they present a single approach for solution and compare their results with known formulae for cylindrical shafts, we consider several non-cylindrical geometries and cross-check our results using several different methods of solution. Lastly, while they further compare their calculations with experimental results for a particular rotor with bearing compliance and damping effects, we restrict our attention to computations for the ideal case but present more detailed numerical comparisons.

In §§6 and 7, we present in detail the two modal projection methods that form the main contribution of this paper. However, to verify our numerical results, we have also conducted alternative, cumbersome, geometrically nonlinear load-stepping calculations in Ansys, and we have further solved the problem using our own geometrically nonlinear finite-element code (which, as explained in §2, need not explicitly account for gyroscopic effects). These two calculations are described briefly in the electronic supplementary material and appendix A, respectively.

## 4. A preliminary calculation

Whirling is, in a sense, like buckling. In Euler buckling of columns (Timoshenko & Young 1935), a linear eigenvalue problem determines the buckling load. Nevertheless, we distinguish between the deformed and undeformed configurations (unlike in linear elasticity); we include in our equilibrium equation a term that is the product of load and displacement (equivalent to stress×strain, which would be treated as second order in linear elasticity); even past the buckling load, the unbuckled solution continues to coexist with the buckled solution (uniqueness results of linear elasticity preclude such solutions). Similarly, we expect the whirling speed to be determined by some linear eigenvalue problem. Yet, we distinguish between the whirling and non-whirling solutions; we retain terms that are linear in the whirling-associated displacements but quadratic in the rotation speed (technically of third order), and even past the whirling speed, the whirling and non-whirling solutions will coexist. Papers on rotors typically do not discuss this nonlinearity (a good discussion in the limited context of Timoshenko rotors is given by Choi *et al*. 1992). To clarify some aspects of this nonlinearity, we begin with a preliminary calculation that, though incomplete and incorrect by itself, carries over to a complete and correct calculation.

We assume that the deformed configuration (or shape) of the shaft can be expressed as a linear combination of a few of its lateral vibration mode shapes, and illustrate the calculation by considering only one mode *ϕ* (the fundamental). In later sections we include more modes as required.

The key advantage of describing the deformed whirling shape in terms of vibration mode shapes is that if the displacement is given by a vibration mode shape *ϕ*, then it involves a stress state *τ*, such that(4.1)where *ρ* is the material density and *ω*_{f} is the natural angular frequency of vibration in that mode. We use this in the following sections.

### (a) Dynamic equilibrium

We start with the dynamic equilibrium equation in reference coordinates (e.g. Gurtin 1981; Jog 2002),Here, ** F** is the deformation gradient;

**is the second Piola–Kirchhoff stress;**

*S**ρ*

_{0}is the undeformed density; and

**is the position vector of the material point of interest. We project this equation onto a single mode, linearize the resulting equation and obtain the incorrect answer that motivates this paper.**

*Χ*Consider a material point initially at position vector ** X** in a rotating frame that spins at the rotor speed. The displacement of this point is taken aswhere

*a*is an infinitesimal coefficient and

*ϕ*is the mass-normalized eigenvector.

In this paper, we adopt the St Venant–Kirchhoff stress–strain relation for nonlinear calculations, although here we linearize immediatelywhere *λ* and *μ* are Lame constants and is the Green strain tensor given by

The deformation gradient is

### (b) Virtual work

Considering a virtual displacement of *δaϕ*, for synchronous whirl, we havewhere *V* is the volume in the reference configuration and the angular velocity ** Ω** is directed along the undeformed centreline of the rotor.

Cancelling *δa* on both sides, we obtain a linear equation in *a*; and non-uniqueness of the whirling solution requires the coefficient of *a* to be zero (in a multi-mode projection, we would look for a singular coefficient matrix). Setting in the zero coefficient condition, we obtain(4.2)The term on the l.h.s., by equation (4.1), givessince the eigenvector is mass-normalized. Substituting this in equation (4.2), we obtain(4.3)Subsequently, the predicted critical speed is(4.4)where the unit vector is along the undeformed rotor centreline.

In equation (4.4), the natural frequency *ω*_{f} and mode shape *ϕ* can be determined using solid elements in a commercial finite-element package (in this paper, we used 10-noded tetrahedral elements in Ansys for all the calculations). It turns out that equation (4.4) is incorrect. For example, for a simply supported Timoshenko rotor of length four times the diameter, this formula captures only approximately half the difference between the bending frequency and the correct whirl frequency. In §6, we identify a single correction term that needs to be included in this naive calculation.

## 5. Simplifying insights

We observe that the Timoshenko rotor theory incorporates gyroscopic moments (i.e. terms of the form ), which are macroscopic (disc elementwise as opposed to continuum pointwise) *inertial* terms. It is possible to model the same whirling rotor in Ansys, and we have also modelled it using our own nonlinear finite-element code (see electronic supplementary material and appendix A, respectively). Both Ansys and our own code, however, use continuum equations and nonlinear displacement and stress terms, not macroscopic inertial terms. Yet, all three approaches will agree on results.

In our search for the bifurcation point (critical speed), since incipient whirling involves truly infinitesimal bending displacements, terms that are quadratic in them may be rigorously dropped. Moreover, terms that are nonlinear purely in the spin-induced displacements are likely to have a negligible physical effect, if the spin-induced geometry changes are small (at any rate, no radial expansion is considered in the Timoshenko theory). Note that this is a genuine physical approximation appropriate for the specific physical problem, although these terms are technically of order of unity, i.e. finite and non-zero. Finally, terms that couple the spin-induced displacements with the bending displacements are technically of first order in infinitesimals and some of them may play a crucial role in determining the whirling speed.

The spin-induced displacements appear to be important not because they are significant compared with the physical dimensions of the rotor, but because they are associated with a significant *stress* state that is in dynamic equilibrium when the rotor is straight. On infinitesimal bending, this stress state is infinitesimally disturbed from dynamic equilibrium and plays an infinitesimal but non-negligible role in the infinitesimal bending dynamics.

Incidentally, since the divergence of the spin-induced stress field is simply a centripetal body force field (countering an inertial force), it is intuitively, if not explicitly, seen how the inertial gyroscopic terms of the Timoshenko rotor theory might be equivalent to our nonlinear displacement and stress-based formulation. Moreover, for those using commercial code, the spin-induced stress field is easy to find by a single axisymmetric analysis, and the effects of this stress field can be largely incorporated by retaining it as a prestress while finding bending modes and frequencies.

These ideas are used to develop our modal projections in the following sections.

## 6. Modal projection method 1

We return to modal projections and again assume that the deformed configuration of the shaft can be expressed as a linear combination of a few lateral vibration mode shapes. But we now consider the effects of spin-induced stresses, and thereby obtain the correct answer.

The calculation is first described in detail for a single-mode projection. The incorporation of several modes is easy, and discussed briefly later.

### (a) Dynamic equilibrium

Again, we start with the dynamic equilibrium equation in reference coordinates,(6.1)As before, ** F** is the deformation gradient;

**is the second Piola–Kirchhoff stress;**

*S**ρ*

_{0}is the density in the undeformed configuration; and

**is the position vector of the material point of interest.**

*Χ*Unlike in §4, we now include the spin-induced displacements (denoted by *u*_{0}), and the displacement of a material point ** X** in the rotating frame is taken as(6.2)where

*ϵ*and

*a*are bookkeeping coefficients and

*ϕ*is a mass-normalized mode of lateral vibrations of the non-spinning rotor.

As mentioned in §4, we are interested in the coefficient matrix of terms that are linear in *a* (when that coefficient matrix is singular, infinitesimal whirling displacements are possible). This idea guides our simplifications below.

Starting again with the St Venant–Kirchhoff stress–strain relation, the second Piola–Kirchhoff stress is written aswhere *λ* and *μ* are Lame constants and is the Green strain tensor given by

However, as discussed in §5, the key nonlinear physical effect that contributes to the whirling speed is that of an infinitesimal disturbance (bending) of a pre-existing significant stress state (spin-induced). This disturbance is accounted for by ** F** in equation (6.1). Strain terms that are nonlinear in the displacement, in our opinion, play an insignificant role, and hence

**in equation (6.1) is here approximated using only linear terms.1 Accordingly, we take , whereWe can then split**

*S***, the second Piola–Kirchhoff stress, into bending and spinning components. The spinning component is given bythe bending component is given byand(6.3)**

*S*We now turn to the deformation gradientAs discussed in §5, the key term of interest involves the bending-induced disturbance of the spin-induced stress state *S*_{0}. This, consistent with neglect of the spin-induced changes in geometry, allows us to ignore *u*_{0} and write(6.4)It has recently come to our attention that substituting equations (6.4) and (6.3) into (6.1) gives equations that match those obtained by Bolotin (1963).

### (b) Virtual work

Considering a virtual displacement , we haveCancelling *δa* on both sides, and retaining the r.h.s. terms that are only linear in *a*, we have(6.5)plus terms that are not linear in *a*. Substituting equations (6.3) and (6.4) into the l.h.s. of (6.5), and retaining only terms that are linear in *a*, we obtain our modal projection(6.6)where *Ω*_{c} is the critical speed and is along the undeformed rotor axis.

Term B in equation (6.6), by equation (4.1), is known in terms of the natural frequency of vibration, as in §4,Such use of mode shapes is not a new idea (e.g. Bolotin 1963).

Finally, we consider term A in equation (6.6). This term, consistent with the qualitative discussion of §5, constitutes the only difference between equation (6.6) and the incorrect prediction of equation (4.4).

Term A involves derivatives of second order, and so direct evaluation would be possible if we used finite elements of sufficiently high-order shape functions (e.g. a 20-noded brick). Here, we transform the term for easier evaluation. For a second-order tensor field and a vector field ** v**, we have the identityUsing this and the divergence theorem, we have A as(6.7)where the unit vector

**(distinct from used previously) is normal to the surface of the rotor.**

*n*The surface *S* bounding the domain *V* can, in some problems, be split into a displacement-specified surface *S*_{u} and a traction-specified surface *S*_{t}. On *S*_{u}, *ϕ* is zero and hence the surface integral term corresponding to *S*_{u} is zero. On the traction-free surface *S*_{t}, we havebecause *S*_{0}** n** is zero on

*S*

_{t}.

Under more general restraints on the rotor, this surface term may turn out to be important and need accurate evaluation.

Now that all the terms in our modal projection (equation (6.6)) are known, *Ω*_{c} can be found. An analytical example from Ewins (2000) is used in §8 to illustrate that equation (6.6) exactly matches the results obtained using explicit gyroscopic terms. Numerical results for a few other rotors are given in §9.

### (c) Multi-mode projections

The previous calculations involved a single-mode projection. However, the method can be easily extended to multiple modes. We write the displacement of the rotor as (corresponding to equation (6.2))Here *ϕ*_{k} is the *k*th retained mode of vibration (lateral or otherwise) and *a*_{k} is an undetermined coefficient. Using *m* different virtual displacements corresponding to each mode *δa*_{k}*ϕ*_{k}, we would obtain *m* equations. Mass orthogonality of the mode shapes would give some simplifications. Eventually, the critical speed *Ω*_{c} would be found by solving an *m*-dimensional eigenvalue problem. Details of this routine calculation are omitted for brevity.

## 7. Modal projection method 2

The foregoing discussion allows us to present an alternative, iterative technique suitable for working with commercial codes such as Ansys. We start with an initial guess for the critical speed. We might, for example, start with *Ω*_{c,0} equal to the fundamental frequency of lateral vibrations of the non-spinning rotor. At the *k*th iteration, let the working guess be *Ω*_{c,k}. The iteration proceeds as follows.

(i) We do an axisymmetric and linear elastic calculation for the rotor in pure spin at a speed *Ω*_{c,k}. This gives a stress state which we call *τ*_{k}. (ii) We specify *τ*_{k} as a prestress (this step finds a wide industrial use) in the non-spinning rotor, and find new bending frequencies *ω*_{f} and mode shapes *ϕ*. (iii) Using these, we conduct the modal projection of §4 (i.e. equation (4.4)) to obtain a new estimate of the whirling speed, . (iv) We stop when is acceptably close to *Ω*_{c,k}.

Note that inclusion of the prestress makes this calculation different from that in §4, although the steps in the calculation are similar. It can be shown that this iterative procedure is in principle equivalent to carrying out the modal projection of §6, except that the boundary term of equation (6.7) is not retained, and that, in this case, the mode shapes correspond to a prestressed (though stationary) rotor, and may differ from those of the unstressed rotor. An interesting point is that this iterative solution will only find real whirling speeds, while the modal projection of §6 finds some imaginary whirling speeds as well (which correspond to imaginable whirling modes that are in fact suppressed due to the gyroscopic effects).

We now consider a detailed analytical example given in §8.

## 8. A detailed analytical example: Ewins's rotor

We consider a rotor from Ewins (2000), as shown in figure 1. Note that this simple system is easily treated using rigid-body mechanics; our formulation is really for finite-element-assisted calculations for deformable bodies with complex shapes. We select Ewins's rotor because it is the simplest system we know on which our ideas can be tested *analytically*.

The massless rigid shaft serves only to kinematically couple the disc's translation and tilt.

### (a) Ewins's solution (including explicit gyroscopic terms)

The equations of motion (Ewins 2000) in laboratory-fixed coordinates are(8.1)where *I*_{0} is the mass moment of inertia of the system about the *X*- or *Y*-axis; *J* is the polar mass moment of inertia of the disc–shaft system; and *M* is the mass of the disc. Matrix *G* in equation (8.1) is the gyroscopic matrix. The natural frequencies of the system at any spin speed, *Ω*, are obtained by Ewins from equation (8.1) aswhereIn particular, as mentioned above, the forward synchronous whirl speed is(8.2)

### (b) Our formulation (no explicit gyroscopic terms)

We start with the governing equation(8.3)Now from dynamic equilibrium of the rotor in pure spin, we have (using linear elasticity)where ** r** is the reference position vector of the point under consideration.

Choosing any virtual displacement *δ*** w**, the virtual work equation is(8.4)where the inconsequential zero subscript in

*ρ*has been dropped.

Following our procedure, we first determine the natural frequency and mode shape of the non-spinning system. For vibrations in the *X*–*Z* plane, the governing equation isThe natural frequency is . The mode shape is shown in figure 2. The natural frequency and mode shape of vibration in the *Y*–*Z* plane is similar to that in the *X*–*Z* plane.

We consider(8.5)where *θ* and *ψ* are small rotations about the *X*- and *Y*-axis, respectively (figure 2).

Next, we need the spin-induced stresses (*S*_{0}) in the disc. The disc is treated as a stiff but isotropic elastic body with Young's modulus *E* and Poisson's ratio *ν*=0 (for simplicity). The plane stress solution for a spinning disc is given in polar coordinates by Timoshenko & Goodier (1970, p. 81). The stress components are rewritten in Cartesian coordinates using formulae given in Timoshenko & Goodier (1970, p. 67). At a point (*x*, *y*) on the disc, we haveThe displacements and stresses are uniform throughout the thickness of the disc, which for integration purposes may be taken as unity.

Since *ϕ* is itself a vibrational mode (being a linear combination of two modes with the same frequency), we have (as discussed previously)

The acceleration iswhere the subscript ‘r.f.’ refers to the derivative taken in the rotating frame that spins at the shaft speed *Ω*. The velocity and acceleration of the point as seen in the rotating frame areUsing the above expressions and substituting angular velocity , we obtainWe consider two independent virtual displacements corresponding to variations *δψ* and *δθ*.

From equation (8.5), these virtual displacements are

Substituting the above relations into equation (8.4) and taking the two virtual displacements separately, carrying out the resulting two integrations, we obtain the following two equations (arranged in the matrix form):(8.6)where and . The equations (8.6) are written in a coordinate system rotating at the shaft spin rate. Substituting and , and transforming *x* and *y* to the stationary coordinates *X* and *Y*, respectively, measured in a laboratory-fixed frame, usingwe obtainDividing each of the above two equations by the constant , we obtain(8.7)Finally, we note that , and, whenceSubstituting the above equation into equation (8.7), we obtain exactly equation (8.1), showing analytically, for this example, that our formulation exactly captures all gyroscopic effects at all speeds and also describes asynchronous whirl and other related motions.

We emphasize that Ewins's problem is easily tackled by simple methods, and our calculation above is intended only to establish the theoretical validity of our approach. We now proceed to numerical results for a few other geometries.

## 9. Results and discussion

We now consider six different rotor geometries, as shown in figure 3. The relevant geometrical dimensions are

cylinder,

*L*=2 m,*D*=0.5 m,truncated cone,

*L*=2 m,*D*=0.5 m,*d*=0.2 m,bottle,

*L*=1.0 m,*h*=0.3 m,*D*_{1}=0.5 m,*D*_{2}=0.45 m,*d*=0.2 m,beam cylinder,

*L*_{1}=0.1 m,*d*=0.003 m,*L*_{2}=0.5 m,*D*=1.0 m,shell,

*L*_{1}=0.1 m,*d*=0.003 m,*L*_{2}=0.3 m,*D*_{1}=1.6 m,*D*_{2}=0.8 m,*D*_{3}=0.8 m,*t*=0.03 m,*t*_{1}=0.025 m, andbell,

*L*_{1}=0.1 m,*d*=0.003 m,*L*_{2}=0.8 m,*D*_{1}=1.8666 m,*D*_{3}=0.8 m,*t*=0.0416 m,*t*_{1}=0.025 m.

The material properties specified are: Young's modulus, 210 GPa; Poisson's ratio, 0.25; density, 7800 kg m^{−3}. The first two rotors are modelled as having simply supported end faces, approximately implemented in the finite-element model by restricting end-face nodes to have only axial displacements; the last four have their left faces fully fixed, with no other restraints.

From the first modal projection method, we recall the first integral in equation (6.7) (a surface term). A rough numerical estimate of this term, for the first rotor, worked out to only approximately 8% of the second term in equation (6.7), and so a more accurate integration was not carried out, and this surface term was entirely dropped for the second rotor. For the remaining four rotors, the surface term is zero.2

Predicted whirling speeds from the two modal projections described previously as well as from verifying calculations done using both Ansys (see electronic supplementary material) as well as our own nonlinear finite-element code (appendix A), for three of the six geometries, are given in table 1.

For the cylinder (simply supported), the Timoshenko theory is applicable, and the formula of Zu & Han (1992) with a shear factor of 0.9 gives a first-mode bending frequency of 1497.52 rad s^{−1}; this matches very well with our *ω*_{f} in table 1 (from ANSYS). The critical speed of the cylinder from the analytical formula (Zu & Han 1992) is 1545.38 rad s^{−1}. This compares well with the critical speeds obtained from all the methods. Among these, the result from our code (1546.46) is here taken as the most accurate one. The result from modal projection method 1, not including the boundary term of equation (6.7), is slightly lower and is close to the result from method 2, to which it is equivalent. By including the surface term in method 1, we obtain the value in parentheses (1546.76), more closely matching our nonlinear code.

Also, for the truncated cone, there are similarly small boundary terms, which we have dropped; results that are obtained resemble those for the cylinder.

For the bottle (as for the remaining three rotors), the boundary terms of equation (6.7) are known to be zero, because all points on the left faces of the rotors are held fixed in all three directions.

The laborious load-stepping calculations of Ansys (see electronic supplementary material) are slightly compromised due to the finite-sized imperfections used, and are presented here only to ensure that the results from other methods are good. These calculations were only carried out for the first three rotors. Also, only these same three rotors were studied using our own code (appendix A), which retains *all* terms while our interest here has been to identify the key terms necessary so that simple modal projections will yield a good answer.

We now turn to the remaining three rotors. These extreme geometries involve dynamics where the bulk of the mass moves like a rigid body, and essentially all deformations are restricted to narrow neck regions of negligible mass. For analytical estimates, these narrow necks may be thought of as massless beams, and the rest of the rotors thought of as rigid bodies. The corresponding rigid-body dynamics analysis is straightforward and not presented here. The final numerical results from those analyses are used here for cross-checking.

The results obtained are given in table 2. All modal projections in this table were done with two modes (actually, four modes in Ansys, with two modes for each frequency due to symmetry).

For each geometry, natural frequencies *ω*_{f} as obtained from Ansys are given in table 2, with the second beam-bending mode frequencies in parentheses. In each case, there are two lateral vibration modes, but only one corresponding whirling mode. It is also interesting to note that the whirl speed is approximately 2, 10 and 3 times the first bending frequencies for the three rotors, respectively; that is, the differences between the non-spinning vibration frequency and the whirl frequencies are large. We note that analytical estimates and modal projections give comparable but not identical results, especially for the last two geometries; this is probably because one of them involves external calculation of deformation gradients and the other does not, and so the calculations differ in detail, and possibly because one involves mode shapes without prestress and the other involves mode shapes with prestress. The prestress-based calculations are limited to seeking only real solutions and fail to capture the imaginary solutions (suppressed modes). Overall, the agreement is good.

## 10. Conclusions

In this paper, whirling speeds of arbitrary axisymmetric rotors have been considered. The aim was to estimate these speeds accurately using modal projections that can be used with routinely available commercial finite-element packages. Two such modal projection methods have been presented, both based on a single and simple new insight, i.e. the gyroscopic terms commonly used in rotordynamics analyses may be viewed as arising out of a state of prestress caused by the non-zero spin rate.

In addition to a simple example considered analytically in full detail, several numerical examples have been presented.

For three rotor geometries, the critical speeds have been computed using the two modal projection methods developed here as well as two different nonlinear finite-element calculations done for cross-checking. Three more rotor geometries have been studied with the same modal projection methods, and the results cross-checked against simple beam-plus-rigid-body models. Overall, the modal projection methods, with distinct strengths and levels of convenience, have given accurate estimates of the whirling speeds even for rotors where the whirling speed differs significantly from the lateral bending frequency.

## Acknowledgments

Sandeep Goyal's ME (Master of Engineering) project identified the incorrectness of the naive modal projection of §4, which led us to this more detailed investigation. A.C. thanks S. B. Koganti for initially bringing this problem to his notice, and A. Ruina and J. M. Papadopoulos for their helpful discussions.

## Footnotes

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2007.0139 or via http://journals.royalsociety.org.

↵Interestingly, the dropped strain terms that are nonlinear in the displacements turn out to be

*identical*with terms representing the effect of spin-induced configuration changes, which we also omitted in equation (6.4), in line with §5.↵Our implementation of simple supports for the first two rotors, through restricting end-face displacements to be purely axial, causes artificial end-face tractions when the rotor spins. An alternative modelling approach might be to allow the rotors to spin up without constraints and to apply the simple-support conditions only for subsequent bending. The surface terms would then be zero for the first two rotors as well.

- Received July 21, 2007.
- Accepted February 22, 2008.

- © 2008 The Royal Society