Tuesday, April 26, 2016

How Would You (mostly) Automate a FEA Problem in MATLAB?

In order to make a MATLAB program to (mostly) automate solving a Finite Element Analysis problem, all you have to do is follow the same steps in the previous post. In summary:

  1. Define Element Type
  2. Define Node Positions/Representations
  3. Define Nodal Connections (i.e. members)
  4. Setup Local Matrices
  5. Determine Positions in Global Matrix
  6. Insert Local Matrices into Global Matrix
  7. Determine Constraints and Trim Global Matrix Accordingly
  8. Apply Forces, and Solve
For this exercise, I solved the example from the last post, which also allowed me to compare the results between the book, the hand calculations, and my Matlab program. The code is commented and clearly shows each step:

Note: If the formatted version is too hard to read, skip to the bottom for the plain text version copied directly from Matlab


FEM Example 2.1

% Sam Ciocco
clear, clc, close all;

Variable Inputs

numElements = 6; %Number of elements
numNodes = 5; %Number of nodes

DOF = 2; %Degrees of Freedom

A = 8; %Area in in^2
E = 1.9E6; %Elastic Modulus, in lb/in^2

forces = [0;0;0;0;0;0;0;-500;0;-500]; %Forces Vector, in order 1x;1y;2x;2y;...

constrained = [1,3]; %List of constrained nodes

Element Table

% For now, going to hardcode the table
% Will add node locations to calculate lengths and angles if I have time

elements = ones(numElements,7);
for i=1:numElements
    elements(i,1) = i;
end

% Declare the first node for each element
elements(:,2) = [1;2;3;2;2;4];
% Declare the second node for each element
elements(:,3) = [2;3;4;4;5;5];
% Declare lengths for each element
elements(:,4) = [36;50.9;36;36;50.9;36];
% Declare angles for each element, in degrees
elements(:,5) = [0;135;0;90;45;0];

%Insert sine and cosine for each element, based on the angle
for j=1:numElements
    elements(j,6) = cosd(elements(j,5));
    elements(j,7) = sind(elements(j,5));
end

%Element table organized as follows:
% Element | Node 1 | Node 2 | Length (in) | Angle (deg) | Cosine | Sine

Set up each local matrix

%Apparently I figured out how to make an array of matrices
for i=1:numElements
    keq(i) = A*E/elements(i,4);
    Klocal(:,:,i) = keq(i).*[elements(i,6)^2,elements(i,6)*elements(i,7),-(elements(i,6)^2),-elements(i,6)*elements(i,7);
                elements(i,6)*elements(i,7),elements(i,7)^2,-elements(i,6)*elements(i,7),-(elements(i,7)^2);
                -(elements(i,6)^2),-elements(i,6)*elements(i,7),elements(i,6)^2,elements(i,6)*elements(i,7);
                -elements(i,6)*elements(i,7),-(elements(i,7)^2),elements(i,6)*elements(i,7),elements(i,7)^2];
end

Set up Global Matrix

% Declare a matrix of correct size for each node to have its number of
% degrees of freedom (DOF)

Kglobal = zeros(numNodes*DOF,numNodes*DOF);

% Now painstakingly fill in each local matrix into the global matrix
% Ensuring each local matrix is placed correctly according to its nodes
for i=1:numElements
    for j=0:1
       for k=0:1
           Kglobal(elements(i,2)*2-1+j,elements(i,2)*2-1+k) = Kglobal(elements(i,2)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+1,k+1,i);
       end
    end

    for j=0:1
       for k=0:1
           Kglobal(elements(i,2)*2-1+j,elements(i,3)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+1,k+3,i);
       end
    end

    for j=0:1
       for k=0:1
           Kglobal(elements(i,3)*2-1+j,elements(i,3)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,3)*2-1+k)+Klocal(j+3,k+3,i);
       end
    end

    for j=0:1
       for k=0:1
           Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+3,k+1,i);
       end
    end

    % Old code, not completely sure why it didn't work

%     Kglobal(elements(i,2)*2-1,elements(i,2)*2-1) = Kglobal(elements(i,2)*2-1,elements(i,2)*2-1)+Klocal(1,1,i);
%     Kglobal(elements(i,2)*2-1,elements(i,2)*2-1+1) = Kglobal(elements(i,2)*2-1,elements(i,2)*2-1+1)+Klocal(1,2,i);
%     Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1) = Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1)+Klocal(2,1,i);
%     Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1+1) = Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1+1)+Klocal(2,2,i);
end

Create list of unconstrained (aka moving) nodes

    % This code assumes each node can either move or not; no partial
    % constraints

nodes = 1:numNodes; %Create list of nodes
movingNum = 1;
moving = zeros(numNodes-length(constrained));
for i=1:numNodes
        if((nodes(i) ~= constrained(1)) && (nodes(i) ~= constrained(2)))
            moving(movingNum) = i;
            movingNum = movingNum+1;
        end
end

Trim Global Matrix in Accordance with Constrained Nodes

% Again assuming any constrained nodes are fully constrained

Kglobal2 = zeros((numNodes-length(constrained))*2,(numNodes-length(constrained))*2);

    for i=1:length(moving)
        for j=1:length(moving)

            for k=0:1
                for m=0:1

                    Kglobal2(i*2-1+k,j*2-1+m) = Kglobal(moving(i)*2-1+k,moving(j)*2-1+m);


                end
            end
        end
    end

% for l=1:length(moving)
% for i=:1
%     for j=0:1
%     Kglobal2(l*2-1+i,l*2-1+j) = Kglobal(moving(l)*2-1+i,moving(l)*2-1+j);
%
%     end
% end
% end

Trim Forces in Accordance with Constrained Nodes

forces2 = zeros(length(moving)*2,1);

for l=1:length(moving)
    forces2(l*2-1) = forces(moving(l)*2-1);
    forces2(l*2) = forces(moving(l)*2);
end

Solve Matrix for movements

movements = Kglobal2\forces2;
%movements is the vector of the first node x, the first node y, the second
%node x, the second node y,...


numMoving = 1;
for i=1:2:length(movements)

            fprintf('\nThe movement of node %.0f in the x-direction is %.3f inches\n',moving(numMoving),movements(i));
            fprintf('The movement of node %.0f in the y-direction is %.3f inches\n',moving(numMoving),movements(i+1));
            numMoving = numMoving+1;

end

The movement of node 2 in the x-direction is -0.004 inches
The movement of node 2 in the y-direction is -0.010 inches

The movement of node 4 in the x-direction is 0.001 inches
The movement of node 4 in the y-direction is -0.011 inches

The movement of node 5 in the x-direction is 0.002 inches
The movement of node 5 in the y-direction is -0.020 inches
Published with MATLAB® R2015a

Plain Text Version:


%% FEM Example 2.1

% Sam Ciocco
clear, clc, close all;

%% Variable Inputs

numElements = 6; %Number of elements
numNodes = 5; %Number of nodes

DOF = 2; %Degrees of Freedom

A = 8; %Area in in^2
E = 1.9E6; %Elastic Modulus, in lb/in^2

forces = [0;0;0;0;0;0;0;-500;0;-500]; %Forces Vector, in order 1x;1y;2x;2y;...

constrained = [1,3]; %List of constrained nodes


%% Element Table

% For now, going to hardcode the table
% Will add node locations to calculate lengths and angles if I have time

elements = ones(numElements,7);
for i=1:numElements
    elements(i,1) = i;
end

% Declare the first node for each element
elements(:,2) = [1;2;3;2;2;4];
% Declare the second node for each element
elements(:,3) = [2;3;4;4;5;5];
% Declare lengths for each element
elements(:,4) = [36;50.9;36;36;50.9;36];
% Declare angles for each element, in degrees
elements(:,5) = [0;135;0;90;45;0];

%Insert sine and cosine for each element, based on the angle
for j=1:numElements
    elements(j,6) = cosd(elements(j,5));
    elements(j,7) = sind(elements(j,5));
end

%Element table organized as follows:
% Element | Node 1 | Node 2 | Length (in) | Angle (deg) | Cosine | Sine

%% Set up each local matrix

%Apparently I figured out how to make an array of matrices
for i=1:numElements
    keq(i) = A*E/elements(i,4);
    Klocal(:,:,i) = keq(i).*[elements(i,6)^2,elements(i,6)*elements(i,7),-(elements(i,6)^2),-elements(i,6)*elements(i,7);
                elements(i,6)*elements(i,7),elements(i,7)^2,-elements(i,6)*elements(i,7),-(elements(i,7)^2);
                -(elements(i,6)^2),-elements(i,6)*elements(i,7),elements(i,6)^2,elements(i,6)*elements(i,7);
                -elements(i,6)*elements(i,7),-(elements(i,7)^2),elements(i,6)*elements(i,7),elements(i,7)^2];
end

%% Set up Global Matrix

% Declare a matrix of correct size for each node to have its number of
% degrees of freedom (DOF)

Kglobal = zeros(numNodes*DOF,numNodes*DOF);

% Now painstakingly fill in each local matrix into the global matrix
% Ensuring each local matrix is placed correctly according to its nodes
for i=1:numElements
    for j=0:1
       for k=0:1
           Kglobal(elements(i,2)*2-1+j,elements(i,2)*2-1+k) = Kglobal(elements(i,2)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+1,k+1,i);
       end
    end
    
    for j=0:1
       for k=0:1
           Kglobal(elements(i,2)*2-1+j,elements(i,3)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+1,k+3,i);
       end
    end
    
    for j=0:1
       for k=0:1
           Kglobal(elements(i,3)*2-1+j,elements(i,3)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,3)*2-1+k)+Klocal(j+3,k+3,i);
       end
    end
    
    for j=0:1
       for k=0:1
           Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+3,k+1,i);
       end
    end
    
    % Old code, not completely sure why it didn't work
    
%     Kglobal(elements(i,2)*2-1,elements(i,2)*2-1) = Kglobal(elements(i,2)*2-1,elements(i,2)*2-1)+Klocal(1,1,i);
%     Kglobal(elements(i,2)*2-1,elements(i,2)*2-1+1) = Kglobal(elements(i,2)*2-1,elements(i,2)*2-1+1)+Klocal(1,2,i);
%     Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1) = Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1)+Klocal(2,1,i);
%     Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1+1) = Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1+1)+Klocal(2,2,i);
end
    
%% Create list of unconstrained (aka moving) nodes

    % This code assumes each node can either move or not; no partial
    % constraints
    
nodes = 1:numNodes; %Create list of nodes
movingNum = 1;
moving = zeros(numNodes-length(constrained));
for i=1:numNodes
        if((nodes(i) ~= constrained(1)) && (nodes(i) ~= constrained(2)))
            moving(movingNum) = i;
            movingNum = movingNum+1;
        end
end

%% Trim Global Matrix in Accordance with Constrained Nodes

% Again assuming any constrained nodes are fully constrained

Kglobal2 = zeros((numNodes-length(constrained))*2,(numNodes-length(constrained))*2);

    for i=1:length(moving)
        for j=1:length(moving)
        
            for k=0:1
                for m=0:1
        
                    Kglobal2(i*2-1+k,j*2-1+m) = Kglobal(moving(i)*2-1+k,moving(j)*2-1+m);
        
        
                end
            end
        end
    end

% for l=1:length(moving)
% for i=:1
%     for j=0:1
%     Kglobal2(l*2-1+i,l*2-1+j) = Kglobal(moving(l)*2-1+i,moving(l)*2-1+j);
%     
%     end
% end
% end

%% Trim Forces in Accordance with Constrained Nodes

forces2 = zeros(length(moving)*2,1);

for l=1:length(moving)
    forces2(l*2-1) = forces(moving(l)*2-1);
    forces2(l*2) = forces(moving(l)*2);
end


%% Solve Matrix for movements

movements = Kglobal2\forces2;
%movements is the vector of the first node x, the first node y, the second
%node x, the second node y,...


numMoving = 1;
for i=1:2:length(movements)
    
            fprintf('\nThe movement of node %.0f in the x-direction is %.3f inches\n',moving(numMoving),movements(i));
            fprintf('The movement of node %.0f in the y-direction is %.3f inches\n',moving(numMoving),movements(i+1));
            numMoving = numMoving+1;
            
end




Posted by at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)