Overview:

The analysis above was done using a homemade FEA solver written in MATLAB.   The user specifies a structural frame, clamps to hold the frame, and the loads to be applied.  Running the code then calculates the resulting strain in each beam and plots the deformed structure. 

(Figure 1)

Theory:

The analysis is based on simple beam theory.   Each element is subject to the following three equations which govern axial loading, bending loads, and bending moments.

\[F_{x}=EAU_{x} \\ F_{y}=-EI\frac{\partial U_{y}}{\partial x} \\ M = EI\frac{\partial^{2} U_{y}}{\partial x^2} \]

Analysis:

The three governing equations from earlier were integrated to find displacements and rotations.  Solving the system for a single element subject to unique loads and moments on each end yields the following equation.  

\[\begin{bmatrix}F_x^1 \\F_y^1 \\ M^1 \\ F_x^2 \\F_y^2 \\ M^2\end{bmatrix} = \begin{bmatrix} \frac{EA}{L} & 0 & 0& -\frac{EA}{L} & 0 & 0 \\0 &\frac{12EI}{L^{3}} & \frac{6EI}{L^2} & 0 &-\frac{12EI}{L^{3}} & \frac{6EI}{L^2} \\ 0 &\frac{6EI}{L^2} & \frac{4EI}{L} & 0 &-\frac{6EI}{L^2} & \frac{2EI}{L} \\-\frac{EA}{L} & 0 & 0& \frac{EA}{L} & 0 & 0 \\0 &-\frac{12EI}{L^{3}} & -\frac{6EI}{L^2} & 0 &\frac{12EI}{L^{3}} & -\frac{6EI}{L^2} \\ 0 &\frac{6EI}{L^2} & \frac{2EI}{L} & 0 &-\frac{6EI}{L^2} & \frac{4EI}{L} \end{bmatrix} \begin{bmatrix}U_x^1 \\U_y^1 \\ \theta^1 \\ U_x^2 \\U_y^2 \\ \theta^2\\ \end{bmatrix}\]

The six by six matrix here is called the local stiffness matrix.  It relates the forces experienced by two nodes to the displacement of those two nodes within the local coordinate system of its respective element.  

In order to preform a multi-element analysis, one of these local stiffness matrices is assigned to each element.  Applying a direction cosine matrix, rotates the local stiffness matrix from the element frame into the global frame.  All local stiffness matrices are them assembled into a global stiffness matrix by mapping each component in the matrix to its specific address in the global stiffness matrix.  For example, if an element spans between node one and node three, each of the four quadrants of the local stiffness matrix are assigned positions (1,1), (1,3), (3, 1), and (3,3) respectively in the global stiffness matrix. This is illustrated in Figure 2.

\[\begin{bmatrix} \\ Forces\\ \\ \end{bmatrix} = \begin{bmatrix} & & \\ & Globa\:Stiffness \\ & & \end{bmatrix} \begin{bmatrix} \\Displacements \\ \ \end{bmatrix}\]

Once the global stiffness matrix is assembled, the rows and columns corresponding to the clamped nodes are deleted.  The force vector is populated with the applied loads and displacements are solved for.  The code ends by backwards calculating strains and stresses from the displacements, plotting the final result.

(Figure 2)

Code:

This is the MATLAB script that preformed the analysis shown in the cover photo.  Feel free to use it and play around with it. 

close all
clear all
%Input Parameters
ult_ten = 300*10^6; %Ultimate Tensile Strength (Pa)
ult_com = 300*10^6; %Ultimate Compression Strength (Pa)
 
%Coordinates of nodes (m)
i = 1;
r = 10;
for theta = 0:3*pi()/15:3*pi()/2
   pointt(i,:)=[r*cos(theta),r*sin(theta)];
   i = i+1;
end

r = 12;
for theta = 0:3*pi()/15:3*pi()/2
   pointt(i,:)=[r*cos(theta),r*sin(theta)];
   i = i+1;
end

pointt = [pointt; 12, -6; 10, -6; 0, -10; 0, -12];
co = pointt;
 
%Clamped Nodes
clamps = [19; 20];
 
%Element Matrix:  
%[Node Connection 1, Node Connection 2,X-section Area, Modulus, Area Moment of Inertia
 
A = 0.009;
EE = 190^(9);
IYY = 4.918*10^(-5);
 
    
e = [17 18 A EE IYY
    9 18 A EE IYY
    1 17 A EE IYY
    1 9 A EE IYY
    1 10 A EE IYY
    2 9 A EE IYY
    2 10 A EE IYY
    2 11 A EE IYY
    3 10 A EE IYY
    3 11 A EE IYY
    9 17 A EE IYY
    1 18 A EE IYY
    1 2 A EE IYY
    9 10 A EE IYY
    10 11 A EE IYY
    2 3 A EE IYY
    3 4 A EE IYY
    11 12 A EE IYY
    3 12 A EE IYY
    4 11 A EE IYY
    4 12 A EE IYY
    4 13 A EE IYY
    5 12 A EE IYY
    4 5 A EE IYY
    12 13 A EE IYY
    5 13 A EE IYY
    5 14 A EE IYY
    6 13 A EE IYY
    5 6 A EE IYY
    13 14 A EE IYY
    6 14 A EE IYY
    6 7 A EE IYY
    6 15 A EE IYY
    7 14 A EE IYY
    14 15 A EE IYY
    7 15 A EE IYY
    7 16 A EE IYY
    8 15 A EE IYY
    15 16 A EE IYY
    7 8 A EE IYY
    8 16 A EE IYY
    8 20 A EE IYY
    16 19 A EE IYY
    8 19 A EE IYY
    16 20 A EE IYY]; 
%Forces [Node, F_X (N), F_Y (N)]
External_Forces = [17 -30000 -30000;18 -20000 -20000];
 
%Displacement Display Exaggeration 
Display_Scale = 0.2;
 
%Output
figure(1)
plot(co(:,1),co(:,2),'.')
title('User Defined Structure')
xlabel( 'x coordinate (meters)')
ylabel( 'y coordinate (meters)')
hold on

%Plot Forces
for pp = 1:length(External_Forces(:,1))
Force_dir = atan2(External_Forces(pp,3),External_Forces(pp,2));

x1 = [co(External_Forces(pp,1),1),6*cos(Force_dir)+co(External_Forces(pp,1),1)];
y1 = [co(External_Forces(pp,1),2),6*sin(Force_dir)+co(External_Forces(pp,1),2)];

quiver( x1(1),y1(1),x1(2)-x1(1),y1(2)-y1(1),'LineWidth',1.1,'Color',[1,0,0])    
end

%Plot Clamps
plot(co(clamps(:),1),co(clamps(:),2),'s','color',[1,0,0])

%Plot Input Structure
for ww = 1:size(e(:,1))
plot([co(e(ww,1),1); co(e(ww,2),1)],[co(e(ww,1),2);co(e(ww,2),2)],'color',[0 0 0],'LineWidth',2)
end
 
%Calculate Element Lengths
for w = 1:size(e(:,1))
    length = sqrt((co(e(w,1),1)-co(e(w,2),1))^2+(co(e(w,1),2)-co(e(w,2),2))^2 );
    e(w,6) = length;
end
 
%Develop Local Stiffness Matrices and Assemble into Global Stiffness
%Matrix
 
Num_el = size(e,1);%number of elements
Num_Nodes = size(co,1); %number of nodes
Num_Nodes_el = 2; %number of nodes per element
Deg_of_Freedom = 3; %degree of freedom per node
 
%Initialize Matrix 
K = zeros(Num_Nodes*Deg_of_Freedom,Num_Nodes*Deg_of_Freedom);
F = zeros(Num_Nodes*Deg_of_Freedom,1);
 
for A = 1:Num_el
 %Create Direction Unit Vector For Each Element
 n = (co(e(A,2),:) - co(e(A,1),:));
 n = n./norm(n);%unit vector
 
 %Calculate DCM for matrix rotation
 theta = atan2(n(2),n(1));
 DCM = [cos(theta) -sin(theta) 0; sin(theta) cos(theta) 0; 0 0 1];
 L = e(A,6);
 Area = e(A,3);
 E = e(A,4);
 Iyy = e(A,5);
 
 %Calculate Stiffness matrices and rotate into proper position
 k11 = DCM*[E*Area/L 0 0; 0 12*E*Iyy/(L^3) 6*E*Iyy/(L^2); 0 6*E*Iyy/(L^2) 4*E*Iyy/L]*(DCM');
 k22 = DCM*[E*Area/L 0 0; 0 12*E*Iyy/(L^3) -6*E*Iyy/(L^2); 0 -6*E*Iyy/(L^2) 4*E*Iyy/L]*(DCM');
 k12 = DCM*[-E*Area/L 0 0; 0 -12*E*Iyy/(L^3) 6*E*Iyy/(L^2); 0 -6*E*Iyy/(L^2) 2*E*Iyy/L]*(DCM');
 
 %local stiffness matrix and force vector
 localstiffness = [k11 k12; k12' k22];
 localforce = zeros(Num_Nodes_el*Deg_of_Freedom,1);
 
 %Assemble Local Stiffness Matrices into Global Stiffness Matrix
 %Global Matrix Assembler Barrowed From AERO 315 Code
 for B = 1: Num_Nodes_el  
     for i = 1: Deg_of_Freedom  
         nK1 = (e(A, B)-1)*Deg_of_Freedom+i;
         nKe1 = (B-1)*Deg_of_Freedom+i;
         F(nK1) = F(nK1) + localforce(nKe1);
         for C = 1: Num_Nodes_el
             for j = 1: Deg_of_Freedom
                nK2 = (e(A, C)-1)*Deg_of_Freedom+j;
                nKe2 = (C-1)*Deg_of_Freedom+j;
                K(nK1, nK2) = K(nK1, nK2) + localstiffness(nKe1,nKe2);
             end
         end
     end
 end
end
 
%Configure Forces
for ww = 1:size(External_Forces(:,1))
    F((Deg_of_Freedom*External_Forces(ww,1)-Deg_of_Freedom+1):(Deg_of_Freedom*External_Forces(ww,1)-Deg_of_Freedom+2)) = External_Forces(ww,2:3);
end
 
%Configure Clamps
%Delete Stiffness Matrix Rows and Columns, and Force Matrix Nodes of Clamped Nodes:  Displacement = 0
del_row = [];
for ww = 1:size(clamps(:))
        del_row = [del_row, (Deg_of_Freedom*clamps(ww)-Deg_of_Freedom+1), (Deg_of_Freedom*clamps(ww)-Deg_of_Freedom+2), (Deg_of_Freedom*clamps(ww)-Deg_of_Freedom+3)];
end
 
%Preform Deletion
K(del_row,:) = [];
K(:,del_row) = [];
F(del_row) = [];
 
%Solve for Vector of Displacements and Rotations
u = K\F;
 
%Exaggerate Displacements and Apply to Structure
% Exaggeration Scale
scale = (max(max(co(:,1)),max(co(:,2)))/max(u))*Display_Scale;
 
%Build co_prime, an Exaggerated Node Coordinate Matrix
%Build co_prime_true, a Non- Exaggeration Node Coordinate Matrix
n = 1;
m = 1;
for ww =1:size(co(:,1))
    %If Not Clamped, Add Exaggeration Displacement
    if clamps(n)~= ww
      co_prime(ww,1) = co(ww,1)+u(m)*scale;
      co_prime_true(ww,1) = co(ww,1)+u(m);
      m = m+1;
      co_prime(ww,2) = co(ww,2)+u(m)*scale;
      co_prime_true(ww,2) = co(ww,2)+u(m);
      m = m+2;
    
    %If Clamped, Use Old Coordinate
    else
        n = n+1;
        co_prime(ww,:) = co(ww,:);
        co_prime_true(ww,:) = co(ww,:);
    end
end
 
%Calculate Stresses in Each Element and Develop Color Coded Stress Display
 color_maps = jet(201);
 stress_val = 0;
for w = 1:size(e(:,1))
    lengths(w) = sqrt((co_prime_true(e(w,1),1)-co_prime_true(e(w,2),1))^2+(co_prime_true(e(w,1),2)-co_prime_true(e(w,2),2))^2 );
    %Axial stress do to tension and compression only, need to add in moment
    %stress.  
    sigma(w) = (lengths(w)-e(w,6))*e(w,4);
    
    if sigma(w)>=0
        if (sigma(w)/ult_ten)>1
           sigma(w) = ult_ten;
        end
          stress_val(w) = sigma(w)/ult_ten*100+100;
           
    else
        if abs(sigma(w)/ult_com) >=1
            sigma(w) = ult_com;
        end   
        stress_val(w) = -sigma(w)/ult_com*100+100;
    end
end
 
%Plot Deformed Bars with Color Coded Stress
figure(2)
plot(co_prime(:,1),co_prime(:,2),'.')
title('Exaggerated Deformation')
xlabel( 'x coordinate (meters)')
ylabel( 'y coordinate (meters)')
colormap(jet(201))
c = colorbar('YTickLabel',{-ult_com/(10^6),'','','','','','','','','',ult_ten/(10^6)});
c.Label.String = 'Stress (MPa)';
hold on
a = [3,5]';
for ww = 1:size(e(:,1))
    plot([co_prime(e(ww,1),1); co_prime(e(ww,2),1)],[co_prime(e(ww,1),2);co_prime(e(ww,2),2)],'color', color_maps(round(stress_val(ww)+1),:),'LineWidth',2)
end