Different graphs returned for the same functions

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.

Answer 1

[Nipun]

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