Different graphs returned for the same functions

6 次查看（过去 30 天）
Janco Tunggal
Janco Tunggal 2021-9-29
回答： Nipun 2024-4-15
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.

请先登录，再回答此问题。

回答（1 个）

Nipun
Nipun 2024-4-15
Hi Janco,
I understand that you intend to know why the two MATLAB codes produce different plots.
After a scrupulous look at both the codes, the second code seems like an iterative solver and very different from the first one.
The calculating logic, the critical speeds and amplitudes all seem to be different. On a high level, even the length of the codes is different - it seems like the second is a subset of first.
These could be the reasons for the dissimilarity in the outputs.
Hope this helps.
Regards,
Nipun

类别

在 Help Center 和 File Exchange 中查找有关 Graph and Network Algorithms 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by