Contents

1 Dynamics of Rotating Machines - Finite Element Modelling, Simulation,

Analysis and Experiments 4
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Mechanical Models of Rotating Machines . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Mathematical Model of Rotating Machine Elements – Kinematics . . . . . . . . . 5
1.3.1 Inertial and Moving Referential Frames and Transformation Matrices . . 5
1.3.2 Absolute Angular Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Mathematical Model of the Rotating Machine Elements – Equations of Motion . 7
1.4.1 Rigid Discs – Gears, Impellers, etc. . . . . . . . . . . . . . . . . . . . . . . 7
1.4.2 Flexible Shaft Element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4.3 Bearings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.4.4 Excitation Forces – Unbalance Mass on the Rigid Disc . . . . . . . . . . . 18
1.4.5 Excitation Forces – Distributed Unbalance along the Shaft Length . . . . 19
1.5 Global Model and Equation of Motion . . . . . . . . . . . . . . . . . . . . . . . . 20
1.6 Example of Rotor-Bearing System Modelling . . . . . . . . . . . . . . . . . . . . 22
1.6.1 Rotor-Bearing Modelling using a MatLab Program . . . . . . . . . . . . 23
1.6.2 Natural Frequencies and Mode Shapes - MatLab Program Results . . . . 30
1.6.3 Validation of Rotor-Bearing System Modelling using Theoretical and Ex-
perimental Natural Frequencies . . . . . . . . . . . . . . . . . . . . . . . . 37
1.6.4 Validation of Rotor-Bearing System Modelling – Experimental Transient
Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
1.6.5 Validation of Rotor-Bearing System Modelling – Theoretical Transient
Analysis and MatLab Program . . . . . . . . . . . . . . . . . . . . . . . . 42
1.6.6 Validation of Rotor-Bearing System Modelling – Theoretical Transient
Analysis and MatLab Program Results in Time and Frequency Domains . 49
1.6.7 Campbell Diagram – MatLab Program and Theoretical Results . . . . . 50
1.7 Basic Phenomenology of Rotor-Bearing Systems . . . . . . . . . . . . . . . . . . 57
1.7.1 Steady-State Response due to Unbalance Excitation - Forward and Back-
ward Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
1.7.2 Cross-Coupling Stiffness and Rotor-Bearing Stability Analysis . . . . . . . 58
1.8 Stability Analysis of Rotor-Journal Bearing Systems . . . . . . . . . . . . . . . . 61
1.9 Further Bibliography related to Journal Bearing Properties – Research and De-
velopment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
1.9.1 On the Adjusting of the Dynamic Coefficients of Tilting-Pad Journal Bear-
ings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
1.9.2 Tilting-Pad Journal Bearings with Electronic Radial Oil Injection . . . . 72
1.10 Exercises using MatLab Program . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
1.11 Project 3 – Rotor-Journal Bearing Stability - Lateral Dynamics . . . . . . . . . . 78

2 Diagnosis of Rotating Machinery Malfunctions & Condition Monitoring 82

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
2.2 Measurement Procedures and Locations . . . . . . . . . . . . . . . . . . . . . . . 82
2.3 Diagnostic Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
2.3.1 Time-Domain or Waveform Analysis . . . . . . . . . . . . . . . . . . . . . 82
2.3.2 Orbital Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
2.3.3 Spectrum Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
2.3.4 Cepstrum Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
2.4 Theoretical and Experimental Example . . . . . . . . . . . . . . . . . . . . . . . 85

2
 2.4.1 Description of the Test Facilities . . . . . . . . . . . . . . . . . . . . . . . 85
2.4.2 Mechanical Model of the Rotor Test Rig . . . . . . . . . . . . . . . . . . . 86
2.4.3 Theoretical Natural Frequencies and Mode Shapes . . . . . . . . . . . . . 87
2.4.4 Experimental Natural Frequencies of the Rotating Components . . . . . . 88
2.4.5 Measuring the Damping Factor (ξ) and the Log Dec (β) of the Rotor Kit. 88
2.4.6 Experimental Natural Frequencies taking into account the Foundation . . 89
2.4.7 Experimental Forward and Backward Orbits . . . . . . . . . . . . . . . . 90
2.4.8 Detecting Experimentally Super-Harmonic Vibrations . . . . . . . . . . . 92
2.5 Identification of Malfunctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
2.6 Condition Monitoring – Steam Turbines1 . . . . . . . . . . . . . . . . . . . . . . 94
2.7 Research & Development – Bearings with Electronic Oil Injection . . . . . . . . . 96

3
1 Dynamics of Rotating Machines - Finite Element Modelling,
Simulation, Analysis and Experiments

1.1 Introduction
The aim of this chapter is to introduce the basic phenomenology related to vibration and stability
in rotating machines. The Finite Element Method is presented in order to achieve mathematical
models for representing rotor-bearing lateral vibration. The dynamics of rigid and flexible rotat-
ing machines is illustrated. Comparison between theoretical and experimental results, obtained
with help of small test rigs, elucidate the most important topics related to rotor dynamics.

1.2 Mechanical Models of Rotating Machines

The basic components of a rotating machines can be seen in figure 1 and listed above:

Figure 1: Basic components of a rotating machine – Discs, shaft, blades, bearing.

• Rigid Discs

• Shaft

4
 • Bearings

• Blades

In the manuscript emphasis will be given to rigid discs, flexible shafts and journal bearings as
you can see in figure 2.

(a) turbine (b) mechanical model

Figure 2: Basic components of a rotating machine – rigid disc, flexible shaft, bearing; (a) steam
turbine and (b) mechanical model composed of flexible shaft and rigid disc.

1.3 Mathematical Model of Rotating Machine Elements – Kinematics

1.3.1 Inertial and Moving Referential Frames and Transformation Matrices
After three consecutive rotations it is possible to define three moving reference frames and three
transformation matrices, which allows to transform the representation of the vectors (velocity,
acceleration, force, moment etc.) from one frame to another without troubles. The three
consecutive rotations are illustrated in figures 3, 4 and 5.

• Transformation matrix TΓ (from the inertial frame I to the moving frame B1 ):

     
 ~i1   ~i  cos Γ sin Γ 0
~j1 = TΓ ~j where TΓ =  − sin Γ cos Γ 0 
 ~   ~ 
k1 k 0 0 1
• Transformation matrix Tβ (from the moving frame B1 to the moving frame B2 ):
     
 ~i2   ~i1  cos β 0 − sin β
~j = Tβ ~j where Tβ =  0 1 0 
 ~2   ~1 
k2 k1 sin β 0 cos β
• Transformation matrix Tφ (from the moving frame B2 to the moving frame B3 ):
     
 ~i3   ~i2  1 0 0
~j = Tφ ~j where Tφ =  0 cos φ sin φ 
 ~3   ~2 
k3 k2 0 − sin φ cos φ
where I is the inertial reference frame, B1 is the moving reference frame (~i1 , ~j1 , ~k1 ) obtained
after the first rotation around Z. B2 is the second moving reference frame (~i2 , ~j2 , ~k2 ) obtained
after the second rotation around Y1 . B3 is the third moving reference frame (~i3 , ~j3 , ~k3 ) obtained
after the third rotation around X2 .

5
1.3.2 Absolute Angular Velocity

Figure 3: Rotation Γ around Z axis of the inertial frame I.

Figure 4: Rotation β around Y1 axis of the moving frame B1 .

Figure 5: Rotation φ around X2 of the moving reference frame B2 .

• Relative angular velocities of the reference frames:

     
 0   0   φ̇ 
I Γ̇ = 0 B1 β̇ = β̇ B2 φ̇ = 0
     
Γ̇ 0 0

Describing the absolute angular velocity (ω) with help of the moving reference frame B3 , which
is attached to the rotor cross section, and represented by ~i3 , ~j3 , ~k3 , one has:

B3 ω = B3 Γ̇ + B3 β̇ + B3 φ̇ (1)

where

6
 B3 Γ̇ = Tφ · Tβ · TΓ · I Γ̇

B3 β̇ = Tφ · Tβ · B1 β̇

B3 φ̇ = Tφ · B2 φ̇

With help of the transformation matrices, one writes:

  
cos β cos Γ cos β sin Γ − sin β  0 
B3 Γ̇ =  sin φ sin β cos Γ sin φ sin β sin Γ + cos φ cos Γ sin φ cos β  0
 
cos φ sin β cos Γ cos φ sin β cos Γ − sin φ cos Γ cos φ cos β Γ̇
 
− sin β 
B3 Γ̇ = Γ̇ sin φ cos β (2)
 
cos φ cos β

    
cos β 0 − sin β  0   0 
B3 β̇ =  sin φ sin β cos φ sin φ cos β  β̇ ⇒ B3 β̇ = β̇ cos φ (3)
   
cos φ sin β − sin φ cos φ cos β 0 − sin φ

    
1 0 0  φ̇   φ̇ 
B3 φ̇ =  0 cos φ sin φ  0 ⇒ B3 φ̇ = 0 (4)
   
0 − sin φ cos φ 0 0

Rewriting equation (1) as a function of equations (2), (3) e (4), one has:
     
 − sin β   0   φ̇ 
B3 ω = Γ̇ sin φ cos β + β̇ cos φ + 0
     
cos φ cos β − sin φ 0

Using matrix notation:

    
 ωa  1 0 − sin β  φ̇ 
ω =  0 cos φ sin φ cos β  β̇
 b   
ωc 0 − sin φ cos φ cos β Γ̇

Setting a constant rotational speed and neglecting torsional vibrations, the spin angle φ can be
written as Ω t where Ω is the constant angular velocity of the rotor.

1.4 Mathematical Model of the Rotating Machine Elements – Equations of

Motion
1.4.1 Rigid Discs – Gears, Impellers, etc.
• Kinetic energy of a rigid disc can be written as

 T   
 T     wa  IP 0 0  wa 
1 V̇ mD 0 V̇ 1
Ekin = + wb  0 ID 0  w (5)
2 Ẇ 0 mD Ẇ 2   b 
| {z } w c 0 0 ID w c
(I)
| {z }
(II)

7
where the first term of the equation (5), term (I), is related to the linear motion of a rotating
disc and the second part, term (II), to its angular motion.

Using equation (??) and the second term (II) of (5) can be written as:
 T   
 wa  IP 0 0  wa 
1
(II) 2  wb  0 ID 0  w =

 b 
wc 0 0 ID wc
 T   
−Γ̇ sin β + φ̇ IP 0 0 −Γ̇ sin β + φ̇
1 
 0
 
Γ̇ cos β sin φ + β̇ cos φ ID 0  Γ̇ cos β sin φ + β̇ cos φ (6)
2   
Γ̇ cos β cos φ − β̇ sin φ 0 0 ID Γ̇ cos β cos φ − β̇ sin φ

Solving eq.(6) one can get the kinetic energy related to the angular motions of the rigid disc:
n   
1 2 sin2 β − 2Γ̇φ̇ sin β + φ̇2 + I 2 2 2
(II) 2 IP Γ̇ D Γ̇ cos βsin φ + 2Γ̇β̇ cos φ sin φ cos β+

  o
β̇ 2 cos2 φ + ID Γ̇2 cos2 βcos2 φ − 2Γ̇β̇ cos φ sin φ cos β + β̇ 2 sin2 φ (7)

For small amplitudes of β one can consider sin β = β and cos β = 1 and the second order terms
can be neglected. In this matter eq.(7) becomes

1n    o
(II) IP −2Γ̇φ̇β + φ̇2 + ID Γ̇2 sin2 φ + Γ̇2 cos2 φ + β̇ 2 cos2 φ + β̇ 2 sin2 φ (8)
2

The kinetic energy, equation (5), is rewritten then as:

 T     T   
1 V̇ mD 0 V̇ 1 β̇ ID 0 β̇ 1
Ekin = + − φ̇Γ̇βIP + IP φ̇2
2 Ẇ 0 mD Ẇ 2 Γ̇ 0 ID Γ̇ 2
(9)

1 2  1  1
Ekin = V̇ mD + Ẇ 2 mD + β̇ 2 ID + Γ̇2 ID − φ̇Γ̇βIP + IP φ̇2 (10)
2 2 2

• Lagrange’s equation:
     
∂ ∂Ekin ∂Ekin ∂Epot
− + = Fi (i = 1, ..., n) (11)
∂t ∂ q˙i ∂qi ∂qi

where qi are the generalized coordinates, i = 1, 2, 3, ..., n is the number of degree of freedom.
For the rigid disc, n = 4 :

q1 = V ; q2 = W ; q3 = β ; q4 = Γ

Assuming that the disc is rigid, no potential energy due to disc deformation is stored, Epot = 0.
Combining equations (10) and (11) one achieves the equations of motion to mathematically
describe the rigid disc movements:

8
 • 1st equation (q1 = V ):
 
d
dt V̇ m D = FV ⇒ V̈ mD = FV

• 2nd equation (q2 = W ):

 
d
dt Ẇ m D = FW ⇒ Ẅ mD = FW

• 3rd equation (q3 = β):

 
d
dt β̇ · ID + φ̇Γ̇IP = Fβ ⇒ β̈ID + φ̇Γ̇IP = Fβ

• 4th equation (q4 = Γ):

 
d
dt Γ̇ · ID − φ̇βIP = FΓ ⇒ Γ̈ID − φ̇β̇IP = FΓ

Rewriting the four equations in matrix form, one has:

       
mD 0 0 0 
 V̈ 
 0 0 0 0 
 V̇ 
 
 FV 

     
 0 mD 0 0 
 Ẅ  0 0 0 0  Ẇ FW

 0 − Ω = (12)
0 ID 0   β̈  0 0 0 −IP   β̇ F
 β
 
  
  
     

0 0 0 ID Γ̈ 0 0 IP 0 Γ̇ FΓ

Eq.(12) can be rewritten as:

Md · q¨d − Ω · Gd · q˙d = Qd (13)

where

M d → mass matrix of the disc considering the linear and angular motions;

Gd → gyroscopic matrix of the disc;

q̈ d → acceleration vector;

q̇ d → velocity vector;

Qd → is a vector with forces and moments acting on the disc, including unbalance for
example.

1.4.2 Flexible Shaft Element

Recalling Beam Theory from section ??, the angular deformation of a beam can be written as

∂W ∂V
β=− ; Γ= (14)
∂s ∂s

The minimal coordinates for representing the flexible shaft movements are: q1 , q2 , q3 , q4 , q5 , q6 ,
q7 and q8 . All them are time depending and describe the linear and angular motion of the shaft
element extremities presented in figure 6.

9
Figure 6: Coordinates of the shaft elements: 4 linear deformations and 4 angular deformations.

  linear and angular motion at s=0   linear and angular motion at s=l

 q1 
 
 q5 


q2
 q1 = linear motion in Y 
q6
 q5 = linear motion in Y
→ q2 = linear motion in Z → q6 = linear motion in Z

 q3 
 
 q7 


q4
 q3 = angular motion around Y 
q8
 q7 = angular motion around Y
q4 = angular motion around Z q8 = angular motion around Z

Aiming at describing the movements of a point inside of the shaft element one can introduce to
form functions of a continuous beam (shaft element). Considering firstly the plane XY, where
translation in Y direction and rotation around Z axis can occur, one obtain four different cases
and four different form functions for representing the beam deformations, as can be seen in
figure 7.

One can write four form functions, where the main characteristic of such functions is that they
are a polynomial functions of 3rd order:

• Form function for describing the linear deformations:

ψ(s) = A + Bs + Cs2 + Ds3

• Form function for describing the angular deformations:

∂ψ(s)
θ= ∂s = B + 2Cs + 3Ds2

10
Figure 7: Form functions for representing the flexible shaft deformations in the plane XY.

• 1st Case – Boundary conditions for the shaft element (plane XY ) are:

q1 = 1 ; q4 = 0 ; q5 = 0 ; q8 = 0

Finding the constants A, B, C and D:

s=0 → ψ(0) = q1 = 1 ⇒ ψ(0) = A + B0 + C02 + D03 = 1 → A=1

s=0 → θ(0) = q4 = 0 ⇒ θ(0) = B + 2C0 + 3D02 =0 → B=0
s=l → ψ(l) = q5 = 0 ⇒ ψ(l) = A + Bl + Cl2 + Dl3 =0 → Cl2 + Dl3 =
−1 (∗)
s = l → θ(l) = q8 = 0 ⇒ θ(l) = B + 2Cl + 3Dl2 = 0 → 2Cl + 3Dl2 =
0 (∗∗)

From (∗) and (∗∗) one gets

    
Cl2 + Dl3 = −1 l2 l3 C −1
⇒ =
2Cl + 3Dl2 = 0 2l 3l2 D 0
        
C 1 3l2 −l3 −1 1 −3l2 −3
l2
= = = 2
D 3l4 −2l4 −2l l2 0 l4 2l l3

−3 2
A=1 B=0 C= l2
D= l3

3 2 2 3
ψ1 (s) = 1 − l2
s + l3
s and θ1 (s) = − l62 s + 6 2
l3
s

11
• 2nd Case – Boundary conditions for the shaft element (plane XY ) are:

q1 = 0 ; q4 = 1 ; q5 = 0 ; q8 = 0

Finding the constants A, B, C and D:

s=0 → ψ(0) = q1 = 0 ⇒ ψ(0) = A = 0 → A=0

s=0 → θ(0) = q4 = 1 ⇒ θ(0) = B = 1 → B=1
s=l → ψ(l) = q5 = 0 ⇒ Cl2 + Dl3 = −l
s=l → θ(l) = q8 = 0 ⇒ 2Cl + 3Dl2 = −1

Solving the linear system as it was done for the 1st case, one gets:
−2 1
A=0 B=1 C= l D= l2

ψ2 (s) = s − 2l s2 + 1 3
l2
s and θ2 (s) = 1 − 4l s + 3 2
l2
s

• 3rd Case – Boundary conditions for the shaft element (plane XY ) are:

q1 = 0 ; q4 = 0 ; q5 = 1 ; q8 = 0

Finding the constants A, B, C and D:

s=0 → ψ(0) = q1 = 0 ⇒ ψ(0) = A = 0 → A=0

s=0 → θ(0) = q4 = 0 ⇒ θ(0) = B = 0 → B=0
s=l → ψ(l) = q5 = 1 ⇒ Cl2 + Dl3 =1
s=l → θ(l) = q8 = 0 ⇒ 2Cl + 3Dl2 =0
3 −2
A=0 B=0 C= l2
D= l3

2 3 6 6 2
ψ3 (s) = 3 sl2 − 2 sl3 and θ3 (s) = l2
s − l3
s

12
 • 4th Case – Boundary conditions for the shaft element (plane XY ) are:

q1 = 0 ; q4 = 0 ; q5 = 0 ; q8 = 1

Finding the constants A, B, C and D:

s=0 → ψ(0) = q1 = 0 ⇒ ψ(0) = A = 0 → A=0

s=0 → θ(0) = q4 = 0 ⇒ θ(0) = B = 0 → B=0
s=l → ψ(l) = q5 = 0 ⇒ Cl2 + Dl3 =0
s=l → θ(l) = q8 = 1 ⇒ 2Cl + 3Dl2 =1
−1 1
A=0 B=0 C= l D= l2

−s2 s3
ψ4 (s) = l + l2
and θ4 (s) = − 2l s + 3 2
l2
s

The linear deformation inside of the shaft element can be described as a linear combination of
the four form functions and the movements of the shaft element extremities (nodes or degree of
freedom of the discrete mathematical model):

V (s, t) = ψ1 (s) · q1 (t) + ψ2 (s) · q4 (t) + ψ3 (s) · q5 (t) + ψ4 (s) · q8 (t) (15)

The angular deformation inside of the shaft element can be described as:

∂V (s, t)
Γ(s, t) = = θ1 (s) · q1 (t) + θ2 (s) · q4 (t) + θ3 (s) · q5 (t) + θ4 (s) · q8 (t) (16)
∂s

Making the same analysis for the plane XZ, the linear deformation inside of the shaft element
can be described as:

W (s, t) = ψ1 (s) · q2 (t) − ψ2 (s) · q3 (t) + ψ3 (s) · q6 (t) − ψ4 (s) · ·q7 (t) (17)

Making the same analysis for the plane XZ, the angular deformation inside of the shaft element
can be described as:

∂W (s, t)
β(s, t) = − = −θ1 (s) · q2 (t) + θ2 (s) · q3 (t) − θ3 (s) · q6 (t) + θ4 (s) · q7 (t) (18)
∂s

Rewriting in matrix form the linear and angular deformations of the flexible shaft in both planes,
XY as well as XZ, one achieves:

13
  
V (s, t)
= Ψ(s)qe (t) (19)
W (s, t)

 
β(s, t)
= Θ(s)qe (t) (20)
Γ(s, t)

where
 
ψ1 (s) 0 0 ψ2 (s) ψ3 (s) 0 0 ψ4 (s)
Ψ(s) =
0 ψ1 (s) −ψ2 (s) 0 0 ψ3 (s) −ψ4 (s) 0
   
Θβ (s) 0 −θ1 (s) θ2 (s) 0 0 −θ3 (s) θ4 (s) 0
Θ(s) = =
ΘΓ (s) θ1 (s) 0 0 θ2 (s) θ3 (s) 0 0 θ4 (s)

and
 T
qe (t) = q1 (t) q2 (t) q3 (t) q4 (t) q5 (t) q6 (t) q7 (t) q8 (t)

• Potential energy related to the bending motion of an infinitesimal disc:

 T   
e 1 V ′′ EI 0 V ′′
dEpot = ds (21)
2 W ′′ 0 EI W ′′

∂2
where E is the elasticity modulus of the material, I area moment of inertia, V ′′ = ∂s2
V (s, t)
∂2
and W ′′ = ∂s 2 W (s, t)

• Kinetic energy related to the linear and angular motions of an infinitesimal disc:

 T     T   
1 V̇ µ 0 V̇ 1 1 β̇ Id 0 β̇
e
dEkin = ds+ φ̇2 Ip ds+ ds− φ̇Γ̇βIp ds
2 Ẇ 0 µ Ẇ 2 2 Γ̇ 0 Id Γ̇
(22)
∂V (s,t) ∂W (s,t)
where µ is distributed mass per length, V̇ = ∂t and Ẇ = ∂t .

Using equations (19) and (20) one achieves:

e 1 T
dEpot = EI · qeT · Ψ′′ · Ψ′′ · qe ds (23)
2

1 T 1 1 T T
e
dEkin = µ· q˙e ·ΨT ·Ψ· q˙e ds+ φ̇2 Ip ds+ Id · q˙e ·ΘT ·Θ· q˙e ds − φ̇Ip · q˙e ·ΘTΓ ·Θβ ·qe ds (24)
2 2 2

14
Aiming at obtaining the total energy related to the movements of the flexible shaft one integrates
eq.(23) and (24) with respect to the shaft length l, resulting in the following matrices:

Rl 
MeT = 0 µ · ΨT · Ψ ds 





Rl 
MeR = 0 Id · ΘT · Θ ds 




Rl (25)
Ne = 0 Ip · ΘTΓ · Θβ ds 




Rl 

EI · Ψ′′ T

KeB = 0 · Ψ′′ ds 




The total potential energy related to the bending deformation of the shaft is rewritten as:

e 1
Epot = qeT KeB qe (26)
2

The total kinetic energy related to the linear and angular movements of the shaft is rewritten
as:

1 T 1 T
e
Ekin = q˙e (MeT − MeR ) q˙e + Ip φ̇2 − φ̇q˙e Ne qe (27)
2 2

Using Lagrange’s equation 11, one can get the equation of motion of the flexible shaft element,
while considering the angular velocity φ̇ = Ω constant, and as a function of its extremities
movements:

(MeT + MeR ) · q¨e − Ω · Ge · q˙e + KeB · qe = Qe (28)

where q¨e is the acceleration vector composed of the extremities of the shaft element, q˙e is
the velocity vector composed of the extremities of the shaft element, qe is the displacement
vector composed of the extremities of the shaft element and Qe excitation vector with forces
and moments acting on the extremities of the shaft element.

The matrices given in equation (28) are obtained from the integration of equation (25), and
Ge = (Ne − NeT ). For the case of a shaft element with constant cross section, one achieves the
matrices following:

15
 • Mass matrix of the shaft element (considering the linear motion)
 
156
 0 156 
 
 0
 −22l 4l2 

µl  22l
 0 0 4l2  T
e
MT = 420   MeT = MeT
 54 0 0 13l 156 

 0 54 −13l 0 0 156 
 
 0 13l −3l2 0 0 22l 4l2 
−13l 0 0 2
−3l −22l 0 0 4l2

• Mass matrix of the shaft element (considering the angular motion)

 
36
 0 36 
 
 0
 −3l 4l2 

µr  3l
2  0 0 4l2  T
e
MR = 120l   MeR = MeR
 −36 0 0 −3l 36 

 0 −36 3l 0 0 36 
 
 0 −3l −l 2 0 0 3l 4l2 
3l 0 2
0 −l −3l 0 0 4l2
• Gyroscopic matrix of the shaft element
 
0
 36 0 
 
 −3l 0 0 
 
2µr 2  0
 −3l 4l2 0  T
e
G = 120l   Ge = −Ge
 0 36 −3l 0 0 

 −36 0 0 −3l 36 0 
 
 −3l 0 0 l2 3l 0 0 
0 −3l −l2 0 0 3l 4l2 0
• Stiffness matrix considering the bending motion of the shaft element
 
12
 0 12 
 
 0
 −6l 4l2 

 6l 0 0 4l2  T
KeB = EI 
l3  −12
 KeB = KeB
 0 0 −6l 12 

 0 −12 6l 0 0 12 
 
 0 −6l 2l 2 0 0 6l 4l2 
6l 0 0 2l2 −6l 0 0 4l2

1.4.3 Bearings
• Mechanical Model – In rotor dynamics modelling the bearings are normally represented by
springer and damper as can be seen in figure 8(a). Such mechanical elements are mathematically
quantified by means of spring and damping coefficients, as can be seen in figure 8(b).
• Mathematical Model – The equation of motion for a linear bearing is represented by the
follow equation:

16
Figure 8: Mechanical model of bearings – Representation of bearing by means of springers and
dampers.

−Cb · q˙b − Kb · qb = Qb (29)

where
 
V
qb = is the linear displacement of the center of shaft, where the bearing is mounted,
W
 
KVb V KVb W
Kb = b b is the bearing stiffness matrix,
KW V KW W
 
CVb V CVb W
Cb = b b is the bearing damping matrix and
CW V CW W

Qb → is the resultant force supported by the bearing.

The elements KVb V , KVb W , KW

b b b b b
V and KW W describe the stiffness and CV V , CV W , CW V and
b
CW W the viscous damping. More about how to achieve such coefficients, will be presented in
section 1.8.

17
1.4.4 Excitation Forces – Unbalance Mass on the Rigid Disc

Figure 9: Concentrated unbalance mass and mathematical representation in real and complex
forms.

For being completed with the notes in the classes!

18
1.4.5 Excitation Forces – Distributed Unbalance along the Shaft Length
Considering a shaft element with eccentricity (η(s); ξ(s)), the unbalance force will be given by:

Z l     
e 2 T η(s) −ξ(s)
Q = µ·Ω ·ψ · · cos Ωt + sin Ωt ds = QeC · cos Ωt + QeS · sin Ωt (30)
0 ξ(s) η(s)

For a linear distributed unbalance along the finite shaft element one has:

 s s  s s
η(s) = ηL 1 − + ηR e ξ(s) = ξL 1 − + ξR (31)
l l l l

where (ηL , ξL ) and (ηR , ξR ) represent the eccentricities of the center of mass in s=0 and s=l
respectively.

With help of eq.(30) and eq.(31), one can get Qec and Qes as:
 7 3   −7 3 

 20 ηL l + 20 ηR l 
 
 20 ξL l − 20 ξR l 

 7 3   7 3 
20 ξL l + 20 ξR l 20 ηL l + 20 ηR l

 
 
 


 −1 2 1 2

 
 −1 2 1 2


20 ξL l − 30 ξR l 20 ηL l − 30 ηR l

 
 
 


 1 2 1 2

 
 −1 2 1 2


Qec = µΩ2 20 ηL l + 30 ηR l Qes = µΩ2 20 ξL l − 30 ξR l
3 7 −3 7


 20 ηL l + 20 ηR l 




 20 ξL l − 20 ξR l 


3 7 3 7
20 ξL l + 20 ξR l 20 ηL l + 20 ηR l

 
 
 


 
 
 

1 2 1 2 1 2 1 2
30 ξL l + 20 ξR l 30 ηL l + 20 ηR l

 
 
 


 
 
 

−1 2 1 2 1 2 1 2
30 ηL l − 20 ηR l 30 ξL l + 20 ξR l

19
1.5 Global Model and Equation of Motion
The equation of motion of the global system will be achieved by combining the equation of motion
of discs, flexible shaft elements and bearings. One introduces here the global displacement vector
(q s ) with contain the local vectors q e , q d e q b , where each qis represents the linear and angular
motion of the nodes of the rotor discrete model. Considering 2 linear motions V e W and 2
angular motions β e Γ, each node has 4 degrees of freedom.

Each rigid disc is placed in a particular node and influences directly this node with 4 degrees
of freedom. That is the same for the bearing element, which is directly attached to a node.
The shaft element has a length l, connecting 2 different nodes. It means the shaft element will
influence 2 nodes(s=0) and (s=l), where each extremity represents one node of the mechanical
model. Thus, the shaft element has 8 degrees of freedom.

Considering the example presented in figure (10) the equation of motion can be structured as
shown in figure 11, where the sub-matrices M1 M2 M3 and M4 correspond to the shaft elements
(A, B, C e D) of fig.(10) including the bearing elements (stiffness and damping) and disc (inertia
and gyroscopic effects) located in the nodes 2, 3 and 4. The gyroscopic and damping effects are
represented in the global matrix Ds .

(a) disc-shaft system (b) Mechanical Model

Figure 10: (a) Rotating machine outside of the housing; (b) Mechanical model with 4 finite shaft
elements.

The global equation of motion can be written as:

Ms · q¨s + Ds · q˙s + Ks ·qs = Qs (32)

where Ds = −(ΩGs − Cs ). The global matrices Ms , Ds , Gs e Ks are mounted from the machine
elements (shaft, disc, bearing) shown in figure 10, 11 and 12.

20
Figure 11: Structure of the global matrices – Equation of motion of the disc-shaft-bearing system
with 5 nodes and 20 degrees of freedom.

Figure 12: Structure of the matrix Ds for a flexible rotor with 5 nodes and 20 degrees of freedom.

21
1.6 Example of Rotor-Bearing System Modelling

Figure 13: Mechanical System and Mechanical Model with 13 nodes – Flexible Rotor with 2 Discs
and 2 Bearings: (1) motor speed controller; (2) motor; (3) rolling bearing housing attached to a
flexible support and two acceleration sensors; (4) rigid disc (5) flexible shaft (6) rigid disc (7)
rolling bearing housing attached to a flexible support (8) flexible beams for generating different
stiffness coefficients in the horizontal and vertical directions.

In figure 13 it is possible to see a flexible rotating machine (laboratory prototype) built by a

flexible shaft, two rigid discs and two rolling bearing housings attached to two flexible supports.
The little test rig was designed aiming at visualizing the natural frequencies and modes shapes by
using the human eyes. The first three natural frequencies are under 40 Hz, and the rotor-bearing
system is very flexible.

22
1.6.1 Rotor-Bearing Modelling using a MatLab Program

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MACHINERY DYNAMICS LECTURES (41614) %
% MEK - DEPARTMENT OF MECHANICAL ENGINEERING %
% DTU - TECHNICAL UNIVERSITY OF DENMARK %
% %
% Copenhagen, February 10th, 2001 %
% %
% Ilmar Ferreira Santos %
% %
% ROTATING MACHINES -- NATURAL FREQUENCIES AND MODES %
% %
% EXPERIMENTAL RESULTS %
% 13.0 (horizontal) %
% 14.9 (vertical) %
% 33.6 (horizontal) %
% 43.0 (horizontal) %
% 46.0 (vertical) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINITION OF THE STRUCTURE OF THE MODEL %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NE=12; % number of shaft elements
GL = (NE+1)*4; % number of degree of freedom
ND=2; % number of discs
NM=2; % number of bearings
CD1=4; % node - disc 1
CD2=10; % node - disc 2}
CMM1=1; % node - bearing 1
CMM2=13; % node - bearing 2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CONSTANTS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = 2.0E11; % {elasticity modulus [N/m^2}
RAco = 7800; % {steel density [kg/m^3]}
RAl = 2770; % {aluminum density [kg/m^3]}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OPERATIONAL CONDITIONS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Omega= 0*2*pi; % angular velocity [rad/s]
Omegarpm = Omega*60/2/pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GEOMETRY OF THE ROTATING MACHINE %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%
%(A) DISCS %
%%%%%%%%%%%%%%%

Rd = 6/100; % disc radius [m]

espD = 1.1/100 ; % disc thickness [m]
MasD = pi*Rd^2*espD*RAl; % disc mass [kg]

23
Id = 1/4*MasD*Rd^2+1/12*MasD*espD^2; % transversal mass moment of inertia of the disc [Kgm^2]
Ip = 1/2*MasD*Rd*Rd; % polar mass moment of inertia of the disc [Kgm^2]

%%%%%%%%%%%%%%%
%(B) BEARINGS %
%%%%%%%%%%%%%%%

MasM = 0.40698; % bearing mass [kg](housing + ball bearings)

h=1/1000; % beam thickness [m]
b=28.5/1000; % beam width [m]
Area=b*h; % beam cross section area [m^2]
I=b*h^3/12; % beam moment of inertia of area [m^4]
lr=7.5/100; % beam length [m]
Kty0=2*12*E*I/(lr^3); % equivalent beam flexural stiffness [N/m]
Ktz0=2*E*Area/lr; % equivalent bar stiffness [N/m]
% Bearing 1 - Damping
Dty1 = 0.0 ;
Dtz1 = 0.0 ;
Dry1 = 0.0 ;
Drz1 = 0.0 ;
% Bearing 2 - Damping
Dty2 = 0.0 ;
Dtz2 = 0.0 ;
Dry2 = 0.0 ;
Drz2 = 0.0 ;
% Bearing 1 - Stiffness
Kty1 = Kty0 ;
Ktz1 = Ktz0 ;
Kry1 = 0.0 ;
Krz1 = 0.0 ;
% Bearing 2 - Stiffness
Kty2 = Kty0 ;
Ktz2 = Ktz0 ;
Kry2 = 0.0 ;
Krz2 = 0.0 ;

%%%%%%%%%%%%%%%
%(C) SHAFT %
%%%%%%%%%%%%%%%

ll = 435/1000; % length of shaft elements [m]

Rext = (6/2)/1000; % shaft external radius [m]
Rint = (0/2)/1000; % shaft internal radius [m]

% length of the shaft elements [m]

l(1) = 0.140/3;
l(2) = 0.140/3;
l(3) = 0.140/3;
l(4) = 0.205/6;
l(5) = 0.205/6;
l(6) = 0.205/6;
l(7) = 0.205/6;
l(8) = 0.205/6;
l(9) = 0.205/6;
l(10) = 0.090/3;
l(11) = 0.090/3;
l(12) = 0.090/3;

% external radius of shaft elements [m]

for i=1:NE,
rx(i)=Rext;
end

24
% internal radius of shaft elements [m]
for i=1:NE,
ri(i)=Rint;
end

% density of shaft elements [kg/m]

for i=1:NE,
ro(i) = RAco;
end

% transversal areal of the shaft elements [m^2]}

for i=1:NE,
St(i) = pi*(rx(i)^2-ri(i)^2);
end

% area moment of inertia of the shaft elements [m^4]}

for i=1:NE,
II(i)=pi*((rx(i)+ri(i))/2)^3*(rx(i)-ri(i));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MOUNTING THE GLOBAL MATRICES %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(’MOUNTING THE GLOBAL MATRICES - WAIT!’)

disp(’ ’)

% Defining the global matrices with zero elements

M=zeros(GL);
G=zeros(GL);
K=zeros(GL);

%%%%%%%%%%%%%%%%%%%%%%
% GLOBAL MASS MATRIX %
%%%%%%%%%%%%%%%%%%%%%%

disp(’MOUNTING THE GLOBAL MASS MATRIX - WAIT!’)

disp(’ ’)

% Mass matrices of shaft elements due to linear and angular movements

a=1; b=8;

for n=1:NE,

MteAux= [156 0 0 22*l(n) 54 0 0 -13*l(n)

0 156 -22*l(n) 0 0 54 13*l(n) 0
0 -22*l(n) 4*l(n)^2 0 0 -13*l(n) -3*l(n)^2 0
22*l(n) 0 0 4*l(n)^2 13*l(n) 0 0 -3*l(n)^2
54 0 0 13*l(n) 156 0 0 -22*l(n)
0 54 -13*l(n) 0 0 156 22*l(n) 0
0 13*l(n) -3*l(n)^2 0 0 22*l(n) 4*l(n)^2 0
-13*l(n) 0 0 -3*l(n)^2 -22*l(n) 0 0 4*l(n)^2];

Mte = ((ro(n)*St(n)*l(n))/420)*MteAux;

MreAux= [36 0 0 3*l(n) -36 0 0 3*l(n)

0 36 -3*l(n) 0 0 -36 -3*l(n) 0
0 -3*l(n) 4*l(n)^2 0 0 3*l(n) -l(n)^2 0

25
 3*l(n) 0 0 4*l(n)^2 -3*l(n) 0 0 -l(n)^2
-36 0 0 -3*l(n) 36 0 0 -3*l(n)
0 -36 3*l(n) 0 0 36 3*l(n) 0
0 -3*l(n) -l(n)^2 0 0 3*l(n) 4*l(n)^2 0
3*l(n) 0 0 -l(n)^2 -3*l(n) 0 0 4*l(n)^2];

Mre = ((ro(n)*II(n))/(30*l(n)))*MreAux;

MauxT=Mte+Mre;

for f=a:b,
for g=a:b,
M(f,g)=M(f,g)+MauxT(f-(n-1)*4,g-(n-1)*4);
end
end

a=a+4; b=b+4;

end

% Adding the mass matrices of the disc elements

M((CD1-1)*4+1,(CD1-1)*4+1)=M((CD1-1)*4+1,(CD1-1)*4+1)+MasD;
M((CD1-1)*4+2,(CD1-1)*4+2)=M((CD1-1)*4+2,(CD1-1)*4+2)+MasD;
M((CD1-1)*4+3,(CD1-1)*4+3)=M((CD1-1)*4+3,(CD1-1)*4+3)+Id;
M((CD1-1)*4+4,(CD1-1)*4+4)=M((CD1-1)*4+4,(CD1-1)*4+4)+Id;
M((CD2-1)*4+1,(CD2-1)*4+1)=M((CD2-1)*4+1,(CD2-1)*4+1)+MasD;
M((CD2-1)*4+2,(CD2-1)*4+2)=M((CD2-1)*4+2,(CD2-1)*4+2)+MasD;
M((CD2-1)*4+3,(CD2-1)*4+3)=M((CD2-1)*4+3,(CD2-1)*4+3)+Id;
M((CD2-1)*4+4,(CD2-1)*4+4)=M((CD2-1)*4+4,(CD2-1)*4+4)+Id;

% Adding the mass matrices of the bearing elements

M((CMM1-1)*4+1,(CMM1-1)*4+1)=M((CMM1-1)*4+1,(CMM1-1)*4+1)+MasM;
M((CMM1-1)*4+2,(CMM1-1)*4+2)=M((CMM1-1)*4+2,(CMM1-1)*4+2)+MasM;
M((CMM2-1)*4+1,(CMM2-1)*4+1)=M((CMM2-1)*4+1,(CMM2-1)*4+1)+MasM;
M((CMM2-1)*4+2,(CMM2-1)*4+2)=M((CMM2-1)*4+2,(CMM2-1)*4+2)+MasM;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLOBAL GYROSCOPIC MATRIX %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(’MOUNTING THE GLOBAL GYROSCOPIC MATRIX - WAIT!’)

disp(’ ’)

% Gyroscopic matrix of shaft elements

a=1; b=8;

for n=1:NE,

GeAux=[0 -36 3*l(n) 0 0 36 3*l(n) 0

36 0 0 3*l(n) -36 0 0 3*l(n)
-3*l(n) 0 0 -4*l(n)^2 3*l(n) 0 0 l(n)^2
0 -3*l(n) 4*l(n)^2 0 0 3*l(n) -l(n)^2 0
0 36 -3*l(n) 0 0 -36 -3*l(n) 0
-36 0 0 -3*l(n) 36 0 0 -3*l(n)
-3*l(n) 0 0 l(n)^2 3*l(n) 0 0 -4*l(n)^2
0 -3*l(n) -l(n)^2 0 0 3*l(n) 4*l(n)^2 0 ];

26
 Ge=Omega*(ro(n)*II(n))/(15*l(n))*GeAux;

for f=a:b,
for g=a:b,
G(f,g)=G(f,g)+Ge(f-(n-1)*4,g-(n-1)*4);
end
end

a=a+4; b=b+4;

end

% Adding the gyroscopic matrices of the disc elements

G((CD1-1)*4+3,(CD1-1)*4+4)=G((CD1-1)*4+3,(CD1-1)*4+4)-Omega*Ip;
G((CD1-1)*4+4,(CD1-1)*4+3)=G((CD1-1)*4+4,(CD1-1)*4+3)+Omega*Ip;
G((CD2-1)*4+3,(CD2-1)*4+4)=G((CD2-1)*4+3,(CD2-1)*4+4)-Omega*Ip;
G((CD2-1)*4+4,(CD2-1)*4+3)=G((CD2-1)*4+4,(CD2-1)*4+3)+Omega*Ip;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLOBAL STIFFNESS MATRIX %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(’MOUNTING THE GLOBAL STIFFNESS MATRIX - WAIT!’)

disp(’ ’)

% Stiffness matrix of shaft elements due to bending

a=1; b=8;

for n=1:NE,

KbeAux= [12 0 0 6*l(n) -12 0 0 6*l(n)

0 12 -6*l(n) 0 0 -12 -6*l(n) 0
0 -6*l(n) 4*l(n)^2 0 0 6*l(n) 2*l(n)^2 0
6*l(n) 0 0 4*l(n)^2 -6*l(n) 0 0 2*l(n)^2
-12 0 0 -6*l(n) 12 0 0 -6*l(n)
0 -12 6*l(n) 0 0 12 6*l(n) 0
0 -6*l(n) 2*l(n)^2 0 0 6*l(n) 4*l(n)^2 0
6*l(n) 0 0 2*l(n)^2 -6*l(n) 0 0 4*l(n)^2];

Kbe = ((E*II(n))/(l(n)^3))*KbeAux;

for f=a:b,
for g=a:b,
K(f,g)=K(f,g)+Kbe(f-(n-1)*4,g-(n-1)*4);
end
end

a=a+4; b=b+4;

end

% Adding the stiffness matrices of the bearing elements

K((CMM1-1)*4+1,(CMM1-1)*4+1)=K((CMM1-1)*4+1,(CMM1-1)*4+1)+Ktz1;
K((CMM1-1)*4+2,(CMM1-1)*4+2)=K((CMM1-1)*4+2,(CMM1-1)*4+2)+Kty1;
K((CMM2-1)*4+1,(CMM2-1)*4+1)=K((CMM2-1)*4+1,(CMM2-1)*4+1)+Ktz2;
K((CMM2-1)*4+2,(CMM2-1)*4+2)=K((CMM2-1)*4+2,(CMM2-1)*4+2)+Kty2;

27
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLOBAL MATHEMATICAL MODEL %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Mglob=[ G M
M zeros(size(M,1))];
Kglob=[ K zeros(size(M,1))
zeros(size(M,1)) -M ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODAL ANALYSIS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(’CALCULATING NATURAL FREQUENCIES AND MODE SHAPES - WAIT!’)

disp(’ ’)

% Calculating Eigenvectors and Eigenvalues

[U,lambda]=eig(-Kglob,Mglob);
[lam,p]=sort(diag(lambda));
U=U(:,p);

% Number of divisions in time for plotting the mode shapes

nn=99;

N=size(U,1);
maximo=num2str((N-2)/2);
ModoVirt=N;

ModoVirt=input([’ Enter the number of the mode shape to be plotted, ...

zero to esc, highest mode ’,maximo,’: ’]);

% For visualizing the mode shapes:

ModoReal=2*ModoVirt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOOP TO PLOT THE MODES SHAPES %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while ModoReal>0,

% Natural frequencies
wn=imag(lam(ModoReal));
ttotal=8/abs(wn);
dt=ttotal/nn;
t=0:dt:ttotal;

% Defining v and w real e imaginary for each node

y=1:4:GL;
z=2:4:GL;

for i=1:(NE+1),
vr(i)=real(U(y(i),ModoReal));
vi(i)=imag(U(y(i),ModoReal));
wr(i)=real(U(z(i),ModoReal));
wi(i)=imag(U(z(i),ModoReal));
end

28
 % Calculation the modal displacements v and w

for i=1:(NE+1),
v(i,:)=vr(i)*cos(wn*t)+vi(i)*sin(wn*t);
w(i,:)=wr(i)*cos(wn*t)+wi(i)*sin(wn*t);
end

Zero=diag(zeros(length(t)))’;
Um=diag(eye(length(t)))’;

for i=1:(NE+1)
pos(i,:)=Zero+(i-1)*Um;
end

clf
hold on
for cont=1:NE+1,
plot3(pos(cont,:),w(cont,:),v(cont,:),’k’,’LineWidth’,2.5);
end

nm=num2str(ModoVirt);
fn=num2str(abs(wn/2/pi));
dfi=num2str(Omegarpm);
title([’Angular Velocity: ’,dfi,’ rpm Mode: ’,nm,’ Nat. Freq.: ’,fn,’ Hz’],’FontSize’,14)
view(-25,20);
grid on;
ModoVirt=input([’ Enter the number of the mode shape to be plotted, ...
zero to esc, highest mode ’,maximo,’: ’]);
ModoReal=2*ModoVirt;
figure(ModoVirt)

end

29
1.6.2 Natural Frequencies and Mode Shapes - MatLab Program Results

Angular Velocity: 0 rpm Mode: 1 Nat. Freq.: 13.8512 Hz

0.8

0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1
1
0.5
0
−3 12
x 10 10
−0.5 8
6
4
−1 2
0

Angular Velocity: 1200 rpm Mode: 1 Nat. Freq.: 13.5682 Hz

−4
x 10

−1

−2

−3
1
0.5
0
−3 12
x 10 10
−0.5 8
6
4
−1 2
0

Figure 14: First mode shape and first natural frequency calculated with the MatLab program, at
zero and 1200 rpm (20Hz).

30
 Angular Velocity: 0 rpm Mode: 2 Nat. Freq.: 14.7288 Hz
−4
x 10

−2

−4

−6
1
0.5
0
12
10
−0.5 8
6
4
−1 2
0

Angular Velocity: 1200 rpm Mode: 2 Nat. Freq.: 14.9916 Hz

−4
x 10

−2

−4

−6
4
2
0
−4 12
x 10 10
−2 8
6
4
−4 2
0

Figure 15: Second mode shape and second natural frequency calculated with the MatLab program,
at zero and 1200 rpm (20Hz).

31
 Angular Velocity: 0 rpm Mode: 3 Nat. Freq.: 33.4345 Hz

0.8

0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1
4
2
0
−4 12
x 10 10
−2 8
6
4
−4 2
0

Angular Velocity: 1200 rpm Mode: 3 Nat. Freq.: 33.4233 Hz

−6
x 10

2.5

1.5

0.5

−0.5

−1

−1.5

−2

−2.5
4
2
0
−4 12
x 10 10
−2 8
6
4
−4 2
0

Figure 16: Third mode shape and third natural frequency calculated with the MatLab program,
at zero and 1200 rpm (20Hz).

32
 Angular Velocity: 0 rpm Mode: 4 Nat. Freq.: 42.256 Hz

0.8

0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1
2
1
0
−4 12
x 10 10
−1 8
6
4
−2 2
0

Angular Velocity: 1200 rpm Mode: 4 Nat. Freq.: 42.1675 Hz

−5
x 10

2.5

1.5

0.5

−0.5

−1

−1.5

−2

−2.5
2
1
0
−4 12
x 10 10
−1 8
6
4
−2 2
0

Figure 17: Fourth mode shape and fourth natural frequency calculated with the MatLab program,
at zero and 1200 rpm (20Hz).

33
 Angular Velocity: 0 rpm Mode: 5 Nat. Freq.: 47.1143 Hz
−4
x 10

0.8

0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1
1
0.5
0
12
10
−0.5 8
6
4
−1 2
0

Angular Velocity: 1200 rpm Mode: 5 Nat. Freq.: 47.1231 Hz

−4
x 10

0.8

0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1
2

0
−5 12
x 10 10
8
6
4
−2 2
0

Figure 18: Fifth mode shape and fifth natural frequency calculated with the MatLab program, at
zero and 1200 rpm (20Hz).

34
 Angular Velocity: 0 rpm Mode: 6 Nat. Freq.: 56.7282 Hz

0.8

0.6

0.4

0.2

−0.2

−0.4

−0.6

−0.8

−1
1
0.5
0
−4 12
x 10 10
−0.5 8
6
4
−1 2
0

Angular Velocity: 1200 rpm Mode: 6 Nat. Freq.: 56.6822 Hz

−6
x 10

−2

−4

−6
1
0.5
0
−4 12
x 10 10
−0.5 8
6
4
−1 2
0

Figure 19: Sixth mode shape and sixth natural frequency calculated with the MatLab program, at
zero and 1200 rpm (20Hz).

35
 Angular Velocity: 0 rpm Mode: 7 Nat. Freq.: 131.2026 Hz
−5
x 10

−1

−2

−3
1
0.5
0
12
10
−0.5 8
6
4
−1 2
0

Angular Velocity: 1200 rpm Mode: 7 Nat. Freq.: 117.5409 Hz

−5
x 10

2.5

1.5

0.5

−0.5

−1

−1.5

−2

−2.5
2

0
−5 12
x 10 10
8
6
4
−2 2
0

Figure 20: Seventh mode shape and sixth natural frequency calculated with the MatLab program,
at zero and 1200 rpm (20Hz).

36
1.6.3 Validation of Rotor-Bearing System Modelling using Theoretical and Exper-
imental Natural Frequencies

MODE SHAPE PLANE Nat. Frequency Nat. Frequency Error

(theoretical) (experimental)
[Hz] [Hz] %
1 horizontal 13.8 13.0 6.1
2 vertical 14.7 14.9 1.3
3 horizontal 33.4 33.6 0.6
4 horizontal 42.3 42.3 0.0
5 vertical 47.1 46.0 2.1
6 horizontal 56.7 55.4 2.1
7 vertical 131. (out of range) –

Table 1: Validation of Rotor-Bearing System Modelling using information about the theoretical
and experimental natural frequencies in the range of frequencies between 0 and 60 Hz, when the
rotor angular velocity is zero.

Figure 21: Acceleration sensors 9 and 10 mounted on the bearing housing, with the goal of
measuring the movements of rotor-bearing system in the horizontal and vertical directions re-
spectively.

37
1.6.4 Validation of Rotor-Bearing System Modelling – Experimental Transient
Analysis

Figure 22: Mechanical System and Mechanical Model with 5 nodes – Flexible Rotor with 2 Discs
and 2 Bearings.

38
 Horizontal Acceleration of Bearing 2 (Excitation on Disc 2)
1
Acc. [m/s2] − Y direction

0.5

−0.5

−1
0 1 2 3 4 5 6 7 8 9 10
Time [s]

40
FFT(Acc.) − Y direction

33.6 Hz

30
13.0 Hz

20

10 42.3 Hz
55.4 Hz

0
0 10 20 30 40 50 60
Frequency [Hz]

Figure 23: Experimental Analysis – Transient lateral vibration of the rotor-bearing system, when
the rotor has no angular velocity, and it is excited by a horizontal perturbation (shock) on the
disc 2 (node 4) and the horizontal acceleration of the bearing 2 (node 5) is measured.

39
 Vertical Acceleration of Disc 2 (Excitation on Disc 2)
1
Acc. [m/s2] − Y direction

0.5

−0.5

−1
0 1 2 3 4 5 6 7 8 9 10
Time [s]

30
FFT(Acc.) − Z direction

46.0 Hz
25

20

15
14.9 Hz
10

0
0 10 20 30 40 50 60
Frequency [Hz]

Figure 24: Experimental Analysis – Transient lateral vibration of the rotor-bearing system, when
the rotor has no angular velocity, and it is excited by a vertical perturbation (shock) on the disc
2 (node 4) and the horizontal acceleration of the disc 2 (node 4) is measured.

40
 FFT of Horizontal and Vertical Acceleration Signals
40
(H) − Horizontal Plane
33.6 Hz (H)
(V) − Vertical Plane
35

30

13.0 Hz (H) 46.0 Hz (V)

25
FFT(Acc.)

20

15
14.9 Hz (V)

10
42.3 Hz (H)
55.4 Hz (H)
5

0
0 10 20 30 40 50 60
Frequency [Hz]

Figure 25: Experimental Analysis - Natural frequencies of the rotor-bearing system: 13.0 Hz and
33.6 Hz in the horizontal plane, and 14.9 Hz and 46 Hz in the vertical place.

41