Dr. R. Tiwari ([email protected].

in) 1
By
Dr. Rajiv Tiwari
Department of Mechanical Engineering
Indian Institute of Technology Guwahati 781039
Under AICTE Sponsored QIP Short Term Course on
Theory & Practice of Rotor Dynamics
(15-19 Dec 2008)
IIT Guwahati
TORSIONAL VIBRATIONS OF ROTOR
SYSTEMS
Dr. R. Tiwari ([email protected]) 2
INTRODUCTION
Torsional vibrations is predominant whenever there is large discs on
relatively thin shafts (e.g. flywheel of a punch press).
Torsional vibrations may originate from the following forcing
Inertia forces of reciprocating mechanisms (such as pistons in IC
engines)
Impulsive loads occurring during a normal machine cycle (e.g. during
operations of a punch press)
Shock loads applied to electrical machinery (such as a generator line
fault followed by fault removal and automatic closure)
Torques related to gear mesh frequencies, turbine blade passing
frequencies, etc.
For machines having massive rotors and flexible shafts (where the
system natural frequencies of torsional vibrations may be close to, or
within, the excitation torque frequency range during normal operation)
torsional vibrations constitute a potential design problem area.
In such cases designers should ensure the accurate prediction of
machine torsional frequencies and frequencies of any torsional load
fluctuations should not coincide with the torsional natural frequencies.
Hence, the determination of torsional natural frequencies of the system is
very important.
Dr. R. Tiwari ([email protected]) 3
Consider a rotor system as shown Fig.1
The shaft is considered as massless and it provides torsional stiffness only.
The disc is considered as rigid and has no flexibility.
If an initial disturbance is given to the disc in the torsional mode and
allow it to oscillate its own, it will execute the free vibrations as shown in
Fig. 2.
It shows that rotor is spinning with a nominal speed of and executing
torsional vibrations, (t), due to this it has actual speed of { + (t)}.
It should be noted that the spinning speed remains same however angular
velocity due to torsion have varying direction over a period.
The oscillation will be simple harmonic motion with a unique frequency,
which is called the torsional natural frequency of the rotor system.
Simple System with Single Rotor Mass
Fig.1a A single-mass cantilever rotor system Fig.1(b) Free body diagram of disc
Fig. 2 Torsional vibrations of a rotor
Dr. R. Tiwari ([email protected]) 4
From the theory of torsion of shaft, we have
where,K
t
is the torsional stiffness of shaft,I
p
is the rotor polar moment of
inertia, kg-m
2
,J is the shaft polar second moment of area, l is the length
of the shaft and is the angular displacement of the rotor. From the
free body diagram of the disc as shown in Fig. 1(b)
Equation (2) is the equation of motion of the disc due to free torsional
vibrations. The free (or natural) vibration has the simple harmonic
motion (SHM). For SHM of the disc, we have
where is the amplitude of the torsional vibration and is the torsional
natural frequency. On substituting Eqs. (3) and (4) into Eq. (2), we get
t
T GJ
K
l
= =
External torque of disc
p t p
I K I = =

(1)
(2)

( ) sin
nf
t t =
2 2

sin
nf nf nf
t = =

(3)
(4)
/
nf t p
K I = (5)

nf

Dr. R. Tiwari ([email protected]) 5

A two-disc torsional system is shown
in Fig. 3. In this case whole of the
rotor is free to rotate as the shaft
being mounted on frictionless
bearings.
From the free body diagram in Fig.
4(a)
For free vibration, we have SHM, so
the solution will take the form
Two-Disc Torsional System
Fig.3 A two-disc torsional system
Fig. 4 Free body diagram of discs
1 2
1 2
External torque and External torque
p p
I I = =


1
1 1 2
0
p t t
I K K + =

2
2 2 1
0
p t t
I K K + =

(6)
(7)
2
1 1 nf
=

2
2 1 nf
=

and
(8)
Substituting equation (4) into equations (6) & (7), it gives
1
2
1 1 2
0
p nf t t
I K K + =
2
2
2 2 1
0
p nf t t
I K K + =
(9)
(10)
Dr. R. Tiwari ([email protected]) 6
In matrix form
The non-trial solution of equation (11) is obtained by taking determinant
of the matrix [K] which gives the frequency equation as
The roots of equation (12) are given as
From equation (11) corresponding to first natural frequency for = 0,
we get
1
=
2
1
2
2
1
2
2
0
0
t p nf t
t t p nf
K I K
K K I

(


=
(
` `

) ( )

[ ]{ } { }
0 K =
or
(11)
( )
1 2
1 2
4 2
0
p p nf p p t nf
I I I I K + =
(12)
( ) ( )
1 2 1 2 1 2
0.5
0 and
nf nf p p t p p
I I K I I
(
= = +

(13)
1
nf

Fig .5 First mode shape

(14)
Dr. R. Tiwari ([email protected]) 7
From Eq. (14) it can be concluded that,
The first root of equation (12) represents the case when both discs
simply rolls together in phase with each other as shown in Fig. 5. It
is the rigid body mode, which is of a little practical significance. This
mode it generally occurs whenever the system has free-free end
conditions (for example aeroplane during flying).
From equation (11), for , we get
The second mode shape (Eq. 15) represents the case when both
masses vibrate in anti-phase with one another.
Fig. 6 shows mode shape of two
-rotor system, showing two discs
vibrating in opposite directions.
2
nf nf
=
(15)
( ) ( )
1 1 2 1 2
1 2

0
t p p p p p t t
K I I I I I K K
(
+ =

2 1
1 2

p p
I I =
or
Fig.6 Second mode shape
Dr. R. Tiwari ([email protected]) 8
From mode shapes, we have
Since both the masses are always vibrating in opposite direction, there
must be a point on the shaft where torsional vibration is not taking place
i.e. a torsional node.
The location of the node may be established by treating each end of the
real system as a separate single-disc cantilever system as shown in Fig.
6.
Since value of natural frequency is known (the frequency of oscillation of
each of the single-disc system must be same), hence we write
where is defined by equation (13), and are torsional stiffness of two
(equivalent) single-rotor system, which can be obtained from equation
(17), as
The length l
1
and l
2
then can be obtained by (from equation 1)
1 2 1 1
1 2 2 2

l
l l l

= =
2 1 1 2 2
2
nf t p t p
K I K I = =
1 2 1 2 2 2
2 2
and
t nf p t nf p
K I K I = =
1 2
1 2
and
t t
l GJ K l GJ K = =
with
1 2
l l l + =
(16)
(17)
(18)
(19)
Dr. R. Tiwari ([email protected]) 9
Fig. 7(a) shows a two-disc stepped shaft. In such cases the actual shaft
should be replaced by an unstepped equivalent shaft for the purpose of
the analysis as shown in Fig. 7(b).
The equivalent torsional spring of shafts connected in series, can be
written as
Substituting equation (1), we have
where are equivalent lengths of shaft segments having equivalent shaft
diameter d
3
and l
e
is the total equivalent length of unstepped shaft having
diameter d
3
as shown in Fig. 7(b).
System with a Stepped Shaft
1 2 3
1 1 1 1
e
t t t t
K K K K
= + +
1 2 3
1 2 3 3 1 2 e e e e e e e
l l J J l J J l J J l l l = + + = + +
1 1 1
1 2 2 3 3
1
/ , / , /
e e e e e e
l l J J l l J J l l J J = = =
with
(20)
(21)
Dr. R. Tiwari ([email protected]) 10
From Fig. 7(b) and noting equations (17) and (19), in equivalent shaft the
node location can be obtained as
Fig. 7(a) Two discs with stepped shaft (b) Equivalent uniform shaft
( )
1
2 1
2
e e nf p
l a GJ I + =
( )
3
2 2
2
e e n p
l b GJ I + =
( )
2
1/2
1 2
1 1 2 2 3 3
1 2
1
and
p p t
e
n t
e
p p
I I K
K
I I l GJ l GJ l GJ

(
+
(
=
=
(
+ +
(

and
where
(22)
(21)
From equation (13), we have
4 4
2 3 2 2
2
2
, ,
64 4
e
e e
J
l l J d J d
J

= = =
Since above equation is for shaft
segment in which node is assumed
to be present, we can write
2 2
and
e e
a a J J b b J J

= =
Above equations can be combined as
a a
b b

(23)
Dr. R. Tiwari ([email protected]) 11
When there are several number of discs in the rotor system it becomes is
multi-DOF system. When the mass of the shaft itself may be significant
then the analysis described in previous sections (i.e. single or two-discs
rotor systems) is inadequate to model such system, however, they could
be extended to allow for more number of lumped masses (i.e. rigid discs)
but resulting mathematics becomes cumbersome.
Alternative methods are:
o Transfer matrix methods
o Methods of mechanical impedance
o Finite element methods
Transfer matrix method:
A multi-disc rotor system, supported on frictionless supports, is shown in
Fig. 7. Fig. 8 shows the free diagram of a shaft and a disc, separately. At
particular station in the system, we have two state variables: the angular
twist and Torque T.
MODF Systems
Dr. R. Tiwari ([email protected]) 12
Governing equations of motion of the whole system.
Point Matrix:
The equation of motion for the disc 2
is given by
For free vibrations, angular
oscillations of the disc is given by
Substituting back into equation (24),
we get
Angular displacements on the either
side of the rotor are equal, hence
Equations (26) and (27) can be
combined as
Fig. 9(a) Free body diameter of shaft section 2
2
2 2 2 p R L
T T I =

2

sin t =
2 2
2 2

sin
nf nf
t = =

2
2
2 2 2 nf p R L
T T I =
2 2 R L
=
Fig. 8 A multi-disc rotor system
(b) Free body diagram of rotor section 2
2
2 2
1 0
1
nf p
R L
I
T T

(

=
` ` (

) )

{ } [ ] { }
R 2 2 2
S
L
P S = or
where {S}
2
is the state vector at station 2 and [P]
2
is the point matrix for station .
(24)
(25)
(26)
(27)
(28)
Dr. R. Tiwari ([email protected]) 13
Field Matrix:
For shaft element 2 as shown in Fig. 9(a), the angle of twist is related to its
torsional stiffness and to the torque, which is transmitted through it, as
Since the torque transmitted is same at either end of the shaft, hence
Combining (29) and (30), we get
where [F]
2
is the field matrix for the shaft element 2.
Now we have
where [U]
2
is the transfer matrix, which relates the state vector at right of station
2 to the state vector at right of station 1.
2 1
2
T
K
=
2 1 L R
T T =
{ } [ ] { }
2 2 1 L R
S F S =
{ } [ ] { } [ ][ ] { } [ ] { }
2
2 2 1 1
2 2 2
R R
R
S P S P F S U S = = ==
(29)
(30)
(31)
(32)
2 2 2
1 1
0 1
R L
k
T T

(
=
` `
(
) )
or
Dr. R. Tiwari ([email protected]) 14
On the same lines, we can write
where [T] is the overall system transfer matrix
The overall transformation can be written as
For free-free boundary conditions,
On using equation (35) into equation(34), the second set of equation gives
On taking
In Eq. (36), t
11
contains
n
so for each value of
n
different value of
R

4
is
obtained and using Eq. (34) relative displacements of all other stations can be
obtained, by which mode shapes can be plotted.
{ } [ ] { }
{ } [ ] { } [ ] [ ] { }
{ } [ ] { } [ ] [ ] [ ] { }
{ } [ ] { } [ ] [ ] [ ] { } [ ]{ }
1 1 0
1 0
1
2 1 2 1 0
2
3 2 3 2 1 0
3
1 0 n n n n n
R
R
R R
R R
R R
S U S
S U S U U S
S U S U U U S
S U S U U U S T S =

=
= =
= =
= =

11 12
21 22
0 n R R
t t
t t T T

(
=
` `
(
) )
0
0
R n R
T T = =
( )
21
0
nf
t =
0
0
R

since
( )
11 0 0
1 we get
nf R R
t = =
(33)
(34)
(35)
(36)
Dr. R. Tiwari ([email protected]) 15
Example 1.
Obtain the torsional natural frequency of the system shown in Fig.10
using the transfer matrix method. Check results with closed form
solution available. Take G = 0.81011 N/m2
Solution:
We have following properties of the rotor
The torsional stiffness is given as
Analytical method: The natural frequencies in the closed form are given
as
Mode shapes are given as
Fig.10
4
;
11 2 4

-6
0.8 10 N/m 0.6 m; (0.1) 9.82 10 m
32
G l J

= = = =
-6
11
6

0.8 10
9.82 10 1.31 10 Nm/rad
0.6
t
GJ
k
l

= = =
( )
1 2
2 2
1 2
6
( )
22.6 5.66 1.31 10
0; and 537.97 rad/sec
22.6 5.66
p p t
n n
p p
I I k
I I

+
+
= = = =

1
0
n
=
{ } { }
2 0
R
=
2
537.77 rad/s
n
= { } { } { }
1
2
2 0 0
4.0
R
p
p
I
I
= =
For
For
Dr. R. Tiwari ([email protected]) 16
Example 1 contd...
Transfer matrix method:
State vectors can be related between stations 0 & 1 and 1 & 2, as
The overall transformation of state vectors between 2 & 0 is given as
On substituting values of various rotor parameters, it gives
1 1 0
2 2 2 1 2 2 1 0
{ } [ ] { }
{ } [ ] [ ] { } [ ] [ ] [ ] { }
R
R R
S P S
S P F S P F P S
=
= =
( )
( ) ( )
2 1 1
2 2
1
2 1 2 2
2 2 2
2 2
2 0 0
2
2 2 2 2
0
1 1
1 0 1 0 1 0
1 1

1 1 1
1 0 1
1 1
1 1
R
t
t
n p n p n p n p n p t
n p t t
n p n p p t n p t
k
k
I I I
I I k T T T
I k k
T
I I I k I k

(
( ( (
(
= = (
` ` ` ( ( (
(


) ) ) (


(

(
=
`
(
)

( )
( )
5 2 7
2 5 4 2 7 2
2 0
1 1.73 10 7.64 10
5.66 9.77 10 22.6 9.77 10 1
R
n
n n n n
T T

(


(
=
` `
(
+ + ) )

Dr. R. Tiwari ([email protected]) 17
Since ends of the rotor are free, the following boundary conditions will
apply
On application of boundary conditions, we get the following condition
Since , we have
which gives the natural frequency as
which are exactly the same as obtained by the closed form solution.
Mode shapes can be obtained by substituting these natural frequencies
one at a time into equation (A), as
Example 1 contd...
0 2
0
R
T T = =
2 5 4
21 0
[ 28.26 9.77 10 ]{ } 0
n n
t

= + =
{ }
0
0
2 5 2
[9.77 10 28.26] 0
n n

=
1 2
0 and 537.77 rad/sec
n n
= =
1
0
n
=
{ } { }
2 0
R
=
2
537.77 rad/s
n
=
{ } { }
2 0
4.0
R
=
which are also exactly the same as obtained by closed form solutions
Dr. R. Tiwari ([email protected]) 18
Example 2.
Find torsional natural frequencies and mode shapes of the rotor system
shown in Figure 1. B is a fixed end and D1 and D2 are rigid discs. The
shaft is made of steel with modulus of rigidity G = 0.8 (10)11 N/m2 and
uniform diameter d = 10 mm. The various shaft lengths are as follows:
BD-1 = 50 mm, and D1D2 = 75 mm. The polar mass moment of inertia of
discs are: Ip1 = 0.08 kg-m2 and Ip2 = 0.2 kg-m2. Consider the shaft as
massless and use (i) the analytical method and (ii) the transfer matrix
method.
Fig. 11
Solution:
Analytical method:
From free body diagrams of discs as
shown in Fig.12, equations of motion can be written as
1
2
1 1 1 2 1 2
2 2 2 1
( - ) 0
( - ) 0
p
p
I k k
I k


+ + =
+ =

The above equations for free vibrations and they are homogeneous second
order differential equations.
In free vibrations discs will execute simple harmonic motions.
Dr. R. Tiwari ([email protected]) 19
For the simple harmonic motion, hence equations of motion take the form
On taking determinant of the above matrix, it gives the frequency
equation as
which can be solved for , as
For the present problem the following properties are gives
1
2
2
1 2 2
1
2
2
2 2
-
0
0
p n
p n
k k I k
k k I

( +

=
(
` `

) ( )

1 2 1 2 2
4 2
2 1 2 1 2
( ) 0
p p n p p p n
I I I k I k I k k k + + + =
( )
1 2 2 1 2 2 1 2
1 2
2
2 1 2 2 1 2 1 2
2
4
2
p p p p p p p p
n
p p
I k I k I k I k I k I k k k I I
I I

+ + + +
=
Fig.12 Free body diagram of discs
1 2
1 2
1 2
1 2
2 2
1878N/m and 523.598N/m
0.08 kgm and 0.2 kgm
p p
GJ GJ
k k
l l
I I
= = = =
= =
Example 2 contd...
Dr. R. Tiwari ([email protected]) 20
Natural frequencies are obtained as
The relative amplitude ratio can be obtained as (Fig.13)
1 2
44.792 rad/s and 175.02 rad/s
n n
= =
2
1 2
2
2
1
2 2
-
0.2336 for and -10.700 for
p n
n n
k I
k

= =
Fig.13 Mode shapes
Example 2 contd...
Dr. R. Tiwari ([email protected]) 21
Transfer matrix method
For Figure 4.14 state vectors can be related as
The above state vector at various stations can be related as
On application of boundary conditions the second equation of equation
(A), we get
Fig.14 Two-discs rotor system with station numbers
{ } [ ] [ ] [ ] [ ] { }
2 2 2 1 1 0
R
P F P F =
1 1
2
2 2
2 1 2 2
2
2 0
1 2
1 1
1- 1
1
n p n p
R
n p
I I
k k k k
T T
I
p
p
k k

(
| | | |
+ ( | |
| |
(

\ . \ .
=
` ` (
| |
) )
(
+
|
(
|
\ .

2
2
0
1 2
0 0 1
n p
R
I
p
p T
k k
| |
= + +
|
|
\ .
0
0
R
T
since
Example 2 contd...
Dr. R. Tiwari ([email protected]) 22
On substituting for p, we get
which can be solved to give
It should be noted that it is same as obtained by the analytical method
Example 2 contd...
2 2
2 1
2 2
2 2
1 2 2
-
1
- 1 1 0
n p n p
n p n p
I I
I I
k k k


(
| |
+ + = ( |
|
(
\ .

2 1 2 1 1 2
2
2 1 2 1 1 2
2 2 1 1
4
2 4
n
p p p p p p
k k k k k k
I I I I I I

| | | |
= + +
| |
| |
\ . \ .
Dr. R. Tiwari ([email protected]) 23
In some machine the shaft may not be continuous from one end of the
machine to the other, but may have a gearbox installed at one or more
locations. So shafts will be having different angular velocities as
shown in Fig.15(a).
For the purpose of analysis the gear system must be reduced to
system with a continuous shaft so that they may be treated as
described in the preceding section.
Geared Systems
Fig. 15(a) Actual system (b) Equivalent system without gearbox
Dr. R. Tiwari ([email protected]) 24
In actual system as shown in Fig. 15(a) the shaft torsional stiffness k
2
and rotor moment of inertia I
p2
.
Let the equivalent system as shown in Fig.15(b) has the shaft torsional
stiffness k
e
and the disc moment of inertia I
pe
.
The strain and kinetic energies must be the same in both the real and
equivalent systems for theoretical model to be valid.
By imagining the rotor I
p2
to be held rigidly whilst shaft 1 is rotated
through some angle
1
at the gearbox. The shaft 2 is rotated through
an angle
1
/n at gearbox, where n is the gear ratio.
The strain energy stored in shaft 2
While applying the same input at the gear location to the equivalent
system results in the stain energy stored in the equivalent shaft and
can be expressed as
( )
1 1 2
2 2 2 1
2 2
2
E k k n = =
1 2
1
2
e
E k =
(37)
(38)
Dr. R. Tiwari ([email protected]) 25
Equating equations (37) and (38) give
(39)
If consideration is now given to the kinetic energies of both the real
and equivalent systems, which must also be equated
where
2
and
e
are the angle of twist of actual shaft 2 and the shaft 2
in equivalent system respectively. It can be seen from Figure 15(b)
that
e
=
1
and
1
and
2
are the angular frequencies of the shafts 1
and 2, respectively. Noting that
Equation (40) can be written as
which simplifies to
2
2
/
e
k k n =
( ) ( )
1
1 1
2 2 1
2 2
2
2 2
p p e
e
I I + = +

(40)
( )
1
1 2
1
2 1
and
e
nT k T k
e
= = =
2
2 2
1
1
2
2
1
1
2
p e
e
d nT d T
I I
n dt k dt k

( (

+ = +
( ` ` (
(
) )

2
1
2
2 2
2 2
1
2
2
p
e
I
d n T d n T
I
dt k dt k
n

( (


( (
+ = +
` `
( (

) )

2
2
/
p e
I I n = (41)
Dr. R. Tiwari ([email protected]) 26
The general rule, for forming the equivalent system for the purpose of
analysis, is to divide all shaft stiffness and rotor polar mass moment of
inertias of the geared system by n
2
(where n is the gear ratio n =
2
/
1
=
geared shaft
/
reference
). When analysis is completed, it should be
remembered that the elastic line of the real system is modified (as
compared to with that of the equivalent system) by dividing the
displacement amplitudes for equivalent shaft by gear ratio n as shown
in Fig.16.
Fig.16 The elastic line in the original system
Dr. R. Tiwari ([email protected]) 27
Example 3
For geared system as shown in Fig.16 find the natural frequency and
mode shapes. Find also the location of nodal point on the shaft (i.e. the
location of the point where the angular twist due to torsional vibration is
zero). The shaft A has 5 cm diameter and 0.75 m length and the shaft
B has 4 cm diameter and 1.0 m length. Take modulus of rigidity of the
shaft G equals to 0.8 1011 N/m2, polar mass moment of inertia of
discs are IA = 24 Nm2 and IB = 10 Nm2. Neglect the inertia of gears.
20 cm diameter
10 cm diameter
G
e
a
r

P
a
i
rB
A
Fig.17
Solution:
On taking shaft B has input shaft (or
reference shaft) as shown in Fig.18 the gear
ratio can be defined as
input speed 20 rpm of reference shaft
Gear ratio 2
output speed 10 rpm of driven shaft
A A B B
B B A A
D T N
n
D T N

| |
= = = = = = = = =
|
\ .
Dr. R. Tiwari ([email protected]) 28
The area moment of inertia and the torsional stiffness can be obtained
as
4 7 4

6.136 10 m ;
32
A A
J d

= =
4
4
6.545 10 Nm/rad
2.011 10 Nm/rad
A
A
A
B
GJ
K
l
K
= =
=
On replacing shaft A with reference to the shaft B by an equivalent
system, the system will look as shown in Figure 4.19. The equivalent
system of the shaft system A has the following torsional stiffness and
mass moment of inertia properties
4 7 4

2.51 10 m
32
B B
J d

= =
4
4 2

2 2 2 2
6.545 10 24
1.6362 10 N m/rad and 6 Nm
2 2
A PA
A PA
e
K I
K I
n n

= = = = = =
which gives the equivalent length as
11 7
4
0.8 10 2.513 10
1.2288m
1.6362 10
B
A
e
A
e
GJ
l
K

= = =

Example 3 contd..
Dr. R. Tiwari ([email protected]) 29
The equivalent stiffness of the full shaft is given as
The equivalent shaft length is given as
The natural frequency of the equivalent two mass rotor system as
shown in Fig.19 is given as
Fig.19 Equivalent single shaft system
4

4 4
1 1 1 1 1
1.1085 10 m/N
2.011 10 1.6362 10
e A B
e
K K K

= + = + =

which gives 9021.2 N/m
e
K =
4
2
e A B A B
4
e
1.2288 1 2.2288 m
B
A
d
l l l n l l
d
= + = + = + =
( )
( )
1/ 2
1/ 2
2
( )
6 10 9021.2/9.81
153.62 rad/sec
( ) 6 10 /9.81
PA PB e
e
n
PA PB
e
I I K
I I

( +
( +
( = = =
(

(


Example 3 contd..
Dr. R. Tiwari ([email protected]) 30
The node location can be
obtained from Fig. 20 as
The negative sign indicates that both discs are either ends of node
location
The absolute location of the node position is given as
Also from Fig. 20 we have which gives
Fig. 20 Mode shape and nodal point location in the equivalent system
1
2
10
1.667
6
A n
e PB
n B PA
e
l
I
l I

= = = =
1 2
1.667
n n
l l =
1 2
2.2288
n n
l l + =
2
0.8358 m
n
l =
Let 1rad then 1.667rad
B A
e
= =
Hence, 0.8333rad
A
A
e
n

= =
Example 3 contd..
Dr. R. Tiwari ([email protected]) 31
The mode shape and node location in the actual system is shown in
Fig. 21
Alternative way to obtain natural frequency is to use the equivalent
two mass rotor (Fig.19) can be considered as two single DOF
systems (Fig. 22).
The stiffness and mass moment of inertia properties of the system is
given as
Fig. 21 Mode shape and nodal point
location in the actual system
Fig. 22 A single DOF system
11
4 2
2
2
0.8 10 2.513 10
2.435 10 N/m and kgm
0.8358 9.81
B
l PB
n
e
GJ
K I
l

= = = =
Example 3 contd..
Dr. R. Tiwari ([email protected]) 32
It gives the natural frequency as
which is same as obtained earlier. The whole analysis can be done by
replacing shaft B with reference to shaft A speed by an equivalent
system. For more clarity some of the basic steps are given as follows.
4
2
2
2.435 10 9
154.62 rad/sec
10
l
PB
n
e
K
I

= = =
Fig.23
Actual and equivalent geared systems
Example 3 contd..
Dr. R. Tiwari ([email protected]) 33
It is assumed here that we are choosing reference shaft as input shaft
(i.e. for present case shaft A is reference shaft hence it is assumed to
be input shaft and according the gear ratio will be obtained).
It is assumed that equivalent shaft (i.e. B) has same diameter as the
reference shaft (i.e. A). The equivalent mass moment of inertia and
stiffness can be written as
The total equivalent length and the torsional stiffness would be
10
0.5
20
B B
A A
D
n
D

= = = =
4
2 4 2

2 2 2
2.011 10
40N m and 8.044 10 N/m
(0.5)
PB PB
PB B
e e
I I
I K
n n

= = = = =
11 7
4
4
0.8 10 6.136 10
8.044 10 0.610 m
8.044 10
A
B B
B
e
e
e
GJ
K l
l

= = = =

A B
11 7
4 2
e
e
0.61 0.75 1.36 m
and
0.8 10 6.136 10
3.61 10 N/m
1.360
A
e
l l l
GJ
K
e
l

= + = + =

= = =
Example 3 contd..
Dr. R. Tiwari ([email protected]) 34
The natural frequencies of two mass rotor system is given by
1
0
n
=
1 1
2 2
2
( )
( )
9.81
PA PB
PA PB
PA PB PA PB
e e
e e
e e
I I
I I
K K
n
I I I I

( (
+
+
( ( = =
( (

and
2
1
2
4
(24 40)
9.81 3.61 10 153.65 rad/s
24 40
n

+
(
= =
(

Fig.24 Equivalent two mass rotor system
The node location can be obtained as
1
2
10
1.667
6
PB n
n PA
e
I l
l I
= = =
we have
1 2 1 2
1.36
e e n n
l l l l + = + =
1 2 1 2
(1.667 ) 1.36 0.85m and 0.51m
n n n n
l l l l + = = =
which gives
Example 3 contd..
Dr. R. Tiwari ([email protected]) 35
The stiffness of will be (equivalent stiffness corresponding to shaft A
speed)
The shaft stiffness corresponding to shaft B speed can be defined in
two ways i.e
On equating above equations the location of the node in the actual
system can be obtained as
which is same as by previous method.
2
2
A
B
e
n
GJ
K
l
=
2
2
B
B
GJ
K
l
=
2 2
2
2
2
A
n
B B
e
GJ
K n K n
l
= =
2
2
2
0.84
B
A
e
J
l n l
J
= =
Example 3 contd..
Dr. R. Tiwari ([email protected]) 36
For the marine vessel power transmission shafts or machine tool
drives, there may be many rotor inertias in the system and gear box
may be a branch point where more than two shafts are attached. In
such cases where there are more than two shafts attached as shown in
Fig. 25 to the gearbox the system is said to be branched
For Fig.25 having branched system, state vectors for different branch
can be written as
For branch A,
Taking
0A
=1 as reference value
for angular displacement
since left hand end of branch A is free end
hence for free vibrations we have T
0A
= 0.
Equation for branch A takes the form
Branched Systems
Fig.25 Branched system
{ } [ ]{ } { } [ ]{ } { } [ ]{ }
;
0 0 0
;
nA nC
A nB B C
S A S S B S S C S = = =
11 12
21 22
1
0
nA
a a
T a a

(

( =
` `
(
) )

(42)
(43)
Dr. R. Tiwari ([email protected]) 37
which can be expanded as
At branch point, between shaft A and B, we have
where n
AB
is the gear ratio between shaft A and B. For branch B, T
nB
= 0,
since end of branch is free.
For branch B from equation (42), and noting the condition above, we have
Above equation can be expanded as
21 11
and
nA nA
T a a = =
nA
oB
AB
n

=
( )
11 12
21 22
0
nA AB
nB
oB
b b n
b b
T


(


(
=
` `
(

(
)
)
( )
12
11
nA AB nB
oB
b n b T = +
( )
( )
( )
22 21 22 21
0
nA AB nA AB oB oB
b n b T T b b n = + =
(44)
(45)
(46)
(47)
and
(48)
Dr. R. Tiwari ([email protected]) 38
At branch C, we have the following condition (nothing equation 45)
Another conditions at branch to be satisfied regarding torque
transmitted at various branches is
Equation (50) can be written as
On substituting equation (48) and (43), we get
Substituting (49) and (51) into equation (42), we get:
11
nA
AC AC
oc
a
n n

= =
1 1 1
2 2 2
nA nA oB nB oC nC
T T T = +

oB oC oB oC
nA nB nA nC AB AC
nA nA
T T T T
T T
n n
= + = +
oB
oC
AB
AC nA
T
T n T
n
(
=
(

21 11
2
22
21 oC
AB
AC
b a
T n a
b n
(
= +
(

(49)
(50)
(51)
(52)
/
11
11 12
21 22
21 21 11 22
2
0
AC
AC AC AB
nc
a n
c c
c c
n a b a n b n

(


(
=
` `
(

(
+ )
)
Dr. R. Tiwari ([email protected]) 39
From equation (52) second equation will give the frequency equation as
The roots of the above equation are system natural frequencies.
As before, these frequencies may then be substituted back into
transfer matrices for each station considered, where upon the state
vector at each station may be evaluated. The plot of angular
displacement against shaft position then indicates the system mode
shapes.
Using this method, there will not be any change in elastic line due to
gear ratio, since these have now already been allowed for in the
analysis.
Moreover, for the present case we have not gone for equivalent
system at all.
For the case when the system can be converted to a single shaft the
equivalent system approach has the advantage.
(53)
22 21 11
22 21
22
21 11
2
0
AC
AB
AC
AC
c b a n c a
n c a
n
b n
+ + =
Dr. R. Tiwari ([email protected]) 40
Example 4
Obtain the torsional critical speeds of the branched system as shown in
Fig. 26. Take polar mass moment of inertia of rotors as: I
PA
= 0.01 kg-
m2, I
PE
= 0.005 kg-m2, I
PF
= 0.006 kg-m2, and I
PB
= I
PC
= I
PD
= 0. Take
gear ratio as: n
BC
= 3 and n
BD
= 4. The shaft lengths are: l
AB
= l
CE
= l
DF
=
1 m and diameters are d
AB
= 0.4 m, d
CE
= 0.2 m and d
DF
= 0.1 m. Take
shaft modulus of rigidity G = 0.8 1011 N/m2
Solution:
The branched system has the following
mass moment of inertias
For branch A the state vector at stations are related as
Fig. 26 A branched rotor system
01 . 0 =
PA
I
0.005
PE
I =
0.006
PF
I =
OA A nA
S U S } { ] [ } { =
[ ]
A AB A
P F U ] [ ] [ =
11 2 9 9
2 2
1 0 1 4.97 10 4.97 10 1 4.97 10
0.01 1 0.01 1 0 1
n
n n

( ( (
= =
( ( (


with
Dr. R. Tiwari ([email protected]) 41
For branch B the state vector at stations are related as
Similarly, for branch C, we have
From equation (53), the frequency equation can be written as
On substitution, we get
{ } [ ] { }
nB B OB
S U S =
[ ]
8
2 10 2
1 7.95 10
[ ] [ ]
0.005 3.97 10 1
B CE
B
n n
U P F

(
= =
(
+

with
{ } [ ] { }
nC C OC
S U S =
[ ]
6
2 9 2
1 1.27 10
[ ] [ ]
0.006 7.64 10 1
B BF
C
n n
U P F

(
= =
(
+

with
22 21 11
22
22
21 11
21
2
0
BD
BC
BD
BD
c a c b a n
n c a
n
b n
+ + =
( )
2 11 2
9 2 2
9 2 2 11 2
10 2
0.006 (1 4.97 10 )
4 1 7.64 10 ( 0.01 )
4
(1 7.64 10 )( 0.005 )(1 4.97 10 ) 4
0
9 (1 3.97 10 )
n n
n n
n n n
n

+

+ =

Example 4 contd..
Dr. R. Tiwari ([email protected]) 42
which can be simplified to
The roots of the polynomial are
Natural frequencies are given as
It can be seen that the rigid body mode exist since ends of the gear
train is free.
2 10 4 19 6
0.04372 3.3921 10 1.21 10 0
n n n


+ =
1,2,3
2 6
0; 135.48 10
n
= and
9
10 64 . 2
0
1
=
n

2
11640
n
=
3
51387
n
= rad/s
Example 4 contd..
Dr. R. Tiwari ([email protected]) 43
Damping may come from
o The shaft material
o Torsional vibration dampers.
The torsional vibration damper is a device which may be used to join
together two-shaft section as shown in Fig. 27. It transmits a torque,
which is dependent upon of the angular velocity on one shaft relative to
the other.
Torsional dampers can be used as
a means of attenuating system vibration
and to tune system resonant frequencies
to suit particular operating conditions.
The damping in the system introduces
phase lag angles to the system
displacement and torque. The displacement and torque parameters must
now be represented mathematically both in-phase and quadrature
components.
Fig. 28 shows a general arrangement of torsional MDOF rotor system
with damping. EOM of the nth rotor from free body diagram (Fig. 29) can
be written as
Damping in Torsional Systems
Fig. 27 Torsional vibration damper
Dr. R. Tiwari ([email protected]) 44
And
Torques
R
T
n
and
L
T
n
may be written in
whilst the angular displacements is
Differentiating equations (56) and (57) with
respect to time and substituting
in equations (54) and (55) leads to
n n R n L n
T T C

L n R n
=
(54)
(55)
Fig. 28 General arrangement of MDOF system with damping.
Fig. 29 Free body diagram of rth rotor.
1 2
sin cos
n n n
T T t T t = +
1 2
sin cos
n n n
t t = +
(56)
(57)
1 1
2 2
1 1
2 2
2
2
1 0 0 0
0 1 0 0
1 0
0 1
n L n
n
T T I C
T T
C I




(

(

(

=
` `
(


(

(
) )

(58)
Dr. R. Tiwari ([email protected]) 45
The previous equation can be simplified as
The characteristics of the shaft element
at station n are represented in the
equation describing the torque applied
to the shaft at the location of rotor n, as
While the torque transmitted through the shaft is the same at each
ends i.e.
Substituting T, and in equation (60) and on separating the in-
phase and quadrature components , we get
{ } [ ] { }
L n
R n n
S P S =
Fig. 30 Free body diagram of nth shaft segment
( ) ( )
( )
1 1 1 R n L n n n n n n n
T T K B

= = +

1 L n R n
T T

=

( ) ( 1) ( 1) ( 1)
1 1 1 1 2 2
R n L n L n n R n n L n n R n
n
T T K K B B

= =
( ) ( 1) ( 1) ( 1)
2 2 2 2 1 1
R n L n n L n n R n n L n n R n
T T K K B B

= = +
(59)
(60)
(61)
(62)
Dr. R. Tiwari ([email protected]) 46
Substituting T, , and in equation (61), we get
On separating in-phase and quadrature components, we get
Combining (62) and (63) we get:
which can be written as
which can be simplified as
From equations (59) and (65), we get

( ) ( )
( 1) ( 1)
1 2 1 2
sin cos sin cos
L n L n R n R n
T t T t T t T t

+ = +
( 1) ( 1)
1 1 2 2
and
L n R n L n R n
T T T T

= =
1 1
2 2
1 1
2 2
1
0 0 1 0
0 0 0 1
0 0 1 0 0 0 1 0
0 0 0 1 0 0 0 1
n n L n R n
K B K B
B K B K
T T
T T

( (
( (

( (
=
` `
( (

( (

) )
1
[ ] { } [ ] { }
n L n n R n
L S M S

=
1
1 1
{ } [ ] [ ] { } [ ] { }
L n n n R n n R n
S L M S F S

= =
1 1
{ } [ ] { } [ ] [ ] { } [ ] { }
R n n L n n n R n n R n
S P S P F S U S

= = =
Remaining analysis will remain same for obtaining natural frequency
& made shapes.
(63)
(64)
(65)
(66)
Dr. R. Tiwari ([email protected]) 47
Equations of motion of continuous system can be derived by (i) force
balance of a differential element (ii) Hamiltons principle (iii) Lagranges
method. While the force balance is a convenient approach for most of the
problems, Hamiltons principle and Lagranges method are applied for a
complex system. These two approaches need the consideration of the
energy of the system.
Hamiltons Principle
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
where L is the Lagrangian, T and V are kinetic and potential energy of
the system and W is the work done. The physical interpretation of
equation (67) is that out of all possible paths of motion of a system during
an interval of time from t
1
to t
2
, the actual path will be that for which the
integral has a stationary value. It can be shown that in fact 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
Torsional Vibration for Continuous Systems
[ ]
2 2
1 1
( ) 0
t t
t t
T V W dt Ldt = + = =
} } (67)
Dr. R. Tiwari ([email protected]) 48
Lagranges Equation
Hamiltons principle is stated as an integral equation where total
energy is integrated over an time interval.
On the other hand, Lagranges equation is differential equations, in
which one considers energies of the system instantaneously in time.
Hamiltons principle can be used to derive the Lagranges equation in
a set of generalized coordinates. Lagranges equation can be written
as
The study of torsional vibration of rotors is very important especially in
applications where high power transmission and high speed are there.
As compared to transverse vibrations governing equations are much
simpler in torsional vibration of shafts and it is identical to the axial
vibration of rods
0 =

|
|
.
|

\
|

i i
q
L
q
L
dt
d
(68)
Dr. R. Tiwari ([email protected]) 49
Governing Differential Equations
To obtain the governing differential equations the
variational principle is used that requires the system
potential energy functional.
P(y, z)
(x, t) (x, t)
M(x, t)
x
0
x y
z
l
0
x
y
z
Figure 31 A rod under time dependent torque
Dr. R. Tiwari ([email protected]) 50
The potential energy functional is obtained by considering the
displacement field of any point P, as shown in Figure 31 and 32.
P (y, z)
P(y, z)
y
z
u
z
u
y
o
(x, t)
P (y, z)
P(y, z)
r
y
z
u
z
u
y
Q
r
z
y

Figure 2 A point under torsional displacement

Dr. R. Tiwari ([email protected]) 51
The displacement field at a point P is
where (x, y, z) is the Cartesian coordinate of the point P, (r, , z) is the
cylindrical coordinate of the point and (x, t) is the angular displacement
and u(x, y, z) is the linear displacement of the point.
The strain field is given by
where etc.
The stress field is given by
) , ( cos ); , ( sin ; 0 t x y r u t x z r u u
z y x
= = = = =
( )
( ) ( )
1
, , , , ,
2
1 1 1 1
, , , , , ,
2 2 2 2
0; 0; 0; 0;
; ;
xx x x yy y y zz z z yz y z z y
xy y z z y x zx z x x z x
u u u u u
u u z u u y


= = = = = = = + =
= + = = + =
x
u
u
x
x x

=
,
(69)
(70)
Dr. R. Tiwari ([email protected]) 52
(71)
The strain energy is given by
(72)
The kinetic energy is given by
(73)
; ;
; 0 ; 0 ; 0
yz
x
Gz G
x
Gy G
G E
xy zx xy zx
yz zz yy xx xx

= =

= =
= = = = = =

1 1 1 1 1 1
2 2 2 2 2 2
2 2 2
1 1 2 2
2 2
0 0
(
1
2
xx xx yy yy zz zz xy xy yz yz zx zx
V
l l
A
U dV
Gz Gy dAdx GJ dx
x x x


(
= + + + + +

(

| | | | | |
= + =
(
| | |

\ . \ . \ .
(

}
}} }
( )
2 2
with
A
J y z dA = +
}
{ }
1 2 2 2 2
2
0 0
1
2
l l
x y z
A
T u u u dAdx J dx
(
= + + =
(

} } }


Dr. R. Tiwari ([email protected]) 53
where , etc. and velocity components, noting the equation (69) can
be written as (74)
The external work done is given by
(75)
The Lagrangian is = T (U+W). The following integration can be
obtained
(76)
From the Hamiltons principle, = 0 , Eq.(76) becomes
which can be written, by performing the integration by parts, as
dt
du
u
x
x
=
) , ( ); , ( ; 0 t x y u t x z u u
z y x

= = =
1 *
0 0
2
0
( )
l
W M M x x dx
(
= +
}
( )
2 2
1 1
1 1 2 2 *
, 0 0
2 2
0
( )
t t l
x
t t
T U W dt J GJ M M x x dxdt
(
= =

} } }

2
1
*
, , 0 0
0
0 ( )
t l
x x
t
J GJ M M x x dxdt
(
= =

} }

Dr. R. Tiwari ([email protected]) 54

(77)
which can be written, again by performing integration by parts, as
(78)
which can be further simplified as
(79)
Since are variations of angular displacement and
they are arbitrary, so
and (80,81)
2 2 2
1 1 1
*
, 0 0
0 0 0
( ) ( ) ( ) 0
l t t l t l
x
t t t
J dtdx GJ dxdt M M x x dxdt
t x


( (
(
+ =
( (


} } } } } }

2
2 2 2
1 1 1
1
2
1
, ,
0 0 0
0
*
0 0
0
( ) ( ) ( ) ( )
( ) 0
l
t
l l t t t l
x xx
t t t
t
t l
t
J dx J dtdx GJ dt GJ dxdt
x
M M x x dxdt

(
( +



(
+ =

} } } } } }
} }

{ }
2
1
*
, 0 0 ,
0
0
( ) ( ) ( ) 0
l t
l
xx x
t
J GJ M M x x dtdx GJ
(
+ =

} }

l
and , ,
0
*
, 0 0
( ) ( ) 0
xx
J GJ M M x x =

0 ) (
0
,
=
l
x
GJ
0
Dr. R. Tiwari ([email protected]) 55
Equations (80) and (81) represent equations of motion and boundary
conditions respectively.
Finite Element Formulation
For the finite element solution let the approximate solution for an element
be, in the form of a general polynomial, defined as
(82)
where , i = 1, 2, , r are called shape functions
Galerkin Method
On substituting equation (82) into equation (80), the residue as,
(83)
{ }
1
( )
2 ( ) 2
1 2
( )
( )
( , ) ( ) ( ) ( ) ( ) ( )
( )
ne
e
r
r
t
t
x t a b cx N x N x N x N x t
t

= + + + = =
( (
`



)

) (x N
i
( ) e
R
( ) ( ) ( ) *
, 0 0
( , ) ( ) ( )
e e e
xx
R J GJ M x t M t x x =

Dr. R. Tiwari ([email protected]) 56

On minimising the residue by using the Galerkin method, we have
(84)
On substituting equation (83) into equation (84), we get
(85)
by integrating by parts of second term in equation (85), so as to give
(86)
The interpolation function of linear form is given by
(87)
where is the dependent variable a and b are constants to be determined
by element end conditions as shown in Fig. 33. (C
0
of continuity and C
1
of
compatibility conditions).
r i dx R N
l
e
i
, , 2 , 1 ; 0
0
) (
= =
}
( ) ( ) *
, 0 0
0
( , ) ( ) ( ) 0 1, 2, ,
l
e e
i xx
N J GJ M x t M t x x dx i r
(
= =

}


( ) ( ) ( ) *
, , , 0 0
0
0 0 0
( , ) ( ) ( ) 0
1, 2, ,
l l l
h
e e e
i i x i x x i
JN dx GJN GJN dx N M x t M t x x dx
i r

(
+ + =

=
} } }

bx a
e
+ =
) (

Dr. R. Tiwari ([email protected]) 57

At any general element (e) the following end condition exists
at i
th
node x = 0 (e) = ; at j
th
node x = l (e) =
j
(88)
On application of equations (88) into equation (87) constants can be
obtained as
(89)
Noting equation (82) the shape function can be written as
(90)
i

( ) l b a
i j i
/ ; = =

{ }
) (
) (
) ( ) (
ne
e
t x N =
i

Fig.33 A beam element with end conditions

Dr. R. Tiwari ([email protected]) 58
with (91)
where N
i
and N
j
are interpolation functions corresponding to i
th
and j
th
nodes,
respectively.
Checking of Compatibility Requirement:
In Figure 34 the node j is the common node of element (e) and (e+1). For
element (e), we can write equation (90), as
(92)
{ } ( )
( )
; ; 1 / ; /
ne i
i j i j
j
N N N N x l N x l

( = = = =
(
`


)
Fig. 34 Two neighbouring elements

)
`

=
j
i
e
l x l x

) / ( ) / 1 (
) (
Dr. R. Tiwari ([email protected]) 59
Similarly for the element (e+1), we can write
(93)
We will check the compatibility requirements at common node j.
For x = l in element (e), equation (92), we get
(94)
For x = 0 in element (e+1), equation (93), we get
(95)
To verify whether the present interpolation function gives compatibility of
higher order, on taking first derivative of the equation (92) with respect to
x, we get
(96)

)
`

=
+
k
j e
l x l x

) / ( ) / 1 (
) 1 (
j
h x
e
=
=
) (
j
x
e
=
=
+
0
) 1 (

{ }
) (
,
) (
,
) (
ne
x
e
x
t N =
Dr. R. Tiwari ([email protected]) 60
with
(97)
For element (e), we can write equation (96), as
(98)
Similarly for the element (e+1), we can write
(99)
For x = l in element (e), equation (98), we get
(100)
For x = 0 in element (e+1), equation (98), we get
(101)
{ } ( )
( )
, , , , ,
; ; 1/ ; 1/
ne
i
x i x j x i x j x
j
N N N N l N l

( ( = = = =
`

)

)
`

=
j
i
e
x
l l

) / 1 ( ) / 1 (
) (
,

)
`

=
+
k
j e
x
l l

) / 1 ( ) / 1 (
) 1 (
,
( ) l
i j
l x
e
x
/
) (
,
=
=
( ) l
j k
x
e
x
/
0
) 1 (
,
=
=
+
Dr. R. Tiwari ([email protected]) 61
The completeness condition ensures the convergence as the element
size is decreased and in the limit it should converge to the exact value
for infinitesimal size element.
Equation (86) can be written in vector notation as
(102)
On substituting equation (90) into equation (102), we get
(103)
{ } { } { } { }
( ) ( ) ( ) *
, , , 0 0
0
0 0 0
( , ) ( ) ( )
l l l
h
e e e
x x x
J N dx GJ N dx GJ N N M x t M t x x dx
(
+ = + +
} } }

{ }
{ } { } { }
{ }
( ) ( )
, ,
0 0
( ) *
, 0 0
0
0
( , ) ( ) ( )
l l
ne ne
x x
h
l
i
e
x
j
J N N dx GJ N N dx
N
GJ N M x t M t x x dx
N


( + =
(



( + +
`

)
} }
}

Dr. R. Tiwari ([email protected]) 62

Noting from equation (91) that N
i
is having 1 and 0 value at node 1 and 2
respectively, and N
j
is having 0 and 1 value at node 1 and 2 respectively.
On substituting these values in equation first term of right hand side of
equation (103), we get
(104)
which can be simplified to standard finite element formulation of
equations of motion as
(105)
with
(106)
{ }
{ } { } { }
{ }
( )
,
0
( ) ( ) *
, , 0 0
( )
0 0 0
,
0
( , ) ( ) ( )
0
e
l l l
x
x
ne ne
x x
e
x
x l
GJ
J N N dx GJ N N dx N M x t M t x x dx
GJ

=
=

( ( + = + +
(
`

)
} } }

[ ] { } [ ] { } { } { }
E R
ne
e
ne
e
T T K M + = +
) (
) (
) (
) (


[ ] { }
2
( )
2
0 0
2 1
1 2 6
l l
e
i i j
j i j
N N N
Jl
M J N N dx J dx
N N N

(
(
= = =
(
(
(
(

} }
Dr. R. Tiwari ([email protected]) 63
(107)
(108)
where is the element mass matrix, is the element stiffness matrix,
is the reaction element torque vector, and is the external
element torque vector. Matrices and are called consistent mass
and stiffness matrices.
From strength of materials pure torsion theory we have
{ }
E
T
[ ] { }
2
( )
, , ,
, ,
2
, , ,
0 0
1 1
1 1
l l
e
i x i x j x
x x
j x i x j x
N N N
GJ
K GJ N N dx GJ dx
N N N l
(

(
( = = =
(
(


(

} }
{ }

=
)
`

=
=
=
l x
e
x
x
e
x
R
R
R
GJ
GJ
T
T
T
) (
,
0
) (
,
2
1

{ } { }
*
0 0
0
( , ) ( ) ( )
l
E
T N M x t M t x x dx ( = +
}
[ ]
( ) e
K
T
GJ l

=
[ ]
( ) e
M
[ ]
( ) e
M [ ]
( ) e
K
{ }
R
T
Dr. R. Tiwari ([email protected]) 64
Rayleigh-Ritz Method
Rayleight-Ritz method can be also used to develop finite element
formulation from the functional i.e. the potential energy, equation (76), of
the system. On substituting equation (82) into equation (76), get
functional for the element
(109)
minimising the functional for element, , is equivalent to
(110)
Hence, for a typical node i, on substituting the equation (109) into the
equation (110), we get
2
1
1 1 ( ) ( )2 ( )2 ( ) * ( )
, 0 0
2 2
0
( )
t l
e e e e e
x
t
J GJ M M x x dxdt
(
=

} }

( ) e

0
) (
=

i
e

r i , , 2 , 1 =
( )
2
1
( )
( ) ( ) ( )
, ( ) ( ) *
, 0 0
0
0 ( )
e
e e e
t l
x e e
x
t
i i i i
J GJ M M x x dxdt

= = +
(


} }

Dr. R. Tiwari ([email protected]) 65

0
On performing integration of parts of first term in left hand side of the
above equation, we get
On substituting equation (14) into the above equation, we get
The above equations can be written in the matrix form as
which gives the standard finite element formulation of equations of
motion, as
( )
2
2
1
1
( )
( ) ( ) ( )
, ( ) ( ) ( ) *
, 0 0
0 0
( ) 0
t
e
e e e
l t l
x e e e
x
t
i i i i
t
J dx J GJ M M x x dxdt

(

+ + + =
(


} } }

{ }
{ }
( )
2
1
( )
( )
*
0 0
0
( ) 0
ne t l
ne
i
i i
t
dN dN
J N N GJ M M x x N dxdt
dx dx

(
(
+ + + =
(
(
(


} }

r i , , 2 , 1 =
{ }
{ }
{ }
( )
{ }
2
1
( )
( )
*
0 0
0 0 0
( ) 0
ne t l l l
ne
t
dN dN
J N N dx GJ dx M M x x N dx dt
dx dx

(
(
+ + + =
(
`
(
(
)

} } } }

Dr. R. Tiwari ([email protected]) 66

(111)
where is the element mass matrix, is the element stiffness matrix
is the elemental reaction torque vector and is the elemental
external torque vector.
Assembled System Equations
Once the elemental equation of the form as equation (105) has been
derived next step is to obtain the system equation by assembling all such
elements of the system as shown in Figure 35,
[ ]
{ }
[ ] { } { } { }
( )
( ) ( ) ( ) ( ) ( )
ne
e e ne ne ne
R E
M K T T + = +

[ ]
( ) e
M [ ]
( ) e
K
{ }
( ) ne
E
T

1

2
h
1
2
x

1

2
1 2
3
(2) (1)
(b) Beam is discretised into two elements
Figure 35. (a) A cantilever beam
{ }
( ) ne
R
T
Dr. R. Tiwari ([email protected]) 67
The element finite element equation for the first element can be written, by
noting equations (105) to (108) and for first element i = 1 and j =2, as
(112)
where and are reaction forces at node 1 and at node 2 respectively.
and represent equivalent external force at node 1 and 2 of element
(1). Similarly, for the second element i = 2 and j = 3, hence the element
finite element equation can be written as by considering the equation (105)
(113)
Now equations (112) and (113) can be assembled in the following form
(114)

+
)
`

=
)
`

+
)
`

) 1 (
) 1 (
2
1
2
1
2
1
2
1
1 1
1 1
2 1
1 2
6
E
E
R
R
T
T
T
T
l
GJ Jl

1
R
T
2
R
T
) 1 (
1
E
T
) 1 (
2
E
T
2
2
3
3
( 2)
2
2
( 2)
3
3
2 1 1 1
1 2 1 1 6
R
E
R
E
T
T
Jl GJ
T l T

( (

+ = +
` ` ` `
( (

)
)
)
)

1 1
2 2 2 2
3 3
(1)
1 1
(1) (2)
2 2
(2)
3 3
2 1 0 1 1
1 2 2 1 1 1 1 1
6
0 1 2 0 0 1
R E
R R E E
R E
T T
Jl GJ
T T T T
l
T T

( (


( (
+ + + = + +
` ` ` `
( (

( (
)
)
) )

Dr. R. Tiwari ([email protected]) 68

Application of Boundary Conditions:
For the present illustration since node 1 is fixed and node 3 is free, hence
the following conditions are specified by the problem
(115,116)
On applying the boundary conditions, i.e. equations (115) and (116) into
equation (114), we get
(117)
On taking last two equations, we get
(118)
The equation (118) can be solved if external forces are specified
1 1
0 and = =

3
0
R
T =

+ +

(
(
(

(
(
(

+
) 2 (
) 2 ( ) 1 (
) 1 (
3
2
3
2
3
2 2
1 1
0
0
0
1 0 0
1 1 1 1
0 1 1 0
2 1 0
1 2 2 1
0 1 2
6
E
E E
E R
T
T T
T T
l
GJ Jl

+
=
)
`

+
)
`

) 2 (
) 2 ( ) 1 (
3
2
3
2
3
2 2
1 1
1 2
2 1
1 4
6
E
E E
T
T T
l
GJ Jl

Dr. R. Tiwari ([email protected]) 69
Free Torsional Vibration
For free vibration, the solution can be assumed of the following form
so that (119,120)
where is the amplitude of vibration, is the natural frequency and
.On substituting equations (119) and (120) into homogeneous
part of equation (118), we get
(121)
Equation (121) is the standard eigen value problem. For non-trivial solution
the following determinant should be zero, i.e.
(122)
t j
n
e
t
t

)
`

=
)
`

3
2
3
2
) (
) (
2 2 2
3
3
( )
( )
n
j t
n
t
e
t

=
` `

)
)

1 j =
)
`

=
)
`

|
|
.
|

\
|
(

+
(

0
0
1 1
1 2
2 1
1 4
6
3
2 2
l
GJ Jl
n

0
3 3
3 3
2 2
2 2
2 2
=
|
.
|

\
|

|
.
|

\
|

|
.
|

\
|

|
.
|

\
|

Jl
l
GJ Jl
l
GJ
Jl
l
GJ Jl
l
GJ
n n
n n

Dr. R. Tiwari ([email protected]) 70

On taking determinate of equation (122), we get a polynomial and it is
called characteristics equation or frequency equation of the following form
(123,124)
Equation (123) is a quadratic in and can be solved as
(125)
Note that only the positive sign is considered in equation (125), since
otherwise the natural frequency will be complex in nature.
Example 5
Determine the natural frequencies and mode shapes for the rotor system
as shown in Figure 5. Neglect the mass of the shaft and assume discs are
having lumped masses.
( ) 0
7 7
10
2
2
2
2
2
= +
b
a
b
a
n n

/ ; / 6 with a GJ l b Jl = =
2
n

|
|
.
|

\
|
|
.
|

\
|
+
|
.
|

\
|
+ =
2
2
2
2
7 7
10
4
7
10
7
10
b
a
b
a
b
a
b
a
n

Dr. R. Tiwari ([email protected]) 71

Figure 36 Two disc rotor system
Since the present problem contains only lumped masses, let us discretise
the shaft into one element as shown in Figure 6. The finite element
equation can be written, noting the equation (38), as
(126)
l
Ip
1
G J
Ip
2
mass less shaft

=
)
`

+
)
`

2
,
1
,
2
1
2
1
1 1
1 1
0
0
2
1
x
x
p
p
GJ
GJ
l
GJ
I
I

Solution:
Example 5 contd..
Dr. R. Tiwari ([email protected]) 72
Figure 36 Discretised single element Two disc rotor system.
Since for the present problem both the nodes are free, the following end
conditions apply
(127)
For simple harmonic motion, we have the following relation
(128)
On substituting equations (127) and (128) into the equation (126), we get
(129)
l
Ip
1
Ip
2

1

2
(1)
1 2
0
2 1
= = T T
1 2
. 0 i e GJ GJ

= =
)
`

=
)
`

2
1 2
2
1

n


)
`

=
)
`

(
(
(

0
0
2
1
2
2
1
1

p n
p n
I
l
GJ
l
GJ
l
GJ
I
l
GJ
Example 5 contd..
Dr. R. Tiwari ([email protected]) 73
For non-trivial solution of the equation (129), the following determinate
should be zero
(130)
;
which gives a frequency equation, it gives two torsional natural frequencies
and are given as
(131)
which are exactly the same as which we obtained by the closed form
solution. Mode shapes can be obtained for these natural frequencies from
equation (129).
0
1
1
2
2
=


p n
p n
I
l
GJ
l
GJ
l
GJ
I
l
GJ

( ) 0
2 1 2 1
2 2
=
(

+
p p p p n n
I I
l
GJ
I I
( )
1 2
1 2
0
p p
p p
k I I
p and p
I I
+
= =
Example 5 contd..
Dr. R. Tiwari ([email protected]) 74
Dr. R. Tiwari ([email protected]) 75
Dr. R. Tiwari ([email protected]) 76
Dr. R. Tiwari ([email protected]) 77
Dr. R. Tiwari ([email protected]) 78
Dr. R. Tiwari ([email protected]) 79
Dr. R. Tiwari ([email protected]) 80
Dr. R. Tiwari ([email protected]) 81
Dr. R. Tiwari ([email protected]) 82
Dr. R. Tiwari ([email protected]) 83
Dr. R. Tiwari ([email protected]) 84
Dr. R. Tiwari ([email protected]) 85
Dr. R. Tiwari ([email protected]) 86
(rad/s)
Dr. R. Tiwari ([email protected]) 87
Dr. R. Tiwari ([email protected]) 88
Geared Element for Branched Systems
Figure 37 Gear Element
Figure 7 shows the example of branched systems, Gear element.
for no slip condition (132,133)
from the state vector of the system.
Gear 1
Gear 2

g1

g2

2
= 1/n
1
1
2
Gear ratio

=
2 1
1
g g
n
=
2 1 2
Since is defined in terms of we can eliminate
g g g

Gear 1
Gear 2
Shaft for gear 2
Dr. R. Tiwari ([email protected]) 89
The equivalent inertia force of the gear pair w. r. t shaft reference 1 is
(134)
The stiffness matrix of shaft element connected to gear 2 will have to be
modified, since the angular deflection of the left hand side of that element is
now (see Figure 38.)
The P.E of the shaft element is
(135)
2
1 1
2
g
g g
I
I
n


| |
+
`
|
\ .
)

1
1
g
n

-1/n
g1
=
1

2
(2) (1)
K
1
Figure 38 Equivalent Gear Element
( )
1
2
1 1 2 1
2
1 2
1
2
1
.
2
g
U K
K
n

=
| |
= +
|
\ .
Dr. R. Tiwari ([email protected]) 90
The external work done is give as
(136)
On applying Lagrangian equation,
(137)
which is equations of motion (stiffness matrix) of the gear element.
{ }
1
1 1
1
1 1
1 1
2
2
1
2
1
1
1 2
2
0
.
0
g
g g
g
U K
K K
n n
n n
F
K
U
K
K
n
n

| |

(
= + =
|
(

\ .
=
(
` `
| |
(
)
= + =
|
(

\ .
)
Example 6.
Obtain the torsional critical speeds of the branched system as shown in
Figure 8. Take polar mass moment of inertia of rotors as: I
PA
= 0.01 kg-
m
2
, I
PE
= 0.005 kg-m
2
, I
PF
= 0.006 kg-m
2
, and I
PB
= I
PC
= I
PD
= 0. Take
gear ratio as: n
BC
= 3 and n
BD
= 4. The shaft lengths are: l
AB
= l
CE
= l
DF
= 1
m and diameters are d
AB
= 0.4 m, d
CE
= 0.2 m and d
DF
= 0.1 m. Take
shaft modulus of rigidity G = 0.8 10
11
N/m
2
.
1
1 1 1 2 2 1 2 2
g
W f f f f
n

= =
Dr. R. Tiwari ([email protected]) 91
A
C
(1)
0
.
4
m
D
B
(3)
0
.
1
m
F
E
0
.
2
m
(2)
Figure 38
Branched power system.
kg-m
2
kg-m
2
kg-m
2
kg-m
2
kg-m
2
kg-m
2
01 . 0 =
PA
I
005 . 0 =
PE
I
006 . 0 =
PF
I
005 . 0 =
PB
I
006 . 0 =
PC
I
Figure 8
006 . 0 =
PD
I
Gear ratio is defined as
Similarly,
/ ; /
BC B C C B BC
n n = =
/
D B BD
n =
Example 6 contd..
Dr. R. Tiwari ([email protected]) 92
Element (1)
(A)
Element (2)
To find stiffness matrix we will consider P.E., from equation (A)
)
`

=
)
`

+
)
`

B
A
B
A
B
A
B
A
T
T
k k
k k
I
I

1 1
1 1
0
0


C C
(2)
2
2
2
2 2
) / ( 2 / 1 ) ( 2 / 1 n k k U
B E C E
+ = =
2
/
BC BC
n I
PE
I
Figure. 39 Element (1)
Figure. 40 Element (2)
Example 6 contd..
Dr. R. Tiwari ([email protected]) 93
2
2, 2 2 2, 2 2
2
(2) 2 2
2
2 2
2 2
2 2
2 2
( / ) ( / ) ; ( / )
/ /
/
/ / / 0
/ 0
B E
BC B BC E Bc B E
BC BC
BC
B C BC
BC BC C BC B
E E BC E E
U k n k n U k n k
k n k n
k
k n k
T n k n k n I n
T
k n k I

= + = +
(
=
(

( (

+ =
` ` ` ( (
) )
)

Now,
Since
Element (3)
BC C BC
n T T / =
D F
(3)
2
/
BC PD
n I
PF
I
Similar to element (2) we can write
)
`

=
)
`

+
)
`

F
BD D
F
B
BD
BD BD
F
B
F
Bd D
T
n T
k n k
n k n k
I
n I
/
/
/ /
0
0 /
3 3
2
2
3
2

Figure. 41 Element (3)
(B)
(C)
Example 6 contd..
Dr. R. Tiwari ([email protected]) 94
1 1
2 2 2 2
1 1 2 3 2 3
2 2
3 3
0 0 0 0 0
/ / / / 0 ( / / ) 0 0 0
0 / 0 0 0 0
0 / 0 0 0 0
A A A
A
BC BD BC BD B C BC D BD B
B
BC E E E
E
BD E F F
F
k k I T
k k k n k n k n k n I I n I n
k n k I T
k n k I T

( (

( (
+ + + +

( (
+ =
` ` `
( (

( (

) )
)

The assembled equations of motion can be written as

Now to have eigen value problem (EVP),
we have
N/m; N/m; N/m
0
2
= M K
6
1
10 062 . 201 = k
6
2
10 566 . 12 = k
3
3
10 398 . 785 = k
6
10
785398 . 0 0 19635 . 0 0
0 566 . 12 18867 . 4 0
19635 . 0 18867 . 4 507 . 202 062 . 201
0 0 062 . 201 062 . 201

(
(
(
(

= k
(
(
(
(

=
006 . 0 0 0 0
0 005 . 0 0 0
0 0 007167 . 0 0
0 0 0 01 . 0
M
(D)
It should be noted that in the assembled form second row in the torque
column is zero, since we have the following condition
/ / 0
B c BC D BD
T T n T n =
(E)
(F)
and
Example 6 contd..
Dr. R. Tiwari ([email protected]) 95
Using equations (F), the natural frequency can be obtained as
rad/s; rad/s;
rad/s rad/s
5
4
10 1472 . 2 =
n

5
3
10 5085 . 0 =
n

5
2
10 1156 . 0 =
n

0
1
=
n

Relative amplitude of angular displacements (i.e. mode shapes)can be

obtained by substituting these natural frequencies one at a time in the
equation (D). It should be noted that angular displacements obtained
from equation (D) are actual displacements and no use of gear ratio is
required.
Example 6 contd..
Dr. R. Tiwari ([email protected]) 96
A multi-cylinder reciprocating machine contains many reciprocating and
rotating parts such as pistons, connecting rods, crankshafts, flywheels,
dampers and cracks.
The system is so complicated that it is difficult, if not impossible, to
undertake an exact analysis of its vibrational characteristics.
The actual system is characterised by the presence of unpredictable
effects like variable inertia, internal dampings, misalignments in the
transmission units, uneven firing order etc.
The analysis can be best carried out, by lumping the inertias of rotating
and reciprocating parts at discrete points on the main shaft. The problem
then reduces to the forced torsional vibration study of an N disc system
subjected to varying torques at different cylinder points.
The crankshaft and the other drive or driven shafts are generally flexible
in torsion, but have low polar moments of inertia, unlike in the case of
some large turbines or compressors.
Modelling of reciprocating machine systems
Dr. R. Tiwari ([email protected]) 97
On the other hand, parts mounted on the shafting, like the damper,
flywheel, generator etc. are rigid and will have very high polar mass
moments of inertia.
The system containing the crankshaft, coupling, generator/ auxiliary drive
shaft/ other driven shaft like pumps, and mounted parts can then be
reduced to a simple system with a series of rigid rotor (representing the
inertias) connected by the massless flexible shafts as shown in Figure 4.42
Figure 4.42
Dr. R. Tiwari ([email protected]) 98
Polar mass moments of inertia
Determination of polar mass moments of inertia is a straightforward
matter for rotating parts, however, it is not quite so simple in the case of
reciprocating parts.
Consider the piston shown in two different positions in Figure 4.43, and
imagine the crankshaft with a polar mass moment of inertia, I
rot
, to be
non-rotating but executing small torsional oscillations.
In first case there is no motion for the piston, with small oscillations of
the crank and hence the equivalent inertia of the piston is zero.
Figure4.43 Equivalent mass moment of inertia of piston
Dr. R. Tiwari ([email protected]) 99
Whereas in second figure, the piston has practically the same
acceleration as that of the crank pin and the equivalent inertia is m
rec
r
2
, where m
rec
is the mass of the reciprocating parts and r is the
crank throw or radius. Hence, the total polar mass moment of inertia
varies from I
rot
to I
rot
+ m
rec
r
2
, when the crankshaft is rotating.
The inertia of connecting rod can be taken care of by considering
dynamic equivalent two masses one at piston & other at crank pin. It
will contribute to both I
rot
and m
rec
by small amount.
We consider as an approximation the system to have an average
inertia given by
where r is crank radius.
1
2
2
rot rec
I I m r = +
Dr. R. Tiwari ([email protected]) 100
Torsional stiffness of shafts connecting the rotors
In determining the torsional stiffness of the shafts connecting the rotors,
the main difficulty arises from the crank webs.
Considering an idealisation of a crank throw into an ordinary shaft having
the same flexibility as shown in Figure 4.44. Through this idealisation is
physically possible, but the calculations involved are extremely difficult.
This is because the crank webs are subjected to bending and the crank
pin to twisting, if the main shaft is subjected to twisting.
Moreover the beam formulae, if used will not very accurate, because of
short stubs involved rather than long beams usually considered.
Further torques applied at the free end also give rise to sidewise
displacement, which is prevented in the machine.
For high-speed lightweight engines, the crank webs are no more
rectangular blocks and application of the theory becomes extremely
difficult.
Dr. R. Tiwari ([email protected]) 101
Because of this uncertainties in analytical calculation to estimate the
torsional stiffness of crank throws, several experiments have been
carried out on a number of crank shafts of large slow speed engines,
In general the procedure that is applied to reduce the reciprocating
machine system to a mathematical model, is to use a basic diameter,
which corresponding to the journal diameter of the crankshaft. The
torsional stiffness are all calculated based on the basic diameter,
irrespective of their actual diameter.
For the end rotors (i.e. the generator rotor) compute the stiffness of the
shaft from the coupling up to the point of rigidity.
In case where one part of the system is connected to the other part
through gears, or other transmission units, it is convenient to reduce all
the inertias and stiffness to one reference speed.
Fig 44 Equivalent length of a crank
Dr. R. Tiwari ([email protected]) 102
Torsional vibration in reciprocating machinery
Finite element methods, which is described in subsequent sections,
could be used for more accurate analysis of such complex system.
Once the mathematical model is developed, it can be used in illustrating
the critical speed calculations, forced vibration response, and the
coefficient of cyclic irregularity
Torsional oscillation in the crankshaft and in the shafting of driven
machinery is vibration phenomenon of practical importance in the design
of reciprocating engines.
The average torque delivered by a cylinder in a reciprocating machine,
is a small fraction of the maximum torque, which occurs during the firing
period. Even though the torque is periodic the fact that it fluctuates so
violently within the period, constitutes one of the inherent disadvantages
of a reciprocating machine, from the dynamics point of view, as
compared with a turbine where the torque is practically uniform.
It is possible to express the torque by a reciprocating engine into its
harmonic components of several orders of the engine speed, and these
harmonic components can be excite the engine driven installations into
forced torsional vibrations.
Dr. R. Tiwari ([email protected]) 103
It is a commonly known fact that failures can occur in reciprocating
machine installations, when the running speed of the engine is at or
near a dominant torsional critical speed of the system.
High dynamic stresses can occur in the main shafting of such engine
installations and to avoid these conditions, it is essential that the
torsional vibration characteristics of the entire installation be analysed
before the unit is put into operation.
Any analysis of torsional vibration characteristics of reciprocating
machinery should finally predict the maximum dynamic stresses or
torque developed in the shafting & couplings of the system, as
accurately as possible, so that they can be compared with the
permissible values, to check the safety of installation.
Dr. R. Tiwari ([email protected]) 104
Thank you Thank you Thank you Thank you