% One dimensional finite element program
clear;
% set number of elements (numele) and number of nodes (numnod) below
numele=50;
numnod=numele+1;
% X is a set of equispaced nodal coordinates
X= 0:100/numele:100;
% Element node numbers for a sequentially numbered set of nodes;
% third row gives the material number
node= [1:numele
2:numele+1
ones(1,numele)];
prop=[30.E6,1];
% prop( gives material and section properties
ifix=[zeros(1,numele),1];
fext=[zeros(1,numnod)];
% fext is external force
fext=[zeros(1,numnod)];
fext(1)=1000;
% above is for specific case where node 1 is loaded by a constant in time
fint=[zeros(1,numnod)];
%----------- the mass -----------------
mass=[zeros(1,numnod)];
for e=1:numele
densityo=.000722;
mat=node(3,e);
young=prop(mat,1);
areao=prop(mat,2);
% compute element length
lengtho=X(node(2,e))-X(node(1,e));
ma=densityo*lengtho*areao/2;
mass(1,node(1,e))=mass(1,node(1,e))+ma;
mass(1,node(2,e))=mass(1,node(2,e))+ma;
end
%---------------------------------------
% loop over elements
% initial conditions
disp=[zeros(1,numnod)];
vel=[zeros(1,numnod)];
% integrate in time
numstep=500;
outvar=[zeros(1,numstep)];
tutvar=[zeros(1,numstep)];
time=[zeros(1,numstep)];
%*******************************************
for nstep=1:numstep;
fint=[zeros(1,numnod)];
%-------------------------------------------
for e=1:numele
mat=node(3,e);
young=prop(mat,1);
areao=prop(mat,2);
cfl=1.046;
% compute element length
lengtho=X(node(2,e))-X(node(1,e));
c=young*areao/lengtho;
delt=0.5*lengtho/sqrt(young/densityo);
%delt=cfl*length/sqrt(young/densityo);
P=young*(disp(node(2,e))-disp(node(1,e)))/lengtho;
fint(1,node(1,e))=fint(1,node(1,e))-areao*P;
fint(1,node(2,e))=fint(1,node(2,e))+areao*P;
end
%--------------------------------------------
a=fext-fint;
a=a./mass;
% support conditions (boundary conditions)
for n=1:numnod
if (ifix(n) == 1);
a(n)=0;
end
end
vel=vel+delt*a;
disp=disp+delt*vel;
outvar(nstep)=disp(26);
tutvar(nstep)=vel(26);
time(nstep)=nstep*delt;
% resp(nstep,1:numnod)=vel(1:numnod);
end
%****************************************
% plot results
figure
subplot(221)
plot(time,outvar)
title('Displt x=50 vs time');xlabel('time');ylabel('displacement')
subplot(222)
plot(time,tutvar)
title('velocity x=50 vs time');xlabel('time');ylabel('velocity')
subplot(224)
plot(X,vel)
title('velocity vs x');xlabel('x');ylabel('velocity')
