&,,,p,,~~,s

d StmctrursVol.40,No. 3.pa. 741-747.19% OMS-?949/91 53.00+O.M

PrinlcdioGreptBritin. 0 1991 F%rgamomRcu plc

FINITE ELEMENT ANALYSIS OF NATURAL WHIRL

SPEEDS OF ROTATING SHAFTS
LIEN-WEN COHENand DFR-MING
Ku
Department of Mechanical Engineering, National Cheng Kung Unive~ity, Tainan, Taiwan 70101,
Republic of China

(Received 3 May 1990)

Ab&mct-A three-nodal, Co Timoshenko beam finite element model is applied to analyze the natural whirl
speeds of a rotating shaft with different end conditions. The effects of translational and rotary inertia,
gyroscopic moments, bending and shear deformation are included in the mathematical model. In order
to facilitate the calculations of the natural whirl speeds, a rotating frame of reference is used in the
formulation. Comparison is made of the present finite etement analysis with other previously publish~
works available in the literature, with results presented for the case of forward and backward precession
to verify the gyroscopic effects. The numerical results show that the present finite element model has high
accuracy and good convergence.

IhTRODUCI’ION ness, many other important effects such as the rotary

inertia, shear deformation and gyroscopic moments
The design trend of modem rotating machinery are all neglected.
is toward higher rotating speeds and lower weight After Ruhl’s works, many investigators studied
by improving the prediction of rotor dynamic similar problems by introducing more different
behaviour. Therefore, the dynamic analysis of rotor- effects. For example, Nelson and McVaugh [7j have
bearing systems has become a technical field of ever developed a Rayleigh beam finite element model
increasing importance and has received considerable including the effects of rotary inertia, gyroscopic
attention by many engineers recently. There have moments and axial load for a flexible rotor system
been a number of studies relating to this field supported on elastic, viscously damped bearings.
in the past decades as indicated in the book by Zorzi and Nelson@] extended the work of [7j by
~imarogonas and Paipetis [I] and in the survey paper inco~rating the effects of both internal viscous
by Rieger[2]. In those published studies, the main and hysteretic damping in the same model.
aspects of rotor-bearing dynamic behavior are the Nelson [9] has also utilized Timoshenko beam theory
rotor vibration due to unbalance forces and different for establishing the shape functions. Based on the
self-excited sources, rotor stability and torsional shape functions, he derived the system matrices of
dynamics of drive train. However, of the many governing equations. The governing equations can
important dynamic characteristics of rotor-bearing be quickly changed into a Rayleigh beam or a
systems, the most commonly predicted for design Euler beam model if the shear parameter in the
purposes are the critical speeds and the unbalance shape functions is set to equal zero. C)zgiiven and
response. &kan [IO] presented the combined effects of shear
In the past, various mathematical methods for deformations and the internal damping in their finite
transverse vibration analysis of rotor-bearing systems element formulation to analyze the natural whirl
have been successfully developed. These methods speeds and unbalance response of multibearing
may be divided into two most prevalent classes. One rotors. Many other works which utilized the finite
is the transfer matrix method and the other is the element technique to the analysis of rotor dynamics
finite element technique (or direct stiffness approach). can be found in Reefs [ 1l-131. Those works showed
The transfer matrix approach allows for a continuous that the finite element technique can provide an
representation of the shaft section and produces accurate modeling of rotor-bearing systems and make
results in good agreement with the experiments [3,4]. it possible to formulate ever increasingly complex
This method has also the advantage of small problems encountered in practical engineering design;
computer memory requirements, but the equations also, the finite element solutions can systematically
of motion are not explicitly written. On the other and easily be obtained and yield high successful
hand, Ruhl[5] and Ruhl and Booker [6] are the first results.
to appIy the finite element technique to study Although the effects of various factors on the
the dynamics of rotor systems. However, the finite dynamics of rotor-bearing systems have been investi-
element fo~ulations contained in their works only gated via the improved finite element model by
included the translational inertia and bending stiff- several authors during the past 20 years, a common

741
142 LIEN-WEN CHEN and DER-MING Ku

feature of the published works is that those finite B(x, t) = -8 W*(x, t)/ax,
elements belong to a C’ class and, also, in order to
take into account the effect of shear deformation, a T(x, t) = 8V,(x, r)/ax. (2)
shear parameter is included in the shape functions.
The objective of the present article is to develop a The potential energy u’ of a rotating shaft element
Co Timoshenko beam finite element model to study of length I, including the elastic bending and shear
the natural whirl speeds of a rotating shaft with deformation energy, is expressed as
different end conditions. In the present Co finite

s‘Ez{(B’)2+(r)21dx
element model, each element has three nodes and
each node has four degrees of freedom: two trans- &;
lations and two rotations. Numerical examples are 0

given to demonstrate the accuracy and convergence

of the present finite element model. +; ‘KGA{(V:)r+(W:)}dx, (3)
I0
FINITE ELEMENT FORMULATION
where the prime denotes differential with respect
A uniform shaft, rotating at a constant speed R, of to axial distance x, E the Young’s modulus, I the
length L is illustrated in Fig. 1. Two reference frames second moment of inertia, K the shear coefficient, G
are employed to describe the system motion. One is the shear modulus, A the cross-sectional area of the
the fixed reference X-Y-Z and the other is a rotating shaft, and
reference x-y-z. The X and x axes are colinear and
coincident with the undeformed rotor center line. The v:=avlax-r, w:=aw/ax+~. (4)
two reference frames are defined by a single rotation
at difference about X axis with o denoting the whirl
The kinetic energy T’ of the rotating shaft element,
speed.
including the translational and rotational forms, is
It is assumed that the axial motion is small and
given by
can be reasonably neglected, therefore, a typical
cross-section of the shaft, located at a distance x
from the left end, in a deformed state is described by
the translations V(x, t) and W(x, t) in the Y- and
Z-directions and small rotations B(x, t) and r(x, t)
+
s
‘PA{(P)‘+(@)2)dx
0

about Y and Z.
+; ‘1d{&2+(f)2}dx
The two translations (V, IV) consist of a contri- s0
bution (V,, IV,) due to bending and a contribution
(V,, IV,) due to shear deformation. The rotations
(B, r) are related to the bending deformations -R ‘l,(i)Bd_x+; ‘1$x, (3
s0 s0
(Vb, IV,). The relationship can be written as
where the dot denotes differential with respect to time
w, r) = Vb(X,t) + Vs(x,t), r, p the mass density of the shaft material, I,, and ZP
the diameter and polar mass moments of inertia per
Wx, r) = W&, r) + Ws(x,r), (1) unit length.

Fig. 1. Displacement variables and coordinate systems of a rotating shaft.

 Natural whirl speeds of rotating shafts 743

In the finite element method, the continuous The performance of the Timoshenko beam element
displacements may be expressed in terms of the is improved if the stiffness terms due to shear defor-
discretized nodal displacements. In the present finite mations are evaluated by using reduced integration.
element model, the displacement field in an element This helps to avoid locking that occurs when the
e is approximated as beam is thin [ 141.Therefore, in the present paper, the
order of integration of the shear stiffness matrix [K,r
{ V(x, t), W(x, t). B(x, r), T(x, r)} is taken as 2 and 3 for the other matrices. The shear
coefficient K is taken as 0.85.
Substituting the element potential energy and the
element kinetic energy, respectively given by eqns (9)
and (IO), into Hamilton’s principle, the following
or matrix equation of motion for the finite rotating shaft
element is produced

(]Mrr + Pf~l%li}~- fJGl%P + ]Kl’{qJ’ = 101,

(11)

where [Gr = [H]' - (H]"is the element gyroscopic

matrix.
where Ni(x) is the one-dimensional quadratic Equation (11) describes the element motion to
Lagrangian shape function, {qje = {V,, IV,, B,, the fixed frame coordinate X-Y-Z and has the
r,,.... V,, IV,, I&, T3jT is the nodal displacement property that the X-Y and X-Z planar motions are
vector, [N,(x)] and [N,(x)] are the translational and coupled due to the skew symmetric gyroscopic
rotational shape function matrices, respectively. matrix, (Cr. However, it will be convenient to
From eqn (7), (4) is related by the nodal displace- utilize a rotating reference frame to analyze the
ment vector as natural whirl speeds or the unbalance response of
a rotating shaft with isotropic supports. To this
end, the equation of transformation between fixed
(41’= b’s(xNq 1’. (8) and rotating frame coordinates, shown in Fig. 1, is
given as
With the aid of eqns (7) and (8), the element
potential energy II’ and the element kinetic energy T’ ISI’ = VW 1’ (12)
can be written in terms of the nodal displacement
vector as, respectively where
Vl PI 101
Lr. = f141’r(]KJ + ]K,l%?j’= ;{4j’r]KY(q1’ (9) [RI = WI Vl WI ’
[ PI PI P-l I
and

1
cos wt -sinot 0 0
T’ = t{4}rr(]Mrr + Pf,J){4}’
sin wt cos cot
PI= o o O O .
- n(qy[Hl’{qj + ~z,n’ (10) cos wt -sin ot
[ 0 0 sin rdt cos 01
where
Introducing the transformation equation (12)
into eqn (11) and premultiplication by [R]r, the
element equation of motion described in rotating
frame coordinates is given

+ (1 - i)[G]3{6}‘+ ([Kr - u’([M,]

+(1-2~)[M,1’)){~}‘={0}, (13)

where [fir)’ = [R]qMt]‘fk]/w, i = n/o is the whirl

ratio.
Assembling the contribution of each element
equation of motion and assuming the supports of the
rotating shaft to be undamped isotropic (in such a
744 LIEN-WEN
CHENand DER-MINGKu

Table 1. Natural frequency parameters (o*, ti* = pAL’02/Ei) for a nonrotating shaft with
various end conditions and slenderness ratios ri2L
C-F H-H C-H c-c
rj2L 1151 1151 1151 1151
0.1 10.42 10.45 70.38 70.35 122.9 125.3 191.7 199.4
0.04 12.00 12.00 91.58 91.52 206.5 207.1 397.1 400.7
0.02 12.27 12.27 95.92 95.86 229.3 229.2 470.8 470.7
0.01 12.34 12.34 97.08 97.02 235.8 235.5 493.9 492.7
0.002 12.37 12.36 97.46 97.40 238.0 237.4 501.9 500.3
C, clamped boundary condition.
H, hinged boundary condition.
F, free boundary condition.

case {S)< equals to a constant), the global finite 0.1 to 0.002 is studied first. A five element idealization
element equation is reduced to of the shaft is used to obtain the numerical results
and they are listed in Table 1, where r is the radius
@I - 4wTl+ (1 - 2~)WA!lN{~)
= w. (14) of the shaft. From Table 1, it can be seen that the
present results agree very well with those of Ref. [15].
It becomes an eigenvalue problem and its non- It is also seen that the natural frequencies of a
trivial solutions are sought if a specified value of the nonrotating shaft is considerably affected by the
whirl ratio I is given. Equation (14) has also the effect of shear deformation and rotary inertia in
property that all the matrices are symmetric and the the case of hinged-clamped or clamped-clamped
motions in x-y and X-Z planes are uncoupled. There- shafts compared to the shafts with other end
fore, for computational purposes, it can simplify conditions.
eqn (14) of a two-planar eigenvalue problem into The second comparison is for the primary forward
an one-planar eigenvalue problem. Then eqn (14) (I = + I) and backward (3. = - 1) critical speeds
becomes for a hinged-hinged rotating shaft with the slender-
ness ratios varying from 0.02 to 0.10. The present
([1y]+- w*([Mr]* + (1 - ~~)[~~I*))~~~* = (O], (15) results are compared to those published by Eshleman
and Eubanks[I6] (equations 41 and 42 in [16])
where the order of the matrices in eqn (15) is reduced and those by Nelson [9]. These comparisons are
to half of those in eqn (14). presented in Table 2. The results from Eshleman and
Eubanks in Table 2 are the approximate solutions of
NCMERICAL EXAMPLES AND DISCUSSION
the differential equations in which the interaction
effects between the transverse shear and rotary inertia
In order to evaluate the accuracy of the present and the transverse shear and gyroscopics are not
finite element model, analysis of the natural frequen- included. In order to demonstrate the good conver-
cies of a uniform nonrotating shaft with various gence of the present finite element model, the critical
end conditions and slenderness ratios varying from speeds of the rotating shaft represented by various

Table 2. Primary critical speeds @,p’= pAL’02/EI) of a hinged-hinged rotating shaft with various slenderness ratio r/2L
i. = + 1 (first mode) i. = + I (second mode)
--
present present
r/2L (3)* (5) (7) [91 [161 (3) (5) (7) [91 [161

0.02 3.1392 3.1360 3.1356 3.1374 3.1373 6.3371 6.249 1 6.2383 6.2532 6.2489
0.04 3.1208 3.1177 3.1173 3.1246 3.1245 6.1816 6.1016 6.0918 6.1551 6.1515
0.06 3.0901 3.0870 3.0867 3.1027 3.1037 5.9373 5.8688 5.8604 5.9873 6.0045
0.08 3.0474 3.0445 3.0442 3.0715 3.0757 5.6415 5.5848 5.5778 5.7623 5.8250
0.10 29939 2.9912 2.9909 3.0311 3.0416 5.3350 5.2881 5.2822 5.5069 5.6287
i = - 1 (first mode) A = - 1 (second mode)
present present
r/2L (3) (5) (7) 191 I161 (3) (5) (7) f91 1161

0.02 3.1271 3.1240 3.1236 3.1253 3.1251 6.2441 6.1603 6.1500 6.1631 6.1560
3.0766 3.0736 3.0733 3.0796 3.0780 5.9107 5.8404 5.8318 5.8748 5.8387
8:: 3.0032 3.0005 3.0001 3.0125 3.0067 5.5351 5.4769 5.4698 5.5380 5.4493
0.08 2.9172 2.9147 2.9144 2.9328 2.9193 5.1847 5.1356 5.1296 5.2146 5.0670
0.10 2.8264 2.8242 2.8239 2.8475 2.8234 4.8763 4.8338 4.8286 4.9238 4.723 1
+ Values in brackets denote the numbers of elements in the finite element formulation.
 Natural whirl speeds of rotating shafts 745

numbers of elements are also indicated in Table 2. As The effects of the end conditions on the primary
compared to Ref. [ 161, the numerical results shows forward and backward critical speeds of a rotat-
that for the first mode, accuracy is well within 2% for ing shaft are graphically illustrated in Fig. 2 for
all cases and under 5% for the second mode. The the first mode and Fig. 3 for the second mode,
accuracy decreases slightly as the slenderness ratio where p = (JTAL~o~/EI)‘~~ is the nondimensional
increases. whirl speed parameter. As compared to the natural

I I I I I

(b)

--- Bernoulli-Euler --- Bernoulll-Euler

a Whirl Ratlo x Whirl Ratlo
1.70 - I I I I I I I I I I
0.00 0.32 0.04 0.06 0.08 0.10 0.12 0.00 0.02 0.04 0.06 0.08 0.10 0.\2

Slenderness Ratlo -& Slenderness Ratlo-&

I I I I

F----q
I

(d)

i -

--- Bernoulli-Euler
1 Whlrl Ratlo
3.3 I I I I I
0.00 0.32 0.04 0.06 0.08 0.10 0.12 0.00 0.02 0.04 0.06 0.08 0.10 0.12

Slenderness Ratlo & Slenderness Ratlo &

Fig. 2. First natural nondimensional whirl speeds (frequencies) of a rotating shaft with various end
conditions. (a) clamped-free; (b) hinged-hinged; (c) clamped-hinged and (d) clampcd-hmped.
 LIEN-WI3 CHENand I&R-MINQKu

6.4

6.2
rr,
Z
g 6.0

a
$ 5.8

2 4.8
=: --- Bernoulli-Euler
x
4.6 --X Whirl Ratio
A Whlrl Ratio
I I I I I 4.4 A
0.00 0.02 0.34 0.06 0.08 0.10 0.12 0.00 0.02 0.04 0.06 0.08 0.10 0.12

Slenderness Ratio -& Slenderness Ratio &

7.2 I 1 I I
Cd) ’
7.0

,” 6.8 -
5
z 6.6 -
z
Lt.
_ 6.4 -
z
.z 6.2 -
E
E” 6.0 -

z
z” 5.8 -

=: 5.t -
3
z” 5.4 -
n
E 5.2
ki
ln --- Bernoulli-Euler
5.0
X Whlrl Ratlo
4.8 1
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.00 0.02 0.04 0.06 0.08 0.10 0.12

Slenderness Ratio -& Slenderness Ratio $

Fig. 3. Second natural nondimensional whirl speeds (frequencies) of a rotating shaft with various end
conditions. (a) clamped-free; (b) hinged-hinged; (c) clamped-hinged and (d) clampedclamped.

frequencies of a nonrotating shaft (1. = 0) from the effects of both shear deformation and rotary
those Figures, it is found that the natural inertia.
whirl speeds are more affected by the hinged- CONCLUSIONS
hinged or clamped-free end conditions than the
others. As the slenderness ratio increases, the A Co Timoshenko beam finite element fo~ulation
natural whirl speeds are significantly affected by is presented to study the natural whirl speeds of
 Natural whirl speeds of rotating shafts 741

a rotating shaft with various end conditions and distributed parameter turborotor systems. J. Engng Ind.
ASME 94, 128-132 (1972).
slenderness ratios. The numerical results indicate that
7. H. D. Nelson and J. M. McVaugh, The dynamics of
the natural whirl speeds are considerably affected by rotor-bearing systems using finite elements. J. Engng
the end conditions of rotating shafts and the effects Znd. 98, 593+500 (1976).
of both shear deformation and rotary inertia. It 8. E. S. Zorzi and H. D. Nelson, Finite element simulation
also demonstrates that the present finite element of rotor-bearing systems with internal damping.
J. Engng Power ASME 99, 71-76 (1977).
model can provide an accurate representation of 9. H. D. Nelson. A finite rotating shaft element using
rotating shaft systems and has good convergence. As Timoshenko beam theory. J. M;?ch. Des. ASME 102;
expected, this model can be further used to study the 793-803 (1980).
unbalance response, stability and other dynamic 10. H. N. Gzgiiven and Z. L. Gzkan, Whirl speeds and
unbalance response of multibearing rotors using finite
problems or rotor-bearing systems.
elements. J. Vibr. Acoust., Stress, Reliab. Des. ASME
106, 72-79 (1984).
REFERENCES 11. K. E. Rouch and J. S. Kao, A tapered beam finite
element for rotor dynamics analysis. J. Sound Vibr. 66,
1.A. D. Dimarogonas and S. A. Paipetis, Analytical 119-140 (1979).
Methods in Rotor Dynamics. Applied Science, New 12. Y. D. Kim and C. W. Lee, Finite element analysis of
York (1983). rotor bearing systems using a modal transformation
2. N. F. Rieger, Rotor-bearing dynamics-state-of-the- matrix. J. Sound Vibr., 111, 441456 (1986).
art. Mech. Mach. Theory 12, 261-270 (1977). 13. G. Sauer and M. Wolf, Finite element analysis of
3. J. W. Lund and F. K. Orcutt, Calculations and exper- gyroscopic effects. Finite Elements in Analysis and
iments on the unbalance response of a flexible rotor. Design 5, 131-140 (1989).
J. Engng Ind. ASME 89, 785-796 (1967). 14. E. Hinton and D. R. J. Gwen, An Introduction to Finite
4. J. Gu, An improved transfer matrix-direct integration Element Computations. Pineridge Press, Swansea (1979).
method for rotor dynamics. J. Vibr. Acoust., Stress, 15. B. P. Shastry and G. V. Rao, Dynamic stability of bars
Rehub. Des. ASME 108, 182-188 (1986). considering shear deformation and rotatory inertia.
5. R. L. Ruhl, Dynamics of distributed parameter rotor Comput. Struct. 19, 823-827 (1984).
systems: transfer matrix and finite element techniques. 16. R. L. Eshleman and R. A. Eubanks, On the critical
Ph.D. dissertation, Cornell University (1970). speeds of a continuous rotor. J. Engng Ind. ASME 91,
6. R. L. Ruhl and J. F. Booker, A finite element model for 1180-1188 (1969).