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.
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:
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:
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
