Free download MATLAB file for Finite Element beam analysis using beam elements + solved examples

Visit http://CivilEngineeringBible.com for newer posts and files

Many structural systems used in practice consist of long slender members that are subjected to loading normal to their longitudinal axis and must resist bending and shear forces. They are called beams. Beams are governed by a fourth-order  differential equation.

beam element
Assumed degrees of freedom for a beam element used in this code.

Example1:
Consider a simple cantilever beam with a circular cross-section of 10 in diameter and a length of 400in. E = 30e6 psi. Load = 1000 lbs in downward direction at the right end of the beam.

Analytical solution for the mximum deflection and slope at the right end of the beam is as follows:
Slope = P(L^2)/2*EI =  0.0054
Deflection = P(L^3)/3*EI = 1.44866

and Finite Element result generated by MATLAB program is shown below:
cantilever

Example2:
Consider a simply supported beam with a circular cross-section of 10 in diameter and a length of 400 in. The Young’s Modulus of the beam is 30e6. Load = 1000 lbs in downward direction at the right end of the beam.

Analytical solution -> Maximum deflection at the center of the beam is:
deflection = Ymax = P(L^3)/(48*EI) = 0.090541
slope = 0.0

and Finite Element results generated by MATLAB program for two element model and 4 element model are shown below:

simply supported

simply supported4element

Feel free to contact me via email (hosseinali.sut@gmail.com) for more details.

MATLAB Code is:

main.m *****************************

%%
% Developed by Massoud Hosseinali (hosseinali.sut@gmail.com) – 2015
% Finite element MATLAB program for beam analysis
%%
clear all;clc;close all;
Nodes = load(‘nodes.txt’); % coordinate of node
Elements = load(‘elements.txt’);    % first node, second node, E, I
numnode = length(Nodes); %Holds number of nodes
numelem = size(Elements,1); %Holds number of elements

K = zeros(2*numnode,2*numnode);
F = zeros(2*numnode, 1);
U = zeros(2*numnode, 1);
act = [ 2 3 4 5 6 7 8 10 ]; %Hold active DOFs
for ie = 1:numelem
DOFs = [2*Elements(ie, 1)-1, 2*Elements(ie, 1), 2*Elements(ie, 2)-1, 2*Elements(ie, 2)]; %Holds element’s DOFs
L = abs(Nodes(Elements(ie,2)) – Nodes(Elements(ie,1)));  %Holds length of element
E = Elements(ie,3); %Holds modolus of elasticiy of element
I = Elements(ie,4); %Holds moment of inertia of element
K(DOFs,DOFs) = K(DOFs,DOFs) + (E*I/(L^3))*[12 6*L -12 6*L; 6*L 4*L^2 -6*L 2*L^2; -12 -6*L 12 -6*L; 6*L 2*L^2 -6*L 4*L^2]; % Calculates the element stiffness matrix and assembles it to the global stiffness matrix
end

F(5) = -1000;
U(act) = K(act,act)\F(act); % Solves linear solutions

icounter= 0;
for ie = 1:numelem
DOFs = [2*Elements(ie, 1)-1, 2*Elements(ie, 1), 2*Elements(ie, 2)-1, 2*Elements(ie, 2)]; %Holds element’s DOFs
L = abs(Nodes(Elements(ie,2)) – Nodes(Elements(ie,1)));  %Holds length of element

for s = 0:(0.1*L):L
icounter = icounter  + 1;
slratio = s/L;
%Shape functions
N1 = (1-3*slratio^2+2*slratio^3);
N2 = L*(slratio-2*slratio^2+slratio^3);
N3 = (3*slratio^2-2*slratio^3);
N4 = L*(-slratio^2+slratio^3);
N = [N1 N2 N3 N4];
x(icounter) = s + Nodes(Elements(ie,1));
v(icounter) = N*U(DOFs);
end
end
hold all
plot(x,v);
plot(Nodes, U(1:2:end,1), ‘o’, ‘MarkerSize’, 10, ‘Color’, ‘Black’)
xlabel(‘X’)
ylabel(‘Deflection’)
grid on

Elements.txt (Sample) *******************

1 2 30e6 490.8739
2 3 30e6 490.8739
3 4 30e6 490.8739
4 5 30e6 490.8739

Nodes.txt (Sample) *******************

0
100
200
300
400

Advertisements