Why is this happening?

Hi Guys,

I firstly want to apologise, because I am going to post a lot of code, but hopefully you won't need to go through all of it.

I am a post grad Civil Engineer and I am using a method called The Direct Stiffness Method to solve for a portal gable frame. I don't know if it will help, but here is what it looks like:

The thing that makes this tricky to solve is that of the non-prismatic eaves haunch (seen in the pink). To model this I have had to divide the section into a number slices that in reality would give a staggered apperance. I have divided it into six equal segements, each with a smaller cross-section in the directon towards the centre of the frame.

Anyway, back to the programming side of things, I am getting some strange results when I do my computations. Please see the following code, but in particular go to where I have lablled 'Relevant Code'.

clear, clc, close all

        % Haunch length
        L = 3.0165;
        % Length of far end of haunch and sharp end
        h1 = 754/1000;
        h5 = 453.4/1000;
        % Haunch angle at column intersection
        theta = atand(L/(h1-h5));
        % Cut haunch into six segments
        d = L/6;
        % increment that the haunch decreases
        hs = d/tand(theta);
        % for loop to calculate height of sections
        for i = 1:6
               H(i) = h1 - hs;
               h1 = H(i);
        endfor
       % properties of rafter sections
        b = 189.8/1000;
        tf = 12.7/1000;
        tw = 8.5/1000;
       % calculating of second moment of area
        for i = 1:6
               Iy(i) = 1/12.*b.*H(i).^3 - 1/12.*(b-tw).*(H(i) - 2*tf).^3;
        endfor
        %calculating cross-sectional area
        for i = 1:6
               Ac(i) = H(i).*b -(b - tw)*(H(i) - 2*tf);
        endfor
        Iy_start = 9.49961E-04;
        Ac_start = 0.01112;
        Iy = [Iy_start Iy(1) Iy(2) Iy(3) Iy(4) Iy(5) Iy(6)];
        Ach = [Ac_start Ac(1) Ac(2) Ac(3) Ac(4) Ac(5) Ac(6)];

        for i = 1:6
        Ix(i) = (Iy(i) + Iy(i+1))/2
        Ach(i) = (Ach(i) + Ach(i+1))/2
        end

        E = 2.1E+08;
        Ic = 8.74999E-04;
        Acr = 0.01442;
        Lc = 8;
        n = 51;
        alpha = 6;
        matrix = zeros(6,n);
        Lr = 15/(cosd(alpha)) - 3.0165;
        Ir = 2.94331E-04;
        Ar = 0.00856;

        % Local stiffness matrix for columns
        Kc = [1 0 0 -1 0 0
        0 12 6*Lc 0 -12 6*Lc
        0 6*Lc 4*Lc^2 0 -6*Lc 2*Lc^2
        -1 0 0 1 0 0
        0 -12 -6*Lc 0 12 -6*Lc
        0 6*Lc 2*Lc^2 0 -6*Lc 4*Lc^2];
        % Local stiffness matrix for rafters
        Kr = [1 0 0 -1 0 0
        0 12 6*Lr 0 -12 6*Lr
        0 6*Lr 4*Lr^2 0 -6*Lr 2*Lr^2
        -1 0 0 1 0 0
        0 -12 -6*Lr 0 12 -6*Lr
        0 6*Lr 2*Lr^2 0 -6*Lr 4*Lr^2];
        % Local stiffness matrix for eaves haunch
        Kh = [1 0 0 -1 0 0
        0 12 6*d 0 -12 6*d
        0 6*d 4*d^2 0 -6*d 2*d^2
        -1 0 0 1 0 0
        0 -12 -6*d 0 12 -6*d
        0 6*d 2*d^2 0 -6*d 4*d^2];

         for i = 1:6
                if i == 1 || i == 4
                Kc(i,:) = Kc(i,:)*E*Acr/Lc;
                Kr(i,:) = Kr(i,:)*E*Ar/Lr;
                Kh1(i,:) = Kh(i,:)*E*Ach(1)/d;
                Kh2(i,:) = Kh(i,:)*E*Ach(2)/d;
                Kh3(i,:) = Kh(i,:)*E*Ach(3)/d;
                Kh4(i,:) = Kh(i,:)*E*Ach(4)/d;
                Kh5(i,:) = Kh(i,:)*E*Ach(5)/d;
                Kh6(i,:) = Kh(i,:)*E*Ach(6)/d;
                Kh7(i,:) = Kh(i,:)*E*Ach(6)/Lr;
         else
                Kc(i,:) = Kc(i,:)*E*Ic/Lc^3;
                Kr(i,:) = Kr(i,:)*E*Ir/Lr^3;
                Kh1(i,:) = Kh(i,:)*E*Ix(1)/d^3;
                Kh2(i,:) = Kh(i,:)*E*Ix(2)/d^3;
                Kh3(i,:) = Kh(i,:)*E*Ix(3)/d^3;
                Kh4(i,:) = Kh(i,:)*E*Ix(4)/d^3;
                Kh5(i,:) = Kh(i,:)*E*Ix(5)/d^3;
                Kh6(i,:) = Kh(i,:)*E*Iy(6)/d^3;
                Kh7(i,:) = Kh(i,:)*E*Iy(6)/Lr^3;
               end
          end

           Kh = cat(6, Kh1, Kh2, Kh3, Kh4, Kh5, Kh6)

          % Transformation matries

          % Left-hand column
          T1 = matrix;
          T1(2,1) = -1;
          T1(1,2) = 1;
          T1(3,3) = 1;
          T1(5,4) = -1;
          T1(4,5) = 1;
          T1(6,6) = 1;

          % Start of haunch left side
          T2 = matrix;
          T2(1,4) = cosd(alpha);
          T2(1,5) = sind(alpha);
          T2(2,4) = -sind(alpha);
          T2(2,5) = cosd(alpha);
          T2(3,6) = 1;
          T2(4,7) = cosd(alpha);
          T2(4,8) = sind(alpha);
          T2(5,7) = -sind(alpha);
          T2(5,8) = cosd(alpha);
          T2(6,9) = 1;

          % Remaing Sections of eaves haunch
          for i = 1:5
                 ThL(:,:,i) = circshift(T2,[0,i*3])
          endfor

          % Rest of Left Rafter
          TrL = matrix;
          TrL(1,22) = cosd(alpha);
          TrL(1,23) = sind(alpha);
          TrL(2,22) = -sind(alpha);
          TrL(2,23) = cosd(alpha);
          TrL(3,24) = 1;
          TrL(4,25) = cosd(alpha);
          TrL(4,26) = sind(alpha);
          TrL(5,25) = -sind(alpha);
          TrL(5,26) = cosd(alpha);
          TRL(6,27) = 1;

          % Start of Right Rafter
          
          TrR = matrix;
          TrR(1,25) = cosd(-alpha);
          TrR(1,26) = sind(-alpha);
          TrR(2,25) = -sind(-alpha);
          TrR(2,26) = cosd(-alpha);
          TrR(3,27) = 1;
          TrR(4,28) = cosd(-alpha);
          TrR(4,29) = sind(-alpha);
          TrR(5,28) = -sind(-alpha);
          TrR(5,29) = cosd(-alpha);
          TrR(6,30) = 1;

          % Start of Eaves Haunch
          
          Th3 = matrix;
          Th3(1,28) = cosd(-alpha);
          Th3(1,29) = sind(-alpha);
          Th3(2,28) = -sind(-alpha);
          Th3(2,29) = cosd(-alpha);
          Th3(3,30) = 1;
          Th3(4,31) = cosd(-alpha);
          Th3(4,32) = sind(-alpha);
          Th3(5,31) = -sind(-alpha);
          Th3(5,32) = cosd(-alpha);
          Th3(6,33) = 1;

          % Remaining sections of eaves haunch
          
          for i = 1:5
                 ThR(:,:,i) = circshift(Th3,[0,i*3])
          endfor

          % Right-hand column
          
          T4 = matrix;
          T4(2,49) = -1;
          T4(1,50) = 1;
          T4(3,51) = 1;
          T4(5,46) = -1;
          T4(4,47) = 1;
          T4(6,48) = 1;
          
          % Transformation from local to global

          Km1 = T1'*Kc*T1;
          Km2 = T2'*Kh(1)*T2;
          Km3 = TrL'*Kh(7)*TrL;
          Km4 = TrR'*Kh(7)*TrR;
          Km5 = Th3'*Kh(6)*Th3;
          Km6 = T4'*Kc*T4;
         
(Relevant Code)  for i = 1:5
                 KmhL(:,:,i) = ThL(:,:,i)'*Kh(i+1)*ThL(:,:,i);
                 KmhR(:,:,i) = ThR(:,:,i)'*Kh(6-i)*ThR(:,:,i);
         endfor   

         % Assembly
         
          Ks = KmhL(:,:,1) + KmhL(:,:,2) + KmhL(:,:,3)+ KmhL(:,:,4)+ KmhL(:,:,5)+ KmhR(:,:,1)+ KmhR(:,:,2)+KmhR(:,:,3)+KmhR(:,:,4)+KmhR(:,:,5) + Km1 + Km2 + Km3 + Km4 + Km5 + Km6

          Ks([1,2,49,50],:)=[];
          Ks(:,[1,2,49,50])= [];

          F = zeros(n,1);
          F(4) = 0.484;
          F(46) = 0.565;
          F([1,2,49,50]) = [];

          % Solving for displacements          
          
          U = inv(Ks)*F
          Ux = [0 0 U(1:46)' 0 0 U(47)]'

For some reason, the KmhL matrix is producing a matrix with all values as zero - like this:

There should be some numbers in this matrix!!

So I really hope you guys don't have to go through the whol;e code, but if you could help, that would be great.

Many Thanks,

Scott

[Walter Roberson]

This is not a matlab question. matlab does not use endfor.

The code appears to be Octave. The purpose of Octave is to force Mathworks out of business (or at the very least to release literally all Mathworks proprietary code as Open Source.)