Finite Element Analysis for an Axial Bar Coding Help

Problem 2: Write the generalized FEM program in Matlab to solve axial bar problems, based on solving P3-12, the output of the program should have displacements (all nodes), stresses (all elements), unknown reactions, internal forces (all elements), and type of internal forces (all elements). Plot the displacement of the tip node (where the force is applied) vs. number of elements and compare with the analytical solution in the plot. (15 points)

I have this code set up from a similar problem but I am having trouble converting it for this problem. All help would be great!

clear all

clc

close all

nodeConn = [3 1; 1 2];

elemProp = [1000; 2000];

disp = [NaN;20e-3; 0.0];

force = [0.0; NaN; NaN];

[disp, force, elemForce, elemForceType] = AEM461_HW1_Spring(nodeConn, elemProp, disp, force)

function [disp, force, elemForce, elemForceType] = AEM461_HW1_Spring(nodeConn, elemProp, disp, force)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%my first FEM program for springs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

node = unique(nodeConn); % total number of nodes

noElem = size(nodeConn,1); % total number of elements

GlobalK = zeros(length(node),length(node));

for elemCount = 1:noElem

node1Pos = find(node == nodeConn(elemCount,1)); % finding the position of node 1 in global array

node2Pos = find(node == nodeConn(elemCount,2)); % finding the position of node 2 in global array

GlobalK(node1Pos, node1Pos) = GlobalK(node1Pos, node1Pos) + elemProp(elemCount);

GlobalK(node2Pos, node2Pos) = GlobalK(node2Pos, node2Pos) + elemProp(elemCount);

GlobalK(node2Pos, node1Pos) = GlobalK(node2Pos, node1Pos) - elemProp(elemCount);

GlobalK(node1Pos, node2Pos) = GlobalK(node1Pos, node2Pos) - elemProp(elemCount);

end

dispUind = find(isnan(disp)==1);%%%% finding the indices of unknown displacements

dispKind = find(isnan(disp)==0);%%%% finding the indices of unknown forces

Kaa = GlobalK(dispUind,dispUind);

Kbb = GlobalK(dispKind,dispKind);

Kab = GlobalK(dispUind,dispKind);

Kba = GlobalK(dispKind,dispUind);

forceK = force(dispUind); %%%% known force vector

dispK = disp(dispKind); %%%% known displacement vector

modForce = forceK - Kab*dispK; %%%%modfied force vector

dispU = Kaa\modForce; %%%% unknown displacement vector

forceU = Kba*dispU + Kbb*dispK; %%%% unknown force vector

disp(dispUind) = dispU; %%%% forming total displacement vector

force(dispKind) = forceU; %%%% forming total force vector

elemForce = zeros(noElem,1); %%% initializing elemental forces

elemForceType = cell(noElem,1); %%%% initializing types of elemental forces

for elemCount = 1:noElem

node1Pos = find(node == nodeConn(elemCount,1)); % finding the position of node 1 in global array

node2Pos = find(node == nodeConn(elemCount,2)); % finding the position of node 2 in global array

dispElemVec = [disp(node1Pos); disp(node2Pos)]; %%% forming elemental displacemental vector

elemForceVec = elemProp(elemCount)*[1 -1; -1 1]*dispElemVec; %%% finding elemental force vector

elemForce(elemCount) = abs(elemForceVec(1)); %%% magnitude of the elemental force

if sign(elemForceVec(1)) == -1

elemForceType{elemCount} = 'T'; %%% elemental force is tensile

else

elemForceType{elemCount} = 'C'; %%% elemental force is compressive

end

end

[Maaz Nayeem]

did you get this code? i am facing difficulty for the same problem