THE UNIVERSITY OF DANANG

UNIVERSITY OF SCIENCE AND TECHNOLOGY

FACULTY OF MECHANICAL ENGINEERING

TECHNICAL REPORT
ON VIBRATION ANALYSIS

Supervisor : Assoc. Prof. Dr. DANG PHUOC VINH

Students : VO MINH NHAT
NGUYEN TANG PHI LONG
DOAN CONG TIN
TRUONG QUOC DAI
Class : 20CDT
Group : 20.04

Da Nang, 9/2024
VIBRATION ANALYSIS TECHNIQUES 20.04

Student's name : Vo Minh Student's name : Nguyen

Nhat Tang Phi Long
Students ID : 101200234 Students ID : 101200231
Class : 20CDT1 Class : 20CDT1

Student's name : Doan

Student's name : Truong
Cong Tin
Quoc Dai
Students ID : 101200293
Students ID : 101200261
Class : 20CDT2
Class : 20CDT2

Group Information

2
VIBRATION ANALYSIS TECHNIQUES 20.04

Table of Contents

CHAPTER 3: LINEAR VIBRATION OF 1 D.O.F SYSTEM...........................................

3.1 Free Undamped vibration of 1 D.O.F system...................................................4
3.1.1 Assignment:.........................................................................................................
3.1.2 Mathematical formula.........................................................................................
3.1.3 Code Matlab by analytical solution:..................................................................
3.1.4 Assignment:.........................................................................................................
3.1.5Mathematical formula:........................................................................................
3.1.6 Code Matlab by analytical solution:..................................................................
3.2. Free damped system...........................................................................................14
3.2.1 Assignment:.......................................................................................................
3.2.2 Mathematical formula:.....................................................................................
3.2.3 Code Matlab.....................................................................................................
3.2.4 Assignment:.......................................................................................................
3.3. Forced vibration of 1 D.O.F system...................................................................22
3.3.1 Assignment:.......................................................................................................
3.3.2 Mathematical formula.......................................................................................
3.3.3 Code matlab......................................................................................................
3.3.4 Topic.................................................................................................................
3.3.5 Assignment 1.....................................................................................................
3.3.6 Assignment 2.....................................................................................................
3.3.7 Asignment 3......................................................................................................
CHAPTER 4: CASE STUDY...........................................................................................
CASE STUDY 1: BEAM..........................................................................................48
4.1 Analytical solution...............................................................................................
4.2 Code matlab :....................................................................................................
4.3 FEM: Numerical model.......................................................................................
4.4 Code matlab :....................................................................................................

3
VIBRATION ANALYSIS TECHNIQUES 20.04

CHAPTER 3: LINEAR VIBRATION OF

1 D.O.F SYSTEM
3.1 Free Undamped vibration of 1 D.O.F system
System model:

3.1.1 Assignment:
1. Calculate the natural frequency ωn of the system
2. Evaluate the free undamped vibrations of the system (displacement and
velocity)
a. Analytical solution
b. Numerical integration (Matlab ODE solver)
3. Play with the coefficients to see the effect on the solution (m , k)
3.1.2 Mathematical formula

Solution:
System: m ẍ + kx = 0
Impose solution form: x(t) = Cest = 0
Characteristic equation: ms2 + k)Cest = 0
Impose non-trivial solution: ms2 + k = 0

√
 S1,2 = ± i k
m
; ωn =
√ k
m

Solution with unknown amplitude/phase: x(t) = C1es1t + C2es2t

Initial conditions: t = 0
t = 0: x(0) = X0  X0 = C1 + C2
t = 0: v(0) = V0  V0 = S1.C1 +S2.C2
x(t) = C1es1t + C2es2t
v(t) = S1.C1es1t + S2.C2est2
ODE:
m ẍ + kx = 0
y (1) = x (t)
y (2) = ẋ (t)

4
VIBRATION ANALYSIS TECHNIQUES 20.04

k k
ẏ (2) = ẍ = - x = 0 y(2) + - y(1)
m m
ẏ (1) = ẋ = y(2) + 0 y(1)

[ ] [
ẏ (2)
 ẏ (1) = 1
0 −k /m
0 ] [ y (2)]
x y (1)  ẏ = A.y
Integrate numerically: ODE ( ẏ = A.y, t, y[0])
Complete time domain solution: y(t)

3.1.3 Code Matlab by analytical solution:

clc
clear
close all
m = 10;
k = 50;
dt = 0.001;
t = 0 : dt : 10;

s1 = sqrt(k/m)*1i;
s2 = -sqrt(k/m)*1i;

a = [1 1 ; s1 s2];
z = [0.01 ;0.2];
c = a\z;
x = c(1)*exp(s1*t) + c(2)*exp(s2*t);
v = s1*c(1)*exp(s1*t) + s2*c(2)*exp(s2*t);
figure('name' , 'Respond of 1 DOF system');
subplot(2,1,1); hold on, title('Free Undamping system'), plot(t,x,'r'), xlabel('Time [s]'),
ylabel('Displacement [m]');
grid on;
subplot(2,1,2);hold on, plot(t,v); xlabel('Time [s]'), ylabel('Velocity [m/s]');
grid on;
hold off;

 Code Matlab by Numerical integration (Matlab ODE solver):

clc
clear
close all

m = 10;
k = 50;
x0 = 0.01;
5
VIBRATION ANALYSIS TECHNIQUES 20.04

v0 = 0.2;

tspan = [0 10];

initial_conditions = [x0; v0];

ode_system = @(t, y) [y(2); -(k/m) * y(1)];

[t, y] = ode45(ode_system, tspan, initial_conditions);

wn = sqrt(k/m);
fn = wn / (2 * pi);

figure('name' , 'Respond of 1 DOF system');

subplot(2,1,1); hold on;

plot(t, y(:,1), 'r', 'LineWidth', 2);
title(['Free Undamping system, \omega_n = ', num2str(fn), ' Hz']);
xlabel('Time [s]');
ylabel('Displacement [m]'); grid on;

subplot(2,1,2);
plot(t, y(:,2), 'b', 'LineWidth', 2);
xlabel('Time [s]');
ylabel('Velocity [m/s]'); grid on;

 Result:
,

Figure 3.1.1 Displacement and Velocity over time for a free undamping system

 Observation:
6
VIBRATION ANALYSIS TECHNIQUES 20.04

- Displacement Graph: The x-axis shows time(s), and the y-axis shows
displacement(m). The waveform indicates harmonic oscillation with constant
amplitude, showing that there is no damping force acting on the system. The
system oscillates around an equilibrium position with a consistent period.
- Velocity Graph: The x-axis also shows time in seconds, and the y-axis shows
velocity in meters per second (m/s). This waveform also exhibits harmonic
oscillation, reflecting how velocity changes during the oscillation of the object.
In summary, both graphs illustrate a free undamped oscillatory system, where
displacement and velocity vary according to harmonic motion.

3.1.4 Assignment:
Using ODE solver, compare responds of system (displacement & velocity)
i. 4 different values of mass (10 – 30 – 50 – 70 kg)
ii. 4 different values of stiffness (50 – 100 – 200 – 300 N/m)

3.1.5Mathematical formula:

x(t) = X1eiωnt + X2e-iωnt = Acos(ωnt + φ c)

v(t) = - ωnAsin(ω nt + φ c)
3.1.6 Code Matlab by analytical solution:
 4 different values of stiffness (50 – 100 – 200 – 300 N/m)

clc
clear
close all

m = 30;
x0 = 1;
v0 = -2;
dt = 0.00001;
t = 0:dt:10;
figure('name','Đáp ứng với các giá trị k');
for k = [50 100 200 300]
wn = sqrt(k/m);
7
VIBRATION ANALYSIS TECHNIQUES 20.04

disp(['Natural frequency: ',num2str(wn/2/pi),' Hz']);

s1 = 1i*sqrt(k/m);
s2 = -1i*sqrt(k/m);
A = [1 1; s1 s2];
b = [x0; v0];
c = A\b;
c1 = c(1);
c2 = c(2);

x = c1*exp(s1*t) + c2*exp(s2*t);
v = s1*c1*exp(s1*t) + s2*c2*exp(s2*t);

subplot(2,1,1);
plot(t, x, 'DisplayName', ['k = ', num2str(k), ' N/m , wn = ',
num2str(wn/2/pi),' Hz'] );
hold on;
xlabel('Timer [s]','FontWeight', 'bold', 'Color', 'blue');
ylabel('Displacement [m]','FontWeight', 'bold', 'Color', 'blue');
title('Effect of stiffness coefficient', 'Color', 'blue');

subplot(2,1,2);
plot(t,v); hold on;
xlabel('Timer [s]','FontWeight', 'bold', 'Color', 'blue');
ylabel('Velocity [m/s]','FontWeight', 'bold', 'Color', 'blue');
end

subplot(2,1,1);
legend( 'Location', 'northeast');
subplot(211), ylim([-2 2]), grid on

subplot(212), ylim([-4 4]), grid on;

set(gca,'FontSize', 10)

 Code Matlab by Numerical integration (Matlab ODE solver):

clc
clear
close all

x0 = 1;
v0 = -2;
tspan = [0 10];
8
VIBRATION ANALYSIS TECHNIQUES 20.04

stiffnesses = [50, 100, 200, 300];

m = 30;

figure('name','Đáp ứng với các giá trị k');

for i = 1:length(stiffnesses)
k = stiffnesses(i);

ode_system = @(t, y) [y(2); -(k/m) * y(1)];

[t, y] = ode45(ode_system, tspan, [x0; v0]);

subplot(2,1,1);
hold on;
plot(t, y(:,1), 'LineWidth', 2);

subplot(2,1,2);
hold on;
plot(t, y(:,2), 'LineWidth', 2);
end

subplot(2,1,1);
title('Effect of Stiffnesses coefficient');
xlabel('Time [s]');
ylabel('Displacement [m]');
legend('k = 50 N/m', 'k = 100 N/m', 'k = 200 N/m', 'k = 300 N/m'); grid on;

subplot(2,1,2);
xlabel('Time [s]');
ylabel('Velocity [m/s]');
legend('k = 50 N/m', 'k = 100 N/m', 'k = 200 N/m', 'k = 300 N/m'); grid on;

 Result:

9
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 3.1.2 Displacement and Velocity over time for a free undamped system with 4 different
value of stiffness (50 – 100 – 200 – 300 N/m)

 Observation:
This graph illustrates the effect of the stiffness coefficient.
Graph Displacement:
- The x-axis represents Time(s) and the y-axis represents
displacement(m).
- There are four curves, each corresponding to different stiffness values:
 k = 50 N/m (blue curve) with ω n = 0.20547 Hz.
 k = 100 N/m (orange curve) with ω n = 0.29058 Hz.
 k = 200 N/m (yellow curve) with ω n = 0.41049 Hz.
 k = 300 N/m (purple curve) with ωn = 0.50329 Hz.
Observation: As the stiffness coefficient k increases, the amplitude of
oscillation decreases and the frequency increases. This shows that the
stiffer the system, the faster it oscillates, but with smaller amplitudes.
Velocity Graph:
- The x-axis again represents Timer(s), while the y-axis represents
velocity(m/s). The curves show the velocity variation corresponding to
different stiffness values.
- Observation: The peak velocity also depends on the stiffness
coefficient. Higher stiffness leads to earlier peak velocities with larger
magnitudes.
10
 VIBRATION ANALYSIS TECHNIQUES 20.04

 Code Matlab by analytical solutioon:

 4 different values of mass (20 – 30 – 50 – 100 kg)

clc
clear
close all

k = 50;
x0 = 1;
v0 = -2;
dt = 0.001;
t = 0:dt:20;

figure('name','Đáp ứng với các giá trị m');

for m = [20 30 50 100]

wn = sqrt(k/m);
disp(['Natural frequency: ',num2str(wn/2/pi),' Hz']);

s1 = 1i*sqrt(k/m);
s2 = -1i*sqrt(k/m);
A = [1 1; s1 s2];
b = [x0; v0];
c = A\b;
c1 = c(1);
c2 = c(2);

x = c1*exp(s1*t) + c2*exp(s2*t);
v = s1*c1*exp(s1*t) + s2*c2*exp(s2*t);

subplot(2,1,1);
plot(t, x, 'DisplayName', ['m = ', num2str(m), ' N/m ; wn = ', num2str(wn/2/pi),'
Hz']);
hold on;
xlabel('Timer [s]','FontWeight', 'bold', 'Color', 'blue');
ylabel('Displacement [m]','FontWeight', 'bold', 'Color', 'blue');
title('Effect of mass coefficient','FontWeight', 'bold', 'Color', 'blue');

subplot(2,1,2);
plot(t,v);
hold on;
xlabel('Timer [s]','FontWeight', 'bold', 'Color', 'blue');
ylabel('Velocity [m/s]','FontWeight', 'bold', 'Color', 'blue');
end
subplot(2,1,1);
legend('Location', 'northeast');
subplot(211), ylim([-3 3]), grid on;

subplot(212), ylim([-3 3]), grid on;

set(gca,'FontSize', 10)
11
VIBRATION ANALYSIS TECHNIQUES 20.04

 Code Matlab by Numerical integration (Matlab ODE solver):

clc
clear
close all

x0 = 1;
v0 = -2;
tspan = [0 10];

mass = [20, 30, 50, 100];

k = 100;

figure('name','Đáp ứng với các giá trị m');

for i = 1:length(masses)
m = mass(i);

ode_system = @(t, y) [y(2); -(k/m) * y(1)];

[t, y] = ode45(ode_system, tspan, [x0; v0]);

subplot(2,1,1);
hold on;
plot(t, y(:,1), 'LineWidth', 2);

subplot(2,1,2);
hold on;
plot(t, y(:,2), 'LineWidth', 2);
end

subplot(2,1,1);
title('Effect of mass coefficient');
xlabel('Time [s]');
ylabel('Displacement [m]');
legend('m = 20 kg', 'm = 30 kg', 'm = 50 kg', 'm = 100 kg');
grid on;

subplot(2,1,2);
xlabel('Time [s]');
ylabel('Velocity [m/s]');
legend('m = 20 kg', 'm = 30 kg', 'm = 50 kg', 'm = 100 kg');
grid on;

 Result:

12
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 3.1.3 Displacement and Velocity over time for a free undamped system with 4 different
value of mass (20 – 30 – 50 – 100 kg)

 Observation:
This graph illustrates the effect of mass, with two plots showing
displacement and velocity over time.
- Displacement Graph:
- The x-axis represents time(s), and the y-axis represents
displacement(m).
- There are four curves, each corresponding to different mass values: m =
20 kg (purple curve), m = 30 kg (red curve), m = 50 kg (yellow curve),
m = 100 kg (blue curve).
- Observation: As the mass m increases, the system oscillates with a
lower frequency and higher amplitude. This indicates that systems with
larger mass oscillate more slowly but exhibit greater displacements.
- Velocity Graph:
- The x-axis represents time(s), and the y-axis represents velocity (m/s).
- The velocity curves correspond to the same mass values as the
displacement graph.
- Observation: As the mass increases, the velocity amplitude decreases.
Heavier masses exhibit smaller velocities, as the system takes longer to
respond and move during each oscillation cycle.

13
VIBRATION ANALYSIS TECHNIQUES 20.04

3.2. Free damped system

m 100 Kg
m ẍ ( t )+ c ẋ ( t ) + kx ( t )=0
x (t=0 )=x 0 k 10000 N/m
ẋ (t=0)=v 0
c 50 Ns/m
x₀ 0.1 m

v₀ -2.0 m/s

3.2.1 Assignment:

1. Calculate the characteristics of the system for different damping coefficients

c:
a) Natural frequency
b) Damping ratio
2. Evaluate the free damped vibrations (displacement and velocity) for different
damping coefficients c using:
a) Analytical solution
b) Numerical integration (Matlab ODE solver)
3. Play with the coefficients to see the effect on the solution
3.2.2 Mathematical formula:
Equation of damped harmonic oscillation
m × ẍ ( t ) +c × ẋ ( t ) +k × x ( t )=0
st
Set x (t )=C ×e
 ẋ (t )=S × C ×e st , ẍ (t )=S 2 ×C ×e st
Substituting into the equation, we get:
( m × S2 +c × S+ k ) × C × e st =0
−c+ √ c 2−4 mk
 S1 ,2 =
2m

√
c
Set w n= k và ξ=
m 2 √ mk
 s1 ,2 =−ξ × wn ± w n × √ξ 2−1
The solution of the characteristic equation:
Case 1: ξ = 1 (critical damping):
S1=S 2=−w n
s ×t s×t
x (t )=C 1 × e + C2 ×e
Case 2: ξ > 1 (overdamping)

14
VIBRATION ANALYSIS TECHNIQUES 20.04

s 1 ×t s 2 ×t
x (t )=C 1 × e +C 2 × e
Case 3: ξ < 1 (underdamping):
s1 ,2 =−ξ × wn ± i× wn × √ 1−ξ 2
Set w d =wn × √1−ξ2
 x (t )=C 1 × e s 1 ×t +C 2 × e s 2 ×t
 x (t )=C 1 × e−ξ ×w + i× w +C 2 × e−ξ× w +i × w
n d n d

To find C1 and C2, we will use two initial conditions:

 Initial displacement: x(0) = x₀

 Initial velocity: v(0) = v₀

Substituting t = 0 into the general equation, we get:
x ( 0 )=C 1+ C2=x 0
Substituting t = 0 into the equation ẋ (t ) , we get:
v ( 0 )=C1 × s1 +C 2 × s 2=v 0
Solving the system of equations, we find C 1 and C 2 in terms of s1, s2, x 0 and v 0.
Case 1: ξ = 1 (critical damping):
C 1=x 0
C 1 × s+ C2=v 0
Case 2: ξ ≠ 1
C 1+C 2=x 0
C 1 × s 1+C 2 × s 2=v 0

3.2.3 Code Matlab

Code Matlab by analytical solution

clc
clear
close all
m = 100;
k = 10000;
c = 50;
dt = 0.001;
t = 0 : dt : 20;
wn = sqrt(k/m);
z= [0.01 ; -2];
h = c/2/sqrt(m*k);
if h == 1
s = - wn;
a = [1 0 ; s 1];
ak = a\z;
x = ak(1).*exp(s*t) + ak(2).*t.*exp(s*t);
v = s.*ak(1).*exp(s*t) + ak(2).*(s*t+1).*exp(s*t);
end
if h > 1
s1 = -h*wn + wn*sqrt(h*h -1);
s2 = -h*wn - wn*sqrt(h*h -1);
b = [1 1 ; s1 s2];
15
VIBRATION ANALYSIS TECHNIQUES 20.04

k = b\z;
x = k(1)*exp(s1*t) + k(2)*exp(s2*t);
v = s1*k(1)*exp(s1*t) + s2*k(2)*exp(s2*t);
end
if h < 1
s1 = -h*wn + 1i*wn*sqrt(1- h*h);
s2 = -h*wn - 1i*wn*sqrt(1-h*h);
e = [1 1 ; s1 s2];
g = e\z;
x = g(1)*exp(s1*t) + g(2)*exp(s2*t);
v = s1*g(1)*exp(s1*t) + s2*g(2)*exp(s2*t);
end
figure('name' , 'Respond of 1 DOF system');
subplot(2,1,1); hold on, title('Effect of damping
coefficient'), plot(t,x,'r'), xlabel('Time [s]'),
ylabel('Displacement [m]');
grid on;
subplot(2,1,2);hold on, plot(t,v); xlabel('Time [s]'),
ylabel('Velocity [m/s]');
grid on;
hold off;

Result :

16
VIBRATION ANALYSIS TECHNIQUES 20.04

Effect of damping coefficient

0.2

Displacement [m]
0.1

-0.1

-0.2
0 2 4 6 8 10 12 14 16 18 20
Time [s]

1
Velocity [m/s]

-1

-2
0 2 4 6 8 10 12 14 16 18 20
Time [s]

Figure 3.2.1: Effect of Damping Coefficient on Displacement and Velocity(Using

analytical solution)

 Code Matlab by Numerical integration (Matlab ODE solver):

clc
clear
close all
m = 100;
k = 10000;
c = 50;
dt = 0.001;
t = 0 : dt : 20;
global A;

x0 = 0.01;
v0 = -2;
z0 = [x0; v0];
A = [0 1; -k/m -c/m];
[t, z] = ode45(@mysys, t, z0);
x = z(:,1);
v = z(:,2);

figure('name' , 'Respond of 1 DOF system');

17
VIBRATION ANALYSIS TECHNIQUES 20.04

subplot(2,1,1); hold on, title('Effect of damping

coefficient'), plot(t,x,'r'), xlabel('Time [s]'),
ylabel('Displacement [m]');
grid on;
subplot(2,1,2);hold on, plot(t,v); xlabel('Time [s]'),
ylabel('Velocity [m/s]');
grid on;
hold off;
function dz = mysys(t, z)
global A;
dz = A * z;
end

Effect of damping coefficient

0.2
Displacement [m]

0.1

-0.1

-0.2
0 2 4 6 8 10 12 14 16 18 20
Time [s]

1
Velocity [m/s]

-1

-2
0 2 4 6 8 10 12 14 16 18 20
Time [s]
Figure 3.2.2: Effect of Damping Coefficient on Displacement and Velocity (Using ODE
solver)
Observation:
The graph clearly demonstrates the impact of the damping coefficient on
the system's oscillation. As the damping coefficient increases, the oscillation
decays more rapidly, as evidenced by the decreasing amplitude of both
displacement and velocity over time. This indicates that the damping coefficient
plays a crucial role in controlling the stability of the system.
 Damping coefficient: A higher damping coefficient leads to faster decay
of oscillations.
18
VIBRATION ANALYSIS TECHNIQUES 20.04

 Oscillation: Both displacement and velocity gradually decrease to zero.

3.2.4 Assignment:

Using ODE solver, compare responds of system (displacement and velocity)

with 3 different values of damping

clc
clear
close all
m = 100;
k = 10000;
dt = 0.001;
t = 0 : dt : 10;
global A;
x0 = 0.01;
v0 = -2;
z0 = [x0; v0];
figure;
for c = [100 2000 5000]

A = [0 1; -k/m -c/m];
[t, z] = ode45(@mysys, t, z0);

x = z(:,1);
v = z(:,2);

subplot(2, 1, 1);
hold on;
plot(t, x, 'DisplayName', ['c = ', num2str(c), ' Ns/m']);
subplot(2, 1, 2);
hold on;
plot(t, v,
'DisplayName', ['c = ', num2str(c), ' Ns/m']);
end
subplot(211), title('Effect of damping coefficient'),
xlabel('Time [s]'), ylabel('Dispalcement [m]');
subplot(211), grid on;
legend('show');
subplot(212), xlabel('Time [s]'), ylabel('Velocity [m/s]');
subplot(212), grid on;
legend('show');
hold off;
function dz = mysys(t, z)
global A;
dz = A * z;
end

 Result :

19
VIBRATION ANALYSIS TECHNIQUES 20.04

Effect of damping coefficient

0.2
c = 100 Ns/m

Dispalcement [m] 0.1 c = 2000 Ns/m

c = 5000 Ns/m

-0.1

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

2
c = 100 Ns/m
1 c = 2000 Ns/m
Velocity [m/s]

c = 5000 Ns/m

-1

-2
0 1 2 3 4 5 6 7 8 9 10
Time [s]
Figure 3.2.3 Displacement and Velocity vs. Time for Varying Damping different values c
(Using ODE solver)

Observation:
The provided graph effectively illustrates the influence of the damping
coefficient on the displacement and velocity of a system over time. Each curve
represents a different value of the damping coefficient (c = 100 Ns/m, 2000
Ns/m, and 5000 Ns/m).
 Effect of damping coefficient:

o Higher damping coefficient: Oscillations decay more rapidly,

indicating a quicker return to the equilibrium position and reduced
vibration.
o Lower damping coefficient: Oscillations persist for a longer
duration with larger amplitudes, potentially leading to prolonged
vibration and affecting system performance.

 Displacement:

20
VIBRATION ANALYSIS TECHNIQUES 20.04

o As the damping coefficient increases, the maximum amplitude of

oscillation decreases.
o The time taken for the system to reach a steady state is reduced.
 Velocity:
o Similar to displacement, the maximum amplitude of velocity also
decreases with increasing damping coefficient.
o The time taken for the velocity to diminish to zero is reduced.
Conclusion:
From the graph, it can be concluded that the damping coefficient plays a
crucial role in controlling the system's oscillations. Selecting an appropriate
damping coefficient value can minimize vibration, enhance stability, and
improve overall system performance.

21
VIBRATION ANALYSIS TECHNIQUES 20.04

3.3. Forced vibration of 1 D.O.F system

Forced vibration
m 100 kg
m ẍ+ c ẋ +kx=F
k 1e5 N/m
x (t=0 )=x 0
c 100 Ns/m
ẋ (t=0 )=v 0
x0 0.02 m
F=F 0 cos ⁡(ωf t +f f )
v0 -3.0 m/s
F0 10 kN
ωf 7 rad/s
3.3.1 Assignment: ff π /8 Rad

Evaluate the forced damped vibrations (velocity and

displacement) for different damping coefficients c using:
a) Analytical solution
b) Numerical integration (Matlab ODE solver)

3.3.2 Mathematical formula

System Differential Equation (ODE)
- The general form of the second-order differential equation for a forced
harmonic oscillator is:
m ẍ+ c ẋ +kx=F 0 cos ( ω f + f )
- This equation models the oscillation of a system with mass m , damping
coefficient c , and stiffness k , under the influence of an external harmonic
force F 0 cos ( ω f + f )
Separation into Homogeneous and Particular Solution
The solution to this differential equation can be split into two parts:
- Homogeneous part (where there is no external force):
m ẍ+ c ẋ +kx=0
- Particular part (related to the external force):
m ẍ+ c ẋ +kx=F 0 cos ( ω f + f )
Solving the Homogeneous Equation
- The homogeneous solution is obtained by solving:
s1 t s2 t
x g ( t )=X 1 e + X 1 e

22
VIBRATION ANALYSIS TECHNIQUES 20.04

- The corresponding characteristic equation is:

2
ms 1, 2 + c s1 , 2+ k=0
Solving this characteristic equation yields the roots s1 and s2, which are
expressed as:
s1 ,2 =−ωn ± ωn √ ξ 2−1
Where ξ is the damping ratio, and ω n is the natural frequency of the system.
Solving the Particular Solution
- Forced Oscillation Equation:

Forced oscillation equation:

m ẍ+ c ẋ +kx=F 0 cos ( ω f + f )
To solve this equation, the particular solution x p is represented as a
harmonic oscillation:
x p= Acos( ωf t )
- Representing the Harmonic Function in the Complex Domain:
The forcing term F 0 cos ( ω f t+ f ) is represented in the complex form using
Euler's formula:
F 0 i(ω + f ) F 0 −i(ω +f )
F 0 cos ( ω f t+ f )=
e + e f f

2 2
Similarly, the particular solution x p is also expressed as the sum of two
complex solutions:
x p=x p + x p 1 2
i ωf t −i ωf t
x p= Acos ( ω f t+ f ) =Y 1 e +Y 2 e
where Y 1 and Y 2 are the coefficients to be determined.
- Solving the Complex Solution:
Next, the solution is split into two parts, corresponding to two complex
solutions x p 1 and x p 2:
F 0 j ω jf
m ẍ p +c ẋ p + k x p = e e f

1 1
2 1

F 0 − j ω − jf
m ẍ p + c ẋ p + k x p = e e f

2 2
2 2

Solving these equations will give the complex solutions Y 1 and Y 2.

- Finding Y 1 and Y 2:
The solutions for Y 1 and Y 2 are:
jf
F0 e
Y 1= 2
2(−m ω f +ic ω f +k )
− jf
F0 e
Y 2= 2
2(−m ω f −ic ω f + k)
These are the complex solutions representing the oscillation of the
system under the influence of the forcing term F 0 cos ( ω f t+ f )
23
VIBRATION ANALYSIS TECHNIQUES 20.04

3.3.3 Code matlab

 Code matlab by Analytical solution :

clc
clear
close all
m = 100;
k = 1e5;
n =100;
Fo = 10000;
wf = 7;
phi = pi/8 ;

dt = 0.001;
t = 0 : dt : 30;
wn = sqrt(k/m);
z= [0.02 ; -3];
c = 100;
h = c/2/sqrt(m*k);
Y1 = ((Fo/2)*exp(1i*phi))/(-m*wf*wf+1i*c*wf+ k);
Y2 = ((Fo/2)*exp(-1i*phi))/(-m*wf*wf-1i*c*wf+ k);
if h == 1
s = - wn;
a = [1 0 ; s 1];
d = a\z;
x = d(1).*exp(s*t) + d(2).*t.*exp(s*t);
v = s.*d(1).*exp(s*t) + d(2).*(s*t+1).*exp(s*t);

end
if h > 1
s1 = -h*wn + wn*sqrt(h*h -1);
s2 = -h*wn - wn*sqrt(h*h -1);
b = [1 1 ; s1 s2];
d = b\z;
x = d(1)*exp(s1*t) + d(2)*exp(s2*t);
v = s1*d(1)*exp(s1*t) + s2*d(2)*exp(s2*t);
end
if h < 1
s1 = -h*wn + 1i*wn*sqrt(1- h*h);
s2 = -h*wn - 1i*wn*sqrt(1-h*h);
b = [1 1 ; s1 s2];
d = b\z;
x = d(1)*exp(s1*t) + d(2)*exp(s2*t);
v = s1*d(1)*exp(s1*t) + s2*d(2)*exp(s2*t);
end

xf = Y1*exp(1i*wf*t) + Y2*exp(-1i*wf*t);
vf = 1i*wf*(Y1*exp(1i*wf*t) - Y2*exp(-1i*wf*t));
xn = x + xf ;
24
VIBRATION ANALYSIS TECHNIQUES 20.04

vn = v + vf ;
subplot(2, 1, 1);
hold on;
plot(t, xn,'r', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);

subplot(2, 1, 2);

hold on;
plot(t, vn,'r', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);

subplot(211), title_str = sprintf('Displacement: wn = %.4f Hz, Damping ratio (h) = %.6f', wn,
h);
title(title_str);
xlabel('Time [s]','Color', 'blu','FontWeight','bold');
ylabel('Dispalcement [m]','Color', 'blu','FontWeight','bold');
subplot(211), grid on;

subplot(212);
title_str = sprintf('Velocity: wn = %.4f Hz, Damping ratio (h) = %.6f', wn, h);
title(title_str);
xlabel('Time [s]','Color', 'blu','FontWeight','bold');
ylabel('Velocity [m/s]','Color', 'blu','FontWeight','bold');
subplot(212), grid on;

hold off;

 Code matlab Numerical integration (Matlab ODE solver)

clc;
clear;
close all;
m = 100;
k = 1e5;
Fo = 10000;
wf = 7;
phi = pi/8;
dt = 0.001;
tspan = 0:dt:30;
z0 = [0.02; -3];
c = 100;
wn = sqrt(k/m);
h = c/(2*sqrt(m*k));
odefun = @(t, z) [z(2); (Fo * cos(wf * t + phi) - c * z(2) - k * z(1)) / m];

[t, z] = ode45(odefun, tspan, z0);

xn = z(:, 1);
vn = z(:, 2);
25
VIBRATION ANALYSIS TECHNIQUES 20.04

figure;

subplot(2, 1, 1);
plot(t, xn, 'r', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);
title_str = sprintf('Displacement: wn = %.4f Hz, Damping ratio (h) = %.6f', wn, h);
title(title_str);
xlabel('Time [s]', 'Color', 'blue', 'FontWeight', 'bold');
ylabel('Displacement [m]', 'Color', 'blue', 'FontWeight', 'bold');
grid on;

subplot(2, 1, 2);
plot(t, vn, 'r', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);
title_str = sprintf('Velocity: wn = %.4f Hz, Damping ratio (h) = %.6f', wn, h);
title(title_str);
xlabel('Time [s]', 'Color', 'blue', 'FontWeight', 'bold');
ylabel('Velocity [m/s]', 'Color', 'blue', 'FontWeight', 'bold');
grid on;

Result :

Figure 3.3.1: The graph corresponding to the displacement and velocity of the oscillating
system under the influence of harmonic force.

Observation:
 Shows the displacement of the system over time. The system oscillates
with large amplitudes, which gradually decrease, before stabilizing with a
nearly constant amplitude.
 Shows the velocity of the system over time. Similar to the displacement
plot, the initial velocity oscillations are large but diminish over time.
26
VIBRATION ANALYSIS TECHNIQUES 20.04

General Observations:
This system has a low damping ratio (h < 1), indicating underdamped
motion. Initially, the system oscillates freely but loses energy due to the
damping force. After some time, the oscillations are fully driven by the external
harmonic force with an angular frequency of wf = 7 rad/s.
The plots clearly show the relationship between displacement and
velocity. When displacement reaches its maximum, velocity crosses zero, and
vice versa—typical behavior of harmonic oscillation.

27
VIBRATION ANALYSIS TECHNIQUES 20.04

3.3.4 Topic

m ẍ+ c ẋ +kx=F
x (t=0 )=x 0
ẋ (t=0 )=v 0
F=F 0 cos ⁡(ωf t +f f )

c=500 Ns /m ; F 0=10000 N
c=50 Ns /m ; F0 =50000 N
c=50 Ns /m ; F0 =10000 N

Compare responds of system (displacement and velocity) with 2 different values

c, F0

 Code matlab

clc;
clear;
close all;

m = 100;
k = 1e5;
wf = 5;
dt = 0.001;
t = 0:dt:30;
wn = sqrt(k/m);
phi = pi/8;
params = [500, 10000;
50, 50000;
50, 10000];

for i = 1:size(params, 1)
c = params(i, 1);
Fo = params(i, 2);
h = c / (2 * sqrt(m * k));

Y1 = (Fo/2 * exp(1i * phi)) / (-m * wf^2 + 1i * c * wf + k);

Y2 = (Fo/2 * exp(-1i * phi)) / (-m * wf^2 - 1i * c * wf + k);
x0 = 0.02 - Y1 - Y2;
v0 = -3 -1i*wf*Y1 + 1i*wf*Y2;
z= [ x0 ; v0];
if h == 1
28
VIBRATION ANALYSIS TECHNIQUES 20.04

s = -wn;
a = [1 0; s 1];
d = a \ z;
x = d(1) * exp(s * t) + d(2) * t .* exp(s * t);
v = s * d(1) * exp(s * t) + d(2) * (s * t + 1) .* exp(s * t);
elseif h > 1
s1 = -h * wn + wn * sqrt(h^2 - 1);
s2 = -h * wn - wn * sqrt(h^2 - 1);
b = [1 1; s1 s2];
d = b \ z;
x = d(1) * exp(s1 * t) + d(2) * exp(s2 * t);
v = s1 * d(1) * exp(s1 * t) + s2 * d(2) * exp(s2 * t);
elseif h < 1
s1 = -h*wn + 1i*wn*sqrt(1- h*h);
s2 = -h*wn - 1i*wn*sqrt(1-h*h);
b = [1 1 ; s1 s2];
d = b\z;
x = d(1)*exp(s1*t) + d(2)*exp(s2*t);
v = s1*d(1)*exp(s1*t) + s2*d(2)*exp(s2*t);

end
xf = Y1 * exp(1i * wf * t) + Y2 * exp(-1i * wf * t);
vf = 1i * wf * (Y1 * exp(1i * wf * t) - Y2 * exp(-1i * wf * t));
xn = x + xf;
vn = v + vf;

subplot(2, 1, 1);
hold on;
plot(t, xn, 'DisplayName', sprintf('c = %d Ns/m, h = %.2f, Fo = %d N, wf =%.1d rad/s', c, h,
Fo, wf));
subplot(2, 1, 2);
hold on;
plot(t, vn, 'DisplayName', sprintf('c = %d Ns/m, h = %.2f, Fo = %d N, wf =%.2 rad/s', c, h, Fo,
wf));
end
subplot(2, 1, 1);
title('Effect of Damping Coefficient','Color','blu');
xlabel('Time [s]','Color', 'blu','FontWeight','bold');
ylabel('Dispalcement [m]','Color', 'blu','FontWeight','bold');
ylim([-1 1]);
grid on;
legend('show');

subplot(2, 1, 2);
xlabel('Time [s]','Color', 'blu','FontWeight','bold');
ylabel('Velocity [m/s]','Color', 'blu','FontWeight','bold');
ylim([-20 15]);
grid on;hold off;
 Results.

29
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 3.3.2: The graph compare responds of system (displacement and velocity) with 2
different values c, F 0

Observation:
Displacement Graph:
 Overview: This plot shows the displacement of the system over time for
three different damping and force combinations.
 Initial oscillations: All three cases start with high-amplitude oscillations,
which then dampen as time progresses. The initial part of the graph shows
noticeable differences in how quickly the oscillations die out, influenced by
the damping coefficient. Higher damping (c = 500 Ns/m) causes the
displacement to reduce more rapidly compared to lower damping (c = 50
Ns/m).
 Steady-state behavior: After around 10 seconds, the system reaches a
nearly steady oscillation in all cases, primarily governed by the external
force's frequency (wf = 5 rad/s), with only slight variations in amplitude
depending on the force's magnitude. Higher force amplitudes, like Fo =
50000 N, result in slightly larger steady-state oscillations.
Velocity Graph:
 Overview: This plot represents the velocity of the system over time for the
same three cases.
 Initial velocity response: In the first few seconds, the velocity oscillations
are very pronounced, especially for lower damping (c = 50 Ns/m). Systems
with higher damping lose velocity faster.
 Steady-state behavior: As time passes (after 10 seconds), the velocity
oscillations settle into regular patterns. The systems with lower damping

30
VIBRATION ANALYSIS TECHNIQUES 20.04

show slightly larger oscillations, though all cases eventually stabilize around
a periodic velocity pattern determined by the external force.
General Observations:
 Damping effect: Systems with higher damping (c = 500 Ns/m) reduce
oscillations faster, both in terms of displacement and velocity, showing a
more rapid decay in initial free vibrations. However, the external force
eventually dominates the response, leading to a steady-state behavior similar
across all cases.
 Force effect: Larger external forces (Fo = 50000 N) produce higher steady-
state displacements and velocities, but they do not significantly alter the time
it takes for the system to reach this steady-state compared to lower forces
(Fo = 10000 N).
 Oscillation pattern: All systems oscillate at the external force frequency
(wf = 5 rad/s), as seen in both displacement and velocity plots after the
transient period.
Summary:
 Higher damping leads to quicker suppression of free oscillations, while the
external force primarily controls the long-term behavior.
 Larger external forces lead to larger steady-state oscillations.
 After an initial transient phase (around 10 seconds), all cases exhibit
periodic oscillations governed by the external force frequency, with slight
differences in amplitude.

3.3.5 Assignment 1

m ẍ+ c ẋ +kx=F
31
VIBRATION ANALYSIS TECHNIQUES 20.04

x (t=0 )=x 0
ẋ (t=0 )=v 0
F=F 0 cos ⁡(ωf t +f f )

Using ODE solver, compare responds of system (displacement and velocity)

with 3 different values F 0,ω f and c
 Code matlab Numerical integration (Matlab ODE solver)
clc;
clear;
close all;

m = 100;
k = 1e5;
dt = 0.001;
tspan = 0:dt:20;
params = [300, 5, 100;
2000, 10, 500;
5000, 10, 0];
results_xn = cell(size(params, 1), 1);
results_vn = cell(size(params, 1), 1);

for i = 1:size(params, 1)
Fo = params(i, 1);
wf = params(i, 2);
c = params(i, 3);
wn = sqrt(k/m);
h = c / (2 * sqrt(m * k));

odefun = @(t, z) [z(2); (Fo * cos(wf * t) - c * z(2) - k * z(1)) / m];

z0 = [0.02; -3];
[t, z] = ode45(odefun, tspan, z0);
xn = z(:, 1);

vn = z(:, 2);
results_xn{i} = xn;
results_vn{i} = vn;

32
VIBRATION ANALYSIS TECHNIQUES 20.04

subplot(2, 1, 1);
hold on;
plot(t, xn, 'DisplayName', sprintf('Fo = %d N, wf = %.1f rad/s, c = %d Ns/m', Fo, wf, c));

subplot(2, 1, 2);
hold on;
plot(t, vn, 'DisplayName', sprintf('Fo = %d N, wf = %.1f rad/s, c = %d Ns/m', Fo, wf, c));
end

subplot(2, 1, 1);
title('Effect of Damping Coefficient', 'Color', 'blue');
xlabel('Time [s]', 'Color', 'blue', 'FontWeight', 'bold');
ylabel('Displacement [m]', 'Color', 'blue', 'FontWeight', 'bold');
ylim([-0.2 0.2]);
grid on;
legend('show');

subplot(2, 1, 2);
xlabel('Time [s]', 'Color', 'blue', 'FontWeight', 'bold');
ylabel('Velocity [m/s]', 'Color', 'blue', 'FontWeight', 'bold');
ylim([-4 4]);
grid on;
hold off;

Results:

Figure 3.3.3: The graph compare responds of system (displacement and velocity) with 3
different values c, F 0 and ω f (Using ODE solver)

Observation:

33
VIBRATION ANALYSIS TECHNIQUES 20.04

The graph displays the simulation results of two quantities:

 The upper graph: Displacement over time.
 The lower graph: Velocity over time.
The graph lines represent three different cases with varying parameters: forcing
force (Fo), forcing frequency (wf), and damping coefficient (c), as indicated in
the legend.
Effects of Forcing Force and Forcing Frequency:
 Fo (Forcing Force): A higher force results in a greater amplitude of
oscillation. In the graph, the curves corresponding to Fo = 5000 N and Fo
= 2000 N show larger amplitudes compared to Fo = 300 N. This is evident
in both the displacement and velocity graphs.
 wf (Forcing Frequency): The forcing frequency also affects the
waveform of the oscillation. In the case of wf = 5 rad/s, the oscillation has
a slower rhythm compared to wf = 10 rad/s. This can be observed through
the number of wave cycles over the plotted time interval.
Effects of Damping Coefficient (c):
 Case of c = 0 (no damping): The oscillation is sustained evenly, with no
reduction in amplitude over time. The oscillatory wave has a large
amplitude and does not diminish.
 Cases of c = 100 Ns/m and c = 500 Ns/m (with damping): Increased
damping leads to a reduction in oscillation amplitude over time.
Particularly at larger values like c = 500 Ns/m, the amplitude decreases
rapidly in the initial seconds and then stabilizes. The decrease in
amplitude during the initial period (0 to 5 seconds) is prominent in both
the displacement and velocity graphs.
 Characteristics of Free Oscillation and Forced Oscillation: Initially,
the system's oscillation is primarily free oscillation, with a rapid decrease
in amplitude due to damping (in cases with damping). After some time,
the forced oscillation from the external force becomes more pronounced,
and the system oscillates with a nearly constant amplitude at the forcing
frequency. This is particularly evident in the case with c = 500 Ns/m: the
beginning of the graph shows a significant decrease in amplitude, but
subsequently, the oscillation stabilizes at the forcing frequency.
Summary:
 The amplitude of oscillation is heavily dependent on the forcing force
(Fo), forcing frequency (wf), and damping coefficient (c).
 A larger forcing force and the absence of damping lead to oscillations
with large amplitudes that do not diminish over time.
34
VIBRATION ANALYSIS TECHNIQUES 20.04

 A large damping coefficient results in a rapid

decrease in oscillation amplitude during the initial m 50 kg
phase, followed by a stable forced oscillation.
k 1e6 N/m
c 50 Ns/m
x0 0.01 m
v0 -2.0 m/s
F0 1e4 kN
f 20 Hz
ff π /8 Rad

3.3.6 Assignment 2

m ẍ+ c ẋ +kx=F
x (t=0 )=x 0
ẋ (t=0 )=v 0
F=F 0 cos ⁡(ωf t +f f )

Assignment:
1. Evaluate the forced damped vibrations of system
a) Analytical solution
b) ODE solver
2. Calculation of Fast Fourier Transform of vibration
and force applied to system

35
VIBRATION ANALYSIS TECHNIQUES 20.04

 Code matlab by Analytical solution :

clc
clear
close all
m = 50;
k = 1e6;
Fo = 1e4;
f = 20;
wf = 2*pi*f;
phi = pi/8 ;

dt = 0.001;
t = 0 : dt : 10;
wn = sqrt(k/m);
z= [0.01 ; -2];
c = 50;
h = c/2/sqrt(m*k);
Y1 = ((Fo/2)*exp(1i*phi))/(-m*wf*wf+1i*c*wf+ k);
Y2 = ((Fo/2)*exp(-1i*phi))/(-m*wf*wf-1i*c*wf+ k);
if h == 1
s = - wn;
a = [1 0 ; s 1];
d = a\z;
x = d(1).*exp(s*t) + d(2).*t.*exp(s*t);
v = s.*d(1).*exp(s*t) + d(2).*(s*t+1).*exp(s*t);

end
if h > 1
s1 = -h*wn + wn*sqrt(h*h -1);
s2 = -h*wn - wn*sqrt(h*h -1);
b = [1 1 ; s1 s2];
d = b\z;
x = d(1)*exp(s1*t) + d(2)*exp(s2*t);
v = s1*d(1)*exp(s1*t) + s2*d(2)*exp(s2*t);
end
if h < 1
s1 = -h*wn + 1i*wn*sqrt(1- h*h);
s2 = -h*wn - 1i*wn*sqrt(1-h*h);
b = [1 1 ; s1 s2];
d = b\z;
x = d(1)*exp(s1*t) + d(2)*exp(s2*t);
v = s1*d(1)*exp(s1*t) + s2*d(2)*exp(s2*t);
end

xf = Y1*exp(1i*wf*t) + Y2*exp(-1i*wf*t);
vf = 1i*wf*(Y1*exp(1i*wf*t) - Y2*exp(-1i*wf*t));
xn = x + xf ;
vn = v + vf ;
subplot(2, 1, 1);
hold on;
plot(t, xn,'b', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);
36
VIBRATION ANALYSIS TECHNIQUES 20.04

subplot(2, 1, 2);
hold on;
plot(t, vn,'b', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);

subplot(211), title_str = sprintf('Displacement: wn = %.4f Hz, Damping ratio (h) = %.6f',

wn/2/pi, h);
title(title_str);
xlabel('Time [s]','Color', 'blu','FontWeight','bold');
ylabel('Dispalcement [m]','Color', 'blu','FontWeight','bold');
ylim([-0.08 0.08]);
grid on;

subplot(212);
title_str = sprintf('Velocity: wn = %.4f Hz, Damping ratio (h) = %.6f', wn/2/pi, h);
title(title_str);
xlabel('Time [s]','Color', 'blu','FontWeight','bold');
ylabel('Velocity [m/s]','Color', 'blu','FontWeight','bold');
ylim([-15 10]); grid on; hold off;

 Code matlab Numerical integration (Matlab ODE solver)

clc;
clear;
close all;

m = 50;
k = 1e6;
Fo = 1e4;
f = 20;
wf = 2*pi*f;
phi = pi/8;

dt = 0.001;
tspan = 0:dt:10;
z0 = [0.01; -2];

c = 50;
wn = sqrt(k/m);
h = c / (2 * sqrt(m * k));

[t, z] = ode45(@(t, z) odefunc(t, z, m, c, k, Fo, wf, phi), tspan, z0);

37
VIBRATION ANALYSIS TECHNIQUES 20.04

xn = z(:, 1);
vn = z(:, 2);
subplot(2, 1, 1);
plot(t, xn, 'b', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);
title_str = sprintf('Displacement: wn = %.4f Hz, Damping ratio (h) = %.6f', wn/2/pi, h);
title(title_str);
xlabel('Time [s]', 'Color', 'blue', 'FontWeight', 'bold');
ylabel('Displacement [m]', 'Color', 'blue', 'FontWeight', 'bold');
ylim([-0.08 0.08]); grid on;

subplot(2, 1, 2);
plot(t, vn, 'b', 'DisplayName', ['c = ', num2str(c), ' Ns/m ']);
title_str = sprintf('Velocity: wn = %.4f Hz, Damping ratio (h) = %.6f', wn/2/pi, h);
title(title_str);
xlabel('Time [s]', 'Color', 'blue', 'FontWeight', 'bold');
ylabel('Velocity [m/s]', 'Color', 'blue', 'FontWeight', 'bold');
ylim([-15 10]);
grid on;

function dzdt = odefunc(t, z, m, c, k, Fo, wf, phi)

F = Fo * cos(wf * t + phi);
dzdt = [z(2); (F - c * z(2) - k * z(1)) / m];
end

 Result :
Displacement: wn = 22.5079 Hz, Damping ratio (h) = 0.003536
0.08

0.06

0.04
Dispalcement [m]

0.02

-0.02

-0.04

-0.06

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

Velocity: wn = 22.5079 Hz, Damping ratio (h) = 0.003536

10

5
Velocity [m/s]

-5

-10

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

38
VIBRATION ANALYSIS TECHNIQUES 20.04

Top plot
(Displacement):
Observing the plot, the displacement gradually decreases over time, which
is consistent with a system experiencing light damping. The displacement
continues oscillating with smaller and smaller amplitudes, showing that the
damping is not strong enough to completely eliminate oscillations.
Bottom plot (Velocity):
The velocity oscillates strongly at the beginning and decreases over time.
Similar to the displacement plot, the velocity becomes more stable and its
amplitude decreases as time approaches 10 seconds.
Observation:
 The system is under the influence of a harmonic forcing function, which
explains the clear harmonic oscillations in both the displacement and
velocity graphs.
 Since the damping ratio is small (h < 1), the system continues oscillating,
but the intensity of the oscillations (amplitude) gradually decreases over
time due to the presence of damping, though the damping is not strong
enough to stop the oscillations immediately.
 Overall, the plots show a gradual decay of oscillations, but with some
sustained smaller oscillations over a long period.

2. Calculation of Fast Fourier Transform of vibration and force applied to

system

 Mathematical formula:
Change signal from time domain to frequency domain and vice versa
 A fast Fourier transform (FFT) is an algorithm to compute the discrete
Fourier transform (DFT) and its inverse
Syntax: fft(x)
 Use the command fft, we have to multiply the result by 2 / N, except for
the value at zero frequency, which is normalized to 1 / N (N is length of
signal)

FFT_X = 2*fft(X)/N
FFT_X(1) = FFT_X(1)/2

FFT is symmetric, throw away second half

X = FFT_X(1:N/2)
Frequency range
f = (0:N/2-1)*Fs/N; with Fs is sampling frequency

39
VIBRATION ANALYSIS TECHNIQUES 20.04

Take the magnitude of fft of x

amplitude = abs(X);
• Take the phase of fft of x

phase = angle(X);
[rad] phase = angle(X)*180/pi; [°]
• Take the mean value of x mean_
value = X(1);

 Code matlab FFT :

n = length(t);
X_f = fft(x);
f = (0:n-1)*(1/(n*dt));
X_f = X_f(1:n/2+1);
X_f(2:end-1) = 2*X_f(2:end-1);
f = f(1:n/2+1);
figure('Name', 'FFT of displacement');
subplot(313),plot(f, abs(X_f)/n, 'bo-', 'LineWidth', 1.5);
xlabel('Frequency [Hz]', 'FontWeight', 'bold');
ylabel('Amplitude [m]', 'FontWeight', 'bold');
xlim([0 40]),yticks(0:5:40),grid on;
ylim([0 0.05]), yticks(-0:0.01:0.05),grid on;
title('FFT of displacement','Color','black','fontsize',15);

 Result :

FFT Plot Analysis:

40
VIBRATION ANALYSIS TECHNIQUES 20.04

The FFT plot shows the amplitude of the displacement response at

different frequencies. The prominent peak on the plot indicates the dominant
frequency component of the displacement.

Observation:

Dominant frequency: The peak occurs at around 20 Hz, which corresponds to

the frequency of the applied harmonic force. This indicates that the system is
resonant at the forcing frequency.

Amplitude: The amplitude of the peak is relatively high, indicating a significant

response of the system at the resonant frequency.

Other frequencies: There are smaller peaks at other frequencies, which may
represent harmonics of the forcing frequency or other components in the system
response.
Explanation:

The FFT plot confirms that the system is experiencing a resonant response
to the applied harmonic force. The displacement amplitude at the resonant
frequency is significantly amplified compared to other frequencies. This is a
characteristic of low damping systems.

3.3.7 Asignment 3

41
VIBRATION ANALYSIS TECHNIQUES 20.04

1. Find respond of system

2. Calculation of Fast Fourier Transform of vibration and force applied to

system

3. Change observation time value with 10 and 10.1 to see effect of leakage error

 Code matlab by Analytical solution :

Clc;
clear;
close all;

m = 50;
k = 1e5;
c = 150;
x0 = 0.01;
v0 = -2;

f1_change = [10, 15, 21];

F0_change = [5000, 7000, 9000];
phi_change = [pi/4, pi/6, pi/8];

line_width = 1;

figure(1);
tiledlayout(2, 1);

figure(2);
hold on;
title('FFT of Displacement', 'FontWeight', 'bold');

for i = 1:length(f1_change)
f1 = f1_change(i);
F0 = F0_change(i);
phi = phi_change(i);
wn = sqrt(k/m);
dt = 0.01;
t = 0 : dt : 10;
wf = 2*pi*f1;
h = c/(2*sqrt(k*m));
Y1 = (F0/2 * exp(1i*phi))/(-m*wf^2 + 1i*c*wf + k);
Y2 = (F0/2 * exp(-1i*phi))/(-m*wf^2 - 1i*c*wf + k);

if h == 1
s = -wn;
A = [ 1 0 ; s 1];
42
VIBRATION ANALYSIS TECHNIQUES 20.04

B = [ x0-Y1-Y2; v0-1i*wf*Y1+1i*wf*Y2];
Y = A\B;
C1 = Y(1);
C2 = Y(2);
xn = C1*exp(s*t) + C2*t.*exp(s*t);
vn = C1*s*exp(s*t) + C2*(s*t+1).*exp(s*t);
elseif h > 1
s1 = -h*wn + wn*sqrt(h^2-1);
s2 = -h*wn - wn*sqrt(h^2-1);
A = [ 1 1 ; s1 s2];
B = [ x0-Y1-Y2; v0-1i*wf*Y1+1i*wf*Y2];
Y = A\B;
C1 = Y(1);
C2 = Y(2);
xn = C1*exp(s1*t) + C2*exp(s2*t);
vn = C1*s1*exp(s1*t) + C2*s2*exp(s2*t);
else
s1 = -h*wn + 1i*wn*sqrt(1-h^2);
s2 = -h*wn - 1i*wn*sqrt(1-h^2);
A = [ 1 1 ; s1 s2];
B = [ x0-Y1-Y2; v0-1i*wf*Y1+1i*wf*Y2];
Y = A\B;
C1 = Y(1);
C2 = Y(2);
xn = C1*exp(s1*t) + C2*exp(s2*t);
vn = C1*s1*exp(s1*t) + C2*s2*exp(s2*t);
end

xf = Y1*exp(1i*wf*t) + Y2*exp(-1i*wf*t);
vf = Y1*1i*wf*exp(1i*wf*t) + Y2*(-1i*wf)*exp(-1i*wf*t);
x = xn + xf;
v = vn + vf;
figure(1);
nexttile(1);
plot(t, x, 'b', 'LineWidth', line_width, 'DisplayName',
['f1 = ', num2str(f1), ' Hz']);
hold on;
xlabel('Time [s]', 'FontWeight', 'bold');
ylabel('Displacement [m]', 'FontWeight', 'bold');
title('Displacement Over Time');
grid on;
nexttile(2);
plot(t, v, 'b', 'LineWidth', line_width, 'DisplayName',
['f1 = ', num2str(f1), ' Hz']);
hold on;
xlabel('Time [s]', 'FontWeight', 'bold');
ylabel('Velocity [m/s]', 'FontWeight', 'bold');
title('Velocity Over Time');
grid on;
figure(2);
n = length(t);
43
VIBRATION ANALYSIS TECHNIQUES 20.04

X_f = fft(x);
f = (0:n-1)*(1/(n*dt));
X_f = X_f(1:n/2+1);
X_f(2:end-1) = 2*X_f(2:end-1);
f = f(1:n/2+1);

plot(f, abs(X_f)/n, 'b-o', 'LineWidth', line_width,

'DisplayName', sprintf('f1 = %d Hz', f1));
hold on;
xlabel('Frequency [Hz]', 'FontWeight', 'bold');
ylabel('Amplitude [m]', 'FontWeight', 'bold');
xlim([0 30]);
grid on;
end

figure(1);
nexttile(1);
legend('show');
nexttile(2);
legend('show');

figure(2);
legend('show');
title('FFT of Displacement', 'FontWeight', 'bold');
hold off;

 Result :

44
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 3.3.4: The graph of Fast Fourier Transform of vibration and force applied to
system

Obeservation:
The graph shows the displacement of an object over time under the action of
forces of different frequencies.
 Multiple frequencies: The graph shows the superposition of multiple
oscillations of different frequencies (10 Hz, 15 Hz, 21 Hz). This creates a
complex signal, which is no longer a simple sine wave.
 Amplitude of oscillation: The amplitude of oscillation varies over time
and is not steady. This indicates the presence of oscillation components of
different amplitudes.
 Angular frequency: The angular frequency of the oscillation components
can be estimated based on the number of times the time axis is crossed in
a unit of time.
 Initial phase: The oscillation components may have different initial
phases, resulting in phase shifts between them.
Physical meaning
45
VIBRATION ANALYSIS TECHNIQUES 20.04

 Multi-degree-of-freedom system: The graph shows a system with multiple

degrees of freedom, i.e. the system can oscillate in many different ways.
 Resonance phenomenon: If there is a frequency among the frequencies
that coincides or is close to the natural frequency of the system, the
oscillation amplitude at that frequency will increase significantly, causing
resonance.
 Effect of forced force: Forces acting on the system with different
frequencies will cause the system to oscillate at corresponding
frequencies.

FFT of Displacement graph analysis

The FFT (Fast Fourier Transform) graph of displacement shows us the
frequency spectrum of the oscillation signal, that is, it shows the frequency
components present in that signal with corresponding amplitudes.
 Horizontal axis: Represents the frequency (Hz).
 Vertical axis: Represents the amplitude of the frequency components (m).
 Peaks: Each peak on the graph represents a frequency component in the
original signal. The height of the peak indicates the amplitude of that
frequency component.
 Three curves: Each curve corresponds to an original signal with
frequencies of 10 Hz, 15 Hz, and 21 Hz.

Observation:
 Main peaks: We can clearly see that the main peaks appear at frequencies
of 10 Hz, 15 Hz, and 21 Hz, corresponding to the frequencies of the
original signals. This proves that the Fourier transform has successfully
separated the different frequency components in the signal.
 Secondary peaks: In addition to the main peaks, there are also secondary
peaks appearing around the main peaks. These secondary peaks can be
caused by many reasons, such as noise, the finite resolution of the Fourier
transform, or harmonic components in the signal.
 Amplitude: The amplitude of the main peaks shows the intensity of the
corresponding frequency components in the original signal.

 Peak width: The width of the peaks shows the purity of the frequency
components. The narrower the peak, the purer the frequency component.
Physical Meaning

46
VIBRATION ANALYSIS TECHNIQUES 20.04

 Frequency Component Analysis: The FFT graph allows us to analyze a

complex signal into simple frequency components, thereby helping us to
better understand the nature of that signal.
 Identifying Resonant Frequency: If the system has natural frequencies, the
peaks in the FFT graph will indicate the resonant frequencies of the
system.
 Error Detection: By comparing the actual frequency spectrum with the
ideal frequency spectrum, we can detect errors or abnormalities in the
system.

CHAPTER 4: CASE STUDY

47
VIBRATION ANALYSIS TECHNIQUES 20.04

CASE STUDY 1: BEAM

4.1 Analytical solution
Natural freq. & Model shapes
Consider a simply supported beam (pinned-pinned)
 Materia: Aluminum
 Length: L = 2 [m]
 Width: b = 100 [mm]
 Thickness: h = 200 [mm]
 Young’s modulus: E = 200e9 [N/m^2]
Mass density: ρ=7850 [kg/m^3]

1. Calculate the first 5 natural frequencies of the system

2. Plot the first 5 mode shapes of the system

Natural frequencies: ω n = A=¿

Where:
 ω n: nth natural angular frequency

 n: Mode number (n = 1, 2, 3, ...)

 L: Length of the beam

 E: Modulus of elasticity of the material

 J: Moment of inertia of the beam's cross-section about the y-axis

 ρ: Mass density of the material

48
VIBRATION ANALYSIS TECHNIQUES 20.04

 A: Cross-sectional area of the beam

Moment of inertia J:
3
bh
 For a rectangular section: J =
12
4
πd
 For a circular section: J =
64
 J represents the resistance of the cross-section to bending.

Mode Shape
 Definition: A mode shape describes the deformed shape of the beam when
it vibrates at a particular natural frequency.
 Formula:
Wn(x) = ¿)
Where:
 Wn(x) : Deflection of the beam at position x for the nth mode

 x: Position along the beam

Free vibration:
w(x,t) = ∑ sin¿ ¿ ).[Ancos ¿ ¿ ) + Bnsin(ω¿¿ n t ¿)¿ ¿]
Where:
 w(x,t): Deflection of the beam at position x and time t

 An, Bn: Constants determined by the initial conditions of the vibration

4.2 Code matlab :
clc;
clear;
close all;
E = 200e9;
rho = 7850;
b = 0.1;
h = 0.2;
length = 2;
J = b*h*h*h/12;
S = b*h;
dx = 0.001;
x = 0:dx:length;
N = 5;
figure('name','Analytical mode shapes')
for n = 1:N
omega(n) = (n*pi/length)^2*sqrt(E*J/rho/S);
omega_hz(n) = omega(n)/(2*pi);
49
VIBRATION ANALYSIS TECHNIQUES 20.04

disp(['omega ',num2str(n),' =
',num2str(omega_hz(n)),' [Hz]'])
mode(n,:) = sin(n*pi*x/length);
subplot(N,1,n)
plot(x,mode(n,:),'linewidth',3),
title(['Mode #',num2str(n),' - Freq. =
',num2str(omega_hz(n)),' Hz'],'FontSize',13)
hold on
plot(0,0,'color','r','marker','.','markersize',20)

plot(length,0,'color','r','marker','.','markersize',20)
end
xlabel('x [m]','FontSize',13)

 Result :
Mode #1 - Freq. = 0.02259Hz
1

0.5

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode #2 - Freq. = 0.090359Hz

1

-1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode #3 - Freq. = 0.20331Hz

1

-1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode #4 - Freq. = 0.36143Hz

1

-1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode #5 - Freq. = 0.56474Hz

1

-1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x [m]

Figure 4.1: Mode shapes of a cantilever beam

Observation:
The vibration modes from 1 to 5 correspond to the natural frequencies:
 Mode 1: This mode has one antinode in the middle of the beam and no
nodes.
 Mode 2: In this mode, the frequency has doubled and a node has appeared
in the middle, where the amplitude is zero. Here, the two halves of the
beam vibrate out of phase.
 Mode 3: In this mode, we can clearly see three antinodes and two nodes,
the vibration amplitude has increased and the shape has become more
complex.

50
VIBRATION ANALYSIS TECHNIQUES 20.04

 Mode 4: The vibration in this mode has three nodes, the system vibrates at
a higher frequency, with four antinodes in this mode and all vibrations are
out of phase with each other.
 Mode 5: This is the highest vibration mode in this graph, with four nodes
and five antinodes, and the highest frequency of all.
The frequency increases with each vibration mode.
 Red dots: represent fixed points.

 Curves: The curves in each graph represent the vibration of the system at
length x, the points with the largest amplitude are maxima and the points
with zero amplitude are nodes.

4.3 FEM: Numerical model

Numerical Approach:
- Number of beam elements : Nel
- Number of nodes : Nnodi = Nel + 1
- 4 DOFs per element
- Total number of DOF : Ngdl = 4·Nnodi
DOF final element of the j-th :

Procedure of FE model :
1. Divide the beam in Nel elements
2. Build the 4 × 4 matrices for each beam element of the rotor
 M = classical mass matrix
 K = stiffness matrix
3. Assembly the total 2·Nnod × 2·Nnod matrices of the beam by means of
the matrices of each beam element
4. Add the contribution of the boundary conditions
Matrices of each beam element

51
VIBRATION ANALYSIS TECHNIQUES 20.04

Assembly global matrix system

By considering all the DOFs of the rotor (2·Nnode ) :

It is necessary to assembly the matrices of each beam element in global matrices

Boundary condition
1. The beam is simple supported pined-pined.
52
VIBRATION ANALYSIS TECHNIQUES 20.04

2. We have to take into account stiffness of this boundary condition

[ K (bc) ] = [ k xx k xy
k yx k yy ]
3. It is necessary to sum the contribution of the boundary condition to
obtain:

[ M (b ) ] ẍ (b ) + [ K (b ) ] x b=W + F(t)
[ K ] = [K(b)] + [K(bc)]
4. For simple, we can neglect damping of boundary condition

Contribution of boundary condition

Sample assembly (stiffness matrix)

The homogeneous equation can be considered: [ M ] ẍ + [ K ]x = 0

4.4Code matlab :
clc
clear
close all
E = 200e9;
rho = 7850;
b = 0.1;
h = 0.2;
La = 2;

53
VIBRATION ANALYSIS TECHNIQUES 20.04

Ndel = 10;
Bn = [1, Ndel + 1];
N_dof = 2;
dx = 0.001;
x =0:dx:La ;
N = 5;

L = ones(1, Ndel) * La/ Ndel;

S = ones(1, Ndel) * b * h;
J = ones(1, Ndel) * (b * h^3) / 12;

Mg = zeros(N_dof * (Ndel + 1));

Kg = zeros(N_dof * (Ndel + 1));

for m = 1:Ndel
l = L(m);
jj = J(m);
s = S(m);
Mj = rho*s * l / 420 * ...
[156 22 * l 54 -13 * l;
22 * l 4 * l^2 13 * l -3 * l^2;
54 13 * l 156 -3 * l^2;
-13 * l -3 * l^2 -3 * l^2 4 * l^2];
Kj = E * jj / l^3 * ...
[12 6 * l -12 6 * l;
6 * l 4 * l^2 -6 * l 2 * l^2;
-12 -6 * l 12 -6 * l;
6 * l 2 * l^2 -6 * l 4 * l^2];
j1 = (m - 1) * N_dof + 1;
j4 = (m + 1) * N_dof;
Mg(j1:j4, j1:j4) = Mg(j1:j4, j1:j4) + Mj;
Kg(j1:j4, j1:j4) = Kg(j1:j4, j1:j4) + Kj;
end

BXstiff = 1e12;
BTstiff = 0;

Bk = zeros(N_dof,N_dof,length(Bn));
for iBn = 1:length(Bn)
Bk(:, :, iBn) = diag([BXstiff, BTstiff]);
end

for n = 1:length(Bn)
j1 = (Bn(n) - 1) * N_dof + 1;
j2 = Bn(n) * N_dof;
Kg(j1:j2, j1:j2) = Kg(j1:j2, j1:j2) + Bk(:, :, n);
end

[modes, lambda] = eig(Mg\Kg);

[wn, indexes] = sort(sqrt(diag(lambda)));
modes = modes(:, indexes);
54
VIBRATION ANALYSIS TECHNIQUES 20.04

vertical = modes(1:2:end, :);

torsional = modes(2:2:end, :);

wn_hz = wn / (2 * pi);

disp(‘The first ten natural oscillation frequencies of the system (Hz):’)

for j = 1:min(10, length(wn_hz))
fprintf('Frequency %d: %.4f Hz\n', j, wn_hz(j));
end

len = La;
beam = Ndel;
xcor = 0:len/beam:len;

interp_factor = 10;
xcor_smooth = linspace(0, len, length(xcor) * interp_factor);
for ii = 1:5
figure(ii)
vertical_smooth = interp1(xcor, vertical(:, ii), xcor_smooth, 'spline');
torsional_smooth = interp1(xcor, torsional(:, ii), xcor_smooth, 'spline');

subplot(211)
plot(xcor_smooth, vertical_smooth, 'r', 'LineWidth', 1.5); %
hold on; grid on
plot(xcor(Bn), zeros(length(Bn), 1), 'rs', 'MarkerEdgeColor', 'k', ...
'MarkerFaceColor', 'g', 'MarkerSize', 10);
ylabel('Vertical Displacement');
title(['Mode Shape ', num2str(ii), ' - Frequency: ', num2str(wn_hz(ii), '%.4f'), ' Hz']);
legend('Vertical');

subplot(212)
plot(xcor_smooth, torsional_smooth, 'b', 'LineWidth', 1.5); %
grid on
xlabel('Length of Beam');
ylabel('Torsional Displacement');
title(['Mode Shape ', num2str(ii), ' - Frequency: ', num2str(wn_hz(ii), '%.4f'), ' Hz']);
legend('Torsional');
end

The first ten natural oscillation frequencies of the system (Hz):

Frequency 1: 114.4373 Hz
Frequency 2: 457.7430 Hz
Frequency 3: 1030.1146 Hz
Frequency 4: 1832.4970 Hz
Frequency 5: 2867.0904 Hz
Frequency 6: 4136.8877 Hz
Frequency 7: 5641.4165 Hz
55
VIBRATION ANALYSIS TECHNIQUES 20.04

Frequency 8: 7357.2675 Hz
Frequency 9: 9139.1221 Hz
Frequency 10: 12909.2714 Hz

 Result :

Mode Shape 1 - Frequency: 114.4373 Hz

0
Vertical Displacement

Vertical
-0.05

-0.1

-0.15

-0.2

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode Shape 1 - Frequency: 114.4373 Hz

0.4
Torsional Displacement

Torsional
0.2

-0.2

-0.4
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length of Beam

56
VIBRATION ANALYSIS TECHNIQUES 20.04

Mode Shape 2 - Frequency: 457.7430 Hz

Vertical Displacement 0.1 Vertical

0.05

-0.05

-0.1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode Shape 2 - Frequency: 457.7430 Hz

0.4
Torsional Displacement

Torsional
0.2

-0.2

-0.4
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length of Beam

Mode Shape 3 - Frequency: 1030.1146 Hz

0.1
Vertical Displacement

Vertical
0.05

-0.05

-0.1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode Shape 3 - Frequency: 1030.1146 Hz

0.5
Torsional Displacement

Torsional

-0.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length of Beam

57
VIBRATION ANALYSIS TECHNIQUES 20.04

Mode Shape 4 - Frequency: 1832.4970 Hz

Vertical Displacement 0.05 Vertical

-0.05

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode Shape 4 - Frequency: 1832.4970 Hz

0.5
Torsional Displacement

Torsional

-0.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length of Beam

58
VIBRATION ANALYSIS TECHNIQUES 20.04

Mode Shape 5 - Frequency: 2867.0904 Hz

0.05
Vertical Displacement Vertical

-0.05
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Mode Shape 5 - Frequency: 2867.0904 Hz

0.5
Torsional Displacement

Torsional

-0.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Length of Beam

• Mode 1:
The graph describes the first mode shape of a beam when it oscillates with a
natural frequency of 114.4373 Hz.
The first graph: Vertical displacement
 Shape: The graph has a concave shape, reaching a minimum value in the
middle of the beam and a maximum value at both ends of the beam.
 Meaning: This shows that when the beam oscillates in this mode shape,
the middle part of the beam will be bent downward, while the two ends of
the beam will be bent upward.
The second graph: Torsional displacement
 Shape: The graph has a nearly straight shape, increasing from one end of
the beam to the other.
 Meaning: This shows that when the beam oscillates in this mode shape,
the entire beam will be twisted gradually from one end to the other.

59
VIBRATION ANALYSIS TECHNIQUES 20.04

General Conclusion
 Complex oscillation: The beam is performing two types of oscillations at
the same time: vertical bending oscillation and torsional oscillation.
 Natural frequency: The natural frequency of this mode shape is 114.4373
Hz, meaning that if the beam is excited by a force with a frequency equal
to or close to this frequency, the beam will oscillate with the largest
amplitude.
• Mode 2:
The graph describes the second oscillation shape (mode shape) of a beam
when oscillating with a natural frequency of 457.7430 Hz.
The first graph: Vertical displacement
 Shape: The graph has the shape of a sine wave, with two maximum points
at the two ends of the beam and one minimum point in the middle of the
beam.
 Meaning: This shows that when the beam oscillates in this mode shape,
the middle of the beam will be bent upwards, while the two ends of the
beam will be bent downwards. Compared to the first mode shape, this
mode shape has a oscillating node in the middle of the beam.
Second graph: Torsional displacement
 Shape: The graph has the shape of a sine wave, with a maximum point
and a minimum point.
 Meaning: This shows that when the beam oscillates in this mode shape,
the entire beam will be twisted, with one part being twisted clockwise and
the other part being twisted counterclockwise.
Compared to mode shape 1
 Number of oscillating nodes: Mode shape 2 has more oscillating nodes
than mode shape 1.
 Shape: The shape of the curves describing the displacement is also more
complex.
General conclusion
 More complex oscillation: Mode shape 2 describes a more complex
oscillation pattern than mode shape 1, with more oscillating nodes and a
more complex curve shape.
60
VIBRATION ANALYSIS TECHNIQUES 20.04

 Higher frequency: The natural frequency of mode shape 2 is higher than

mode shape 1, which shows that the beam will oscillate faster in this
mode shape.
Mode 3:
The graph describes the third mode shape of a beam when oscillating with
a natural frequency of 1030.1146 Hz.
The first graph: Vertical displacement
 Shape: The graph has the shape of two consecutive sine waves, with three
maximum points and two minimum points.
 Meaning: This shows that when the beam oscillates in this mode shape,
the beam will bend up and down many times, with two maximum bending
positions and one minimum bending position in the middle. Compared to
the previous mode shapes, this mode shape has more oscillation nodes.
The second graph: Torsional displacement
 Shape: The graph has the shape of two consecutive sine waves.
 Significance: This shows that when the beam oscillates in this mode
shape, the entire beam will be twisted many times, with parts twisted in
opposite directions.
Comparison with previous mode shapes
 Number of oscillating nodes: Mode shape 3 has more oscillating nodes
than mode shapes 1 and 2.
 Shape: The shape of the curves describing the displacement is also more
complex.
 Frequency: The natural frequency of mode shape 3 is higher than that of
mode shapes 1 and 2, which indicates that the beam will oscillate faster in
this mode shape.
General conclusion
 Complex oscillation: Mode shape 3 describes a more complex oscillation
pattern than the previous mode shapes, with more oscillating nodes and a
more complex curve shape.
 High frequency: The natural frequency of mode shape 3 is the highest
among the considered mode shapes, which indicates that the beam will
oscillate fastest in this mode shape.
61
VIBRATION ANALYSIS TECHNIQUES 20.04

• Mode 4:
The graph describes the fourth mode shape of a beam when oscillating at a
natural frequency of 1832.4970 Hz.
The first graph: Vertical displacement
 Shape: The graph has the shape of three consecutive sine waves, with four
maximum points and three minimum points.
 Meaning: This shows that when the beam oscillates in this mode shape,
the beam will bend up and down many times, with three maximum
bending positions and two minimum bending positions. Compared to the
previous mode shapes, this mode shape has more oscillation nodes.
The second graph: Torsional displacement
 Shape: The graph has the shape of three consecutive sine waves.
 Meaning: This shows that when the beam oscillates in this mode shape,
the entire beam will twist many times, with parts twisted in opposite
directions.
Compared to previous mode shapes
 Number of oscillation nodes: Mode shape 4 has more oscillation nodes
than the previous mode shapes.
 Shape: The shape of the curves describing the displacement is also more
complex.
 Frequency: The natural frequency of mode shape 4 is higher than the
previous mode shapes, which indicates that the beam will oscillate faster
in this mode shape.
General Conclusion
 Complex oscillation: Mode shape 4 describes a more complex oscillation
pattern
General Observations of Vibration Modes:
 Increasing frequency: The natural frequency of the modes increases from
mode 1 to mode 5. This is reasonable as modes with higher frequencies
typically have more complex deformation shapes and involve more
localized vibrations.
 Increasing number of nodes: The number of nodes (points where the
amplitude of vibration is zero) increases with the mode number. This
62
VIBRATION ANALYSIS TECHNIQUES 20.04

indicates that the structure is divided into more segments vibrating out of
phase with each other as the frequency increases.
 Mode shapes: The mode shapes become more complex as the frequency
increases. Mode 1 usually has the simplest shape, while higher modes
may have more inflection points and more complex deformation patterns.
 Vibration amplitude: The vibration amplitude at antinodes (points where
the vibration amplitude is maximum) can vary depending on the mode
and the location of the antinode.

63
VIBRATION ANALYSIS TECHNIQUES 20.04

CASE STUDY 2: ROTORKIT

PART 1: IMPACT TEST
1. Data acquisition:
• Kind of test: Test A and B.
• P = number of tests (p = 1, ..., P).
• N = number of accelerometers or measuring point (N = 3).
• i = position of the accelerometer.
Acquire the signals of 03 accelerometers and the 01 force of the hammer for:
• P = 10 impact tests (used for the average process)
• Time acquisition: 15s
• Sampling frequency: Fs = 10240 Hz
• Filename example: A1.mat , B1.mat , ...

INPUT (X) xp(t) signal of the hammer for the p-test

OUTPUT (Y) yp(t) signal of the i-th accelerometer for the p-test

2 Load signal:
Load data from file “Group 1” to collect signal.

3. Trigger and window:

Eliminate signals that have double hits and only signal without double hits.
1
Frequency resolution: ∆f=
Tw

Figure 4.7 analys the raw signal and discard signal

The "double hit" phenomenon in signal analysis and processing occurs when a
peak or event in the signal is detected multiple times instead of once. This can
lead to several issues, especially when analyzing parameters such as frequency,
amplitude, or pulse width of the signal.
When the "double hit" phenomenon occurs, the characteristics of the signal are
affected in various ways, potentially leading to incorrect analysis results,
including:
- Errors in calculating frequency and amplitude.
64
VIBRATION ANALYSIS TECHNIQUES 20.04

- Inaccurate pulse width.

- Incorrect classification results.

Figure 4.8 without double hit and double hit

Figure 4.9 double hit

In the field of signal analysis and processing, a "double hit" can refer to the
duplication or doubling of an event within the signal data.
Noise or errors in peak detection: In signal processing, while identifying peaks
in a signal, a "double hit" can occur if a peak is detected multiple times due to
noise or error. This can lead to inaccuracies in calculating frequency, amplitude,
or other parameters of the signal.
In filtering or signal transformation techniques: The filtering or transformation
process (such as FFT) can create "double hits" in the form of spurious signals,
which appear as real signals but are actually artifacts of the processing.

65
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.10 double hit

4.Check saturation:
In step, you can set if it is necessary.
5..FFT of signal:
Evaluate the frequency spectrum of the hammer and each accelerometer,
considering the sensitivity of each sensor.
yp(t)  Yp(f) Accelerometer signal (output).
x (t)  X (f)
p p
Hammer signal (input)
The DC component in an FFT (Fast Fourier Transform) refers to the zero-
frequency component of the signal. It represents the average or constant part of
the signal, essentially the signal's mean value over time.
 Identifying Signal Offset: The DC component indicates whether there is
an average offset in the oscillating signal. If the DC component is large, it
suggests that the signal does not oscillate around the zero axis but has a
fixed part. This could be due to a shift in the measurement system or
because the signal has an added constant value (offset).
 Removing Unnecessary Components: In many cases, when analyzing
oscillating signals - especially periodic oscillating signals like sine or
square waves—the DC component does not provide useful information
about the oscillations. Therefore, removing the DC component (by
subtracting the signal's mean value) helps focus on analyzing the actual
oscillation frequencies.

6. Cross/Auto-Spectrum
66
VIBRATION ANALYSIS TECHNIQUES 20.04

Evaluate the cross-spectrum of each test for the accelerometer signal and the
auto-spectrum for the hammer signal.
Gyx(f) = Yp(f) *. Xp(f) cross - spectrum for the pth – test.
Gxx(f) = p(f) *. Xp(f) auto - spectrum for the pth – test (Hammer Signal).
Complex Vector YP(f) = a + ib
 This is a complex vector in the complex plane, consisting of:
 a: the real part of the complex vector.
 ib: the imaginary part, where b is the imaginary component and i is the
imaginary unit (the square root of -1).
 This complex vector represents a point in the complex plane, including
both the real and imaginary parts.
 The complex vector represents the signal in the frequency domain.

Complex Conjugate YP(f)* = a - ib

 This is the complex conjugate of vector YP(f). The complex conjugate of a
complex number is created by reversing the sign of the imaginary part,
while the real part remains the same. Specifically:
 The real part a remains unchanged.
 The imaginary part ib becomes -ib.
 The complex conjugate is crucial in many fields, such as signal
processing and Fourier analysis. To calculate real quantities like the
magnitude or power of a signal, the complex vector must be multiplied by
its conjugate.
 Calculate real quantities such as amplitude or signal energy.

7. Average value:
Evaluate the average values of the cross-spectrum and of the autospectrum for
p-th tests.
P
1
G yx = P ∑ Y P (f)* XP (f)
p=1
P
1
G xx = P ∑ X P (f)* XP (f)
p=1

8. Estimate of the FRF:

The i-th transfer function (excitation at k-th node) can be evaluated by means of
the H1 estimator:
G yx ( f ) Cross−spectrum(input /output )
H1 = G xx (f ) Auto−spectrum(input=hammer )
This estimator gives the better estimate of transfer function when the signal of
the
response (output) is affected by noise.

67
VIBRATION ANALYSIS TECHNIQUES 20.04

G yy ( f )
H2 estimator H2 = G yx ( f )
Auto−spectrum(output=accelerometer)
Cross−spectrum(input /output)
It gives the better estimate of transfer function when the signal of the input
(excitation / hammer) is affected by noise.
Evaluate transfer function by means of the H1 estimator and coherence
2
G yx ( f ) |G AB (f )| X (f )
H1 = G xx (f )
γ AB(f) = G AA ( f ) .G BB (f )
X(f) = ω2

Figure 4.11 estimate H1

9. Coherence:
For the analysis of systems in which one input and one output are measured, it is
necessary to evaluate the linearity between the output and the input.
This is important, due to the fact that the frequency analysis is based on the
linearity of the system for which the superposition principle could be applied.
In frequency domain, the coherence function is a index of the linearity between
two signals a(t) and b(t) for a given frequency:
2

γ AB(f)
|G AB (f )|
= G AA ( f ) .G BB ( f )
This function assumes values from 0 to 1. The value 1 means a perfect
correlation
between a(t) and b(t), that is the output is given by the properties of the system
and
is generated only by the input and not by other sources.

68
VIBRATION ANALYSIS TECHNIQUES 20.04

10. Code matlab:

clc
clear
close all

Fs = 10240;
Twindow = 10;
P = 10;
t = 15;
freq = 1/Twindow;
N1 = 200;
Hammer_chanel = 4;
ACC_chanel = 1:3;
folderPath = 'GROUP 1';
Gyx_mean = zeros(51200,3);
Gyy_mean = zeros(51200,3);
Gxx_mean = zeros(51200,1);
Gyx = zeros(51200,3);
Gyy = zeros(51200,3);
H1 = zeros(51200,3);
for i = 1:10
if (i ~= 8)
if (i ~= 10)
signal = load(fullfile('Chapter4', sprintf('A%d.mat', i)));
y = signal.signal(:, 1:3);
x = signal.signal(:, 4);
y = y/10e-3;
x = x/2.25e-3;
N = Fs*Twindow;
[max_hammer, index] = max(x);
first_point = index - 0.1*Fs;
last_point = first_point + Twindow*Fs -1;
t = linspace(0,10,Twindow*Fs+1);
x = x(first_point:last_point,:);
y = y(first_point:last_point,:);

n = length(y(:,1));
FFT_Y = 2*fft(y)/n;
FFT_Y(1) = FFT_Y(1)/2;
Y = FFT_Y(1:n/2,:);

FFT_X = 2*fft(x)/n;
FFT_X(1) = FFT_X(1)/2;
X = FFT_X(1:n/2);
Gxx = conj(X).*X
69
VIBRATION ANALYSIS TECHNIQUES 20.04

for k = 1:3
Gyx(:,k) = conj(Y(:,k)).*X;
Gyy(:,k) = conj(Y(:,k)).*Y(:,k);
Gyx_mean(:,k) = Gyx_mean(:,k) + Gyx(:,k);
Gyy_mean(:,k) = Gyy_mean(:,k) + Gyy(:,k);
end
Gxx_mean = Gxx_mean + Gxx;
end
end
end

Gyx_mean = Gyx_mean/8;
Gyy_mean = Gyy_mean/8;
Gxx_mean = Gxx_mean/8;
f = (0:(N/2-1))/N*Fs;
w = 2*pi*freq;
for k = 1:3
H1(:,k) = Gyx_mean(:,k)./Gxx_mean;
displacement = Y(:,k)./-w^2;
num = Gyx_mean(:,k).*conj(Gyx_mean(:,k));
den = Gxx_mean.*Gyy_mean(:,k);
coherence(:,k) = num./den;
end
for i=1:3
figure;
subplot(211);
plot(f,abs(H1(:,i)),'LineWidth', 2); ylabel(sprintf('A_{%d} [m/s^2/N]',
i),'FontWeight', 'bold');
xlim([5 200]);grid on;
subplot(212);
plot(f,angle(H1(:,i))*180/pi,'LineWidth', 2);
xlabel('Tan so [Hz]','FontWeight', 'bold'),ylabel('Pha [^o]','FontWeight',
'bold')
xlim([5 200]);grid on;
end
figure;
subplot(311);
plot(f, abs(H1(:,1)), 'DisplayName', 'A1', 'LineWidth', 2); hold on;
plot(f, coherence(:,1), '--', 'DisplayName', 'Coherence', 'LineWidth', 2);
ylabel(sprintf('A_{%d} [m/s^2/N]', 1), 'FontWeight', 'bold');
xlim([5 200]); grid on; legend show;
subplot(312);
plot(f, abs(H1(:,2)), 'DisplayName', 'A2', 'LineWidth', 2); hold on;
plot(f, coherence(:,2), '--', 'DisplayName', 'Coherence', 'LineWidth', 2);
ylabel('Pha [^o]', 'FontWeight', 'bold');
70
VIBRATION ANALYSIS TECHNIQUES 20.04

xlim([5 200]); grid on; legend show;

subplot(313);
plot(f, abs(H1(:,3)), 'DisplayName', 'A3', 'LineWidth', 2); hold on;
plot(f, coherence(:,3), '--', 'DisplayName', 'Coherence', 'LineWidth', 2);
xlabel('Frequency [Hz]', 'FontWeight', 'bold');
ylabel('Pha [^o]', 'FontWeight', 'bold');
xlim([5 200]); grid on; legend show;

71
VIBRATION ANALYSIS TECHNIQUES 20.04

 Result:
0.3

A1 [m/s /N]
2 0.2

0.1

0
20 40 60 80 100 120 140 160 180 200

200

100
Pha [ ]
o

-100

-200
20 40 60 80 100 120 140 160 180 200
Tan so [Hz]
Figure 4.12 amplitude A1
Observation:
Top Graph (Magnitude Response):
- The y-axis (labeled A in units of m/s^2, N^-1) likely represents the
amplitude or acceleration response per unit force, indicating how much the
structure vibrates at different frequencies.
- The x-axis shows the frequency in Hertz (Hz), likely up to 200 Hz.
- Peaks in the magnitude graph suggest resonance frequencies, where the
structure responds strongly to the applied impact force. Notable peaks occur
around 40 Hz, 120 Hz, and 160 Hz.
- These peaks reveal the natural frequencies of the structure, where the energy
input from the impact force aligns with the natural vibration modes.

Bottom Graph (Phase Response):

- The phase graph shows the phase shift in degrees as a function of frequency.
- The phase response undergoes significant changes, particularly near the
resonance frequencies, where it jumps sharply. This is typical, as the structure's
response often lags or leads to the input force by larger amounts at resonance
frequencies.
- Rapid phase changes near certain frequencies (e.g., around 40 Hz and 120
Hz) suggest a shift in the system’s dynamic behavior, indicating critical points
of resonance or damping effects.

72
VIBRATION ANALYSIS TECHNIQUES 20.04

10

A2 [m/s /N]
2
5

0
20 40 60 80 100 120 140 160 180 200

200

150
Pha [ ]
o

100

50

0
20 40 60 80 100 120 140 160 180 200
Tan so [Hz]
Figure 4.13 amplitude A2
Observation:
Top Graph (Magnitude Response A2):
- The y-axis represents A2 in units of m/s^2, N^-1), which again likely
indicates the amplitude response per unit force.
- There is a prominent peak at approximately 25 Hz, with a very high
magnitude compared to other frequencies. This indicates a strong resonance at
this frequency, where the structure has a very high response to the applied
impact force.
- The graph shows a smaller rise in magnitude around 110 Hz, though this
peak is far less pronounced than the one at 25 Hz.
- Beyond 110 Hz, the magnitude response becomes relatively stable and low,
indicating that the structure does not respond strongly at higher frequencies.

Bottom Graph (Phase Response - Phase Angle):

- The phase response starts near zero at low frequencies, gradually increases,
and shows a steep phase shift around 25 Hz, corresponding to the strong
resonance peak in the magnitude plot.
- There is another significant change in phase around 110 Hz, associated with
the smaller peak in the magnitude response.
- The phase response eventually stabilizes at around 150 degrees for higher
frequencies (above ~120 Hz), indicating a consistent phase lag in the response
relative to the input force.
Comparison with Previous Graph:
Compared to the first set of plots you shared, this structure shows a more
pronounced primary resonance (at 25 Hz), while the first set of data had multiple
73
VIBRATION ANALYSIS TECHNIQUES 20.04

peaks across the frequency range. This could imply that this structure (or
configuration) has a simpler or less complex vibrational profile, dominated by a
single mode at low frequencies.
4

3
A3 [m/s /N]
2

0
20 40 60 80 100 120 140 160 180 200

200

100
Pha [o ]

-100

-200
20 40 60 80 100 120 140 160 180 200
Tan so [Hz]
Figure 4.14 amplitude A3
Observation:
Top Graph (Magnitude Response A2):
Shows amplitude in units of m/s²/N
Has two main peaks:
- First peak at around 25-30 Hz, with amplitude of about 3.5 m/s²/N
- Second and highest peak at around 110-120 Hz, reaching approximately 4
m/s²/N
The curve generally decreases after the second peak
Has a stable region between 40-80 Hz with low amplitude of about 0.7 m/s²/N

Bottom Graph (Phase Response - Phase Angle):

Shows phase angle in degrees (°)
Has a dramatic shift at around 35-40 Hz, dropping from +180° to -180°
After this point, the phase angle gradually increases from -180° to about +30° at
200 Hz
The initial portion (0-35 Hz) shows slight oscillations around 0°

74
VIBRATION ANALYSIS TECHNIQUES 20.04

1
A1
0.8 Coherence

A1 [m/s2/N]
0.6

0.4

0.2

0
20 40 60 80 100 120 140 160 180 200

10
A2
8 Coherence
Pha [o]

0
20 40 60 80 100 120 140 160 180 200

5
A3
4 Coherence
Pha [o ]

0
20 40 60 80 100 120 140 160 180 200
Frequency [Hz]

Figure 4.15 amplitude A1 A2 A3

Observation:
Top Graph (A1):
- Blue line: Shows relatively low amplitude oscillation, around 0.1-0.2
- Red dashed line (Coherence): Fluctuates significantly between 0.2-1.0
- Multiple sharp drops in coherence, particularly noticeable in the 20-60 Hz
range
- Overall coherence becomes more stable after 80 Hz

Middle Graph (A2):

- Features a prominent peak at around 30 Hz with an amplitude of about 8
- Red dashed coherence line remains relatively stable at a low level
- Amplitude decreases rapidly after the peak and maintains low levels
- Shows a smaller peak around 120 Hz
- Generally, follows a resonance curve pattern

Bottom Graph (A3):

- First peak at around 30 Hz with amplitude of about 3.5
- Second, higher peak at around 110-120 Hz reaching amplitude of about 4
- Coherence line (red dashed) stays stable around 1
- Shows good agreement with expected system behavior

General observations:
- All three graphs share the same frequency range (0-200 Hz) on x-axis
- Correlation exists between peaks at around 30 Hz and 120 Hz across graphs
- Coherence indicates measurement reliability, with values closer to 1 being
more reliable
- A2 and A3 show similar variation patterns but with different amplitudes

75
VIBRATION ANALYSIS TECHNIQUES 20.04

- These graphs likely represent vibration measurements or modal analysis of a

mechanical system

76
VIBRATION ANALYSIS TECHNIQUES 20.04

PART2 : RUNNING
1. Experimental Vibration Analysis
Experimental vibration analysis is conducted through three types of tests:
run-up, constant speed, and run-down. Each of these processes is described as
follows:
 Run-up: In this process, the system starts from a low speed (e.g., 250
RPM) and gradually increases to maximum speed. Measuring vibrations
during this stage helps to determine how the system reacts during
acceleration and can clarify phenomena such as resonance or excessive
oscillations when reaching the system's natural frequencies.
 Constant Speed: After reaching maximum speed, the system maintains a
fixed speed for a duration (typically 10 seconds). This process is crucial
for observing the characteristics of stable vibrations and detecting
continuous oscillations at the main operating speed.
 Run-down: The system decelerates from maximum speed to a lower level.
This deceleration process allows for the observation of vibrations as the
system slows down and may highlight any abnormal signs that appear
during the stopping phase.
 Testing Parameters:
 Sampling frequency F s=20 kHz
 Signal recording time for constant speed is T w =10 s
 Maximum speed depends on the machine group configuration, for
example, 8000 RPM for Group 1 and 9000 RPM for Group 2.

2. Measurement and Signal Processing

To measure the vibrations of the system, proximity sensors are installed at
various positions on the rotating shaft, including:
 NDE (Non-Driven End): The unloaded end, i.e., the end without the
motor attached.
 DE (Driven End): The loaded end, i.e., the end with the motor attached.

At each position, the sensors can be mounted in two orthogonal directions:

 H (Horizontal): Horizontal direction.
 V (Vertical): Vertical direction.

Thus, there are four main signals from the sensors:

 NDEH: Non-driven end in the horizontal direction.
 NDEV: Non-driven end in the vertical direction.
 DEH: Driven end in the horizontal direction.
 DEV: Driven end in the vertical direction.
77
VIBRATION ANALYSIS TECHNIQUES 20.04

Tacho Signal: In addition to the signals from proximity sensors, the tacho
(tachometer) signal is used to measure and monitor the rotational speed of the
shaft. This signal provides information about the number of revolutions of the
system and serves as the basis for synchronizing vibration signals with the actual
rotational speed.
Each vibration signal from the sensors is recorded on a separate channel in
the measurement system, for example:
 Channel 1: NDEH
 Channel 2: NDEV
 Channel 3: DEH
 Channel 4: DEV
 Channel 5: Tacho

Figure 4.16.. Experimental Setup

3. Frequency Spectrum Analysis (Spectrogram)
Frequency spectrum analysis is a method for converting vibration signals
from the time domain to the frequency domain, helping to identify the frequency
components of the oscillation. The spectrogram is an important tool in this
analysis, especially when observing the changes in the frequency spectrum over
time.
Spectrogram: This is a chart that represents the frequency spectrum of the
signal over short time intervals, with overlapping time windows to ensure
continuity. The result is a 3D graph with:
 The x-axis representing time.
 The y-axis representing frequency.
78
VIBRATION ANALYSIS TECHNIQUES 20.04

 The z-axis (or color) indicating the intensity or amplitude of the signal at
each frequency at every moment.

Spectrogram Calculation Equation: The spectrogram of a signal can be

calculated by dividing the signal into small time windows and computing the
Fourier spectrum for each window:
N −1
S ( f , t )=STFT {x ( t ) }= ∑ x ( n ) ⋅w (t−n ) e− j 2 πfn
n=0

where:
 S ( f , t ) is the time-frequency spectrum at frequency f and time t ,
 x ( n ) is the time-domain signal,
 w (t−n ) is the window function (e.g., Hamming, Hanning),
 N is the window length.

Important Parameters in the Spectrogram:

 window: Window size.
 noverlap: The number of overlapping samples between windows to ensure
continuity.
 nfft: The length of the FFT for each window.
 Fs: The sampling frequency of the signal.

This frequency spectrum can be displayed in 3D using mesh, surf, or

waterfall commands in MATLAB for easy observation of changes in frequency
components over time.
4. Tacho Signal
The tacho signal (tachometer) is a type of signal used to measure the
rotational speed of a shaft or rotating machinery. In vibration analysis,
processing the tacho signal is crucial for obtaining characteristics related to
speed and angular position over time.

Processing the tachometer

 Set threshold value of tacho signal :

79
VIBRATION ANALYSIS TECHNIQUES 20.04

 First, an appropriate threshold value needs to be determined for the

tachometer signal.
 This threshold value serves as a reference level to distinguish the valid
tachometer signal from noise/unwanted signals.
 Selecting the right threshold value is important to ensure that only the
crucial tachometer signal pulses are considered for further analysis.

Substraction tacho signal from threshold value :

 After setting the threshold value, the original tachometer signal is
subtracted from this threshold.
 This creates a new signal that has clear peaks whenever the
tachometer signal exceeds the established threshold.
 Subtracting the threshold removes unwanted noise and fluctuations,
and highlights the important pulses/peaks in the tachometer signal.

Figure 4.17 : Tacho Signal

80
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.18. Tacho Signal Analysis

Determining Signal Samples: To identify the point at which the tacho signal
reaches a threshold value s, two consecutive samples that satisfy the following
condition are selected:
x [iΔt ]< s
x [ ( i+1 ) Δt ]> s
where x [iΔt ] is the value of the signal at time t=iΔt , and Δt is the time
interval between consecutive samples.
Linear Interpolation: After finding the two points t [iΔt ] and t [ (i+1 ) Δt ] at
which the signal crosses the threshold s, linear interpolation is used to determine
the time t k when the signal equals s. The interpolation formula is expressed as
follows:
( t( i+1)−ti )
t k =t i+ ( x k −x i )
( x ( i+1)− xi )
where:
 t i and t (i+1 ) are the times at the signal samples before and after the signal
crosses the threshold,
 x i and x (i+ 1) are the values of the signal at the corresponding times,
 t k is the time when x (t )=0 (or reaches the desired threshold).

Application in Vibration Analysis: Accurately identifying the crossing points

through the threshold enables tracking the speed and phase of the rotating
system. This is important for detecting abnormal vibrations or diagnosing faults
in the system.
Processing signals from the tachometer to determine the shaft rotation speed:

81
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.19. Tacho Signal Analysis

Identifying Time Points t k in a Cycle: The signal from the tachometer
(shown in blue) is represented as a periodic waveform, with sudden increases
and decreases.
The red points mark the time points t k , which are the intersections of the
signal with a certain threshold (horizontal axis), indicating the start of each
rotation cycle of the shaft.
Determining Rotational Speed: The rotational speed of the shaft is
calculated based on the time between two consecutive points t k and t (k +1). Speed
Calculation Formula:
1
speed=
( t (k +1)−t k )
This formula helps to determine the shaft's rotational speed in revolutions
per minute or other time units, depending on the units of t .
Number of Revolutions: The number of revolutions is calculated by
subtracting 1 from the length of the vector t k :
N rev =length ( t k )−1
This indicates the total number of cycles (revolutions) recorded during the
analysis period.
5. Order Tracking
Order tracking is an important technique in vibration analysis that helps
convert vibration signals from the time domain to the order domain, allowing for
observation of vibration components corresponding to the rotational speed of the
shaft.
Order Concept: The order of an oscillation is determined by comparing its
frequency to the rotational speed of the shaft. For example, 1X oscillation has
the same frequency as the shaft's rotational speed, 2X oscillation has a frequency
double that of the rotational speed, and so on.

82
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.20. Resampling of the times

Time Resampling Process: Starting from the time vector {t1, ..., tk}, we
need to create a time vector {t*[j/θ]} such that:
¿ ¿ 1
t [ ( j−1 ) Δθ ]−t [ jΔθ]=Δt , with Δt = as a constant .
Fs
¿ ¿ 2π
θ ( t [ ( j−1 ) Δθ] ) −θ ( t [ jΔθ] )=Δθ , with Δθ= as a constant .
N
This divides the time period in the angular domain from 0 to 2 π into N
2π
equal parts, with Δθ= N .

Figure 4.21.Dividing [0 - ɷ.2π] in N points

From there, we can calculate the time vector {t*[j/θ]} corresponding to each
division of the angular domain.
Angular Interpolation Formula: If there are N points in each rotation, with T k
being the times defining each rotation of the shaft, we have:
83
VIBRATION ANALYSIS TECHNIQUES 20.04

θ=0 , Δθ , 2 Δθ ,… , ( N −1 ) Δθ
where Δθ is the angular distance between interpolated points.
Using spline or cubic interpolation, the vibration signal can be interpolated into
a series of values according to the rotation angle points, forming a signal in the
order domain.
 Process of Angular Domain Division and Signal Interpolation

Figure 4.22. Resampling signal

 Dividing the Angular Domain:

o The angular domain from 0 to 2 π is divided into N equal points, with

each point separated by an interval of Δθ .

o The angular interval is calculated as:

2π
Δθ=
N
o For example, if N=100 , then:
2π
Δθ=
100
 Creating the Time Vector:

o From the time vector t [iΔt ], create a new time vector t ¿ [ jΔθ]
corresponding to the divided angular intervals.

 Interpolating the Signal:

o Use interpolation techniques, such as spline interpolation, to calculate

the values of the signal at the new time points t ¿ [ jΔθ].

o Spline interpolation helps to produce a smoother signal, improving

accuracy in analysis.

84
VIBRATION ANALYSIS TECHNIQUES 20.04

 Results:

o The final result is a new time vector and the resampled signal values for
each angular point.

o This signal can be used for deeper analysis, such as analyzing vibration
components at various speeds and conditions, thereby helping to detect
potential issues in the rotating system early.

6. Synchronous Averaging
The objective of this section is to obtain a time vector corresponding to
constant angular increments Δθ for each rotation of the shaft. This helps
standardize the vibration signal according to a uniform time frame, thereby
allowing for a more accurate analysis of the vibration components.

Figure 4.23. Synchronous Average (OT-SA)

 Synchronous Averaging Process:

 Synchronous Averaging:

o Synchronous averaging is a method that helps reduce noise and increase

the accuracy of the main vibration components in the signal. This
technique is carried out by synchronizing and averaging the vibration
cycles over fixed time intervals (each rotation) of the shaft.

o Operating Principle: When the vibration signal is divided into

subsequences corresponding to each rotation, these sequences are

85
VIBRATION ANALYSIS TECHNIQUES 20.04

synchronized and averaged to minimize random noise, highlighting the

specific vibrations.

o Synchronous Averaging Calculation Formula: Assuming there are P

rotations and x i [k ] is the signal of the i-th rotation with k being the
sample index, synchronous averaging is calculated as follows:
ˉ P
1
x ┴ [k ]= ∑ x [k ]
P i=1 i
ˉ
where x ┴ [k ] is the synchronized average value at sample k .

Result: After synchronous averaging, the random components (noise) of

the signal will be significantly reduced, and the defined vibrations will become
clearer. This result can be used to accurately analyze vibration components at
different orders.
7. Deterministic-Random Separation (DRS):
Deterministic-random separation is a technique that helps to separate
deterministic vibrations (originating from rotating parts and harmonic motion)
from random components (due to noise and non-synchronous vibrations). This
method is often combined with the synchronous averaging technique to optimize
the vibration signal.

86
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.24. Deterministic-Random Separation (DRS)

 Operating Principle:

o The original signal is divided into segments of equal length, each

segment corresponding to one rotation cycle. Then, these segments are
averaged to obtain the deterministic vibration component, while the
random noise component will be minimized.

 Component Separation Formula:

o Assuming the signal x (t ) includes a deterministic component x d ( t ) and a

random component x r ( t ), we have:
x (t )=x d ( t ) + x r ( t )
o When synchronous averaging is performed:
P
1
x´d ( t )= ∑ x i ( t )
P i=1
o The averaged signal x´d ( t ) will primarily contain deterministic
components, while the random component will diminish over cycles.
87
VIBRATION ANALYSIS TECHNIQUES 20.04

 Result:

o The DRS method helps clarify the characteristics of deterministic

vibrations, laying the foundation for further in-depth analyses, such as
Fourier spectrum analysis, in subsequent sections.

8. Average Synchronous
 Introduction:

o This is a method used to average vibration signals over time, especially

when machinery operates at a constant speed. This method helps
eliminate noise and highlight useful components in the signal.

Figure 4.25. Synchronous Averaging of Vibration Signals.

 Constant Speed:

o When machinery operates at a constant speed, the vibration signal can

be analyzed synchronously. This means that the vibration signals can be
divided into equal time intervals and averaged over these intervals.

 Averaging Formula:

o The signal is divided into P parts, each with equal length. During each
cycle, the averaging method will be applied to minimize the impact of
noise. The averaging formula typically takes the form:

88
VIBRATION ANALYSIS TECHNIQUES 20.04

N
1
M = ∑ xk
N k=1
o Where M is the average value, N is the number of points in each cycle,
and x k is the value of the signal at point k .

After the signal has been phase-averaged, the next step is to apply the Fast
Fourier Transform (FFT) to analyze the frequency components present in the
signal. This helps to identify the main frequencies contained in the signal, as
well as their amplitudes.

Figure 4.26. Signal Analysis in Angular Domain and Temporal Domain

FFT Formula
When applying FFT to the phase-averaged signal, we can use the following
formula:
2⋅ FFT ( X syn )
FF T proxi=
avg

N¿
Where:
 FF T proxi is the result of the FFT analysis for the phase-averaged signal.

 X syn is the phase-averaged signal.

avg

 N ¿ is the number of division points per revolution.

Objectives of FFT Analysis

The objective of using FFT after averaging the signal is to:
 Clarify the important frequency components, including the fundamental
frequency (1X) and harmonics (2X, 3X, etc.).

89
VIBRATION ANALYSIS TECHNIQUES 20.04

 Facilitate the identification of issues related to machinery operation, such as

abnormalities during operation.

Images and Results

The results of the FFT analysis will be presented in graphical form,
allowing for a comparison of the frequency spectrum of the signal before and
after applying the phase-averaging method. This will highlight the
improvements in identifying frequency components and help engineers gain a
clearer understanding of the operational condition of the machinery.
 Not constant speed (run up and run down)

For the case of "run up and run down" speed variations, the Average
Synchronous method can also be applied, with some adjustments:
 The vibration signal is divided into equal time intervals, not just for the
constant speed case.
 The signal is averaged within each time interval to eliminate the effect of
noise.
 Fast Fourier Transform (FFT) is then applied to the averaged signal to
analyze the key frequency components.

Specifically:
 During the "run up" phase, the speed increases gradually, and the
vibration signal will change accordingly.
 During the "run down" phase, the speed decreases gradually, and the
vibration signal will also change.
 The Average Synchronous method will handle these acceleration and
deceleration phases by dividing the signal into equal time intervals and
taking the average.

90
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.27 : Synchronous Averaging of Vibration Signals

(Run up and run down)

o . The averaging formula typically takes the form:

N
1
Average= ∑ x k
M k=1

91
VIBRATION ANALYSIS TECHNIQUES 20.04

 Fast Fourier Transform (FFT)

Figure 4.28. Signal Analysis in Angular Domain and Temporal Domain

(Run up and run down)
 Input signal: The yellow dots represent the vibration signal data collected
during the acceleration ("run up") and deceleration ("run down") phases
of the machinery.
 FFT: The Fast Fourier Transform is performed on the vibration signal,
converting it from the time domain to the frequency domain.
 Frequency spectrum: The output of the FFT displays the key frequency
components present in the vibration signal during the "run up" and "run
down" phases.
 Peaks in the spectrum: The peaks in the frequency spectrum represent the
dominant frequencies and their corresponding amplitudes.

92
VIBRATION ANALYSIS TECHNIQUES 20.04

 Code Matlab
 CONSTANT_SPEED_20KHZ
%%
clc
clear
close all
%%
signal = load('CONSTANT_SPEED_20KHZ.mat');
proximity = signal.signal(:,1:4)/7.88e-3;
XF = signal.signal(:,5)/7.88e-3;
N = length(XF);
Fs = 20000;
t = 0:1/Fs:(N-1)/Fs;
a = {'NDEH','NDEV','DEH','DEV'};
window = round(length(XF)/100);
noverlap = round(length(XF)/200);
nfft = 2:200;
figure('name','Analytical')
for i=1:4
subplot(4,1,i)
plot(t,proximity(:,i));
grid on
title(sprintf('%s', a{i}));
ylabel('Amp. [\mum]')
end
%% 2D
figure('Name', 'Spectrogram of Proximity Probes');
for i = 1:4
subplot(2,2,i);
spectrogram(proximity(:,i), window, noverlap, nfft, Fs,
'xaxis');
title(sprintf('%s', a{i}));
xlabel('Frequency (Hz)');
ylabel('Time (s)');
yyaxis right
ylabel('Power/frequency (dB/Hz)');
colorbar('eastoutside');
end

f = (0:N/2-1)*Fs/N;
figure;
for i=1:4
subplot(2,2,i);
FFT_X = 2*fft(proximity(:,i))/N;
93
VIBRATION ANALYSIS TECHNIQUES 20.04

FFT_X(1) = FFT_X(1)/2;
X = FFT_X(1:floor(N/2));
stem(f, abs(X), 'Marker','o','LineWidth', 1);
grid on
title(sprintf('%s', a{i}));
xlabel('Frequency [Hz]'), ylabel('Amp. [\mum]');
xlim([0 200])
ylim([0 10]);
end

figure;
THR = (max(XF) + min(XF)) / 2;
XF_mean = XF - THR;
Tk = []; Nrev = 0; Tk1 = 0; speed = [];
plot(t, XF_mean);
hold on;
xlabel('Time [s]');
ylabel('XF');
title('Tachometer Signal');
xlim([0 10]);

for k = 1:length(XF_mean)-1
if(XF_mean(k)< 0 && XF_mean(k+1)>0)
Tk_pre = Tk1;
Tk1 = k + (k+1 - k)/(XF_mean(k+1)-XF_mean(k))*(-
XF_mean(k));
plot(Tk1/Fs, 0, 'color', 'k', 'marker', '.',
'markersize', 10);
plot(k/Fs, XF_mean(k), 'color', 'r', 'marker', '+',
'markersize', 5, 'linewidth', 1.4);
plot((k + 1)/Fs, XF_mean(k + 1), 'color', 'm',
'marker', '+', 'markersize', 5, 'linewidth', 1.4);
if Tk_pre >0
Nrev = Nrev+1;
Tk(Nrev)=Tk1;
speed(Nrev) = (1 / ((Tk1 - Tk_pre) / Fs)) * 60;
end
end
end

figure;
plot(Tk, speed, 'linewidth', 1.5);
xlabel('Time [s]');
ylabel('Speed [rpm]');
title('Shaft Speed');
grid on;
94
VIBRATION ANALYSIS TECHNIQUES 20.04

ylim([9100 9300])
%% Ham noi suy
Ndiv = 100;
teta = 0:2*pi:Nrev*2*pi-1;
delta_teta = 0:2*pi/Ndiv:Nrev*2*pi;
tw = interp1(teta,Tk,delta_teta,'*spline')/Fs;
x_w = interp1(t,proximity, tw ,'*Pchip');
P = Nrev;
X_syn_avg = zeros(Ndiv,4);
X_syn_avg_mean = zeros(Ndiv,4);
for j = 1:4
for i = 1:P
X_syn_avg(:,j) = X_syn_avg(:,j) + x_w(1+(i-
1)*Ndiv:i*Ndiv,j);
end
X_syn_avg(:,j) = X_syn_avg(:,j)/P;
end

for j = 1:4
X_syn_avg_mean(:,j) = X_syn_avg(:,j) -
mean(X_syn_avg(:,j));
end

freq_fft_syn = 0:Ndiv-1;
FFT_proxi = fft(X_syn_avg_mean)*2/(length(X_syn_avg_mean));
FFT_proxi(1,:) = FFT_proxi(1,:)/2;

DE = ["NDEH", "NDEV", "DEH", "DEV"];

figure
for i = 1:4
subplot(2,2,i);

stem(freq_fft_syn ,abs(FFT_proxi(:,i)),'Marker','none','linew
idth',2),grid on
xlim([0 10])
xlabel('nX [Hz]') , ylabel('[µm]') , title(DE(i))
end
theta = linspace(0,360,100);
subplot(211)
plot(theta, X_syn_avg_mean(:,1), 'b', 'linewidth', 2);
hold on, grid on,
plot(theta, X_syn_avg_mean(:,2), 'r--' , 'linewidth', 2);
ylabel('Amp [µm]'), xlabel('\Theta [\circ]')
grid on,
xlim([0 360]);
95
VIBRATION ANALYSIS TECHNIQUES 20.04

title('Non-driven Bearing')
legend('NDEH', 'NDEV')

subplot(212)
plot(theta, X_syn_avg_mean(:,3), 'b', 'linewidth', 2);
ylabel('Amp [µm]'), xlabel('\Theta [\circ]')
hold on, grid on,
plot(theta, X_syn_avg_mean(:,4), 'r--', 'linewidth', 2 );
grid on,
xlim([0 360]);
title('Driven Bearing')
legend('DEH', 'DEV')

figure('name','Direct Orbit')

subplot(121),
plot(X_syn_avg_mean(:,1),X_syn_avg_mean(:,2),'linewidth',2),
grid on, hold on
plot(X_syn_avg_mean(1,1),X_syn_avg_mean(1,2),'ro','markersize
',15,'linewidth',2)
xlabel('X [µm]'),
ylabel('Y [µm]'),
title('Non-driven Bearing')

subplot(122)
plot(X_syn_avg_mean(:,3),X_syn_avg_mean(:,4),'linewidth',2),
grid on, hold on
plot(X_syn_avg_mean(1,3),X_syn_avg_mean(1,4),'ro','markersize
',15,'linewidth',2)
xlabel('X [µm]'),
ylabel('Y [µm]'),
title('Driven Bearing')

ang = linspace(0,0.97*2*pi,1000);
for k = 2:4
x1 = abs(FFT_proxi(k,1)).*cos(ang +
angle(FFT_proxi(k,1)));
y1 = abs(FFT_proxi(k,2)).*cos(ang +
angle(FFT_proxi(k,2)));
x2 = abs(FFT_proxi(k,3)).*cos(ang +
angle(FFT_proxi(k,3)));
y2 = abs(FFT_proxi(k,4)).*cos(ang +
angle(FFT_proxi(k,4)));

figure('name',['Orbit ',num2str(k-1),'X'])

96
VIBRATION ANALYSIS TECHNIQUES 20.04

subplot(121)
plot(x1,y1,'linewidth',2),grid on, hold on
plot(x1(1),y1(1),'ro','markersize',15)
xlabel('X [µm]'),
ylabel('Y [µm]'),
title('Non-driven Bearing')

subplot(122)
plot(x2,y2,'linewidth',2),grid on, hold on
plot(x2(1),y2(1),'ro','markersize',15)
xlabel('X [µm]'),
ylabel('Y [µm]'),
title('Driven Bearing')
end

 Results

NDEH
Amp. [ m]

-460

-470
0 1 2 3 4 5 6 7 8 9 10
NDEV
-440
Amp. [ m]

-445
-450

0 1 2 3 4 5 6 7 8 9 10
DEH
Amp. [ m]

-450

-460
0 1 2 3 4 5 6 7 8 9 10
DEV
Amp. [ m]

-710

-720

-730
0 1 2 3 4 5 6 7 8 9 10

Figure 4.29. Proximity in time domain

Observation:
The amplitude signal chart from the proximity sensors in the time domain
indicates that the motor is operating at a constant speed, reflected by the regular
repeating waveforms. There are four curves representing signals from different

97
VIBRATION ANALYSIS TECHNIQUES 20.04

proximity sensors (NDEH, NDEV, DEH, DEV), showing that the shaft's
vibration levels range from 440 to 470 μm, which is considered acceptable.
However, a significant observation is the asynchrony between the signals
from the four sensors, suggesting possible imbalance or misalignment within the
system. This asynchrony requires further investigation through additional
analyses such as tachometer signals, spectral analysis, and orbit plotting to gain
a clearer understanding of the operational condition and the underlying causes of
this phenomenon. Overall, while the system appears to be stable, the detection
of asynchrony indicates a need for further investigation to ensure optimal system
performance.

Figure 4.30. Spectrogram

Observation:
 For the NDEH and DEH Sensors:

Energy is primarily concentrated in low frequencies below 30 Hz, with

relatively high amplitude vibrations. This suggests the possibility of issues such
as misalignment, dynamic imbalance, or damage to certain rotating components.
The concentration of energy in the low-frequency range may be associated with
the fundamental vibrations of the system.
 For the NDEV and DEV Sensors:

Energy is more widely distributed across higher frequencies, reaching up

to approximately 150 Hz. This may reflect the presence of more complex
vibration components, possibly due to faults or deformation in certain parts. The

98
VIBRATION ANALYSIS TECHNIQUES 20.04

broader energy distribution towards higher frequencies suggests potential issues

such as surface damage, looseness, or misalignment.

 Significant Differences in Energy Levels and Distribution Among Sensors:

The notable differences in energy levels and distribution among the

sensors indicate significant lack of synchronization within the system. This
requires thorough investigation to determine the specific causes, which may
relate to imbalance, misalignment, or damage in certain components.
 Summary:

The spectrogram charts provide valuable insights into the kinematic

condition of the system. The inconsistencies among the sensors serve as a
warning sign of potential issues that need timely investigation and resolution.
Additional analyses, such as orbit analysis and vibration analysis, should be
conducted to accurately diagnose the causes and propose appropriate corrective
measures.

Figure 4.31. FFT before Average Synchronous

Observation:
All FFT charts demonstrate a concentration of energy primarily in low
frequencies, below 50 Hz. This indicates that the system operates with vibration
components mainly in the low-frequency range.
There are significant differences in energy levels and distribution among
the sensors. Some sensors, such as NDEH and DEH, exhibit higher energy
levels compared to NDEV and DEV. For NDEH and DEH, energy is primarily
99
VIBRATION ANALYSIS TECHNIQUES 20.04

concentrated in low frequencies below 30 Hz, while NDEV and DEV show a
broader energy distribution towards higher frequencies.
The inconsistencies in energy levels and distribution among the sensors
suggest the possibility of dynamic imbalance issues within the system.
To provide more accurate assessments and recommendations, further
analyses should be conducted, such as orbit analysis and tacho signal analysis, to
gain clearer insights into the vibration characteristics of the system.

Figure 4.32. FFT before Average Synchronous

Observation:
 Stability of Shaft Speed:

The chart shows that the Tacho signal maintains a stable level, indicating
that the rotating shaft operates at a constant speed. This is important in
confirming that the system is functioning normally without unexpected speed
variations, which could lead to undesirable oscillations or equipment damage.
 Frequency of the Tacho Signal:

A shaft rotating at a constant speed generates a stable series of Tacho

pulses. The frequency of this signal is directly proportional to the shaft’s
rotational speed; therefore, a constant frequency indicates a stable shaft speed
during the measurement period. If this frequency fluctuates, it may indicate
irregularities in the rotation speed or undesirable vibration phenomena.
 Accuracy of Measurement Process:

100
VIBRATION ANALYSIS TECHNIQUES 20.04

The Tacho signal also reflects the accuracy in tracking the shaft's
rotational speed. Regular amplitude and spacing between pulses show that the
measuring device is well-calibrated and free from interference or external
influences.
 Significance in Vibration Analysis:

The Tacho signal is a key tool in vibration analysis and monitoring the
condition of rotating equipment, as it provides a reference for rotational speed.
During vibration analysis, this signal can be used to synchronize other vibration
signals with the rotational speed, thereby helping to identify and distinguish
characteristic frequencies (1X, 2X, …) of the system.

Figure 4.33. Shaft speed

Observation:
 Stability of Shaft Speed:

The chart shows that the shaft speed remains nearly constant, fluctuating
slightly around 9120 to 9160 RPM. This stability reflects high reliability in the
rotating system’s performance, with no significant speed variations. Any major
fluctuations in speed could indicate mechanical or electronic issues, such as
shaft misalignment, imbalance, or speed control errors.
 Ability to Detect Potential Anomalies:

If sudden or irregular speed variations appear, the chart will display peaks
or large oscillations. This can provide an early warning of potential issues,
helping operators detect and address them with maintenance or adjustments
before severe damage occurs. Minor speed fluctuations are generally not
101
VIBRATION ANALYSIS TECHNIQUES 20.04

concerning if they fall within the allowable range but should still be monitored
closely to identify any unusual trends.
 Efficiency of the Speed Control System:

The chart demonstrates that the speed control system is effectively

maintaining the shaft speed within a narrow range. This indicates a well-
performing control system capable of quick responses to adjust speed in the
presence of slight disturbances or other influencing factors. For rotating systems
that require constant speed (such as pumps, generators, or other industrial
equipment), maintaining this stable speed is essential to ensure the equipment's
performance and longevity.
 Impact of External Factors:

The constant shaft speed also indicates that the system is not significantly
affected by external factors, such as fluctuating loads or interference. In cases of
load changes or external disturbances (like temperature variations or vibrations
from other machines), the speed may vary, and this would be visible on the
chart. This chart helps operators monitor the impact of external factors and make
adjustments if necessary.

Figure 4.34. Oder of harmonic

Observation:

102
VIBRATION ANALYSIS TECHNIQUES 20.04

 Harmonic Analysis (nX):

The chart displays the frequency components of harmonic orders (1X, 2X,
3X, ...) corresponding to the rotational frequency of the shaft. These harmonics
represent periodic vibrations occurring at multiples of the shaft's fundamental
frequency (1X). This chart helps identify cyclic vibrations, such as harmonics
that may arise from imbalance (1X), misalignment (2X), or more complex faults
in the system (3X and above).
 Comparison of Signals from Different Sensors:

The chart shows data from NDEH, NDEV, DEH, and DEV sensors,
which are positioned at the non-drive and drive ends of the shaft in both
horizontal and vertical directions. Comparing these sensors allows for
pinpointing the origin of the vibrations. Differences in harmonic amplitudes
between sensors may indicate that the vibrations originate from specific areas on
the shaft, such as from the drive end or the non-drive end.
 Identifying Potential Mechanical Issues:

Prominent harmonics, such as 1X and 2X, are often indicators of

fundamental issues within the rotating system. Vibrations at the 1X harmonic
are typically associated with imbalance, while 2X may relate to misalignment or
axial displacement issues. High amplitudes at these harmonics may signal
mechanical problems that require attention. Higher harmonics, such as 3X or
4X, may indicate the presence of more complex vibration phenomena, such as
bearing impacts or irregular vibrations from other components.
 Assessing the Stability of the Rotating System:

Higher amplitudes of harmonics may indicate strong vibrations within the

rotating system, impacting the performance and lifespan of the equipment. This
chart can serve as an index to measure the system's stability. If harmonic
amplitudes remain low, it indicates that the system is operating stably without
severe issues.

103
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.35. Proximity signal in one revolution

Observation:
 Amplitude Oscillation Analysis of the Signals:

The graph shows the variation of the proximity signal throughout one
revolution, with distinct oscillation amplitudes at each sensor. This oscillation
amplitude indicates the displacement of the shaft at each measurement position.
For the signals from the NDEH and NDEV sensors (non-driven positions) and
DEH and DEV sensors (driven positions), the differences in amplitude between
positions indicate an uneven oscillation along the shaft, which may reflect
phenomena such as imbalance or uneven vibrations.
 Comparison of Oscillation Between Non-driven and Driven Positions:

The marked difference in amplitude and waveform between the sensors at

the driven positions (DEH, DEV) and non-driven positions (NDEH, NDEV)
may indicate the varying impact of force and torque factors along the shaft.
Greater oscillation at the driven positions could be related to issues such as
uneven loading or shaft misalignment. This may be a sign that the driven areas
experience more stress than the non-driven areas.
 Phase Analysis and Phase Shift Between Signals:

The graph reveals the phase differences between signals from the sensors,
especially between horizontal (H) and vertical (V) sensors at the same position.
This phase difference provides information about the primary direction of
oscillation of the shaft. If the signals at the sensors exhibit excessive phase shift,

104
VIBRATION ANALYSIS TECHNIQUES 20.04

it may indicate more complex issues such as shaft twisting or misalignment,

which need thorough investigation.
 Signal Variation Over One Revolution:

The graph allows for observation of signal changes through each degree
of the revolution (from 0° to 360°). By analyzing this detail, we can detect
positions with the greatest deviations during the revolution—often the points
where the shaft experiences the highest load or torque. If there is a point in the
revolution where the amplitude suddenly spikes, this may be a sign of imbalance
or internal impact within the rotating shaft.

Figure 4.36. Oder of harmonic

Observation:
 Shaft Orbit Shape and Vibration State:

The orbit chart shows the shaft's movement in a two-dimensional space,

with the horizontal axis (X) and vertical axis (Y) representing the shaft's lateral
and vertical vibrations. An elliptical or ideally circular orbit indicates stable
shaft rotation. However, if the orbit appears irregular or elongated to one side,
this may signal misalignment, imbalance, or uneven forces acting on the shaft.
 Comparing Orbit Between Non-driven and Driven Sides:

On the non-driven bearing side, the orbit typically has smaller amplitude
and less distortion, indicating that this area is less affected by the torque and
forces from the drive source. On the driven bearing side, the orbit generally
exhibits larger amplitude and may be elongated or distorted, reflecting that this
side endures higher loads. Any imbalance or fault in the drive components can
105
VIBRATION ANALYSIS TECHNIQUES 20.04

cause stronger vibrations in this area. Differences between the two orbits can
provide clues about the location and extent of factors like imbalance or
eccentricity.
 Assessing Shaft Stability and Balance:

A stable, balanced orbit (close to a circular or evenly elliptical shape)

suggests that the shaft operates under stable conditions without major balance or
alignment issues. If the orbit is elongated in a specific direction, this may
indicate imbalance. An imbalanced shaft usually creates an orbit stretched in the
direction of the weight offset, which can place significant stress on the bearings
and reduce equipment longevity.
 Indicating Complex Phenomena Like Shaft Twist or Impact:

If the orbit shows a “butterfly” shape or other asymmetric patterns, there

may be twisting or impact forces on the shaft during rotation. These are critical
warning signs, as they suggest the shaft is experiencing uneven forces or
structural deformation. This can indicate bearing damage or issues within the
drive system and requires thorough inspection to prevent further damage.

Figure 4.37. 1X Orbit

Observation:
 1X Orbit Shape and Fundamental Vibration State:

The 1X orbit shows the shaft’s movement at its fundamental frequency,

reflecting the shaft's displacement with each complete rotation. The ideal shape
for the 1X orbit is elliptical or circular, indicating balanced vibration without
106
VIBRATION ANALYSIS TECHNIQUES 20.04

uneven forces. If the 1X orbit is a regular ellipse, this suggests stable shaft
rotation without significant imbalance or misalignment. However, if the orbit is
elongated in one direction or appears as a “butterfly” shape, this may indicate
mechanical issues such as imbalance or uneven force application.
 Comparing Orbits Between Non-driven and Driven Positions:

The chart reveals clear differences in amplitude and orbit shape between
the non-driven bearing and driven bearing positions. At the driven bearing, the
orbit typically has a larger amplitude and may tend to be more distorted or
elongated. This reflects that the driven side endures higher loads and pressures,
making it more prone to issues like misalignment or stronger vibrations
compared to the non-driven side. This distinction allows operators to pinpoint
the origin of vibrations and identify areas needing more thorough maintenance
or inspection.
 Identifying Imbalance and Misalignment:

1X vibration is closely associated with imbalance, as it is the primary

vibration component reflecting the shaft's rotational frequency. If the 1X orbit is
elongated in one direction, this may signal imbalance due to an unbalanced mass
in the rotating shaft. Imbalance will cause the shaft to exhibit strong vibrations
at the 1X frequency, leading to a skewed orbit. This is an important indicator for
operators to identify and rebalance the shaft.
 Importance of 1X Vibration in Preventive Maintenance:

1X vibration is a critical parameter in preventive maintenance because it

helps detect fundamental issues like imbalance before they become severe.
Monitoring the 1X orbit can aid in early detection and planning maintenance
before the issues cause equipment damage. In industrial rotating systems, 1X
vibration is often the primary factor contributing to wear on components like
bearings and drive elements, making it essential to monitor the 1X orbit to
minimize repair costs and extend equipment life.

107
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.38. 2X Orbit

Observation:
 2X Orbit Shape and Misalignment Abnormalities:

The 2X orbit represents shaft vibrations at twice the fundamental

rotational frequency, providing important information on irregular vibrations.
The ideal shape for the 2X orbit is a smaller ellipse or circle, indicating stable
vibration at this frequency. If the 2X orbit is elongated or distorted, this can
indicate misalignment, meaning the shaft is not completely straight or properly
centered. Misalignment can create unbalanced forces in the system, increasing
stress on bearings and other drive components.
 Comparing Vibration Between Non-driven and Driven Positions:

The vibration amplitude at the driven bearing position is usually higher

than at the non-driven bearing, indicating that the driven side endures greater
loads and is more susceptible to mechanical issues. Differences in the 2X orbit
between the two positions can help pinpoint the source of irregular vibrations. If
the orbit on the driven side is more elongated, this may reflect uneven loading or
other issues in the drive area.
 Identifying Bearing Issues and Shaft Misalignment:

2X vibrations are often related to bearing issues and shaft misalignment.

If the 2X orbit has an unusual shape, such as sharp edges or asymmetry, this
may be a sign of wear or damage in the bearings. Such problems can lead to
greater vibrations at the 2X frequency as the shaft no longer rotates evenly

108
VIBRATION ANALYSIS TECHNIQUES 20.04

around the center. Notably, if a bearing fails, 2X vibrations tend to increase

suddenly, causing uneven pressure on the associated parts.
 Impact of Uneven Loading:

The 2X orbit also shows the influence of uneven loading along the shaft.
If the shaft is subjected to an eccentric force or varying load, the 2X orbit will
become distorted, with higher amplitude at certain rotational angles. Uneven
loading may stem from factors like material inconsistencies, design flaws, or
non-ideal operating conditions. These factors not only reduce equipment
lifespan but can also lead to undesirable vibrations, affecting the system's
performance.

Figure 4.39. 3X Orbit

Observation:
 3X Orbit Shape and Mechanical Issues:

The 3X orbit reflects the vibrations of the shaft at three times the
fundamental rotational frequency. The ideal shape of the 3X orbit is a small
ellipse or circle, indicating stable vibrations at this frequency. If the 3X orbit is
elongated or distorted, this may indicate issues such as misalignment, bearing
damage, or uneven loading along the shaft.
 Comparing 3X Vibrations Between Non-driven and Driven Positions:

Similar to the 2X orbit, the amplitude of 3X vibrations at the driven

bearing position is usually greater than at the non-driven bearing position. This

109
VIBRATION ANALYSIS TECHNIQUES 20.04

suggests that the driven side experiences higher loads and mechanical stresses.
Differences in the 3X orbit between these two positions can help identify the
source of mechanical issues. If the orbit on the driven side is more distorted, it
indicates greater loads and stresses in that area.
 Analyzing Changes in the 3X Orbit Over Time:

Monitoring the evolution of the 3X orbit over time is crucial. A sudden

change or an increasing trend in the amplitude of 3X vibrations may signal
impending failure in the bearings or other components. This monitoring helps to
detect potential mechanical issues early, allowing for timely preventive
maintenance actions.
 The Connection Between 2X and 3X Orbits:

Analyzing the 3X orbit in conjunction with the 2X orbit provides a more

comprehensive view of the operational condition of the rotating shaft.
Differences and correlations between the orbits at different frequencies can help
pinpoint the origin and severity of mechanical problems.

 RUN_UP_20KHZ

clear
close all
clc
Fs = 20000;
load('RUN_UP_20KHZ.mat');
proximity_data = signal(:,1:4) / 7.88e-3;
key_phasor_data = signal(:,5) / 7.88e-3;
data_length = length(key_phasor_data);
time_vector = 0:1/Fs:(data_length-1)/Fs;
titles = {'NDEH','NDEV','DEH','DEV'};

figure;
for idx = 1:4
subplot(2, 2, idx);
A1 = fft(proximity_data(:,idx));
A2 = length(A1);
B1 = abs(A1/A2);
B2 = B1(1:A2/2+1);
B2(2:end-1) = 2*B2(2:end-1);
frequency_vector = Fs*(0:(A2/2))/A2;
plot(frequency_vector, B2);
title(titles{idx});
xlabel('Frequency [Hz]')
ylabel('Amplitude')
xlim([0 200]);
ylim([0 5]);
grid on;
110
VIBRATION ANALYSIS TECHNIQUES 20.04

end

figure;
for idx = 1:4
subplot(4,1,idx)
plot(time_vector, proximity_data(:,idx))
title(titles{idx});
xlabel('Time [s]')
ylabel('Amplitude [um]')
end

figure;
titles_spectrogram = {'NDEH','NDEV','DEH','DEV'};
for idx = 1:4

subplot(2,2,idx);
window_size = round(length(key_phasor_data) / 50);
overlap_size = round(length(key_phasor_data) / 200);
nfft_range = 2:200;
spectrogram(proximity_data(:,idx), window_size, overlap_size, nfft_range, Fs, 'xaxis')
title(titles_spectrogram{idx});
ylabel('Time [s]')
xlabel('Frequency [Hz]')
end

for idx = 1:4

figure;
[S, F, T] = spectrogram(proximity_data(:, idx), window_size, overlap_size, nfft_range, Fs);
mesh(T,F,abs(S))
surf(T,F,abs(S))
waterfall(T,F,abs(S))
view(45, 30);
title(titles_spectrogram{idx}, 'FontSize', 14, 'FontWeight', 'bold', 'Color', 'red');
xlabel('Time[s]', 'FontSize', 12);
ylabel('Freq(Hz)', 'FontSize', 12);
zlabel('Power/Freq.(dB/Hz)', 'FontSize', 12);
xlim([0 50]);
ylim([0 200]);
zlim([0 12e5]);
colorbar;
hold on;
end

threshold_value = (max(key_phasor_data) + min(key_phasor_data)) / 2;

key_phasor_adjusted = key_phasor_data - threshold_value;
figure;
plot(time_vector, key_phasor_adjusted, 'b'), grid on, hold on
ylim([-300 300]), yticks(-300:100:300), grid on;
xlim([0 10]), xticks(0:1:10), grid on;
threshold_mean = (max(key_phasor_adjusted) + min(key_phasor_adjusted))/2;
positives = [];
negatives = [];
zero = [];
tk = [];

for i = 1:length(key_phasor_adjusted) - 1
if key_phasor_adjusted(i) < threshold_mean && key_phasor_adjusted(i + 1) > threshold_mean
tk = time_vector(i) + (time_vector(i + 1) - time_vector(i)) / (key_phasor_adjusted(i + 1) -
key_phasor_adjusted(i)) * (threshold_mean - key_phasor_adjusted(i));
positives = [positives; time_vector(i + 1), key_phasor_adjusted(i + 1)];
negatives = [negatives; time_vector(i), key_phasor_adjusted(i)];
zero = [zero; tk];
111
VIBRATION ANALYSIS TECHNIQUES 20.04

tk = [tk; tk];
end
end

plot(positives(:, 1), positives(:, 2), 'mo','MarkerSize', 3);

plot(negatives(:, 1), negatives(:, 2), 'ro','MarkerSize', 3);

t_k = [];
for k = 2:length(key_phasor_adjusted)
if key_phasor_adjusted(k-1) < 0 && key_phasor_adjusted(k) >= 0
ti = time_vector(k-1);
ti1 = time_vector(k);
xi = key_phasor_adjusted(k-1);
xi1 = key_phasor_adjusted(k);
tk = ti + (ti1 - ti) / (xi1 - xi) * (0 - xi);
t_k = [t_k, tk];
end
end

plot(t_k, zeros(size(t_k)), 'k*','MarkerSize',5 )

speed = diff(t_k) .^ -1 * 60;
Nrev = length(t_k);
figure;
plot(t_k(1:end-1), speed, 'b'), grid on
xlabel('Time[s]')
ylabel('Speed(RPM)')
title('Shaft Speed')

 Results

112
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.40 : Proximity in time domain

Observation:
Overall, the signals from the proximity sensors (NDEH, NDEV, DEH, DEV)
clearly demonstrate the system's operation in the run-up mode. This is reflected
by the synchronized increase in the vibration amplitudes of the signals over
time.
While there is a certain degree of disparity between the signals, the system
maintains relatively synchronized behavior throughout the acceleration process.
This indicates the system is operating in a stable manner, without any significant
signs of imbalance or misalignment

113
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.41. Spectrogram

Observation:
Overall, the spectrograms from the four sensors (NDEH, NDEV, DEH, DEV)
all indicate the system is operating in an acceleration mode. This is evidenced by
the synchronized increase in the vibration energy levels within the key
frequency bands across the spectrograms over time.

While there are some minor differences between the sensors regarding the
precise locations of the high-energy frequency bands, the signals generally
maintain relatively synchronized behavior throughout the acceleration process.
This suggests the system is operating in a stable manner, without any significant
signs of imbalance or misalignment.

Collectively, the spectra evolve in a pattern consistent with the acceleration

process, with increasing vibration energy levels across the frequency range over
time. Further in-depth analyses of the signals could provide additional detailed
information about the system's operating condition.

114
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.42 : 3D Spectrogram

115
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.43: FFT before Average Synchronous

Observation:
The FFT plots clearly exhibit the presence of prominent frequency peaks
concentrated within a specific frequency range. This suggests that the source
signals all contain a dominant frequency component residing within this
frequency band.
The synchronized positioning and magnitudes of the main frequency peaks
across the four plots indicate that the signals are generated from closely related
sources. This may be due to these signal sources being associated with the same
underlying system or operational process.
In addition to the main frequency peaks, the plots also display the emergence of
several secondary peaks at other frequencies. The appearance of these secondary
peaks could be related to other frequency components within the signals,
reflecting additional dynamic characteristics of the system.
The amplitudes of the frequency peaks are not entirely uniform across the plots,
implying that the signal sources may have differences in the relative magnitudes
of the key frequency components.
The overall similarity in the shape of the FFT plots also suggests that the source
signals share common spectral characteristics, despite some variations in the
peak amplitudes.
In summary, the FFT plots provide detailed information about the frequency
structure of the signals, revealing the synchronized and closely related nature of
the signal sources. Further analysis of these frequency characteristics could
provide additional insights into the related system or operational process.

116
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.44 : Tacho signal

Observation:
"Run-up" mode: The waveform indicates that this is the "run-up" or acceleration
phase of the system, where the rotational speed is being gradually increased.
This explains the unsteady changes observed in the Tacho signal.
• Irregular waveform pattern: The fluctuations in the frequency and amplitude of
the Tacho pulses are normal during this run-up phase. These variations reflect
the continuous changes in the rotational speed.
• Characteristics of "run-up" mode: The instabilities in the Tacho signal are
signatures of a controlled acceleration process, aimed at safely and smoothly
bringing the system up to the desired operating speed.
• Operational assessment: While the Tacho signal exhibits changes, this is
consistent with the run-up phase. The key is to closely monitor the process to
ensure the acceleration is progressing smoothly, without any abnormal
indications.
In summary, the variations observed in the Tacho signal waveform are
appropriate for the "run-up" mode you have described. This is a normal
operating state as the system is being brought up to speed in a controlled
manner. Continued monitoring is needed to ensure the system transitions safely
and stably.

117
VIBRATION ANALYSIS TECHNIQUES 20.04

Figure 4.45 : Speed

Observation:
Consistent Acceleration Trend: The graph displays a clear linear increase in
speed over time, indicating that the acceleration of the system is being executed
in a steady and controlled manner. This consistent and continuous acceleration
pattern is characteristic of a well-executed "run-up" process.
Lack of Fluctuations: Unlike the previous Tacho signal, the speed graph does
not exhibit any significant fluctuations or irregularities. The speed increases in a
smooth and linear fashion, suggesting that the rotational speed is being
effectively managed during this "run-up" phase.

118
 VIBRATION ANALYSIS TECHNIQUES 20.04

PART 3: FEM
I. Theoretical Basis
1. Physical Model

 Main Components of the Rotor System

 Shaft: This is the primary element of the rotor system, responsible for
transmitting rotational motion from the power source (such as a motor) to
other components. The shaft is designed to withstand axial and radial
forces. Its physical properties, including density, stiffness, and Young's
modulus, significantly impact the system's vibration frequencies and
modes.
 Disk: These are significant mass elements mounted on the shaft. They play
a role in increasing the system's inertia and typically experience centrifugal
forces when the rotor rotates. The disk's mass characteristics affect the
stability of the system and contribute to the mass matrix when creating the
FEM model. The mass of the disks is often uneven, which may lead to
irregular vibrations in the system.
 Coupling: The coupling connects parts of the shaft or links the shaft with
other system components (e.g., the motor). It is crucial for transmitting
rotational motion and mitigating shaft misalignment between rotor
elements. Coupling characteristics, such as stiffness and damping, affect the
system's overall vibration, especially under different operating conditions.

119
VIBRATION ANALYSIS TECHNIQUES 20.04

 Bearing: Bearings support the shaft, maintain its stability, and reduce
friction during rotation. They are positioned at specific points along the
shaft to ensure the required stability and stiffness, helping to minimize
vibrations and prevent unwanted oscillations. Each bearing is modeled
with stiffness and damping parameters, contributing to system balance and
effective vibration reduction.
 Physical Parameters of the Rotor Model
 Mass: The mass of each component (shaft, disk, and coupling) affects the
system's inertia and influences its natural vibration frequencies. Mass is
used to construct the mass matrix in the FEM model, allowing for accurate
simulation of rotor oscillations under external forces.
 Stiffness: Each element's stiffness describes its resistance to deformation
when subjected to force. The shaft, disk, and bearing each have distinct
stiffness values, which are used in the system's stiffness matrix. Stiffness
directly impacts rotor stability and vibration characteristics.
 Damping: Damping prevents oscillation build-up, helping the system reach
equilibrium after external force impact. Bearings typically have built-in
damping functions that help control vibration during operation.
 Rotational Inertia: This parameter is critical for both the disk and shaft,
affecting how easily they rotate around the axis. Higher rotational inertia
stabilizes the system but also increases sensitivity to natural oscillations.
 Interaction Between Rotor System Components
 Centrifugal Effect: When the rotor rotates, centrifugal force acts on the
disks and other shaft components, causing oscillations due to changes in
inertia. This effect is typically considered in the gyroscopic matrix when
building the FEM model.
 Gyroscopic Effect: For high-speed rotating systems, the gyroscopic effect
plays a vital role in adjusting the shaft's direction when rotational
moments change. This effect is represented through the gyroscopic matrix,
helping model vibration changes due to rotor rotation.
 Component Linkage: Couplings are used to connect the shaft with other
rotating components, allowing them to partake in the system's rotation.
Any misalignment or imbalance in coupling can affect rotor vibration,
leading to irregular and unstable oscillations.
 Vibration Analysis Impact of Components
 In the FEM model, each rotor component (shaft, disk, coupling, bearing) is
represented as a finite element with its characteristic matrices. These
matrices reflect each component's physical properties and their
interactions.

120
 VIBRATION ANALYSIS TECHNIQUES 20.04

 This model allows analysis of natural vibration frequencies and modes,

helping identify susceptible areas to damage or resonance-prone
frequencies, which require adjustments in design and maintenance.
 In summary, the physical model of the rotor system forms the foundation
for establishing finite element matrices in the FEM approach, enabling
accurate simulation of rotor dynamics and determining factors that affect
system stability, durability, and performance.
2. Finite Element Model (FEM)
Element Representative Matrices
 Mass Matrix: The mass matrix MMM represents the mass of each element in the
system. It describes the impact of mass on vibration and is determined based on
the mass and mass distribution of each element. This matrix influences the
inertia properties of the rotor.

Where:
Sm: Cross section of the j-th element
L: Length of the j-th element
ρ : Density j-th element
Dm: Diameter of the mass of the j-th element
 Rotational Inertia Matrix: When the rotor rotates, the elements experience
rotational inertia not only along the axis but also in different directions. The
rotational inertia matrix M s describes the effect of rotational inertia on the
system's vibrations, including additional impacts due to centrifugal forces.

121
 VIBRATION ANALYSIS TECHNIQUES 20.04

Where:
I m: Moment of inertia of the cross section of the j-th element
L: Length of the j-th element
ρ : Density j-th element
Dm: Diameter of the mass of the j-th element
 Gyroscopic Matrix: The gyroscopic matrix GGG models the
impact of the gyroscopic effect due to the rotor's rotational motion. The
gyroscopic effect is the change in the rotor's rotational direction caused by
external forces or centrifugal forces, and it plays a crucial role in maintaining the
system's dynamic stability. In high-speed rotating systems, the gyroscopic
matrix can significantly affect vibration frequencies and vibration modes.

Where:
I m: Moment of inertia of the cross section of the j-th element
L: Length of the j-th element
ρ : Density j-th element
Dm: Diameter of the mass of the j-th element
 Stiffness Matrix: The stiffness matrix KKK describes an
element's resistance to deformation under applied forces. The stiffness of an
element determines its ability to resist deformation and directly influences the
system's vibration frequency. In the FEM model, stiffness is distributed across
each element to reflect how each contributes to the overall durability of the

122
VIBRATION ANALYSIS TECHNIQUES 20.04

system.

Where:
I k : Moment of inertia of the cross section of the j-th element
L: Length of the j-th element
E : Young's modulus of the element j-th element
Dk : Diameter of stiffness of the j-th elemen
a : Shear coefficient of the j-th element
Gk : Shear modulus of the j-th element
v : Poisson's ratio of the element j-th element
Sk : Reduced area of the cross section of the j-th
element

Global Matrix Assembly

 Element Matrix Synthesis: After defining the mass, stiffness, damping, and
gyroscopic matrices for each element, these are assembled into global matrices.
Each global matrix includes component matrices from all elements, connecting
them to form a complete system. These matrices reflect the vibrational behavior
of the entire rotor system.

By considering all the DOFs of the rotor (4·Nnode )

 System Equation of Motion: After assembly, the rotor system's equation of

motion takes the form:

123
VIBRATION ANALYSIS TECHNIQUES 20.04

Where:

 M: Global mass matrix of the system, representing the rotor's total inertia.
 C: Global damping matrix, representing the system's energy dissipation capability.
 G: Global gyroscopic matrix, reflecting gyroscopic effects from rotational motion.
 K: Global stiffness matrix, representing the system's overall resistance to deformation.
 F: Force acting on the system (may include centrifugal, frictional, or external forces).

Bearing Contribution

 Bearing Matrix: Bearings play a critical role in maintaining rotor stability

and reducing vibration. Bearings are also modeled through stiffness and
damping matrices, which are added to the system's global matrices.

124
 VIBRATION ANALYSIS TECHNIQUES 20.04

Bearing Stiffness Matrix K Bearing:

The bearing stiffness matrix typically takes the form:

Where:

 K xx: Stiffness of the bearing along the x-axis.

 K zz: Stiffness of the bearing along the z-axis.
 K xz và K zx: Coupled stiffness between axes, describing interactions between
axes when the bearing is under load.

Bearing Damping Matrix C Bearing:

The bearing damping matrix generally has the form:

Trong đó:

 C xx : Damping of the bearing along the x-axis.

 C zz : Damping of the bearing along the z-axis.
 C xz và C zx : Coupled damping between axes, describing interactions between
axes when the bearing is under load.

Incorporation into the Global Matrix

125
VIBRATION ANALYSIS TECHNIQUES 20.04

 Code Matlab

close all
clear
clc
rho = 7850;
E = 206e9;
v = 0.3;

N_dof = 4;
P = 52;
Lt = 0.570
N = 53;
y = linspace(0, Lt, N

L = [ones(1,2)*12.5, ones(1,5)*9, ones(1,10)*8.5, ones(1,2)*12.5, ...

ones(1,10)*12, ones(1,2)*12.5, ones(1,20)*11.35, ones(1,1)*18] * 1e-
3;
Dm = [ones(1,2)*25, ones(1,5)*10, ones(1,10)*10, ones(1,2)*75, ...
ones(1,10)*10, ones(1,2)*75, ones(1,20)*10, ones(1,1)*10] * 1e-3
Dk = [ones(1,2)*10, ones(1,5)*10, ones(1,10)*10, ones(1,2)*10, ...
ones(1,10)*10, ones(1,2)*10, ones(1,20)*10, ones(1,1)*10] * 1e-3;

Bn = [8, 52];
figure();
hold on;

for j = 1:P
rectangle('Position', [y(j), -Dk(j)/2, L(j), Dk(j)], ...
'FaceColor', 'r', 'EdgeColor', 'k'
rectangle('Position', [y(j), -Dm(j)/2, L(j), Dm(j)], ...
'FaceColor', 'y', 'EdgeColor', 'k'
plot(y(j), 0, 'k.', 'MarkerSize', 15);
end

for i = 1:length(Bn)
fill(y(Bn(i)) + [0 0.01 -0.01], [0 -0.01 -0.01], 'b');
end

line([0 y(end)], [0 0], 'LineStyle', '--', 'Color', 'b');

daspect([1, 1, 1]);
xlabel('Position along the shaft (m)');
ylabel('Diameter (m)');
title('Shaft Visualization with Bearings');
grid on;

hold off?

Sm = zeros(1, length(Dm));
Im = zeros(1, length(Dm));
Ik = zeros(1, length(Dm));
Sk = zeros(1, length(Dm));
a = zeros(1, length(Dm));
Gk = E / (2 + 2*v);

for i = 1:length(Dm)
Sm(i) = Dm(i)^2 * pi / 4;

126
VIBRATION ANALYSIS TECHNIQUES 20.04

Im(i) = Dm(i)^4 * pi / 64;

Ik(i) = Dk(i)^4 * pi / 64;
Sk(i) = pi * pi * 0.25 * Dk(i);
a(i) = (12 * E * Ik(i)) / (Gk * Sk(i) * (L(i)^2));
end

Kg = zeros(N_dof*N);
Gg = zeros(N_dof*N);
Mg = zeros(N_dof*N);

for k = 1:P
M = rho*Sm(k)*L(k)/420*...
[ 156 0 0 -22*L(k) 54 0 0
13*L(k) ;
0 156 22*L(k) 0 0 54 -
13*L(k) 0 ;
0 22*L(k) 4*L(k)^2 0 0 13*L(k) -
3*L(k)^2 0 ;
-22*L(k) 0 0 4*L(k)^2 -13*L(k) 0 0
-3*L(k)^2 ;
54 0 0 -13*L(k) 156 0 0
22*L(k) ;
0 54 13*L(k) 0 0 156 -
22*L(k) 0 ;
0 -13*L(k) -3*L(k)^2 0 0 -22*L(k)
4*L(k)^2 0 ;
13*L(k) 0 0 -3*L(k)^2 22*L(k) 0 0
4*L(k)^2 ];

Ms = rho * Im(k) / (L(k) * 30) * ...

[ 36 0 0 -3*L(k) -36 0 0 -3*L(k);
0 36 3*L(k) 0 0 -36 3*L(k) 0;
0 3*L(k) 4*L(k)^2 0 0 -3*L(k) -L(k)^2 0;
-3*L(k) 0 0 4*L(k)^2 3*L(k) 0 0 -L(k)^2;
-36 0 0 3*L(k) 36 0 0 3*L(k);
0 -36 -3*L(k) 0 0 36 -3*L(k) 0;
0 3*L(k) -L(k)^2 0 0 -3*L(k) 4*L(k)^2 0;
-3*L(k) 0 0 -L(k)^2 3*L(k) 0 0
4*L(k)^2];

G = rho*Im(k)/(15*L(k))*...
[ 0 -36 -3*L(k) 0 0 36 -3*L(k) 0
;
36 0 0 -3*L(k) -36 0 0 -3*L(k)
;
3*L(k) 0 0 -4*L(k)^2 -3*L(k) 0 0
L(k)^2 ;
0 3*L(k) 4*L(k)^2 0 0 -3*L(k) -L(k)^2
0 ;
0 36 3*L(k) 0 0 -36 3*L(k) 0
;
-36 0 0 3*L(k) 36 0 0 3*L(k)
;
3*L(k) 0 0 L(k)^2 -3*L(k) 0 0 -
4*L(k)^2 ;
0 3*L(k) -L(k)^2 0 0 -3*L(k) 4*L(k)^2
0 ];

127
VIBRATION ANALYSIS TECHNIQUES 20.04

K = E * Ik(k) / ((1 + a(k)) * L(k)^3) * ...

[ 12 0 0 -6*L(k) -12 0 0
-6*L(k);
0 12 6*L(k) 0 0 -12 6*L(k)
0;
0 6*L(k) (4 + a(k)) * L(k)^2 0 0 -6*L(k)
(2 - a(k)) * L(k)^2 0;
-6*L(k) 0 0 (4 + a(k)) * L(k)^2 6*L(k) 0
0 (2 - a(k)) * L(k)^2;
-12 0 0 6*L(k) 12 0 0
6*L(k);
0 -12 -6*L(k) 0 0 12 -6*L(k)
0;
0 6*L(k) (2 - a(k))*L(k)^2 0 0 -6*L(k)
(4 + a(k)) * L(k)^2 0;
-6*L(k) 0 0 (2 - a(k)) * L(k)^2 6*L(k) 0
0 (4 + a(k)) * L(k)^2 ];

j1 = ((k-1)*N_dof + 1);
j8 = ((k+1)* N_dof);
Mg(j1:j8 , j1:j8) = Mg(j1:j8 , j1:j8) + M + Ms;
Gg(j1:j8 , j1:j8) = Gg(j1:j8 , j1:j8) + G;
Kg(j1:j8 , j1:j8) = Kg(j1:j8 , j1:j8) + K;
end
Bk = zeros(N_dof,N_dof,length(Bn));
Bc = zeros(N_dof,N_dof,length(Bn));
Kxx = 0.5e8 ; Kxz = 0 ; Kzx = 0 ; Kzz = 0.09e8;
Cxx = 4.46102e7 ; Cxz = 0 ; Czx = 0 ; Czz = 4.46102e7;
K_bearing = [Kxx Kxz 0 0 ;
Kzx Kzz 0 0 ;
0 0 0 0 ;
0 0 0 0 ];
C_bearing = [Cxx Cxz 0 0 ;
Czx Czz 0 0 ;
0 0 0 0 ;
0 0 0 0 ];
for iBn = 1:length(Bn)
Bk(:,:,iBn) = K_bearing;
Bc(:,:,iBn) = C_bearing;
end
alpha = 0;
beta = 0;
Cg = alpha*Mg + beta*Kg;
for n = 1:length(Bn)
j1 = (Bn(n) - 1) * N_dof + 1;
j2 = Bn(n)*N_dof;
Kg(j1:j2 ,j1:j2) = Kg(j1:j2 , j1:j2) + Bk(:,:,n);
Cg(j1:j2 ,j1:j2) = Cg(j1:j2 , j1:j2) + Bc(:,:,n);
end
[V, D] = eig(M \ K);

[XX0, lam0] = eig(Mg\ Kg);

128
 VIBRATION ANALYSIS TECHNIQUES 20.04

[wn0, indexes] = sort(sqrt(diag(lam0)));

XX0 = XX0(:, indexes);

disp('Global system mode shapes and natural frequencies (Hz):');

for ii = 1:10 ;
disp(['Mode ', num2str(ii), ': Frequency = ',
num2str(wn0(ii)/(2*pi)), ' Hz']);

end

N_mode = 10;

x = linspace(0, Lt, N);

figure();
for j = 1:N_mode
subplot(ceil(N_mode/2), 2, j)
zz = real(XX0(1:N_dof:end, j) ./ max(XX0(:, j)));
xx = real(XX0(2:N_dof:end, j) ./ max(XX0(:, j)));

plot(x(1:length(zz)), zz,'-*', 'LineWidth', 1);

hold on;
plot(x(1:length(xx)), xx,'-o' ,'LineWidth', 1

ylabel('Amp');
title(['Mode ' num2str(j),'-', 'Freq = ', num2str(wn0(j)/(2*pi))]);
xlim([0 Lt]);
grid on;
end

 Result :
Global system mode shapes and natural frequencies (Hz):
Mode 1: Frequency = 28.8182 Hz
Mode 2: Frequency = 28.8563 Hz
Mode 3: Frequency = 125.6697 Hz
Mode 4: Frequency = 126.0075 Hz
Mode 5: Frequency = 348.3833 Hz
Mode 6: Frequency = 358.8712 Hz
Mode 7: Frequency = 413.0562 Hz
Mode 8: Frequency = 432.3652 Hz
Mode 9: Frequency = 671.3488 Hz
Mode 10: Frequency = 696.1613 Hz

129
 VIBRATION ANALYSIS TECHNIQUES 20.04

Shaft Visualization with Bearings

0.04

0.02
Diameter (m)

-0.02

-0.04
0 0.1 0.2 0.3 0.4 0.5 0.6
Position along the shaft (m)

Figure 4.46: Diagram of the Shaft Element Shapes

Figure 4.47: Vibration Diagram of the Shaft

Observation:

Mode #1 and Mode #2 - Frequency 28.8 Hz:

 Both modes have the same frequency of 28.8 Hz, but they exhibit different
oscillation patterns.
 This is the most basic mode of vibration, likely representing the bending mode
of the shaft.
 The shaft has a slight curvature from end to end. This is the fundamental mode
and is also the most susceptible to resonance when excited at the frequency of
28.8 Hz.
 The red and blue points illustrate the out-of-phase oscillation along the length
of the shaft. The points oscillate with the highest amplitude at both ends and
decrease towards the center of the shaft.

Mode #3 and Mode #4 - Frequency 125.8 Hz:

130
VIBRATION ANALYSIS TECHNIQUES 20.04

 Modes 3 and 4 have a higher frequency of 125.8 Hz, corresponding to the

second mode of vibration of the shaft.
 This mode introduces a node (a point of no motion) in the middle, creating a
region of out-of-phase oscillation.
 This vibration pattern is more complex than Modes 1 and 2, as the amplitude at
the nodes is lower, but the amplitude at other points is greater.

Mode #5 and Mode #6 - Frequencies 347.1 Hz and 357.5 Hz:

 In these modes, the shaft exhibits more nodes (at least two nodes at different
positions along the shaft).
 The higher frequency leads to a more complex oscillation pattern, while the
amplitudes at the node points remain close to zero.
 This indicates that the vibrations in these modes have a phase shift in the
opposite direction across different regions of the shaft.

Mode #7 and Mode #8 - Frequencies 410.8 Hz and 429.8 Hz:

 As the frequency increases, the number of nodes continues to rise.

 The shaft has additional nodal points, leading to a more complex vibration
pattern. This means that the shaft's oscillation has more phase transition areas,
resulting in a mixed oscillation mode along its length.

Mode #9 and Mode #10 - Frequencies 662.4 Hz and 686.3 Hz:

 These are the highest vibration modes on the graph, with very high frequencies
and numerous nodes.
 The significant increase in the number of nodes indicates that the vibrations of
the shaft in these modes are complex and segmented into various regions, each
exhibiting opposite phase oscillations.
 In these modes, the oscillation amplitude of the shaft is highly sensitive to the
node points, and the shaft appears to be divided into many small oscillating
segments.

131