Different results compared to hand calculation?

Sorry for asking a question as simple as this but can anybody tell me why the program returns different values of T (torque) from what I get from a calculator?

This is my take on Holzer's method used to iterate natural frequencies of torsional vibration in engine shaft.

clear all
%Parameters [Inertia (kgm^2) and Flexibility (rad/Nm)]
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; 

%Vibration analysis
for i=1:1000
 	%w(i)= 0.2*(i-1);
	w(i)= sqrt(580.6) %Trial frequency
    %Initial conditions
    a(1,i) = 1;                           %Amplitude (assume)
    T(1,i) = J(1)*a(1,i)*(w(i)^2);        %Torque
    S(1,i) = T(1,i);                      %Residual Torque on 1st member
    for n=2:11                            %Members 2 to 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);    %a(n) = (sum of preceding res torques)*Flexibility
    end
end

for w^2 = 580.6,

If i was to evaluate T by hand or directly from the command window:

T(3) = J(3)*a(3,i)*w(i)^2 = 5.5381e+06

T(4) = J(4)*a(4,i)*w(i)^2 = 5.5461e+06 and these are the correct values of T.

but if I typed T to evaluate all the torques in members n = 1:11, I get different results from T(3) onwards i.e. T(3) = 5.565e+06, T(4) = 5.5132e+06 even though the values of 'a' are all correct (referring to the manufacturer's calculation booklet). This causes a huge discrepancy in the final calculations, how should I correct this?

[KSSV]

Check are you using correct indices while checking.

[N/A]

Do you mean checking with the find() function?Sorry I'm not really good at this.

Answer 1

[Mike Croucher]

I don't know the algorithm you are implementing but I can see some things in your code that look strange to me and can explain why the code gives different behaviour to what you see at the command line.

First, the loop over the variable i. You seem to be creating the exact same results 1000 times. We can see this with the variable T. Once you have run your code, it is a matrix with 11 rows and 1000 columns:

size(T)
ans =
          11        1000
          

The results in each column are identical

Here are the first few coolumns of row 1.

          >> T(1,1)
          ans =
          6009210
          >> T(1,2)
          ans =
          6009210
          >> T(1,3)
          ans =
          6009210
          

If you look at all of the columns of this row you'll see they are identical. I'll only show the first few columns of the output to save space here

          >> T(1,:)
          ans =
          Columns 1 through 7
          6009210     6009210     6009210     6009210     6009210     6009210     6009210         

Similarly for all of the other rows. Without seeing the formula you are trying to implement, I can't guess what your intent was with the loop over i but as the code currently stands, i is redundant.

On to why the code gives different results to the what you find when you type at the command line. Inside your loop you have this. I added a couple of extra comments that hopefully explains whats going on

a(n,i) = a(n-1,i) - (J(n-1)*E(n-1)*a(n-1,i)*w(i)^2);  %OLD a
T(n,i) = J(n)*a(n,i)*w(i)^2;  % T is computed from the a(n,i) given in the line above.
S(n,i) = S(n-1,i) + T(n,i);
a(n,i) = a(n-1,i) - S(n-1,i)*E(n-1); % You compute a NEW value for a(n,i)

When at the command line you ask for,e.g. T(3,i), you get the value that was computed from what I've called OLD a(3,i). But when you subsequently perform this computation at the command line:

T(3,i) = J(3)*a(3,i)*w(i)^2

You are using the NEW value of a(3,i) which was computed in the final line of your loop.

Hope this helps

[N/A]

Thank you for your answer. I understand that i is redundant but i'm leaving it as is for now because I put a trial frequency "580.6" instead of using w(i)=0.2*(i-1). The program was initially made to plot (w,S(11)) for different values of w, so that I can infer the natural frequency where the curve intersects the X axis.

If that's the case, how should I calculate T from the new value of a? I tried a few other ways but then there's "index exceeds array bounds" errors.

[Mike Croucher]

I guess you could do this in your main loop:

   for n=2:11                            %Members 2 to 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);    %a(n) = (sum of preceding res torques)*Flexibility
        T(n,i) = J(n)*a(n,i)*w(i)^2;
   end
   

The resulting T are what you say is correct when done by hand. Does that help?

[N/A]

Wow

Yes it does, thank you very much!!