FEM – Bar Problem 1

% Clear previous data and graph
clc; clear all; cla

%Input
E=200e6; % Young's modulus (kN/m)
A=0.04; % Area (m^2)
L1=2.5; % Length 1st Element (m)
L2=2.5; % Length 1st Element (m)
P=10; % Point Load (kN)
q=2; % Uniform load (kN/m)

%Create local stifness n forces
klocal1=[(A*E)/L1 -(A*E)/L1;
-(A*E)/L1 (A*E)/L1];

klocal2=[(A*E)/L2 -(A*E)/L2;
-(A*E)/L2 (A*E)/L2];

flocal1=[q*L1/2;q*L1/2];
flocal2=[q*L2/2;q*L2/2];

%Assemble
Kglobal=zeros(3,3);
Kglobal(1:2,1:2)=Kglobal(1:2,1:2)+klocal1;
Kglobal(2:3,2:3)=Kglobal(2:3,2:3)+klocal2;

Fglobal=zeros(3,1);
Fglobal(1:2,1)=Fglobal(1:2,1)+flocal1;
Fglobal(2:3,1)=Fglobal(2:3,1)+flocal2;
FG=Fglobal;
Fglobal(3,1) =Fglobal(3,1)+P; % Point Load

KG=Kglobal;

%Impose BC
Kglobal(1,:)=[];
Kglobal(:,1)=[];

Fglobal(1,:)=[];

%Solve for global displacement
U=Kglobal\Fglobal;

% local result Element 1
U1=[0;U(1,1)];
f1=(klocal1*U1)-flocal1;
stress1=f1/A;

% local result Element 2
U2=[U(1,1);U(2,1)];
f2=(klocal2*U2)-flocal2;
stress2=f2/A;

% Reacttion
UU=[0;U(1,1);U(2,1)];
R=KG*UU-FG;

% Plot Result
L = L1 + L2;
X = 0:L/20:L;
Uexact = -q/(E*A) * X.^2/2 + (P+q*L)/(E*A)*X; % Exact

FEM2Nodex = [0;L1;L];
FEM2Nodey = [0;U(1,1);U(2,1)];
hold on
plot(FEM2Nodex,FEM2Nodey,'--*b')
plot(X,Uexact,'-k')
hold on

title('Displacements u(x) over the distance, x');
xlabel('x');
ylabel('u (x) ');
legend('Ulinear ','Uexact')