FINAL INDEPENDENT PROJECT

Aim:

   To Calculate the Nodal displacement of a beam using Octave.

Description:

•     The Nodal displacements are usually calculated using FEM principles. For doing this the model is first discredited into smaller parts called elements.
•    The Required values are calculated at some finite points called nodes which generally bound the elements.
•    These values are later interpolated with some shape function(approximate function assumed to represent the distribution of the property).

Procedure:

• Beam properties such as Area, length of each element, Young's modulus is collected from the user.
• The stiffness matrix of each element is calculated.

• Then the stiffness matrices are joined by considering common nodes.
• This process is called stitching.
• After sticking the final matrix obtained is called a global stiffness matrix.
• Then the load data is obtained from the user.
• From this dataset the spatial and angular displacement of each node is obtained.
• These displacements define the final state of the beam after the application of load.

Plot:

 INPUTS

MODEL BEAM

(here * represents Fixed support and circle represents roller support)

Code:

Main code:

%setting up the workspace
clc;
close all;
clear all;

%Reading and setting up the initial conditions
elems=input('no of elements :');
nodes=elems+1;
A(elems)=0;
In(elems)=0;
E(elems)=0;
L(elems)=0;
F(2*nodes)=0;
for i=1:elems
fprintf('\nenter the details %d element\n',i);
A(i)=input('enter the area of element(mm^2) : ');
In(i)=input('enter the moment of inertia about bending axis(mm^4) : ');
E(i)=input('enter the youngs modulus of element(MPa) : ');
L(i)=input('enter the length of element(mm) : ');
const(i)=E(i)*In(i)/L(i)^3;
end
%finding the stiffness matrices of individual elements
km=zeros(4,4,elems);
k_tmp=zeros(4,4,elems);
for i=1:elems
k_tmp(1,1,i)=12*const(i);
k_tmp(1,2,i)=const(i)*6*L(i);
k_tmp(1,3,i)=-12*const(i);
k_tmp(1,4,i)=const(i)*6*L(i);
k_tmp(2,1,i)=const(i)*6*L(i);
k_tmp(2,2,i)=4*const(i)*L(i)^2;
k_tmp(2,3,i)=const(i)*-6*L(i);
k_tmp(2,4,i)=const(i)*2*L(i)^2;
k_tmp(3,1,i)=-12*const(i);
k_tmp(3,2,i)=-const(i)*6*L(i);
k_tmp(3,3,i)=12*const(i);
k_tmp(3,4,i)=-const(i)*6*L(i);
k_tmp(4,1,i)=const(i)*6*L(i);
k_tmp(4,2,i)=2*const(i)*L(i)^2;
k_tmp(4,3,i)=const(i)*-6*L(i);
k_tmp(4,4,i)=const(i)*4*L(i)^2;
end
for i=1:elems
km(1,1,i)=km(1,1,i)+k_tmp(1,1,i);
km(1,2,i)=km(1,2,i)+k_tmp(1,2,i);
km(1,3,i)=km(1,3,i)+k_tmp(1,3,i);
km(1,4,i)=km(1,4,i)+k_tmp(1,4,i);
km(2,1,i)=km(2,1,i)+k_tmp(2,1,i);
km(2,2,i)=km(2,2,i)+k_tmp(2,2,i);
km(2,3,i)=km(2,3,i)+k_tmp(2,3,i);
km(2,4,i)=km(2,4,i)+k_tmp(2,4,i);
km(3,1,i)=km(3,1,i)+k_tmp(3,1,i);
km(3,2,i)=km(3,2,i)+k_tmp(3,2,i);
km(3,3,i)=km(3,3,i)+k_tmp(3,3,i);
km(3,4,i)=km(3,4,i)+k_tmp(3,4,i);
km(4,1,i)=km(4,1,i)+k_tmp(4,1,i);
km(4,2,i)=km(4,2,i)+k_tmp(4,2,i);
km(4,3,i)=km(4,3,i)+k_tmp(4,3,i);
km(4,4,i)=km(4,4,i)+k_tmp(4,4,i);
end
%stiching the individual matrices to find global stiffness matrix
K=zeros(2*nodes,2*nodes);
for i=1:elems
od_cnt=2*i-1;
K(od_cnt:od_cnt+3,od_cnt:od_cnt+3)=K(od_cnt:od_cnt+3,od_cnt:od_cnt+3)+km(1:end,1:end,i);
end
k1=K(:,:);
a=[];
for j=1:nodes
i=j;
fprintf('enter the type of support at node %d :',i);
Support_type(i)=input('(2 if fixed supprot ,1- if roller support,0-no support)');
switch(Support_type(i))
case 2
a=[a,2*i,2*i-1];
case 1
a=[a,2*i-1];
case 0
otherwise
fprintf('Invalid!!');
endswitch
end
F_tmp=0;
for i=1:nodes
fprintf('Is there a point load at node %d(press 1 if there is a load and 0 othrwise)',i);
load_type=input(':');
switch(load_type)
case 1
fprintf('Enter the load value(N)');
loc=2*i-1;
F_tmp=input(':');
F(loc)=F(loc)-F_tmp;
case 0
F(2*i-1)=0;
otherwise
fprintf('Invalid!');
endswitch
end
for i=1:nodes
fprintf('Is there a moment at node %d(press 1 if there is a moment and 0 othrwise)',i);
moment_type=input('(clock wise "+" ancti clockwise "-"):');
switch(moment_type)
case 1
fprintf('Enter the moment value(N-mm)');
mom_tmp=input(':');
F(2*i)=F(2*i)+mom_tmp;
case 0
F(2*i)=0;
otherwise
fprintf('Invalid!');
endswitch
end
for i=1:elems
loc=2*i-1;
fprintf('If there is a UDL on element %d press 1 , if its a UVL press 2 or press 0 otherwise',i);
dist_load(i)=input(':');
switch(dist_load(i))
case 1
p=input('enter the load per unit length(N/mm):');
F(loc)=F(loc)-(p*L(i)/2);
F(loc+1)=F(loc+1)-(p*(L(i))^2/12);
F(loc+2)=F(loc+2)-(p*L(i)/2);
F(loc+3)=F(loc+3)+(p*(L(i))^2/12);
case 2
p=input('enter the load per unit length(N/mm) :');
strt=input('enter starting node of of UVL');
if(strt==i)
F(loc)=F(loc)-(p*L(i)/6);
F(loc+1)=F(loc+1)-(p*(L(i))^2/3);
F(loc+2)=F(loc+2)-(p*L(i)/3);
F(loc+3)=F(loc+3)+(p*(L(i))^2/6);
else
F(loc+2)=F(loc+2)-(p*L(i)/6);
F(loc+3)=F(loc+3)+(p*(L(i))^2/3);
F(loc)=F(loc)-(p*L(i)/3);
F(loc+1)=F(loc+1)-(p*(L(i))^2/6);
endif
endswitch
end
n=sort(a);
f1=F(:);
for i=length(a):-1:1
k1(n(i),:)=[];
k1(:,n(i))=[];
f1(n(i))=[];
end
if(det(k1)==0)
x='error:: System improperly constrained!:stiffiness matrix is singular';
disp(x);

else
Q=k1^-1*f1;
end
k=1;

for j=1:2*nodes
b=1;
for i=1:length(n)
if(j==n(i))
b=0;
break;
endif
end
if(b)
q(j)=Q(k);
k=k+1;
else
q(j)=0;
endif
end
q=q';
Overall_dispalcement=q
visualize(A,L,q,nodes,dist_load,Support_type);

Function:

function visualize(A,L,q,nodes,dist_load,Support_type)
elems=nodes-1;
hold on;
x1=0;
for i=1:elems
x2=x1+L(i)/100;
plot([x1 x2],[0 0],'linewidth',A(i)/10);
x1=x2;
end
x1=0;
for i=1:elems
if(Support_type(i)==2)
plot(x1,0,'marker','*','markersize',A(i)/9.9);
elseif(Support_type(i)==1)
plot(x1,0,'marker','o','markersize',A(i)/9.9);
end
x1=x1+L(i)/100;
end
if(Support_type(nodes)==2)
plot(x1,0,'marker','*','markersize',A(elems)/9.9);
elseif(Support_type(nodes)==1)
plot(x1,0,'marker','o','markersize',A(elems)/9.9);
endif
endfunction

 Results:

•  The overall displacement matrix is obtained and from this the stresses can be obtained.