CHAPTER 7

TORSIONAL VIBRATIONS OF ROTORS-II:

THE CONTINUOUS SYSTEM AND FINITE ELEMENT METHODS

In previous chapter, we considered torsional vibrations of rotor systems by the Newtons second law
of motion and by the transfer matrix method. Especially the transfer matrix method proved its
applicability for analysing even large rotor systems in more systematic (algorithmic) fashion. The
advantage of the method is that the size of the matrices does not change with the degrees of freedom
of the system. The only disadvantage of the method is that natural frequencies have to be obtained for
large systems using root searching numerical technique. Moreover, it has some other short comings
that we will observe especially when we will be analysing transverse vibrations in the subsequent
chapters. This method requires special treatment when rotor systems have multiple supports, and
when the shaft is treated as a continuous (i.e., when the mass and the stiffness of shaft are distributed
throughout its span). In the present chapter, we will firstly be analysing torsional vibrations using the
analytical approach by treating the shaft as a continuous, which has infinite DOFs. For this case the
governing equation becomes partial differential equation (i.e., identical to the wave equation). For
simple boundary conditions, the governing equations are solved in the closed form. However, for
complicated rotor-support systems the continuous system is very difficult to analyse by analytical
method. For analysing complex systems by an approximate method, we have chosen the procedure of
the finite element method (FEM). The FEM is well known for its versatility and flexibility to analyse
the complex rotor and support systems. The basic steps involved in the solution of torsional vibrations
by FEM are presented. Variety of cases (a simple to multi-disc systems, geared systems, branched
systems, reciprocating engines, etc.) is illustrated through simple numerical examples. The only
limitation of the FEM as compared to TMM, as we will see, that with the degrees of freedom of the
system the size of the matrices to be handled increases. Large size matrices require larger computer
memory, higher computational cost and inherent ill-conditioning problem (i.e., possibility of
becoming nearly singular matrices). The present chapter will pave the way for analysing more
complex transverse vibration problems in rotor-bearing systems in subsequent chapters.

7.1 Torsional Vibrations of Continuous Shaft Systems

Till now a flexible shaft is considered as massless (Fig. 7.1a) and a long heavy shaft is treated as rigid
(Fig. 7.1b) for the transverse and torsional vibrations. However, there are cases where the mass (i.e.,
the mass, the diametral and polar mass moments of inertia) and the stiffness (e.g., the transverse, the
longitudinal or the torsional) both are distributed along the shaft, and such flexible heavy rotors (Fig.
7.1c) falls under the category of continuous systems. Figure 7.1 is for transverse vibrations since it is

362

easy to visualise the motion, however, the shaft models are still valid for other types of vibrations
(i.e., the torsional and axial vibrations).

Fig. 7.1 Different shaft models in transverse vibrations of a rotor system

(a) a flexible shaft (b) a rigid shaft (c) a continuous shaft

In continuous systems, equations of motion takes the form of a partial differential equation with
spatial (e.g., z in the present one-dimensional model; however, for the general solid model it can have
x, y, and z) and temporal (i.e., t, time) derivatives. Equations of motion of a continuous system can be
derived by (i) the force and the moment balance of a differential element (ii) Extended Hamiltons
principle (iii) Lagranges method. While the force and the movement balance is a convenient
approach for simple rotor systems, the Hamiltons principle and the Lagranges method are applied
for complex rotor systems. These two approaches need the consideration of the potential and kinetic
energies of the system. The Hamiltons principle and the Lagranges method will be briefly discussed
here, for more detailed information readers are referred to excellent texts by Meirovitch (1986), Rao
(1992), and Thomson and Dahleh (1998).

7.1.1 The Hamiltons principle

The extended Hamiltons principle is stated as an integral equation in which the energy is integrated
over an interval of time. Mathematically, the principle can be stated as

t2

( T + W ) dt = 0

(7.1)

t1

where is the variational operator (it is similar to commonly used the differential operator in most of
mathematical operations), T is the variation of the kinetic energy, and W is the virtual work that
includes contributions from both the conservative and non-conservative (e.g., dissipative) forces, and
can be expressed as

W = Wc + Wnc

(7.2)

363

where subscripts c and nc represent the conservative (e.g., the strain (potential) energy in elastic
bodies, gravitational potential energy, etc.) and non-conservative (the energy to due damping forces,
virtual work due to external forces, etc.) parts. The extended Hamiltons principle is quite general and
can be used to derive governing dynamic equations for a system of particles, for a system of rigid
bodies, and continuous (distributed parameter) systems. In deriving the Hamiltons principle (it is
valid for conservative dynamic systems only), it is assumed that virtual displacements must be
reversible. It put restriction that the constraint forces must not perform any work, which means the
principle can not be used for systems with friction forces. The virtual work performed by the
conservative forces can be expressed as the negative of variation of the potential energy, i.e.

Wc = U

(7.3)

where U are the potential energy of the system. Defining the Lagrangian, L, as
L = T U

(7.4)

Equation (7.1) takes the following form

t2

( L + W )dt = 0
nc

with

L = T U

(7.5)

t1

Fig. 7.2 The actual and variational paths of the motion

In a special case, when non-conservative forces are zero. Equation (7.5) gives

t2

Ldt = 0
t1

(7.6)

364

which is the Hamiltons principle for conservative systems. The physical interpretation of equation
(7.6) is that out of all possible paths of motion of a system during an interval of time from t1 to t2 (Fig.
7.2), the actual path will be that for which the integral has a stationary value.

It can be shown that the stationary value will be, in fact, the minimum value of the integral. The
Hamiltons principle can yield governing differential equations as well as boundary conditions. For
illustration of the method, a very simple, however, representative case of a single-DOF spring-massdamper rotor system is considered through an example.

Example 7.1 For a single-DOF spring-mass-damping rotor system (Fig. 7.3) subjected to an external
force; obtain the equation of motion using the extended Hamiltons principle. The shaft stiffness and
damping is represented by kt and ct, respectively; and the polar mass moment of inertia of the disc is
I p.

Fig. 7.3 (a) A single-DOF spring-mass-damping rotor system (b) A free body diagram

Solution: For illustrating the extended Hamiltons principle a simplest single-DOF spring-massdamping rotor system is considered, subsequently which will help in the formulation of governing
dynamic equations for continuous systems.

For such a system, the kinetic, T, and strain, U, energies can be written as

T = 12 I p
z2

and

U = Wc = 12 kt z2

so that

T = I p
z (
z )

and

U = Wc = kt z ( z )

(a)

It should be noted that the variation operator is similar in operation to the differential operator. And
the non-conservative virtual work done is (Fig. 7.3b)

Wnc = TE (t ) z (ct
z ) z

(b)

365

where TE(t) is an external torque. The virtual work done by the external torque will be a positive,
however, the virtual work due to the damping force will be a negative. The extended Hamiltons
principle is written as

t2

( T + W

+ Wnc ) dt = 0

(c)

t1

or
t2

( T U + W ) dt = 0

(d)

nc

t1

On substituting equations (a) and (b) into equation (d), we get

t2

{I
(
) k ( ) + T
p

(t ) z (ct
z ) z } dt = 0

(e)

t1

On integrating by parts of only the first term of equation (e), we get

t2

t2

I p
z ( z ) {I p

z + kt z TE (t ) + ct
z } ( z )dt = 0

t1

(f)

t1

The first term of equation (f) will vanish, since the variation is not defined in temporal domain; and
the second term will give (since the variation z is arbitrary)

I p

z + ct
z + kt z = TE (t )

(g)

which is a standard equation of motion of the single-DOF spring-mass-damper rotor system.

7.1.2 Lagranges equation

The Hamiltons principle is stated as an integral equation, where the total energy is integrated over a
time interval. On the other hand, the Lagranges equation is differential equations, in which one
considers the energies of the system instantaneously in time. The Hamiltons principle can be used to
derive the Lagranges equation in a set of generalized coordinates, i. The Lagranges equation can be
written as

d T

dt
i

T
= qi (t ) ,

i = 1, 2, , N

(7.7)

366

with

qi = qci + qnci

(7.8)

where q(t) is the torque and subscripts: c and nc represent the conservative and the non-conservative,
respectively. Equation (7.7) is Lagranges equations of motion in their most general form. The
conservative generalised torques have the following form

qci =

U
,
i

i = 1, 2, , N

(7.9)

where U is the potential energy. Substituting equation (7.9) into equation (7.7), we get

d T

dt
i

T U
+
= qnci (t )

i i

i = 1, 2, , N

(7.10)

Since the potential energy does not depend upon the velocity, equation (7.10) can be written as

d L L
= qnci (t )

dt
i i

i = 1, 2, , N

(7.11)

where L = (T U ) is the Lagrangian. Lagranges equations are mainly used for discrete systems
including the rigid body. However, for continuous systems the Hamiltons principle is preferred over
Lagranges equations.

Example 7.2: For a single-DOF spring-mass-damping rotor system subjected to an external torque,
obtain the equation of motion using the Lagranges equation.

Solution: From previous example, we have

L = T U = 12 I p
z2 12 kt z2

and

qnci (t ) = TE (t ) + (ct
z )

(a)

From equation (a), we have

d L

dt
z

d
= ( I p
z ) = I p

z
dt

and

L
= kt z
z

On substituting equations (a) and (b) into the Lagranges equation (7.11), we get

(b)

367

I p

z + kt z = TE (t ) ct
z
or

I p

z + ct
z + kt z = TE (t )

(c)

which is a standard equation of motion. This example had motive to clear the concept of application
of the Lagranges equation in the simplest form.

7.1.3 Governing differential equations

Governing equations for torsional vibrations of continuous shafts are much simpler as compared
transverse vibrations. It should be noted that the former is identical to axial vibrations of rods (by
replacing the torsional angular displacement with the axial linear displacement, and the polar mass
moment of inertia with the mass). To obtain the governing differential equations, the variational
principle (i.e., the Hamiltons principle) is used that requires the system potential energy functional.
The potential energy functional is obtained by considering the displacement field of a point P in the
shaft cross section when the external torque is applied onto the shaft, as shown in Figure 7.4. Let T(z,

t) be the distributed external toque applied in some span of the rod, and T0(t) be the concentrated
external torque applied at z0. Let z(z, t) is the angular displacement of any plane at location z about zaxis and L is the total length of the rod. Consider a point P (Fig. 7.4) in the cross section at an axial
position of z from the origin of the coordinate system x-y-z. From Fig. 7.5, we have following
geometric relations

sin =

y
r

and

cos =

x
r

(7.12)

where (x, y, z) is the Cartesian coordinate of the point P, whereas (r, , z) is the cylindrical coordinate
of the same point.

Figure 7.4 A continuous rod under time dependent external torques

368

Figure 7.5 (a) Displacement of a point on the cross-section (b) its geometrical details

Derivation of displacement, strain and stress fields: Due to pure twisting of a plane the point P moves
to P1 (Fig. 7.5(a)) so that PP1 = rz for angular displacement. The displacement field of the point P,
located on the cross-section at an axial distance of z from the origin, can be expressed as

u x = PQ = r z sin = y z ( z , t ) ;

u y = PQ
= r z cos = x z ( z, t ) ;
1

uy = 0

(7.13)

where ui(x, y, z) is the linear displacement of the point with i = x, y, z. The strain field is given by

xx = u x , x = 0 ;

yy = u y , y = 0 ;

yz = 12 ( u y , z + u z , y ) = 12 x z , z

where u x , x =

zz = u z , z = 0 ;
and

xy = 12 ( u x , y + u y , x ) = 0

zx = 12 ( u z , x + u x , z ) = 12 y z , z

(7.14)

u x
u
, u x , y = x , etc. The stress field is given by
y
x

xx = E xx = 0 ;
yz = G yz = 12 Gx z , z

yy = zz = 0 ;
and

xx = G xy = 0 ;
zx = G zx = 12 Gy z , z

Derivation of various energy terms: The strain energy is given by

(7.15)

369

U = ( 12 xx xx + 12 yy yy + 21 zz zz + xy xy + yz yz + zx zx dV
V

(7.16)
L

1
2

Gy 2 z2, z + 12 Gx 2 z2, z dAdz =

1
2

0 A

with

GJ

2
z,z

dz

J = x 2 + y 2 dA
A

where J is the polar moment of inertia of the shaft cross section, V is the volume, and A is the area of
shaft cross section. The velocity components, noting equation (7.13), can be written as

u
x = y
z ( z, t );

where u
x = u x ,t =

u
y = x
z ( z , t );

u
z = 0;

(7.17)

u x
, etc. The kinetic energy is given by
t
L

T = 12 u
x2 + u
y2 + u
z2 dAdz =

1
2

0 A

J
dz
2
z

(7.18)

where is the density of the shaft material. The external work done (i.e., the non-conservative energy)
for a distributed external torque, T(z, t), and for a concentrated torque, T0(t) , as shown in Figure 7.4, is
given by

Wnc = T + T0 * ( z z0 ) z dz

(7.19)

where * is the Direc delta function, which is defined as

+

( z z0 )dz = 1 for z = z0, and

= 0 for z z0

(7.20)

Derivation of equations of motion and boundary conditions: On substituting various energy terms
from equations (7.21), (7.22), and (7.23) in the extended Hamiltons principle (equation (7.24)), it
gives

370

t2

t2

( T U + W ) dt =
nc

t1

t1

J
z
z GJ z , z z , z + T ( z , t ) z + T0 * ( z z0 ) z dzdt = 0 (7.25)

which can be written, in the following form

t2

t1

t2 L
t2 L

T ( z, t ) z + T0 * ( z z0 ) z dzdt = 0

)
dtdz

GJ

)
dzdt
+
,
z
z
z
z
z

0
0
t
t
1
1
t
z

(7.26)

On performing the integration by parts of the first and second terms with respect to t and z,
respectively; equation (7.26) gives
0

t2

J
z ( z ) dz

t2

t2

[ J

z ( z )]dtdz t
0 t

t1

GJ z , z ( z ) dt
0

(7.27)
+

t2

t1

t2

GJ z , zz ( z ) dzdt +
t1

T ( z , t ) z + T0 * ( z z0 ) z dzdt = 0

The first term is zero since the variation of the angular displacement is not defined at two time
instants t1 and t2. Equation (7.27) can be written in a compact form as

{ J

( ) GJ
t2

t1

Since z , z 0 , and z

z , zz

L
T ( z , t ) T0 (t ) * ( z z0 )} z dzdt GJ z , z ( z ) 0 = 0

(7.28)

are variations of the angular displacement and they are arbitrary, so we

have

J

z GJ z , zz T ( z , t ) T0 (t ) * ( z z0 ) = 0

(7.29)

and
L

GJ z , z ( z ) 0 = 0

(7.30)

Equations (7.29) and (7.30) represent equations of motion and boundary conditions, respectively, for
torsional vibrations of continuous shaft. It should be noted that GJ z , z represents the internal reaction
toque, TR; and z is the angular displacement of the shaft about its longitudinal axis. Equation (7.29)
is similar to the wave equation when external torques are absent and can be solve exactly by using the
standard separation of variables method (Kreyszig, 2006) and it is briefly described for completeness
as follows.

371

Natural frequencies and mode shapes for a cantilever rod: Let us obtain torsional natural frequencies
(or eigen values) and mode shapes (or eigen functions) for a continuous shaft with cantilever
boundary conditions as shown in Figure 7.6.

Fig 7.6 A shaft with cantilever boundary conditions

The torsional rigidity of a circular shaft is GJ and the mass density of the shaft material is . At fixed
end there is no angular displacement of the shaft and at the free end the reaction toque is zero. Hence,
boundary conditions for the present case are (noting equation (7.30))

z ( z, t ) z =0 = 0

(7.31)

and

GJ

z ( z , t )
z

=0

(7.32)

z=L

Let us assume the solution of governing equation (7.33) in the following form

z ( z, t ) = ( z ) (t )

(7.34)

where ( z ) is the eigen function, and for the free vibration the function (t ) is a harmonic function
with the frequency equal to the natural frequency, nf , of the system and has the following form

(t ) = A cos nf t + B sin nf t

(7.35)

In view of equation (7.35), on taking double derivatives with respect to time, t, of equation (7.34), we
get

z ( z , t ) = nf2 ( z ) (t ) = nf2 z ( z , t )

(7.36)

372

where the dot, (
), represents the derivative with respect to the time. Now, on taking double
derivatives with respect to z of equation (7.34), we get

z( z, t ) = ( z ) (t )

(7.37)

where the prime, ( ), represents the derivative with respect to special coordinate, z. On substituting
equations (7.36) and (7.37) into equation of motion (7.29), and for free vibrations taking external
torque terms equal to zero, we get

2
nf

( z ) + G ( z )} (t ) = 0

(7.38)

In equation (7.38), since in general the function (t ) cannot be zero, hence we have

d 2 (z)
= nf2 ( z )
2
dz

(7.39)

The above equation has to be satisfied throughout the domain 0< z <L, and at boundaries we need to
satisfy equations (7.31) and (7.32). Equation (7.39) can be written as

d 2 ( z)
+ 2 ( z) = 0
2
dz

(7.40)

with

2=

nf2

(7.41)

which has the solution of the following form

( z ) = C cos z + D sin z

(7.42)

On taking the first derivative with respect to z of equation (7.42), we get

( z ) = C sin z + D cos z
(7.43)
On applying boundary conditions (7.31) and (7.32), respectively, in equations (7.42) and (7.43), we
get

(0) (t ) = 0 (0) = 0 C = 0

(7.44)

( L) (t ) = 0 ( L) = 0

(7.45)

and
D cos L = 0

373

where D and cannot be zero in a general case, otherwise whole response will be zero. Hence, the
frequency equation is obtained as
cos L = 0 ,

(7.46)

with the solution as

i L =

i
2

i = 1,3,5,...

(7.47)

On substituting equation (7.41) into equation (7.47), we get the natural frequency as

nf =
i

i G
,
2L

i = 1,3,5,...

(7.48)

Figure 7.7 Torsional mode shapes for the first three modes of a cantilever rod

On substituting equations (7.47) and (7.44) into equation (7.42), the mode shape is given by

( z ) = D sin

i
2L

z,

i = 1,3,5,...

(7.49)

The first three modes are plotted in Figure 7.7. It could be observed that the second and third modes
have one and two nodes (a zero angular displacement), respectively. Hence, on substituting equations
(7.35), (7.48) and (7.49) into equation (7.34), the torsional free vibration response of the cantilever
rod is obtained as

374

z ( z, t ) =

sin

i =1,3 ,...

i z
i G
i G
t + Bi sin
t
Ai cos
2L
2L
2 L

(7.50)

where Ai and Bi are constants, which are determined by initial conditions of the problem. It should be
noted that the actual response is summation of all the modes of the system and contributions of each
mode (i.e., constants Ai and Bi) will depend upon the initial conditions. For example for zero initial
conditions, i.e., z ( z , 0) = 0 ( z ) and
z ( z , 0) =
0 ( z ) , from equation (7.50) we get

Ai =

2 L
i z
dz
0 ( z )sin

L 0
2L

(7.51)

Bi =

2 L
i z
dz

0 ( z )sin

0
L
2L

(7.52)

and

which are obtained by the use of the orthogonality of mode shapes (Thomson and Dahleh, 1998).

A ( z) ( z)dz = 0
i

for i = j

(7.53)

for i j

(7.54)

and
L

A ( z) ( z)dz 0
i

It can be checked that above orthogonality conditions is valid for mode shapes of the present problem,
i.e., i ( z ) = sin

i
2L

z and j ( z ) = sin

j
2L

z . Table 7.1 summarizes natural frequencies and mode

shapes for some common type of boundary conditions.

Table 7.1 Natural frequency and mode shapes for torsional vibrations of rods
S.N.

Boundary conditions

Natural frequency

Mode shape

(rad/s)
1

Fixed-free

Fixed-fixed

Free-free

i G
, i = 1,3,5,...
2L

sin

i
L

, i = 1, 2,3,...

sin

i
L

, i = 0,1, 2,...

cos

i
2L

i
L

i
L

375

Axial Vibrations: The detailed treatment on axial (longitudinal) vibrations would not be done here.
Since there is a one-to-one analogy could be made between the torsional and axial vibrations by
replacing z with u z , J with A, and G with E. Here u z is the axial displacement, A is the area of the
shaft cross-section, and E is the Youngs modulus of the shaft material. It should be noted that for
axial vibrations the strain and kinetic energy can be written as

U=

1
2

EAu

2
z,z

dz

and

T=

1
2

Au
dz
2
z

(7.55)

Hence, the equations of motion for free vibrations would be

2uz E 2uz
=
t 2
z 2
(7.56)

Hence all the analysis described in this section for free vibrations is equally valid for the longitudinal
vibrations as well and now the natural frequency would be called as the longitudinal natural
frequency. Now through a simple example, a torsional free vibration analysis would be done for a
special case when the boundary condition of the problem depends upon the eigen value (i.e., the
natural frequency).

Example 7.3 Consider a circular shaft (Fig. 7.8) with rigidly fixed at one end, and carries a rigid disc
of the polar mass moment of inetia, IP, at the free end. The torsional rigidity of the shaft is GJ and the
mass density of the shaft material is . Obtain torsional natural frequencies and mode shapes.

Fig 7.8 (a) A cantilever shaft with a disc at free end (b) a free body diagram of the disc
Solution: The equation of motion of the continuous shaft remains the same as equation (7.29).
However, now boundary conditions are (refer Fig. 7.8)

376

z ( z , t ) z =0 = 0

(a)

and

GJ

z ( z, t )
z

= Ip
z =L

2 z ( z , t )
t 2

GJ ( z ) z = L = nf2 I p ( z )

(b)

z=L

z=L

Hence, the reactive torque on the shaft at the free end is balanced by the inertia of the disc; this can be
obtained by using the Newtons second law of motion from the free diagram of the disc as shown in
Fig. 7.8(b). From equation (7.42), we have a general mode shape of the following form

( z ) = C cos z + D sin z

(c)

With

2=

nf2

(d)

On taking the first derivative with respect to z of equation (d), we get

( z ) = C sin z + D cos z

(e)

On applying boundary conditions of equations (a) and (b), respectively, in equations (c) and (e), we
get

(0) = 0 = C

C=0

(f)

and
GJ ( D cos L ) = nf2 I p ( D sin L )

D GJ cos L nf2 I p sin L = 0

(g)

where D cannot be zero for a general case, otherwise the whole response would be zero. Hence, from
equation (g), noting equation (d), the frequency equation is obtained as

tan L =

2
GJ
GJ 2 1
GJ nf 1 GL 1

=
=
=
nf2 I p nf2 I p nf2 I p GJ
Ip L

(h)

which is a transdental equation to be solved for L by numerical methods (roots searching methods,
Newton-Raphson method, etc.) and the first three roots are given as

1L = 0.8605 ,

2 L = 3.4256 ,

and

3 L = 6.4373 ,

(i)

377

Hence, from equation (d), the natural frequency can be expressed as

nf = i L
i

G
,
L2

i = 1, 2,3,...

(j)

where values of i L for the first three modes are given by equation (i). The mode shape is expressed
as i(z) = D sin i z.

7.2 Applications of Finite Element Methods

The method described in the previous section for the analysis of the continuous system having infinite
degrees of freedom is feasible for only simple boundary conditions. For analyses of practical rotors,
the continuous system model is discretized intoa finite-DOF as an approximation. One such
discretization method is called the finite element method. One of the popular method of developing
finite element equations from governing equations in the differential form is the weighted-residual
method (the Galerkin method as a special case). However, by using the Rayleigh-Ritz method the
functional (various energies of the system) can be used directly to derive the finite element
formulation of the system without deriving the governing equations in the differential form. For the
present case, both methods will be presented and it will be shown that they give the same finite
element formulation. However, the former has advantage that boundary conditions come out from
finite element formulations itself. The finite element method will be briefly introduced in the present
section, however, for more detailed information on finite element methods readers are referred to
excellent text by Reddy (1993), Cook et al. (2002), Huebner et al. (2001), Bathe (1989), Zienkiewicz
and Taylor (1989), and Dixit (2009).

Figure 7.9 shows a cantilever rod as a whole system (Fig. 7.9a) and when it is divided into to several
parts (e.g., six in numbers in Fig. 7.9b) that are called elements. The total length of the rod is L and the
element length is l. For each element, for the present case, two nodes are present and a free body
diagram of one such element under torsional loadings is shown in Fig. 7.9c. An axis system is
attached to the element and it is called the local coordinate system. Apart from this we can have
global coordinate system for the whole rod. The element has two node numbers i and j. The direction
of nodal angular displacements ( zi (t ) and z j (t ) ) and nodal reaction torques ( TRi (t ) and TR j (t ) ) at
both ends are taken counter clockwise as positive (while looking from the free end along the polar
axis, i.e. the positive z-axis direction). Hence, for the element shown in Fig. 7.9c the DOF is two, i.e.,

z (t ) and z (t ) .
i

378

Fig. 7.9 A cantilver rod (a) actual system (b) a discretization of the rod into several elements
(c) a free body diagram of an element

In the finite element method, we express an approximate solution of field variable, z ( e ) ( z , t ) , within
an element in terms of nodal variables as a polynomial, and is in general defined as

z1 (t )

( ne )
z2 (t )
(e)
2
z ( z , t ) = a + bz + cz +  = N1 ( z ) N 2 ( z )  N r ( z )
= N ( z ) { z (t )}

(t )
zr

(7.57)

where N i ( z ) , i = 1, 2, , r, are called the shape function (or the approximating function or the
interpolation function) and they are a function of spatial coordinates, z, only; r is the number of DOFs
on the element, zi (t ) is the field variable (i.e., the angular displacement) value at ith node.
Superscripts: (e) and (ne) represent the element and the node of the element. The row and column
vectors are represented by

and

{ } , respectively. Equation (7.57) represents a general form the

approximate solution of equations of motion (7.30) within an element.

7.2.1 Galerkin method

On substituting the assumed approximate solution from equation (7.57) into the equation of motion
(7.29), it will give some residue, R(e), and could be expressed as

R ( e ) ( z , t ) = J

z( e ) ( z , t ) GJ z( e, zz) ( z, t ) T ( z , t ) T0 (t ) * ( z z0 )

(7.58)

379

Here it is assumed that z0 is the location of the concentrated external torque on an element in local
coordinate system also for simplicity. On minimising the residue over the element by using the
Galerkin method, we have

N R

(e)

dz = 0;

i = 1, 2,  , r

(7.59)

where l is the element length and Ni(z) is the shape function. On substituting equation (7.58) into
equation (7.59), we get

( e)
(e)
*
J

z GJ z , zz T ( z , t ) T0 (t ) ( z z0 ) dz = 0

i = 1, 2,  , r

(7.60)

At this stage the choice of the interpolation could be done, however, since the highest derivative of
field variable with respect to the spatial variable is two, we have to choose an interpolation function of
at least quadratic form so as to satisfy the completeness condition (i.e., the field variable value should
not vanish). The completeness requirement can be weaken by performing the integration by parts with
respect to z of the second term in equation (7.60), so as to give

JN


i

(e)

dz GJN i

(e) l
z,z 0

+ GJN i , z
0

(e)
z,z

dz N i T ( z , t ) + T0 (t ) * ( z z0 ) dz = 0
0

i = 1, 2,  , r
(7.61)
While choosing the interpolation function we need to satisfy two conditions, i.e., the completeness
and compatibility requirements. The compatibility requirement is to be satisfied at the inter-element
boundaries, whereas the completeness requirement to be satisfied inside of the element. In integral
equation (7.61), the highest derivative of the field variable with respect to the spatial parameter, z, is
now reduced to one. Hence, the interpolation function completeness requirement is up to the first
derivative with respect to the spatial variable, z, i.e. of z and z , z . Moreover, the compatibility
requirement will be one order less of highest differentiation with respect to z in the integral only, i.e.
up to the field variable, z, only. Hence, the interpolation function of a linear form would be sufficient
and could be expressed as

z ( e ) ( z, t ) = a + bz

(7.62)

380

where a and b are constants to be determined by field variables specified at boundaries of the element
(i.e., at nodes i and j) as shown in Figure 7.10. The z is the local coordinate system. At any general
element (e), the following boundary conditions exist (Fig. 7.10)
at ith node:

z = 0,

at jth node:

z = l,

z( e ) (0, t ) = z (t ) ;

(7.63)

and

z( e ) (l , t ) = z (t ) z( e ) = z
j

(7.64)

Figure 7.10 A beam element with field variables at boundaries

On application of equations (7.63) and (7.64) into equation (7.62), constants of the interpolation
function can be obtained as

a = zi (t )

b=

and

z (t ) z (t )
j

(7.65)

Hence, equation (7.62) takes the following form

z( e ) ( z , t ) = z (t ) +
i

z (t ) z (t )
j

or

z( e ) ( z, t ) = 1 z (t ) + z (t )
l
l

(7.66)

In view of equation (7.66), the shape function of the form as equation (7.57) can be written as
( ne )

z( e ) ( z , t ) = N ( z ) { z (t )}

(7.67)

with

N ( z ) = Ni ( z ) N j ( z ) ;

( ne )

{ z }

zi
= ;

z j

z
Ni = 1 ;
l

Nj =

z
l

(7.68)

where Ni and Nj are interpolation functions corresponding to ith and jth nodes, respectively. These have
characteristics that they have unit value at the node they correspond and have value zero at other node
(e.g., Ni = 1 at z = 0, Ni = 0 at z = l; and Nj = 0 at z = 0, Nj = 1 at z = l).

381

Checking of Compatibility Requirements:

In Figure 7.11, the jth node is the common node of (e)th and (e+1)th elements. For (e)th element, we can
write equation (7.67), as

z zi (t )

l z j (t )

z ( e ) ( z , t ) = 1
l

(7.69)

Figure 7.11 Two neighbouring elements with field variables at boundaries

Similarly for the (e+1)th element, we can write

z z j (t )

l zk (t )

z ( e +1) ( z , t ) = 1
l

(7.70)

We will check the compatibility requirements at common node, i.e. the jth node. Hence, we have
following conditions:
For z = l in (e)th element, from equation (7.69) we get

z (t )
z ( e ) ( z , t ) z =l = 0 1
= z (t )
(t )
z
i

(7.71)

For z = 0 in (e+1)th element, from equation (7.70) we get

z (t )
z ( e +1) ( z , t ) z = 0 = 1 0
= z (t )
z (t )
j

(7.72)

Hence, by using interpolations of two different elements, nodal values of the field variable at the
common node are same, as can be seen from equations (7.71) and (7.72). This ensures the
compatibility of the field variable, z ( z , t ) between two elements. Since in the present analysis, all

382

other elements will have similar interpolation functions, hence the above compatibility will be
ensured for other common nodes as well. Now, to verify whether the present interpolation function
gives compatibility of higher order also, on taking the first derivative of equation (7.69) with respect
to z, we get

( ne )

z( e, z) ( z, t ) = N, z ( z ) { z (t )}

(7.73)

with

N , z ( z ) = N i , z ( z )
N i , z = 1/ l ;

N j , z ( z ) ;

( ne )

{ z (t )}

N j , z = 1/ l

zi (t )
=
;
z j (t )

(7.74)

(7.75)

For (e)th element, we can write equation (7.73), as

z (t )
z( e, z) ( z , t ) = (1/ l ) (1/ l )

(t )
z
i

(7.76)

Similarly for (e+1)th element, we can write

z (t )
z ,(ze +1) = (1/ l ) (1/ l )

z (t )
j

(7.77)

For z = l in (e)th element, equation (7.76) gives

z ,(ze )

x =l

z z
j

(7.78)

For z = 0 in (e+1)th element, equation (7.77) gives

z ,(ze +1)

x =0

z z
k

(7.79)

As it can be seen from equations (7.78) and (7.79), the nodal values of the first spatial-derivative of
the field variable at common node of two neighbouring elements are not same. Thus, the linear
interpolation function chosen does not give compatibility of the first spatial-derivative of the field
variable, which is also not needed for the present case.

383

Checking of Completeness Requirements:

The completeness condition ensures the convergence as the size of elements is decreased and in the
limit it should converge to the exact value for elements of the infinitesimal size. The completeness
requirement of the interpolation function, for the present problem, is that it should not vanish up to the
first spatial-derivative. Since the interpolation function is a linear, hence, it satisfies the completeness
condition.

Now, equation (7.61) can be written in a vector form, as

l

J { N ( z )}

(e)
z

dz +

GJ { N }
,z

(e)
z,z

dz = GJ { N ( z )} z( e, z) +
0

{ N ( z )} T ( z, t ) + T (t )
0

( z z0 ) dz
(7.80)

From the assumed solution of the form of equation (7.67), we have

( ne )

z( e ) ( z , t ) = N ( z ) {

z (t )}

z( e, z) ( z , t ) = N, z ( z ) { z (t )}

and

( ne )

(7.81)

On substituting above equations into equation (7.80), we get

J { N } N dz {

z ( ne ) } + GJ { N , z } N , z dz { z ( ne ) } =
0

(7.82)
l

N i
GJ z , z ( e ) +
N j
0

{N } T ( z, t ) + T (t )
0

( z z0 ) dz

From equation (7.68) it should be noted that Ni gives values of 1 and 0 at nodes ith and jth,
respectively, and Nj gives values of 0 and 1 at nodes ith and jth, respectively. On substituting these
values in the first term of the right hand side of equation (7.82), we get

J { N } N dz {

z ( ne ) } + GJ { N , z } N , z dz { z ( ne ) } =
0

0 GJ ,(ze )

z =0

+
(e)
GJ , z z = l 0

{ N } T ( z, t ) + T (t )
0

( z z0 ) dz

(7.83)

The above equation can be simplified to a standard form of equations of motion of a rod element, as

[ M ] {

z ( ne ) } + [ K ] { z ( ne ) } = {TR } + {TE }

(e)

(e)

(7.84)

384

with
l

N i2
j N i

[ M ]( e) = J { N } N dz = J N
(e)

[K ]

Ni N j
Jl 2 1
dz =

2
N j
6 1 2

l
l
N i2, z
= GJ { N , z } N , z dz = GJ
N j , z N i , z
0
0

(e)

TR1 GJ z , z z = 0
T
=
=
{ R}
;
(e)
TR2 GJ z , z z =l

Ni,z N j, z
GJ
dz =
2
N j , z
l

1 1
1 1

(7.85)

(7.86)

{TE } = 0 { N } T ( z, t ) + T0 (t ) * ( z z0 ) dz

(7.87)

where [M](e) is the element mass matrix, [K](e) is the element stiffness matrix, {TR} is the element
reaction torque vector (since from the strength of materials for the pure torsion theory we have T/GJ
= z / l , hence T = GJ z( e, z) ) and {TE} is the element external torque vector. Matrices [M](e) and [K](e)
are called the consistent mass and stiffness matrices. The mass (or the mass moment of inertia) of an
element is some times lumped at its nodes based on some physical reasoning and this simplifies the
analysis drastically since the mass matrix becomes diagonal in nature. Such mass matrix is called the
lumped mass matrix. If there is a rigid disc in the rotor system then generally we put a node at the disc
and the mass (or the mass moment of inertia) of the disc is considered to be lumped at that node and
the lumped mass appears in the diagonal of the mass matrix corresponding to that node (i.e., the nodal
variable) at which it is lumped. The element external torque vector is equivalent torques acting at
nodes corresponding to the actual loading inside the element. For the present case two types of
loadings, i.e. the distributed and concentrated toques have been considered. The distributed loading
can have different forms (e.g., the linear, quadric, cubic, etc.) and may act inside the element. The
concentrated torque may act at nodes directly (more often a node is chosen in such locations) or if it is
acting inside the element then its equivalence has to be obtained at element nodes. Now, to have a
basic understanding of the Rayleigh-Ritz method the above finite element formulation is repeated in
the following section.

7.2.2 Rayleigh-Ritz Method

The Rayleight-Ritz method can also be used to develop the finite element formulation from the
functional (conservative and non-conservative energies of the system, i.e. from equation (7.25)
without taking variational operator). On substituting approximate solution of equation (7.57) into the
functional, we get the functional for the element, as

385

(e) =

t2

t1

1
2

J
z( e )2 12 GJ z( ,ez)2 T z( e ) T0 * ( z z0 ) z( e ) dzdt

(7.88)

The minimisation of the functional for an element, (e ) is equivalent to

( e )
=0
zi

i = 1,2,  , r

(7.89)

Hence, for a typical ith node, on substituting equation (7.88) into equation (7.89), we get
( e )
=0=
zi

t2

z( e, z)

z( e )
z( e )
( e )
(e)
*

GJ

T
+
T
(
z

z
)

dzdt
( 0
z
z,z
0 )
0
zi
zi
zi

t1

(7.90)
i = 1, 2, , r

On performing the integration by parts with respect to time of the first term in the left hand side of
equation (7.90), we get
0

(e)
J
z( e ) z
0
zi

dz

t1
t2

t2

(e)
(e)

z( e )
( e ) z
( e ) z , z
*

J
+
GJ
+
T
+
T
(
z

z
)

dzdt = 0
( 0
0 )
z
z,z
0
zi
zi
zi

t1

(7.91)

i = 1, 2, , r

The first term is zero for all boundary conditions. Now, by noting equation (7.57), we have following
derivatives

z( e )

( ne )
=
=
N1 z(1ne ) + N 2 z(2ne ) +  + N i z(ine ) +  + N r z(rne ) = Ni
N ( z ) { z (t )}
zi zi
zi

z( ,ez) =

z( ,ez)
zi

z( e )

( ne )
( ne )
=
=
N1, z z(1ne ) (t ) + N 2, z z(2ne ) (t ) +  + N r , z z(rne ) (t ) = N , z ( z ) { z (t )}
N ( z ) { z (t )}
z
z
z

zi

( N

( ne )

,z

( z ) { z (t )}

) = ( N

1, z

z( ne ) + N 2, z z( ne ) +  + Ni , z z( ne ) +  + N r , z z( ne ) ) = Ni , z
1

zi

and

z( e )

( ne )
( ne )
=
=
N1
z(1ne ) (t ) + N 2
z(2ne ) (t ) +  + N r
z(rne ) (t ) = N ( z ) {
z (t )} (7.92)
N ( z ) { z (t )}
t
t
t

386

On substituting equations (7.57) and (7.92) into equation (7.91), we get

t2

t1

JN N {

}

( ne )

+ GJ

dN i
( ne )

N , z { z } + N i (T + T0 * ( z z0 ) ) dzdt = 0 ;
dz

i = 1,2,  , r

Now above equations can be combined as

t2

l J { N } N dz {

}( ne ) +
z

0

t1

GJ { N , z } N , z dz { z }

( ne )

(T + T
0

( z z0 ) ){ N } dz dt = 0

which gives the standard finite element formulation of equations of motion, as

(e)
(e)
[ M ] {

z }( ne) + [ K ] { z }( ne) = {TE }( ne)

(7.93)

where [M](e) is the element mass matrix, [K](e) is the element stiffness matrix, and {TE} is the element
external torque vector as defined earlier. Matrices [M](e) and [K](e) here also are the consistent mass
and stiffness matrices, respectively. It should be noted that the element reaction torque vector, {TR},
does not appear automatically from the formulation, and it has to be added separately.

This element equation is of the general form for the torsional vibration of rod. It could be used to
obtain the element equation for all other elements depending upon the geometrical and physical
parameters of elements. However, the question now remains as, How to bring these equations
together?, so that it represents the system as a whole.

7.2.3 Assembled System Equations

Once the elemental equation of the form as equation (7.84) has been derived for all elements, the next
step would be to obtain the system equation by assembling all elements of the system. This is
illustrated for a simple cantilever beam, which for simplicity divided into only two equal length
elements as shown in Figure 7.12. By noting equations (7.84) to (7.87) and for the first element we
have i = 1 and j = 2; the finite element equation for the first element can be written as

(1)
Jl 2 1

z GJ 1 1 z TR TE
+
=
+
6 1 2

z l 1 1 z TR TE(1)
1

(7.94)

387

where TR1 and TR2 are reaction torques at node 1 and at node 2, respectively. The negative sign
represents that torque is in negative direction relative to the positive sign convention. TE(11) and TE(12 )
represent equivalent external toques, respectively, at nodes 1 and 2 of element (1) due to actual
external loadings on element (1). When a concentrated load (torque) is acting at one of the node, say 2
which is a common node for elements (1) and (2), then that torque can be considered to be acted at
node 2 of element (1) or of element (2).

Figure 7.12 (a) A cantilever beam discretised into two elements (b) element 1 (c) element 2

Similarly, for the second element we have i = 2 and j = 3, hence in view of equation (7.84) the finite
element equation for the second element can be written as

Jl 2 1

z GJ

+
6 1 2

z l
2

(2)
1 1 z2 GJ z , z 2 TE2
+
(2)
1 1

z3 GJ z , z 3 TE3

(7.95)

Since, the beam is discretised into elements of equal lengths, the stiffness and mass matrices for
elements 1 and 2 are same. However, the element can be of different lengths, cross-section and
materials; and accordingly l, J and and G would be changed in the elemental equation. Now
equations (7.94) and (7.95) can be assembled in the following form

1
0

z1
2
Jl


+ GJ
1
2
+
2
1
z2 l
6
0
1
2

z
3

GJ z , z 1
1 1
z1

1 1 + 1 1 = GJ
z , z 2 GJ z , z

z2

0 1 1
GJ z , z 3
z3

T (1)
E1

(1)
(2)
+
T
+
T

E
E
2
2
2

(2)
TE3

(7.96)

388

For the clarity of assembly procedure the matrix elements are kept in expanded form. For example, in
the assembled mass and stiffness matrices the contribution from corresponding matrices of the first
element is in the first two rows and columns, whereas, from the second element it is in the second and
third rows and columns. Similarly, for torque vectors, the contribution from the first element is in the
first two rows of assembled vectors, and the contribution from the second element is in last two rows.
Now it is clear that if we divide the rod in three elements or more how the form of the assembled
governing equation would be. For details readers are encouraged to refer to basic books on finite
element methods. Till now we have not considered the boundary condition of the problem, e.g. for the
present case boundary conditions of a cantilever rod.

7.2.4 Application of Boundary Conditions

For the present illustration (Fig. 7.12), since node 1 is fixed and node 3 is free, hence following
boundary conditions are specified by the problem

z =

z = 0
1

TR3 = 0

and

GJ z , z 3 = 0

(7.97)

On applying above boundary conditions into equation (7.96), we get

1
0 0
2
Jl


+ GJ
1
2
+
2
1
z2 l
6
0
1
2

z
3

(1)

1 1 0 0 TR1 TE1

1 1 + 1 1 = 0 + T (1) + T (2)
E2
E2

z2
0 1 1 0 T (2)
E3

z3

(7.98)

Equation (7.98) is basically three equations and unknowns are angular twists z2 and z3 (and its
second derivatives with respect to time) and the reaction torque at node 1, TR1 . The first equation
contains the reaction force, TR1 , along with the angular twist z2 (and its time derivative) as
unknowns. Whereas, last two equations contain angular twists z2 and z3 (and its time derivatives)
as unknown, however, with no reaction torque as unknown; and these can be solved if external
torques are specified. Once angular twists (and its time derivatives) are obtained from last two
equations, it can be used in the first equation to obtain the unknown reaction toque, TR1 , also. Hence,
on taking last two equations of (7.98), we get

(1)
(2)
Jl 4 1

z GJ 2 1 z TE + TE
+
=

6 1 2

z l 1 1 z TE(2)
2

(7.99)

389

Equation (7.99) can be solved if external torques are specified, this analysis is known as forced
vibrations. If external torque is zero, the problem is called the free vibration.

7.2.5 Free Torsional Vibrations

For free vibrations, the solution of equation (7.99) can be assumed of the following form

z2 (t ) z2 jnf t

=
e
z3 (t ) z3

z2
2
= nf

z3

so that

z2 jnf t

e
z3

(7.100)

where is the amplitude of torsional vibrations, nf is the torsional natural frequency, and j = 1 .
On substituting equation (7.100) into homogeneous part of equation (7.99), we get

2 Jl 4 1 GJ 2 1 z2 0
nf
=

6 1 2
l 1 1 z3 0

(7.101)

Equation (7.101) is a standard eigen value problem and for a non-trivial solution the following
determinant should be zero, i.e.

2 Jl GJ
Jl
2GJ
nf2
nf2

3 l
6
l
=0
GJ
GJ
2 Jl
2 Jl
nf
nf

6 l
3
l

(7.102)

On taking the determinate of equation (7.102), we get a polynomial in terms of nf. It is called the
characteristics equation or the frequency equation and has the following form

nf4

60G 2 36G 2
nf + 2 4 = 0
7 l 2
7 l

(7.103)

Equation (7.103) is a quadratic polynomial of the square of nf, and can be solved as

30G
30G
36G 2
=

2
2 4
7l 2
7l 7 l
2
nf

or

nf2 =

6G 5.367G 0.6492G

=
l 2
l 2
l 2

Hence, the torsional natural frequencies are (on substituting l = 0.5 L)

and

7.923G
l 2

390

nf = 0.8057
1

G
G
= 1.6114
2
l
L2

and

nf = 2.8146
2

G
G
= 5.6292
2
l
L2

(7.104)

which are more than that obtained by the analytical method with nf1 = 1.57 G / ( L2 ) and

nf = 4.71 G / ( L2 ) . The finite element always gives upper bounds of natural frequencies since it
2

over-predicts the stiffness. With increased number of element, it is expected that it should gradually
converged to the analytically calculated natural frequencies. For a single element on equating the first
diagonal element of matrix in equation (7.102) to zero, we get nf1 = 3.464 G / ( L2 ) , which is still
higher as given by equation (7.104) for two elements. The procedure of the finite element method has

been introduced in the present section with the help of a simple cantilever rod. Natural frequencies
and mode shapes could also be obtained by using the eigen value. Equation (7.101) has the following
standard form of the eigen value problem

([ A] [ I ]){} = {0}

(7.105)

with

GJ
[ A] = [ K ] [ M ] =
l
1

2 1 Jl 4 1 l 2
1 1 6 1 2 = 6G

1 1 4 1 L2 5 3
1 2 1 2 = 24G 6 5

and

= 1 / nf2

Eigen values and eigen vectors of the matrix, [A], is given by

2
1 1 / nf1 9.2426 / 24 L2
{} = = 2 =

2 1 / nf2 0.7574 / 24 G

nf 1.6114
1=

nf2 5.6292

G
L2

and
0.5774

[ ] = 0.8165

0.5774

0.8165 2

which gives the same natural frequencies as obtained in Eq. (7.104) along with the mode shapes in the
form of eigen vectors that could be normalised according to the convenience/requirement.

391

Now through more examples the application of FEM will be illustrated.

Example 7.4 Determine natural frequencies and mode shapes for a rotor system as shown in Figure
7.13. Neglect the mass of the shaft and assume that discs as lumped masses.

Figure 7.13 A two-disc rotor system

Solution: Since the present problem contains only lumped masses, let us discretise the shaft into only
one element as shown in Figure 7.14. Noting equation (7.95) the finite element equation of the oneelement rotor system can be written as

I p1

0

z1 GJ
+
I p2

z2 L

1 1 z1 GJ z , z 1

1 1 =

z2 GJ z , z 2

(a)

It should be noted that mass of the shaft is neglected hence consistent mass matrix is not present, and
only lumped masses as diagonal elements are present.

Figure 7.14 A two-disc rotor system discretised as a single element

Since for the present problem both nodes are free, hence, the following boundary conditions will
apply

T1 = T2 = 0

GJ z 1 = GJ 2 = 0

(b)

392

For the simple harmonic motion, we have the following relation

z1
2
= nf


z2

z1

z2

(c)

On substituting equations (b) and (c) into equation (a), we get

GJ
GJ

L nf I p1
0
L

z1 =

0
GJ
GJ

nf2 I p2 z2

L
L

(d)

For the non-trivial solution of equation (d), the following determinate should be zero
GJ
GJ

nf2 I p1

L
L

=0
GJ
GJ

nf I p2

L
L

GJ

I p1 + I p2 = 0
nf2 nf2 I p1 I p2
L

(e)

which gives a frequency equation. It gives torsional natural frequencies as

nf = 0

and

nf =
2

GJ I p1 + I p2
L
I p1 I p2

(f)

which are exactly the same as obtained by the closed form solution with tosional stiffness of the shaft
kt = GJ/L.
On substituting nf1 = 0 into equation (d), we get

z = z
1

Hence, the zero natural frequency is corresponding to the rigid body mode of the shaft, in this case
both rigid discs will have same angular displacements and shaft will not experience any twisting
moment which is of little practical importance. On substituting nf 2 into equation (d), we get

GJ GJ I p + I p

GJ
1
2

I p1 z1
z = 0
L

L
I p1 I p2
L 2

which gives

z
Ip
=
Ip
z
1

393

Hence, mode shapes are also exactly same as solution obtained by analytical method in the closed
form.

Example 7.5 A cantilever continuous rod has a rigid disc at the free end and also it is supported by a
torsional spring of stiffness kt at the free end as shown in Fig. 7.15. Use the following parameters: the
polar mass moment of inertia of the disc Ip = 0.02 kg-m2, torsional stiffness of the spring at the free
end kt = 100 Nm/rad, the length of rod L = 0.4 m, the diameter of the rod d = 0.015 m, the mass
density of the rod material = 7800 kg/m3, and the modulus of rigidity of the rod material G =
0.81011 N/m2. Obtain first two torsional natural frequencies of the system.

Fig. 7.15 A cantilever rod with a rigid disc and a toriosnal spring at the free end.

Solution: We have the following data of the rotor system when the shaft is discretised into three
elements as shown in Fig 7.16. The shaft is divided into two equal length elements.
Ip = 0.02 kg-m2,

kt = 100 N-m/rad,

J = 32
d 4 = 4.9710-9 m4,

GJ/l =1988.04 N-m,

L = 0.4 m,

d = 0.015 m,

G = 0.81011 N/m2,

l = L/2 = 0.4/2 = 0.2 m

Jl/6 = 1.29210-6 kg-m2

Fig. 7.16 Discretisation of the rotor system into three elements

The governing equation of 1st element can be written as

2 1

z1
1 1 z1 TR1
1.292 106
+ 1988.04

=
1 2

z2
1 1 z2 TR2

(a)

394

Let us first consider that there is no disc and no spring attached to the shaft at node (3) then the 2nd
element equation becomes

2 1

z2
1 1 z2 TR2
1.292 106
+ 1988.04

=
1 2

z3
1 1 z3 TR3

(b)

However, 3rd element contains a rigid disc at node 3, and it will contribute towards only the inertia

I p

z3 . Also a torsional spring is attached between nodes 3 and 4, hence we can write

TR3 kt z3 z4 = I p

z3

and

TR4 + kt z4 z3 = 0

(c)

which can be alternatively written in the following matrix form

1 0

z3
1 1 z3 TR3
0.02
+
100

1 1 = T
0 0

z4

z4 R4

(d)

It should be noted that equations (c) and (d) are one and same, only the form of the presentation is
different. Now, on assembling equations (a), (b) and (d), we get

z TR

1
0
1
1
0
z1

1 1

1.29 106 1 (2 + 2)
1
1


z2 + 1988.04 1 (1 + 1)
z2 = 0

0.02

z3
100 z3 0

0
0

1
2
1
1
+

1.29 106

1988.04

(e)
On simplification, we get

0

z1
2 1
1 1 0 z1 TR1

1.29 106 1 4
1

z2 + 1988.04 1 2
1 z2 = 0
0 1 15480

z
0 1 1.05 z 0
3
3

(f)

On application of boundary condition at node 1, that is z1 =

z1 = 0 , the above equation will reduce
to (on considering last two equations, i.e., after eliminating the first row and the first column)

1

z2
4
2 1 z2 0
1.29 106
+
1988.04

1 1.05 = 0
1 15480

z3

z3

(g)

395

This is a homogeneous equation. For a simple harmonic motion, we can write above equation as

3976.08 5.168 106 nf2

6 2
1988.04 1.292 10 nf

1988.04 1.292 106 nf2

2087.44 0.02nf2

(h)

For a non-trivial solution, the determinant of the above matrix should be equal to zero, i.e.

(3976.08 5.168 10

nf2

)( 2087.44 0.02 ) ( 1988.04 1.292 10

2
nf

nf2 ) = 0

(i)

which gives a polynomial of the following form

nf4 7.695 108 nf2 + 4.2047 1013 = 0

(j)

Hence, natural frequencies become

nf = 233.76 rad/s (37.20 Hz);

1

nf = 27730 rad/s (4414.7 Hz).

2

(k)

For comparison, by considering the shaft as massless and with no spring at free end, the natural
frequency is given as (for this simple case only one natural frequency exist)

nf1 =

GJ
0.8 1011 4.97 10 9
=
= 222.93 rad/s = 35.48 Hz
LI d
0.4 0.02

(l)

which is far less than the fundamental natural frequency obtained by FEM (i.e., 37.20 Hz) for the
actual rotor system. It is expected as (i) the additional spring at free end provides more stiffness, so
the natural frequency is expected to be more in actual case as compared simplified rotor system with
massless shaft and without spring (ii) FEM always give upper bound and for the present analysis
number of elements (only two) are too less. However, the polar mass moment of inertia of the shaft
has very little effect as it can be seen below.

The polar mass moment of inertia of the shaft has effect of decreasing the natural frequency, when it
is considered in the analysis. For a continuous cantilever shaft with a disc at free end, the closed form
solution is given as (see Example 7.3)

396

nf1 = 0.8605

GJ
0.8 1011 4.97 10 9
= 0.8605
= 191.83 rad/s =30.53 Hz
LI d
0.4 0.02

(m)

and

nf2 = 3.4256

GJ
0.8 1011 4.97 10 9
= 3.4256
= 763.68 rad/s =121.54 Hz
LI d
0.4 0.02

(n)

On comparing equation (l) and (m), it can be seen that while considering the shaft inertia there is a
decrease in the first natural frequency. On comparing equations (k) and (n), it can be observed that the
second mode natural frequency has large difference for the cases of with and without spring at free
end. These comparisons have been made to highlight effectiveness of the FEM and possible trend of
natural frequencies due to various complexities in the system (i.e., additional spring, additional rigid
disc, addition of shaft inertia, etc.). It should be noted that in the finite element analysis, only two
elements were considered, if we increase the number of elements natural frequencies of the actual
system is expected to be refined.

Example 7.6 A marine reciprocating engine, flywheel and propeller are approximately equivalent to
the following three-rotor system. The engine has a crack 50 cm long and a connecting rod 250 cm
long. The engine revolving parts are equivalent to 50 kg at crank radius, and the piston and pin masses
are 41 kg. The connecting rod mass is 52 kg and its center of gravity is 26 cm from the crankpin
center. The mass of the flywheel is 200 kg with the radius of gyration of 25 cm. The propeller has the
polar mass moment of inertia of 6 kg-m2. The equivalent shaft between the engine masses and the
flywheel is 38 cm diameter and 5.3 m long and that between the flywheel and the propeller is 36 cm
diameter and 11.5 m long. Find the torsional natural frequencies and mode shapes of the system by
the finite element method.

Solution: Using the TMM the same problem was solved in the previous chapter (see example 6.11).
Hence, from the previous chapter; we have the following data of the three-disc rotor system
corresponding to the reciprocating engine, flywheel, and propeller as shown in Fig. 7.17.

Figure 7.17 A three-disc rotor model

397

The equivalent polar moment of inertia of the engine is: I pe = 29.95 kg-m2. The polar mass moment
of inertia of the flywheel is: I p f = 12.5 kg-m2. For the propeller, the polar mass moment of inertia is
given as I p p = 6 kg-m2. The torsional stiffness of shaft segments (1) and (2) are given as: kt1 =3107
N-m/rad, and kt1 =8.67107 N-m/rad.

Elements (1) and (2) are shown in Figure 7.18.

(a)

(b)
Figure 7.18 (a) 1st element and (b) 2nd element

The elemental equation for element (1) is

1 z1 T1
29.95 0

z1
7 1
+ 3 10
0

=
0

z2

1 1 z2 T2

(a)

Similarly, on considering element (2) consists of two discs, one at either end of the shaft, we have

1 z2 T2
12.5 0

z2
7 1
0 6

+ 8.67 10 1 1 = T

z3

z3 3

(b)

On assembling equations (a) and (b), we get

0

z1
3
0 z1 T1
29.95 0
3

0
12.5 0

z2 + 107 3 3 + 8.67 8.67 z2 = 0

0
0
0
6

z
8.67
8.67 z T3
3
3

(c)

For the present case boundary conditions are: T1 = T3 = 0 , and for the simple harmonic motion of free
vibration with frequency nf , we have
0 z1
3
0 z1 0
29.95 0
3

nf2 0
12.5 0 z2 + 107 3 11.67 8.67 z2 = 0
0
0 8.67 8.67 z 0
0
6 z
3
3

(d)

398

or
3 107 29.95nf2

3 107

3 107
11.67 107 12.5nf2
8.67 107

z1 0
0

8.67 107 z2 = 0
8.67 107 6nf2 z 0
3

(e)

On taking determinant of the above matrix to zero, we get the natural frequency of the system as

nf = 1.58 103 rad/s

nf = 0

nf = 4.72 103 rad/s

(f)

in which one of the natural frequency is zero because of the free-free end conditions. The same can
also be obtained by using eigen value analysis. For this let us pre-multiplying equation (d) by inverse
of the mass matrix [M] to get

([ D] [ I ]){ } = {0}
2
nf

(g)

with
1.0 106

[ D] = [ M ]1[ K ] = 2.4 106

1.0 106
9.34 106
14.45 106

0
6
6.94 10
14.45 106

(h)

The eigen value of matrix [D] would be square of the natural frequencies and, columns of eigen
matrix are the corresponding mode shapes. These are given as

nf2 1
0

2
6
Eigen values = nf2 = 2.496 10
2 2.228 107

nf3

nf1
0

nf2 = 1579.87 rad/s

4.720.17

nf3

and
-0.577

Eigenvector matrix = -0.577

-0.577
1

0.387

-0.587
-0.711

0.022

-0.476
0.879

It should be noted that natural frequencies obtained are close to what we obtained in example 6.14.
From eigen vector matrix mode shapes could be plotted after normalisation (e.g., by assigning the first
row elements as 1 and correspondingly other two row can be scaled as follows

399

1

Normalised eigenvector matrix = 1
1
1

-1.517
-1.837

-21.636
39.954

Example 7.7 For a continuous shaft of 3 m length and 0.03 m diameter, obtain torsional natural
frequencies up to the fifth mode and plot corresponding mode shapes for the free-free boundary
conditions. Perform the converge study by increasing the number of elements (i.e., 5, 10, 20, 50, 100,
and 500) and compare natural frequencies with the closed form analytical solutions and discuss the
results. The following properties of the shaft should be taken: = 7800 kg/m3 and G = 0.81011 N/m2.

Solution: The shaft is discretised into five equal numbers of elements as shown in Figure 7.19 for
illustration; however, torsional natural frequencies are then presented with higher numbers of
elements also.

Figure 7.19 A continuous rod discretised into five elements with six nodes

The properties of the elements are:

l = L / 5 = 3 / 5 = 0.6 m, J =

Jl
6

32

(0.03)4 = 7.95 108 m4,

7800 7.95 108 0.6

= 6.20 105 kg-m2,
6

GJ 0.8 1011 7.95 108

=
= 1.06 104 N-m
l
0.6

The elemental equation of jth element for torsional vibrations is given as

1 z j T j
2 1

z j
4 1
6.20 105

+ 1.06 10
=

1 2

z j +1
1 1 z j+1 T j +1

j = 1, 2,,5

(a)

The entire element has same size and properties, for example, for third element the above equation has
the following form

400

1 z3 T3
2 1

z3
4 1
6.20 105
+
1.06

10

1 1 = T
1 2

z4

z4 4

(b)

On assembling equations of all five elements, we get

1 0 0 0 0 z1
1 1 0 0 0 0 z1 T1

1 2 1 0 0 0 z 0
4 1 0 0 0

z2

2

1 4 1 0 0 z3
0 1 2 1 0 0 z3 0
4
+ 1.06 10
=
0 1 4 1 0

z4
0 0 1 2 1 0 z4 0
0 0 0 1 2 1 0
0 0 1 4 1

z
5

z5
0 0 0 1 2


0 0 0 0 1 1 T6
z6
z6

2
1

5 0
6.20 10
0
0

(c)

which has the standard form as

I p {

z } + [ Kt ]{ z } = {TR }

(d)

Boundary conditions are free-free ends of the rod, hence, we have T1 = T6 = 0 in equation (d). Hence,
after application boundary conditions also the size of matrices remains the same, i.e. 66. For simple
harmonic oscillation with frequency nf , equation (d) takes the following standard eigen value
problem form

([ K ]
t

2
nf

I p { z } = {0}

(e)

or

([ D] [ I ]){ } = {0}
2
nf

(e)

[ D ] = I p [ K t ]

(f)

with
1

On obtaining the eigen value of the matrix [D], we get the following natural frequencies after taking
the square root of eigen values

n = 0;
f1

n = 3398; n = 7129;
f2

f3

n = 11467
f4

and n f = 16052 (all are in rad/s)

5

401

Table 7.2 Convergence of torsional natural frequencies with elements number

S.
N.

Natural
frequency

Number of rod elements

(5)

nf

(10)

(20)

(50)

Exact
value
(100)

(500)

3 409.1

3 367.5

3 357.2

3 354.3

3353.9

3 353.7

3 353.7

7 152.2

6 818.2

6 735.0

6 711.8

6708.5

6 707.4

6 707.4

11 503.0

10 436.0

10 154.0

10 076.0

10 065.0

10 061.0

10 061.0

16 114.0

14 304.0

13 636.0

13 450.0

13 424.0

13 415.0

13 415.0

18 490.0

18 490.0

17 202.0

16 838.0

16 786.0

16 769.0

16 769.0

3
4
5
6

nf

nf

nf

nf

nf

0.1
(2)

0.08

(1)

0.06

Relative amplitudes

(3)

0.04
0.02

(4)
(5)

0
-0.02
-0.04
-0.06
-0.08
-0.1

0.5

1
1.5
2
Length of the beam (m)

2.5

Figure 7.20 Mode shapes (1) first mode, (2) second mode, (3) third mode, (4) fourth mode, and
(5) fifth mode.
The corresponding mode shape can be obtained by the corresponding eigen vectors from the eigen
matrix. Figure 7.20 shows mode shapes up to the fifth mode. It should be noted that flexible modes
have nodes (e.g., 2nd mode has one node, 3rd node has two nodes, etc.). For comparison with exact
solution, we have

nf =

n
L

(g)

402

with L = 3m, G = 8.0 1011 N/m 2 , = 7850 kg/m . We have the following exact natural frequencies of
the rotor system

nf = 0; nf = 3353.7; nf = 6707.4; nf = 10061,

1

and nf = 16769.0

(all are in rad/s)

On comparison with the FEM solution, it shows that lower natural frequencies are quite close to the
exact values; however, for higher modes the natural frequencies by FEM are quite high as compared
to the exact values. For getting more accuracy for higher modes by FEM more number of elements
need to be considered. Table 7.2 shows such convergence study up to 500 elements.

7.3 Development of the Finite Element for a Gear-pair

Figure 7.21 shows a gear-pair element. The counter clockwise direction is taken as the positive
direction for all angular displacements. It is the practice and convenient in finite element
formulations. The gear ratio, n, (which is an inverse of the speed ratio) is defined as

n=

1
2

(7.106)

where 1 and 2 are the angular speed of the driver (pinion or gear 1) and driven (gear 2) shafts.

Figure 7.21 A gear-pair element

For no slip and no backlash condition, and for perfectly rigid gear teeth (it is assumed that gear tooth
are not providing any torsional flexibility; however, in the real case bending of the tooth and the
deformation at the contact zone due to the Hertzian contact stress might give rise to appreciable

403

amount of torsional flexibility, especially for the case when the power-transmission is very high) from
Figure 7.21, we have

1
n

z = z
g2

(7.107)

g1

where zg and zg are angular displacements of gears 1 and 2, respectively. Since zg is defined
1

in terms of zg , we can eliminate zg from the state vector of the system.

1

The Mass Matrix: The equivalent inertia force of the gear-pair with respect to the reference shaft 1 is
given as (refer chapter 6 for more details regarding equivalent inertia force of discs in two shaft
systems with different speeds)

I pg 2

I pg 1 + 2

z g 1
n

where I pg 1 and I pg 1 are polar mass moment of inertia of gears 1 and 2, respectively. Above expression
could be written in the mass matrix form as

I pg 2
I pg 1 + 2 0

zg1 0
n =

z2 0
0
0

(7.108)

The Stiffness Matrix: The stiffness matrix of the shaft element connected to gear 2 will have to be
modified, since the angular deflection of the left hand side of that shaft element is now equal to (

zg / n ) instead of zg 2 (see Figure 7.22).

1

Figure 7.22 An equivalent shaft element connected to gear 2

The potential energy of the shaft element in Figure 7.22 is

404

1
2

U1 = kt1 z2 z1

= kt1 g 1 + z1
n

1
2

(7.109)

The work done is give as

W1 = T1 z1 T2 z2 = T1

g1

T2 z2

(7.110)

On applying the Raleigh-Ritz method, we get

(U1 + W1 )
zg 1

T
kt1 zg 1
+ z2 + 1 = 0

n
n n

= kt1 g 1 + z2 T2 = 0
n

kt1
n

z +
g1

kt1
n

z =
2

T1
n

(7.111)

and

(U1 + W1 )
zg 2

kt1
n

z + kt z = T2
g1

(7.112)

Gear-pair Element Equations: Equations (7.111) and (7.112) can be combined in the stiffness matrix
form as

kt1 / n 2

kt1 / n

kt1 / n zg1 T1 / n
=

kt1 z2 T2

(7.113)

On combining equations (7.114) and (7.115), we get

I g2
2
I g1 + 2 0

zg1 kt1 / n
n +

z kt / n

0
0 2 1

kt1 / n zg1 T1 / n
=

kt1 z2 T2

which is the equation of motion of the gear-pair element. Following examples will illustrate the use of
the gear-pair element for the geared and branched systems by using FEM.

Example 7.8 For a geared system as shown in Figure 7.23, find torsional natural frequencies. The
shaft A has the diameter of 5 cm and the length of 0.75 m, and the shaft B has the diameter of 4 cm
and the length of 1.0 m. Take the modulus of rigidity of the shaft G equals to 0.8 1011 N/m2, the
polar mass moment of inertia of discs and gears are I pA = 24 Nm2, I pB = 10 Nm2, I pgA = 5 Nm2, and
I pgB = 3 Nm2.

405

Figure 7.23 Two-discs with a geared system

Solution: Now this problem will be solved by using the FEM for illustration of the method to geared
system. The pinion and the gear have appreciable amount of the polar mass moment of inertia. Let us
divide the geared system in two elements (Fig. 7.23). Denote the node number of the disc on branch A
as 1, the gear as 2, the gear as 3 and the node number of the disc on branch B as 4 (Fig. 7.23).

The following data could be obtained for the present rotor system:

JA =

4
d A = 6.136 107 m 4 ;
32

JB =

4
d B = 2.51 107 m 4 ;
32

and

k1 = k A =

GJ A 0.8 1011 6.136 10-7

=
= 6.545 104 Nm/rad;
lA
0.75

k2 = k B = 2.011 104 Nm/rad .

Now elemental equations will be written one by one for each element, and the same would be then assembled to
get global equation of motion of the geared system.

Fig. 7.24 Elements with nodal variables (a) 1st element (b) 2nd element

Element (1): Elemental equation can be written as (neglecting the polar mass moment of inertia of the shaft)

406

6.545 z1 T1
24 0

z1
4 6.545
0 5

+ 10 6.545 6.545 = T

z2

z2 2

(a)

Element (2): Elemental equation can be written as

2.011 z3 T3
3 0

z3
4 2.011
0 10

+ 10 2.011 2.011 = T

z4

z4 4

Now from above equation

z3 has to be replaced by z2 , since we have, z3 = z2 / n

(b)

where n is the gear

ratio (n =2). Hence, it takes the form

2
3 / 22 0

z2
2.011/ 2 z2 T3 / n
4 2.011/ 2
+
10

=
10

z4
2.011 z4 T4
0
2.011 / 2

(c)

Now equation (a) and (c) can be assembled to get the global governing equation

0
0

z1
6.545
0 z1
T1
24
6.545

0 5 + 0.75 0

+ 104 6.545 6.545 + 0.503 1.006 = T T / n

z2

z2 2 ( 3 )

0
10

z
1.006
2.011 z
T4
0
0

4
4

(d)

At junction, we have T2 (T3/n) = 0; and discs at ends are free, hence, T1 = T4 = 0. Thus, on
application of boundary conditions equation (d) takes the following form for free vibrations

2
nf

0
0
0 z1 0
24
6.545 6.545
0 5.75 0 + 10 4 6.545 7.048 1.006 = 0

z2
0
3
1.006 2.011 z 0
0
0
4

(e)

From which we can write

1

0
0
0
24
6.545 6.545
1

4
[ D ] = [ M ] [ K ] = 0 5.75 0 10 6.545 7.048 1.006
0
0
0
3
1.006 2.011

(f)

From the eigen value analysis, natural frequencies and mode shapes can be obtained from the matrix
[D] and are given as

407

nf = 0,
1

z 1 = 1 1 0.5 ;

nf = 123.66 rad/s,
3

nf = 79.39 rad/s,
2

z 2 = 1 1.31 10.9 ;

z 2 = 1 4.64 1.79 .

Example 7.9 Obtain torsional natural frequencies of a branched rotor system as shown in Figure 7.25.
Take the polar mass moment of inertia of rotors and gears as: I p A = 0.01 kg-m2, I pE = I pB = 0.005
kg-m2 and I pF = I pC = I pD = 0.006 kg-m2. Take gear ratio between various gear pairs as: nBC = 3
and nBD = 4. Shaft lengths are: lAB = lCE = lDF = 0.25 m and its diameters are dAB = 0.03 m, dCE = 0.02
m, and dDF = 0.02 m. Take the shaft modulus of rigidity G = 0.8 1011 N/m2.

Figure 7.25 A branched system with gears

Solution: The branched system is divided into three elements. Figure 7.25 shows various element
numbers of the branched rotor system. Now various element equations can be written as follows:

Element (1): Figure 7.26 shows element (1) with nodal variables. Since this shaft element is the input
shaft, hence, there is no change in the inertia terms as well as in the stiffness terms.

Figure 7.26 Element (1)

Equations of motion for element (1) can be written as (on neglecting the inertia of the shaft)

408

I pA

0

zA kt1
+
I pB

zB kt1

kt1 z A TA
=

kt1 z B TB

(a)

Element (2): Figure 7.27 shows element (2) with nodal variables. The actual nodal variable at gear C
is related with the nodal variable of gear B, hence nodal variable of gear B is used. This affects both
the mass and stiffness matrices as well as the internal torque vector.

Figure 7.27 Element (2)

Equations of motion for element (2) can be written as

2
I pC / nBC

2
0

zB kt2 / nBC
+
I pE

zE kt2 / nBC

kt2 / nBC zB TC / nBC

=

kt2 zE TE

(b)

Element (3): Figure 7.28 shows element (3) with nodal variables. Here also instead of actual angular
displacement of gear D the angular displacement of gear B is used.

Figure 7.28 Element (3)

For element (3), we can write the equation of motion as

2
I pD / nBD

2
0

zB kt3 / nBD
+
I pF

zF kt3 / nBD

kt3 / nBD zB TD / nBD

=

kt3 zF TF

(c)

409

In the present problem we have discretised the system into three elements. However, because at the
branched point angular displacements of three gears (i.e., B, C, D) are related, and this leads to a finite
element system model of only four-DOF system. Hence, on assembling equations (a-c), we get

I pA

{I + ( I
pB

pC

) (

2
2
+ I pD / nBD
/ nBC

kt1

kt
+ 1
0

)}

I pE

0


z
A
0

zB

0

zE

I pFD

zF

kt1

{k + ( k
t1

t2

) (

2
2
+ kt3 / nBD
/ nBC

)}

kt2 / nBC

kt2 / nBC

kt2

kt3 / nBD

z
A TA
kt3 / nBD zB 0
=

0 zE TE

kt3 zF TF
0

(d)

It should be noted that in the assembled form, the second row in the torque column is zero, since we
have the following condition at the branch point

TB

Tc
T
D =0
nBC nBD

(e)

From free boundary conditions at shaft ends A, E and F, we have

TA = TE = TF = 0

(f)

Hence, on application of boundary conditions in equation (d), all the terms in the right hand side
vector are zero. Hence, equation (d) forms an eigen value problem with the matrix size of 44. From
equation (d), now to have torsional natural frequencies, we have

[ K ] nf2 [ M ] = 0

(g)

where [K] and [M] are the assembled stiffness and mass matrices. Equation (g) gives frequency
equation as a polynomial and roots of the polynomial gives torsional natural frequency of the system.
However, for a matrix of size 44 it is quite cumbersome. Alternatively, equation (g) could be written
as

[ D] nf2 [ I ] = 0

with

[ D ] = [ M ]1[ K ]

(h)

410

In which eigen value of matrix [D] is related with natural frequencies (square root of the eigen value
is the natural frequency) and eigen vector is related with mode shape. Hence, this method requires
eigen value analysis and it is preferred for large size matrices. It should be noted that the stiffness
matrix, [K], is a singular matrix since we have free-free boundary conditions. For the present case, we
have following data

I pA = 0.01 kg-m2;

I pE = 0.005 kg-m2;

I pF = 0.006 kg-m2

I pB = 0.005 kg-m2;

I pC = 0.006 kg-m2;

I pD = 0.006 kg-m2

nBC = 3 ;

nBD = 4;

kt AB =

GJ AB
= 2.54 104 N-m/rad;
l AB

ktCE = ktDF = 0.50 104 N-m/rad

The assembled stiffness and mass matrices will have the following form

0
0
25.4 25.4
25.4 26.32 1.68 1.26
103 N-m/rad
[K ] =
0
1.68 5.02
0

1.26
0
5.02
0
and
0
0
0
0.01
0
0.006
0
0

[M ] =
kg-m2
0
0
0.005
0

0
0
0.006
0
Using equations (g) or (h), natural frequencies can be obtained as

nf = 0 rad/s ;

nf = 922.22 rad/s

nf = 1015.70 rad/s;

nf = 2614.50 rad/s

In the former method, relative amplitudes of angular displacements (i.e. mode shapes) are obtained by
substituting these natural frequencies one at a time in equation (d). In latter method, mode shapes
1

(eigen vectors) are obtained by the eigen value analysis of the matrix [ D ] = [ M ] [ K ] . The eigen
value will give the square of natural frequencies and eigen vectors will give mode shapes as described
above. Mode shapes are shown in Fig. 7.29 for first three modes only. It should be noted that relative
angular displacements obtained from equation (d) are actual displacements and no further scaling by

411

the gear ratio is required as we did for the equivalent geared system. The angular displacements of
discs C and D are related (i.e., these are not independent) with the angular displacement of disc B with
respective gear ratio. These angular displacements have been obtained with these relations only for
the present case for completeness of mode shapes.

Relative Linear displacements

1.5

0.5

-0.5

0.1

0.2

E
0.3

0.4

0.5

Position on discs and gears

Figure 7.29(a) Mode shape of the geared system for nf1 = 0

Relative linear displacements

12

10
8
6
4
2

0
-2
0

C,D
0.1

0.2

E
0.3

0.4

Position on discs and gears

Figure 7.29(b) Mode shape of the geared system for nf2 = 922.22

0.5

412

Relative linear displacements

2
A
0

-2

B
C,D

0.1

0.2

0.3

0.4

Position on discs and gears

Figure 7.29(c) Mode shape of the geared system for nf3 = 1015.70

Fig. 7.29(d) Mode shape of the geared system for nf 4 =2623.8

0.5

413

Concluding Remarks:
To summarise, in the present chapter torsional vibrations of rods by the continuous system approach
is initially presented. The Hamiltons principle is used to obtain the governing equations, and
associated geometric (i.e., related to the angular displacement) and natural boundary (i.e., related to
the torque) conditions. The form of the governing equation is similar to the wave equation. The closed
form solutions for simple boundary conditions are presented by method of separation of variables, and
corresponding natural frequencies (eigen frequencies) and mode shapes (eigen functions) are
tabulated and plotted, respectively. Since for more complex rotor-support systems obtaining the
closed form solution is very difficult and some times it is impractical, hence a need of an approximate
method (e.g., FEM) is felt. The finite element analysis for the torsional vibration is presented by the
Galerkin and Raleigh-Ritz methods. The analysis presented for the torsional case is also of the similar
form as that of axial vibrations of rods. In that case the axial displacement will replace the angular
displacement, and the axial force will replace the torque. The axial stiffness of the rod is given as
ka = EA / l (where A is the cross sectional area of the rod), whereas the torsional stiffness is given as kt

= GJ/l. With these analogies all the analysis of the torsional vibration can be used for axial vibrations
also. The finite element analysis is extended to rotors with branched systems also by developing
elemental equations for the gear-pair element. However, the flexibility of the gear teeth is not
considered and they are assumed to be rigid. There is another coupling effect is present between the
transverse and axial directions, when a rod has initial constant tension or compression. Due to initial
constant tension, it is expected that the effective transverse stiffness increases and due to this the
natural frequency also increases. However, for initial compression of the rod the transverse stiffness
decreases and correspondingly the natural frequency also decreases. In the limit when this
compressive load is equal to or greater than the buckling load the stiffness and so natural frequency of
the system becomes zero. This case will be considered while discussing the transverse vibrations in
subsequent chapters. Another issue we have not touched in the present chapter that is how to tackle
the increasing number of matrix size while using FEM. There are condensation schemes (or the
reduction schemes) that allow elimination of the unwanted DOFs to be eliminated from the assembled
system matrices and it will be discussed during the application of FEM in transverse vibrations.

414

Exercise Problems

Exercise 7.1 Obtain the following equation of motion by using the Newtons second law of motion of
a rod subjected to a uniform torque T ( z , t ) for the pure torsional vibration. The notations have the
following nomenclature: GJ as the modulus of rigidity, L is the length, is the mass density, and J is
the polar area moment of inertia.

J

z GJ z , zz Tz ( z , t ) = 0
e

[Hint: The free body diagram is shown in Fig. E7.1, where Tze ( z , t ) is the distributed external torque
on the rod,

Tz ( z , t ) = GJ z , z is the reactive torque at the axial location of z, and

Tz
= GJ z , zz is
z

the increment in the reactive torque from the axial location of z to (z + dz).

Tz

Tze

Tz +

Tz
z

Fig. E7.1 A free body diagram of an element with a distributed torque

With the help of above free body diagram the equation of motion could be obtained].

Exercise 7.2 A cantilever rod has a thin rigid disc (of mass M with r as the radius of gyration) at free
end and also it is supported by a torsional spring of stiffness kt at the free end (See Fig. E7.2). Let us
taje GJ as modulus of rigidity, L is the length, Ip(z) = J as the mass polar moment of inertia per unit
length, is the mass density, and J is the polar area moment of inertia. Write down equations of
motion and boundary conditions for torsional vibrations. Obtain the expressions of frequency
equations and mode shapes in the closed form.

Fig. E7.2 A cantilever rod with a rigid disc and a toriosnal spring at free end.

415

[Hint: Equations of motion remains the same as that of continuous rod, however, boundary condition
would change due to a disc and a spring at the free end].

Exercise 7.3 Consider a rod with E as the Youngs modulus, A is the area of the cross section, L is the
length, m(z) = A as the mass per unit length, and uz(z) is the axial displacement at a position z. Write
down equations of motion and boundary conditions for axial vibrations. Obtain the expressions of
frequency equations and mode shapes in closed form for the following boundary conditions (i) fixedL

1
u
free (ii) free-free, and (ii) fixed-fixed. [Answer: U =
EA z dz ; T =
20
z

1
2

Au
dz ;
2
z

Au

z EAu z , zz Fz ( z , t ) Fz 0 (t ) * ( z z0 ) = 0 ; EAu z , z ( u z ) 0 = 0 ; refer Table E7.3].

Table E7.3 Natural frequency and mode shapes for axial vibrations of rods
S.N.

Boundary conditions

Natural frequency

Mode shape

(rad/s)
1

Fixed-free

Fixed-fixed

Free-free

i E
, i = 1,3,5,...
2L

sin

i
L

, i = 1, 2,3,...

sin

i
L

, i = 0,1, 2,...

cos

i
2L

i
L

i
L

Exercise 7.4 Formulate the equation of motion and the eigenvalue problem of a non-uniform rod for
torsional vibrations, which is clamped at one end and free at the other end. Consider only the torsional
vibrations. Let GJ ( z ) = 2GJ o {1 ( z / L)} and I p ( z ) = 2 I po {1 ( z / L)} , where J o and I po are constant
quantities. Obtain expressions of natural frequencies and mode shapes. [Hint: equation of motion will
remain the same, however, now geometrical parameters of the rod are no more constants and vary
with spatial coordinate, z].

Exercise 7.5 Obtain natural frequencies and mode shapes of the rotor system shown in Figure E7.5
for the following parameters: I p1 = 0.02 kg-m2, I p2 = 0.08 kg-m2, l = 1 m, d = 0.01 m, = 7800
kg/m3 and G = 0.81011 N/m2. Use FEM (at least two elements) and compare the results by
considering with and without mass of the shaft.

416

Figure E7.5 Two disc rotor system

Exercise 7.6 All the exercise problems provided for the TMM in previous chapter can be attempted
with the FEM. Compare your results obtained by TMM and FEM and interpret physically
observations made. Whether TMM under-estimate and FEM over-estimate natural frequencies?

Exercise 7.7 For a continuous shaft of 3 m length and 0.03 m diameter, obtain torsional natural
frequencies up to fifth mode and plot corresponding mode shapes for the following boundary
conditions (i) fixed-free and (ii) fixed-fixed. Use finite element method. Compare natural frequencies
with closed form analytical solutions, and perform a converge study by increasing the number of
elements (i.e., 5, 10, 50, and 100) and discuss the result. The following properties of the shaft should
be taken: = 7800 kg/m3 and G = 0.81011 N/m2.

Exercise 7.8 For a continuous shaft of 3 m length and 0.03 m diameter, obtain torsional natural
frequencies up to fifth mode and plot corresponding mode shapes for the following boundary
conditions (i) a cantilevered shaft with a disc at the free end (ii) a fixed-fixed shaft with a disc at the
midspan,. The polar mass moment of inertia of the disc Ip = 0.02 kg-m2. Use finite element method.
Compare natural frequencies with closed form analytical solutions. The following properties of the
shaft should be taken: = 7800 kg/m3 and G = 0.81011 N/m2.

Exercise 7.9 Obtain torsional natural frequencies and mode shapes of an epi-cyclic gear train as
shown in Figure E7.9. Find also the location of nodal point on the shaft. The gear mounted on shaft
B is a planetary gear and the gear on shaft A is a sun gear. Consider the polar mass moment of
inertia of the shaft, the arm, and gears as negligible. Shaft A has 5 cm of diameter and 0.75 m of
length and shaft B has 4 cm of diameter and 1.0 m of length. Take the modulus of rigidity of the
shaft G equals to 0.8 1011 N/m2, the polar mass moment of inertia of discs are IA = 24 Nm2 and IB =
10 Nm2. Use FEM. [Hint: Obtain the gear ratio nAB between shaft A and B using epi-cyclic gear train,
now obtain the equivalent two mass rotor system for which closed form expressions of natural
frequencies are available]

417

Figure E7.9 An epi-cyclic geared system

Exercise 7.10 A rod has two rigid discs at either ends and it is supported by two torsional springs of
stiffness kt1 = kt2 = kt at each end of the rod. Use the following parameters: the mass polar moment of
inertia of the disc I p = 0.002 kg-m2, torsional stiffness of the spring at free end kt = 2 102 Nm/rad,
length of the rod l = 0.5 m, diameter of the rod d = 0.01 m, mass density of the rod material = 7800
kg/m3 and modulus of rigidity of the rod material G = 0.81011 N/m2. Obtain natural frequencies of
the system by considering rod as two elements only.

Exercise 7.11 Find torsional natural frequencies of an overhung rotor system as shown in Figure
E7.11. Consider the shaft as massless and is made of steel with the modulus of rigidity of 0.8(10)11
N/m2. A disc is mounted at the free end of the shaft with the polar mass moment of inertia 0.01 kg-m2.
In the diagram all dimensions are in cm. Use the FEM and compare results with the TMM.

Figure E7.11

Exercise 7.12 Find torsional natural frequencies and mode shapes of the rotor system shown in Figure
E7.12. B1 and B2 are frictionless bearings, which provide free-free end condition; and D1, D2, D3 and
D4 are rigid discs. The shaft is made of the steel with the modulus of rigidity G = 0.8 (10)11 N/m2 and
a uniform diameter d = 20 mm. Various shaft lengths are as follows: B1D1 = 150 mm, D1D2 = 50 mm,
D2D3 = 50 mm, D3D4 = 50 mm and D4B2 = 150 mm. The mass of discs are: m1 = 4 kg, m2 = 5 kg, m3 =

418

6 kg and m4 = 7 kg. Consider the shaft as mass-less. Consider discs as thin and take diameter of discs
as d1 = 8 cm, d 2 = 10 cm, d 3 = 12 cm, and d 4 = 14 cm.

Figure E7.12 A multi-disc rotor system

419

References:
1. Dixit, U.S., 2009, Finite Element Methods for Engineers, Cengage Learning, New Delhi.
2. Reddy, J.N., 1993, An Introduction to the Finite Element Method, McGraw-Hill, New York.
3. Huebner, K.H., Dewhirst, D.L., Smith, D.E., and Byrom, T.G., 2001, The Finite Element
Method for Engineers, Fourth Edition, Wiley-Interscience Publication, John Wiley & Sons,
Inc., New York.
4. Cook, R.D., Malkus, D.S., Plesha, M.E., and Witt, R.J., 2002, Concepts and Applications of
Finite Element Analysis, 4th edition, John Wiley & Sons, New York, 2002.
5. Bathe, K.J., 1982, Finite Element Procedures in Engineering Analysis, Prentice-Hall,
Englewood Cliffs, NJ.
6. Kreyszig, E, 2006, Advanced Engineering Mathematics, John Wiley & Sons, Hoboken.
7. Meirovitch, L., 1986, Elements of Vibration Analysis, McGraw Hill Book Co., NY.
8. Rao, J.S., 1992, Advanced Theory of Vibration, Wiley Eastern Limited, New Delhi.
9. Hughes, T.J.T., 1986, The Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ.
10. Thomson, W.T. and Dahleh, M.D., 1998, Theory of Vibration with Applications, Fifth
Edition, Pearson Education Inc., New Delhi.
11. Zienkiewicz, O.C., and Taylor, R.L., 1989, The Finite Element Method, 3rd ed. McGraw-Hill,
New York.
----------------------------&&&&&&&&&&&&&&-----------------------------