Different graphs returned for the same functions

1 view (last 30 days)
This program is my take on Holzer's method used to find the natural frequency of torsional vibration, as well as plotting the torsional stress in a shaft with 11 rotating masses.
Using the torsional stress formula (16*T/pi*D^3) where T is the accumulated torque at each member plotted against the critical speed converted from the natural frequency, I found that the graph returned looks very different from how it should be, so I made another script file which returns the graph i'm looking for. Given that the codes for the "Torsional stress" section are exactly the same, i'm wondering how the program returns different values of stress for the same inputs.
The complete version:
clear all
%Required Parameters=======================================================
% Inertia (kg.m^2)
J(1) = 10350;
J(2) = 9668;
J(3) = 9668;
J(4) = 9668;
J(5) = 9668;
J(6) = 9668;
J(7) = 9668;
J(8) = 2525;
J(9) = 20190;
J(10)= 399;
J(11)= 51800;
%Flexibility (Rad/Nm)
E(1)= 0.6560*10^-9;
E(2)= 0.8140*10^-9;
E(3)= 0.8020*10^-9;
E(4)= 0.8300*10^-9;
E(5)= 0.8050*10^-9;
E(6)= 0.7670*10^-9;
E(7)= 0.5680*10^-9;
E(8)= 0.3650*10^-9;
E(9)= 40.680*10^-9;
E(10)= 9.927*10^-9;
E(11)= 0;
%Shaft diameters(mm)
Dout(1)= 720;
Dout(2)= 720;
Dout(3)= 720;
Dout(4)= 720;
Dout(5)= 720;
Dout(6)= 720;
Dout(7)= 720;
Dout(8)= 720;
Dout(9)= 395;
Dout(10)=550;
Din(1)=150;
Din(2)=150;
Din(3)=150;
Din(4)=150;
Din(5)=150;
Din(6)=150;
Din(7)=150;
Din(8)=150;
Din(9)= 0;
Din(10)= 0;
%Vibration Analysis========================================================
for i=1:1000 % Initial conditions
w(i)= 0.15*(i-1); % Frequency
a(1,i) = 1; % Amplitude
T(1,i) = J(1)*a(1,i)*(w(i)^2); % Torque
S(1,i) = T(1,i); % Sum of T
for n=2:11
a(n,i) = a(n-1,i) - (J(n-1)*E(n-1)*a(n-1,i)*w(i)^2);
T(n,i) = J(n)*a(n,i)*w(i)^2;
S(n,i) = S(n-1,i) + T(n,i);
a(n,i) = a(n-1,i) - S(n-1,i)*E(n-1);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Natural frequency of 2 nodes must be found to advance %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Deflections===============================================================
%1st Node
for wn1 = 26.04 %1st node natural freq (from above)
a1(1) = 1; %Initial conditions
T(1) = J(1)*a1(1)*(wn1^2);
S(1) = T(1);
%c1(1)= ((wn1*60)/(2*pi)); %crit. speed conversion (rad/s to rpm)
for n=2:11 %rot. members 2 to 11
a1(n) = a1(n-1) - (J(n-1)*E(n-1)*a1(n-1)*wn1^2);
T(n) = J(n)*a1(n)*wn1^2;
S(n) = S(n-1) + T(n);
a1(n) = a1(n-1) - S(n-1)*E(n-1);
end
end
%2nd Node
for wn2 = 130.5 %2nd node natural freq
a2(1) = 1; %initial conditions
T(1) = J(1)*a2(1)*(wn2^2);
S(1) = T(1);
%c2(1)= ((wn2*60)/(2*pi)); %crit. speed conversion (rad/s to rpm)
for n=2:11 %rot. members 2 to 11
a2(n) = a2(n-1) - (J(n-1)*E(n-1)*a2(n-1)*wn2^2);
T(n) = J(n)*a2(n)*wn2^2;
S(n) = S(n-1) + T(n);
a2(n) = a2(n-1) - S(n-1)*E(n-1);
end
end
%Torsional Stress==========================================================
%(to be plotted against crit. speeds of each vib. node)
%1st Node
for w1 = 26.04 %1st node vib nat freq (wn1)
c1(1)=248.66; %1st order crit. speed
amp1(1) = 1;
T(1) = J(1)*amp1(1)*(w1^2);
S(1) = T(1);
s1(1) = ((16*T(1))/(pi*(Dout(1)^3 - Din(1)^3))); %Stress @ member 1
%c1(1)= ((wn1*60)/(2*pi)); %crit. speed conversion (rad/s to rpm)
for n=2:10 %11 members: 10 "springs"
amp1(n) = amp1(n-1) - (J(n-1)*E(n-1)*amp1(n-1)*w1^2);
T(n) = J(n)*amp1(n)*w1^2;
S(n) = S(n-1) + T(n);
amp1(n) = amp1(n-1) - S(n-1)*E(n-1);
s1(n) = ((16*S(n))/(pi*(Dout(n)^3 - Din(n)^3)));%Stress @ member n
for z = 2:20 %crit. speeds up to 20th order
c1(z)= c1(z-1) - (c1(z-1)/z);
T(n,z) = J(n)*c1(z)*amp1(n);
S(n,z) = S(n-1,z-1)+T(n,z);
s1(n,z)=((16*(S(n,z)*c1(z)*amp1(n)))/(pi*(Dout(n)^3 - Din(n)^3)));
end
end
end
%2nd Node
for w2 = 130.5 %2nd node vib nat freq (wn2)
c2(1)=1246.18; %2nd order crit. speed
amp2(1) = 1;
T(1) = J(1)*amp2(1)*(w2^2);
S(1) = T(1);
%c2(1)= ((w1*60)/(2*pi)); %crit. speed conversion (rad/s to rpm)
s2(1) = ((16*T(1))/(pi*(Dout(1)^3 - Din(1)^3)));
for n=2:10
amp2(n) = amp2(n-1) - (J(n-1)*E(n-1)*amp2(n-1)*w2^2);
T(n) = J(n)*amp2(n)*w2^2;
S(n) = S(n-1) + T(n);
amp2(n) = amp2(n-1) - S(n-1)*E(n-1);
s2(n) = ((16*S(n))/(pi*(Dout(n)^3 - Din(n)^3)));
for z = 2:20 %crit. speeds up to 20th order
c2(z)= c2(z-1) - (c2(z-1)/z);
T(n,z) = J(n)*c2(z)*amp2(n);
S(n,z) = S(n-1,z-1)+T(n,z);
s2(n,z)=((16*(S(n,z)*c2(z)*amp2(n)))/(pi*(Dout(n)^3 - Din(n)^3)));
end
end
end
%PLOTS=====================================================================
figure(1) %natural frequency @S(11)=0
plot (w,S(11,:));
grid on
title('Natural Frequency');
xlabel('Frequency (Rad/s)');
ylabel('Residual Torque');
figure(2) %amplitudes
plot(a1) %first node deflection
hold on
plot(a2,'--') %second node deflection
grid on
title('Amplitude (deflection)');
xlabel('Member(n)');
ylabel('Rad(s)');
figure(3) %first node torsional stress
plot(c1,s1,'-')
hold on
grid on
title('Stresses of Node 1 Vibration');
xlabel('Critical Speed (RPM)');
ylabel('Stress (MPa)');
figure(4) %second node torsional stress
plot(c2,s2,'-')
hold on
grid on
title('Stresses of Node 2 Vibration');
xlabel('Critical Speed (RPM)');
ylabel('Stress (MPa)');
here, the graph returned for the 1st node stress looks like this:
whereas for the trial script i made (exact same code):
clear all
%parameters
J(1) = 10350;
J(2) = 9668;
J(3) = 9668;
J(4) = 9668;
J(5) = 9668;
J(6) = 9668;
J(7) = 9668;
J(8) = 2525;
J(9) = 20190;
J(10)= 399;
J(11)= 51800;
E(1)= 0.6560*10^-9;
E(2)= 0.8140*10^-9;
E(3)= 0.8020*10^-9;
E(4)= 0.8300*10^-9;
E(5)= 0.8050*10^-9;
E(6)= 0.7670*10^-9;
E(7)= 0.5680*10^-9;
E(8)= 0.3650*10^-9;
E(9)= 40.680*10^-9;
E(10)= 9.927*10^-9;
E(11)= 0;
Dout(1)= 720;
Dout(2)= 720;
Dout(3)= 720;
Dout(4)= 720;
Dout(5)= 720;
Dout(6)= 720;
Dout(7)= 720;
Dout(8)= 720;
Dout(9)= 395;
Dout(10)=550;
Din(1)=150;
Din(2)=150;
Din(3)=150;
Din(4)=150;
Din(5)=150;
Din(6)=150;
Din(7)=150;
Din(8)=150;
Din(9)= 0;
Din(10)= 0;
for w1 = 26.04 %1st node vib nat freq (trial)
c1(1)=248.66; %1st order crit. speed
amp1(1) = 1; %initial amplitude/deflection
T(1) = J(1)*amp1(1)*(w1^2);
S(1) = T(1);
%c1(1)= ((w1*60)/(2*pi));
s1(1) = ((16*T(1))/(pi*(Dout(1)^3 - Din(1)^3))); %stress gov. eqn (node1)
for n=2:10
amp1(n) = amp1(n-1) - (J(n-1)*E(n-1)*amp1(n-1)*w1^2);
T(n) = J(n)*amp1(n)*w1^2;
S(n) = S(n-1) + T(n);
amp1(n) = amp1(n-1) - S(n-1)*E(n-1);
s1(n) = ((16*S(n))/(pi*(Dout(n)^3 - Din(n)^3)));
for z = 2:20
c1(z)= c1(z-1) - (c1(z-1)/z);
T(n,z) = J(n)*c1(z)*amp1(n);
S(n,z) = S(n-1,z-1)+T(n,z);
s1(n,z)=((16*(S(n,z)*c1(z)*amp1(n)))/(pi*(Dout(n)^3 - Din(n)^3)));
end
end
end
for w2 = 122.76 %2nd node vib nat freq (trial)
c2(1)=1172.3; %1st order crit. speed
a2(1) = 1; %initial amplitude/deflection
T(1) = J(1)*a2(1)*(w2^2);
S(1) = T(1);
%c2(1)= ((w2*60)/(2*pi));
s2(1) = ((16*T(1))/(pi*(Dout(1)^3 - Din(1)^3))); %stress gov. eqn (node2)
for n=2:10
a2(n) = a2(n-1) - (J(n-1)*E(n-1)*a2(n-1)*w2^2);
T(n) = J(n)*a2(n)*w2^2;
S(n) = S(n-1) + T(n);
a2(n) = a2(n-1) - S(n-1)*E(n-1);
s2(n) = ((16*S(n))/(pi*(Dout(n)^3 - Din(n)^3)));
for z = 2:20 %crit. speed from 2nd to 20th order
c2(z)= c2(z-1) - (c2(z-1)/z);
T(n,z) = J(n)*c2(z)*a2(n);
S(n,z) = S(n-1,z-1)+T(n,z);
s2(n,z)=((16*(S(n,z)*c2(z)*a2(n)))/(pi*(Dout(n)^3 - Din(n)^3)));
end
end
end
figure(5)
plot(c1,s1,'-')
title('Stresses of Node 1 Vibration');
xlabel('Critical Speed (RPM)');
ylabel('Stress (MPa)');
hold on
figure(6)
plot(c1,s2,'-')
title('Stresses of Node 2 Vibration');
xlabel('Critical Speed (RPM)');
ylabel('Stress (MPa)');
hold on
grid on
The plot looks like this (what i'm looking for):
Sorry i'm not really good at this.

Answers (0)

Categories

Find more on Vibration Analysis in Help Center and File Exchange

Products


Release

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!