Matlab Code for solution by Torsional equivalence method
%% THIS PROGRAM DETERMINES THE NATURAL FREQUENCIES AND
LOCATIONS OF NODE POINTS FOR ALL
%% MODES FOR THREE ROTOR GEARED TORSIONAL VIBRATING
SYSTEM(WITH GEAR DRIVE)
clc
%% ENTER THE VALUES OF INERTIAS
Ia=input ('Enter the first inertia value in kg-m^2. ');
Ic=input ('Enter the third inertia value in kg-m^2. ');
Ip=input ('Enter the inertia value of pinion in kg-m^2. ');
Ig=input ('Enter the inertia value of Gear in kg-m^2. ');
% Enter value of equivalent lengths
L1=input ('Enter length between 1st roter and 2nd roter on actual system in m- ');
L2=input ('Enter length of first part of stepped shaft between 2nd roter and 3rd roter on actual
system in m- ');
L3=input ('Enter length of second part of stepped shaft between 2nd roter and 3rd roter on
actual system in m- ');
L4=input ('Enter length of third part of stepped shaft between 2nd roter and 3rd roter on
actual system in m- ');
% ENTER THE VALUES OF DIAMETERS
d1=input ('Enter the diameter of the part between 1st roter and 2nd roter on actual system in
m. ');
d2=input ('Enter the diameter of the first part of the stepped shaft between 2nd roter and 3rd
roter in m. ');
d3=input ('Enter the diameter of the second part of the stepped shaft between 2nd roter and
3rd roter in m. ');
d4=input ('Enter the diameter of the third part of the stepped shaft between 2nd roter and 3rd
roter in m. ');
% d1=input ('Enter; the equivalent shaft dia in m. ');
N1=input ('Enter the speed of the driver machine in rpm. ');
N2=input ('Enter the speed of the driven machine in rpm. ');
g=N1/N2
LE1=L1
LE2=(g^2*L2*(d1/d2)^4)+(g^2*L3*(d1/d3)^4)+(g^2*L4*(d1/d4)^4)
LE=LE1+LE2
Ice=Ic/g^2
Ibe=Ip+(Ig/g^2)
x=Ice/Ia;
A=Ice*(1+x)+Ibe*x
B=-Ibe*LE1-Ibe*x*LE2-Ice*LE
C=Ibe*LE1*LE2
Lc1=(-B+sqrt(B*B-4*A*C))/(2*A)
Lc2=(-B-sqrt(B*B-4*A*C))/(2*A)
La1=x*Lc1
La2=x*Lc2
G=84*10^9;
J=(pi/32)*d1^4;
wn2=sqrt((G*J)/(Lc1*Ice));
wn3=sqrt((G*J)/(Lc2*Ice));
disp('wn2=');
disp(wn2)
disp('wn3=');
disp(wn3)
fn2=wn2/(2*pi);
fn3=wn3/(2*pi);
% disp('fn2=');
% disp(fn2)
NPA1= sprintf('The distance of first node point form end A in m = %f m\n\n',La1);
disp(NPA1);
NPA2= sprintf('The distance of second node point form end A in m = %f m\n\n',La2);
disp(NPA2);
NPC1= sprintf('The distance of first node point form end C in m = %f m\n\n',Lc1);
disp(NPC1);
NPC2= sprintf('The distance of second node point form end C in m = %f m\n\n',Lc2);
disp(NPC2);
OP2= sprintf('The second torsional Natural frequency in Hz = %f Hz\n\n',fn2);
disp(OP2);
% disp('fn3=');
% disp(fn3)
OP3= sprintf('The third torsional Natural frequency in Hz = %f Hz\n\n',fn3);
disp(OP3);
Matlab code for solution by Eigen value method
clc
%% Enter Values of Mass
I1=input('Enter the first mass moment of inertia value in kg-m^2. ');
I2=input('Enter the second mass moment of inertia value in kg-m^2. ');
I3=input('Enter the third mass moment of inertia value in kg-m^2. ');
%% Enter Values of Stiffness
kt1=input('Enter the first Torsional stiffness value in N-m/radian. ');
kt2=input('Enter the second Torsional stiffness value in N-m/radian. ');
kt3=input('Enter the third Torsional stiffness value in N-m/radian. ‘)
%% Code
I=[1 0 0; 0 1 0; 0 0 1]; % IDENTITY MATRIX
I1=[I1 0 0; 0 I2 0; 0 0 I3] %MASS MOMENT OF INERTIA MATRIX
KT=[(kt1+kt2) -kt2 0; -kt2 (kt2+kt3) -kt3; 0 -kt3 kt3]% TORSIONAL STIFFNESS
MATRIX
II=inv(I1) % INVERSE OF INERTIA MATRIX
D= II*KT % DYNAMIC MATRIX
[V E] =eig(D) % Finding Eigen Values and Eigen Vectors
Omega3= sqrt(E(1,1)) % First Natural Frequency
Omega1= sqrt(E(2,2)) % Scond Natural Frequency
Omega2= sqrt(E(3,3)) % Third Natural Frequency
fn1=Omega1/(2*pi);
fn2=Omega2/(2*pi);
fn3=Omega3/(2*pi);
% disp('fn2=');
% disp(fn2)
OP1= sprintf('The first torsional Natural frequency in Hz = %f Hz',fn1);
disp(OP1);
OP2= sprintf('The second torsional Natural frequency in Hz = %f Hz',fn2);
disp(OP2);
OP3= sprintf('The third torsional Natural frequency in Hz = %f Hz',fn3);
disp(OP3);
Mode1=[((V(1,1)/V(1,1)));((V(2,1)/V(1,1)));((V(3,1)/V(1,1)))] %First Mode
Mode2=[((V(1,2)/V(1,2)));((V(2,2)/V(1,2)));((V(3,2)/V(1,2)))] %Second Mode
Mode3=[((V(1,3)/V(1,3)));((V(2,3)/V(1,3)));((V(3,3)/V(1,3)))] %Third Mode
%% Plots
X=[1;2;3];
figure(1)
plot(X,Mode1,'-o',X,Mode2,'-o',X,Mode3,'-o');grid
x2= 2;
y2= Mode1(2,1);
x3=3;
y3= Mode1(3,1);
txt1=['',num2str(y2)];
txt2=['',num2str(y3)];
text(x2,y2,txt1)
text(x3,y3,txt2)
x2= 2;
y2= Mode2(2,1);
x3=3;
y3= Mode2(3,1);
txt3=['',num2str(y2)];
txt4=['',num2str(y3)];
text(x2,y2,txt3)
text(x3,y3,txt4)
x2= 2;
y2= Mode3(2,1);
x3=3;
y3= Mode3(3,1);
txt5=['',num2str(y2)];
txt6=['',num2str(y3)];
text(x2,y2,txt5)
text(x3,y3,txt6)
title('Mode shapes')
xlabel('Distances between mass stations')
ylabel('Modal vectors')
hold off
