To visualize Nodal deflections of a beam under different Load conditions using Octave. 

Aim:

To calculate and visualize Nodal deflections of a beam applied with different Load conditions using Octave.

Description:

• For any continuous body, there are so many degrees of freedom (i.e, any property of the body can be calculated and infinitely many points). So, in order to calculate or approximate the behavior of the body "Finite Element Approach" is followed.
• In a finite element approach, the body is split into discrete elements and the properties are calculated only at nodes that are finite in number. This process is called 'Discretisation'.
• And these properties are interpolated so as to calculate the property at some intermittent point.

Methodology

To solve any finite element problem the following steps are followed :

1. First, the body under consideration is split into smaller elements based on load, supports, and geometrical constraints.
2. Here each element is considered analogous to be a stiff spring connected between the nodes.
3. So the total force is given by the relation :
4. Here the stiffness value is given by- (where L is the length of each element, E is youngs modulus and I is the moment of the inertia)
5. Likewise, the stiffness matrix for each element is calculated.
6. All these mattresses are joined to obtain a "Global stiffness matrix".
7. Using this matrix the slope and deflection are calculated at each node.
8. Since we have slope and displacement at each node the displacement at any point in the element is modeled using a "Bezier curve".
9. Now the displacement at each point is calculated using the bezier approximation and the element is colored accordingly.

Steps to execute the program:

INPUTS:

After executing the program all the data need to be entered as prompted :

Results:-

Code :

Main Program :

%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 elementn',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);
bezier=[2 -2 1 1;-3 3 -2 -1;0 0 1 0;1 0 0 0];
x1=0;
mi=Inf;
mx=-Inf;
partition=100;
for i=1:elems
k=linspace(x1,x1+L(i),partition);
for x=1:length(k)
u=((k(x)-x1)/L(i));
nod_dis(x,i)=[u^3 u^2 u 1]*bezier*[q(2*i-1) q(2*i+1) q(2*i) q(2*(i+1))]';
end
for x=1:length(k)
if(nod_dis(x,i)<mi)
mi=nod_dis(x,i);
endif
if(nod_dis(x,i)>mx)
mx=nod_dis(x,i);
endif
end
x1=x1+L(i);
end
me=(mx+mi)/2;
x1=0;
figure(2)
hold on;
for i=1:elems
k=linspace(x1,x1+(L(i)/100),partition);
for x=1:length(k)
u=(k(x)-x1)/(L(i)/100);
if(u==0)
tmp_x=x1;
endif
x2=tmp_x+((L(i)/100)/(partition-1));
Qe_1(x,i)=(1/((me-mi)*(mx-me)))*((nod_dis(x,i))^2-(mx+mi)*(nod_dis(x,i))+me*(mx+mi-me));
Qe_2(x,i)=(1/((me-mi)*(me-mx)))*((nod_dis(x,i))^2-(mx+mi)*(nod_dis(x,i))+mi*mx);
if(nod_dis(x,i)>me)
c_r=0;
c_g=Qe_2(x,i);
c_b=Qe_1(x,i);
else
c_r=Qe_1(x,i);
c_b=0;
c_g=Qe_2(x,i);
endif
c=[c_r c_g c_b];
plot([tmp_x x2],[0 0],'color',c,'linewidth',A(i)/5);
tmp_x=x2;
end
x1=x1+(L(i)/100);
end

Visualize 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