ENS 5147: MATLAB Code Elaboration for 2D Plane Stress Bar Analysis

This document presents a solution to an assignment involving the finite element analysis of a 2D plane stress bar using MATLAB. The solution elaborates on several MATLAB code functions, including Jacobian4N.m, GaussWeights.m, GaussWeights2D.m, and FEAnalysis4N.m. The code snippets and their functionalities are explained, covering aspects such as nodal coordinate calculation, Gauss point conversion, and stress recovery techniques. The assignment involves modifying the provided code to analyze the bar with specific element configurations and extracting results like nodal displacements and stresses. This document serves as a comprehensive guide to understanding the implementation of FEM in MATLAB for structural analysis, contributed by a student and available on Desklib, a platform providing study resources for students.

ENGINEERING MATHEMATICS
MATLAB CODE ELABORATION
From the Matlab files, the files contains some functions which when executed, the matlab code will
execute with preferred results
Jacobian4N.m
function J = Jacobian4N(st,nc)
s = st(1); t = st(2);
dNds = 0.25*[-1+t 1-t 1+t -1-t];
dNdt = 0.25*[-1+s -1-s 1+s 1-s];
J = [dNds*nc(:,1) dNds*nc(:,2);dNdt*nc(:,1) dNdt*nc(:,2)];
From this code snippet, this is a function which is given the name Jacobian4N, the function has
arguments st of parameters nc.
The function has variable s which is assigned the values of st(1) and another variable t which is assigned
value st(2) ,
DNds are the nodes at s, this is a code formula to calculate the node displacement at state s and
DNdt are the nodes at t . this is a code formulae to calculate the node displacement which is at state t
J = [dNds*nc(:,1) dNds*nc(:,2);dNdt*nc(:,1) dNdt*nc(:,2)]; This code executes the final results when the
function of variable j is implemented.
GaussWeights.m
function weight = GaussWeights(n)
switch n
case 1
weight(1) = 2.0;
case 2
weight(1) = 1.0;
weight(2) = 1.0;
case 3
weight(1) = 5.0/9.0;
weight(3) = 5.0/9.0;
weight(2) = 8.0/9.0;
case 4

weight(1) = 0.347854845137454;
weight(2) = 0.652145154862546;
weight(4) = 0.347854845137454;
weight(3) = 0.652145154862546;
…………………………………………………..
weight(1) = 0.129484966168870;
weight(2) = 0.279705391489277;
weight(3) = 0.381830050505119;
weight(7) = 0.129484966168870;
weight(6) = 0.279705391489277;
weight(5) = 0.381830050505119;
weight(4) = 0.417959183673469;
end
This is a class with a conditional statement switch, the syntax is purely a switch statement which in this
case does the following.
When the file is executed, the program will run from top – down checking if the entries meet the cases
stated above;
For case 1, when the input resembles the values of the case, the program will execute at that particular
point and execute the results at case 1. The same will happen to the rest of the switch cases.
GaussWeights2D
function wt = GaussWeights2D(ngp1D)
wt1D = GaussWeights(ngp1D);
wt = [];
for i=1:ngp1D
for j=1:ngp1D
wt = [wt;wt1D(i)*wt1D(j)];
end
end
The class contains the function named GaussWeights2D with variable wt and argument ngp1D. The
initial variable wt1D is allocated a value GaussWeights(ngp1D) and the wt variable is given an empty
array[].
A for loop is introduced in the function , the loop is to replicate the value of wt in accordance with the
values of wt in the program.

FEAnalysis4N.m
Perform FE analysis of a plane stress bar with varying cross section. Use symmetry and only model the
upper half of the bar.
The code has a functional error handler “function error = FEAnalysis4N(nElementsX,nElementsY)”.
The below code is a snippet of the class code.
length = 0.6;
E = 10000; % Young's modulus
v = 0.3; % Poisson's ratio
syms x;
b = 0.05 + 0.1*x/length; % half width of bar
t = 0.02; % thickness of bar
%q = -10/(2*b) + 20*heaviside(x-0.5*length)/(2*b); % distributed force applied as body load
q = 10 - 20*heaviside(x-0.5*length);
endTraction = -1/subs(b,length); % force at right end of bar applied as a traction in x direction
D = (E/(1-v^2))*[1 v 0;v 1 0;0 0 0.5*(1-v)];
nElements = nElementsX*nElementsY;
nNodesX = nElementsX + 1;
nNodesY = nElementsY + 1;
nNodes = nNodesX*nNodesY;
% work out nodal coordinates
nodeCoord = zeros(nNodes,2); % each row represents one node, col 1 is x coord, col 2 is y coord
k = 0;
The symbol percentage (%) represents the start of the comment which gives more explanation on the
operation of the line of code. The comment starts immediately after the % to the end of the line.
The code for this function is for working out the nodal coordinates.
k = 0;
for i=1:nNodesX
xc = (i-1)*length/nElementsX;
ymax = subs(b,xc);
for j=1:nNodesY
k = k + 1;
nodeCoord(k,1) = xc;
nodeCoord(k,2) = ymax*(j-1)/nElementsY;
end

end
The for loop above which is part of the code in the same function is used for auto incrementing the
nodes coordinates.
% apply superconvergent patch recovery to find stresses at nodes
rStress = zeros(nNodes,3);
count = zeros(nNodes,3);
gp = GaussPoints2D(2);
for i=2:nNodesX-1 % patches are centred on internal nodes
for j=2:nNodesY-1 % i.e. not the nodes on the boundary
xgp = []; scStresses = [];
for k1=1:4 % elments around the patch node, anticlockwise starting from top left
elnum = (i-2)*nElementsY+j-1;
if k1==3 || k1==4
elnum = elnum + nElementsY;
end
if k1==1 || k1==4
elnum = elnum + 1;
end
xegp = Element4NInterpolate(gp,nodeCoord(element(elnum,:),:)); % convert gauss points from
local to global coords
for j1=1:4
indexArray([2*j1-1 2*j1])= [2*element(elnum,j1)-1 2*element(elnum,j1)];
end
scStresses = [scStresses;Element4NStress(gp,nodeCoord(element(elnum,:),:),u(indexArray),D)];
xgp = [xgp;xegp];
end
A = [ones(size(xgp,1),1) xgp(:,1) xgp(:,2) xgp(:,1).*xgp(:,2)];
pcoeffs = inv(A'*A)*A'*scStresses; % fit curve to stresses at superconvergent points
indexArray = [(i-1)*nNodesY+j]; % node in the centre of the patch
for k1=1:4 % elments around the patch node, anticlockwise starting from top left
elnum = (i-2)*nElementsY+j-1;
if k1==3 || k1==4
elnum = elnum + nElementsY;
end
if k1==1 || k1==4
elnum = elnum + 1;
end
% handle edges
if i==2 && (k1==1 || k1==2)
indexArray = [indexArray element(elnum,[4 1])];
end
if i==nElementsX && (k1==3 || k1==4)
indexArray = [indexArray element(elnum,[2 3])];
end
if j==2 && (k1==2 || k1==3)

indexArray = [indexArray element(elnum,[1 2])];
end
if j==nElementsY && (k1==1 || k1 ==4)
indexArray = [indexArray element(elnum,[3 4])];
end
end
indexArray = unique(indexArray); % remove repeats
rStress(indexArray,:) = rStress(indexArray,:) + [ones(size(indexArray,2),1) nodeCoord(indexArray,1)
nodeCoord(indexArray,2) nodeCoord(indexArray,1).*nodeCoord(indexArray,2)]*pcoeffs;
count(indexArray,:) = count(indexArray,:) + 1;
end
end
rStress = rStress./count; % actually in this case there is no need to average - only one recovery for each
node
The code snippet above is used to change or convert the gauss points from local to global. More
explanation of the code can be traced from the convert % statement which indicates the comment
of the code line