%MATLAB PROGRAM FOR 2D FRAME
clear all;clc;close all;%clear all existing variables
%Obtain input file(Excel sheet)
[num]=xlsread('C:\Users\Roopak\Desktop\SEM-1\ASA\project\stiffness_input_l_i.xlsx')
xy_node=num(:,1:2);%coordinates vector
xy_node(any(isnan(xy_node),2),:)=[];
mem=num(:,4:5);
mem(any(isnan(mem),2),:)=[];
E=num(:,6);
E(any(isnan(E),2),:)=[]
A=num(:,7);
A(any(isnan(A),2),:)=[];
I=num(:,8);
I(any(isnan(I),2),:)=[];
dof_r=num(:,10:12);
dof_r(any(isnan(dof_r),2),:)=[];
ext_load=num(:,14:16);
ext_load(any(isnan(ext_load),2),:)=[];
no_elem=size(mem,1);% Found number of elements
no_node=size(xy_node,1);% founf number of nodes
dof_all=1:3*no_node; %martix with all dofs
sms=zeros(3*no_node); %initialised global stiffness matrix with zeros
d=zeros(3*no_node,1); %initialised displacement matrix with zeros
a=zeros(3*no_node,1); %initialised force vector with zeros
dof_res=[]; %initialised an empty matrix
for j=1:size(dof_r,1) %to find the dof corr. to restrains
dof_flag=3*(dof_r(j,1)-1)+dof_r(j,2);
dof_res=[dof_res dof_flag]; %concatinating restrained dof
d(dof_flag)=dof_r(j,3); %setting the input displacement corr. to restrains
end
dof_free=dof_all;%duplicating all dof matrix
dof_free(dof_res)=[];%deletinf the restained dofs
for j=1:size(ext_load,1)% Extracting loads from input file
a(3*(ext_load(j,1)-1)+ext_load(j,2))=ext_load(j,3);
end
for i=1:no_elem %For loop to obtain element stiffness matrix
elem_nodes=mem(i,1:2); %Present Element Node
elem_cord=xy_node(elem_nodes,:); % Current element coordinate
E1=[(elem_cord(2,1)-elem_cord(1,1)) (elem_cord(2,2)-elem_cord(1,2))];
len_elem=sqrt((E1(1)^2)+((E1(2))^2)); %length of element
E1=E1/len_elem; % Direction Cosine
E2=[-E1(2) E1(1)]; % Rotational matrix
sm_2356=[(((12*E(i)*I(i)))/(len_elem^3)) (((6*E(i)*I(i)))/(len_elem^2)) -
(((12*E(i)*I(i)))/(len_elem^3))
(((6*E(i)*I(i)))/(len_elem^2));(((6*E(i)*I(i)))/(len_elem^2))
(((4*E(i)*I(i))/len_elem)) -(((6*E(i)*I(i)))/(len_elem^2))
(((2*E(i)*I(i))/len_elem));-(((12*E(i)*I(i)))/(len_elem^3))
-(((6*E(i)*I(i)))/(len_elem^2)) (((12*E(i)*I(i)))/(len_elem^3))
-(((6*E(i)*I(i)))/(len_elem^2));(((6*E(i)*I(i)))/(len_elem^2))
(((2*E(i)*I(i))/len_elem)) -(((6*E(i)*I(i)))/(len_elem^2))
(((4*E(i)*I(i))/len_elem))];
sm_14=(((E(i)*A(i))/len_elem))*[1 -1;-1 1];
sm_loc([1,4],[1,4])=sm_14;
sm_loc([2,3,5,6],[2,3,5,6])=sm_2356;
trans=[E1; E2];
trans(3,3)=1;
R=[trans zeros(3);zeros(3) trans]; %Rotational matrix
k_elem=R'*sm_loc*R;% element Stiffness matrix obtained
dof_elem=3*(elem_nodes(1)-1)+1:3*(elem_nodes(1));
dof_elem=[dof_elem 3*(elem_nodes(2)-1)+1:3*(elem_nodes(2))];% Obtained dofs of
element
sms(dof_elem,dof_elem)=sms(dof_elem,dof_elem)+k_elem; %Concatinated element
stiffness with Global stiffness matrix
end
d(dof_free)=((sms(dof_free,dof_free))^-1)*(a(dof_free)-
(sms(dof_free,dof_res)*d(dof_res)))% obtained Displacements
a(dof_res)=sms(dof_res,:)*d%obtained actions
for i=1:no_elem % For loop to obtain member forces
elem_nodes=mem(i,1:2);
elem_cord=xy_node(elem_nodes,:);
E1=[(elem_cord(2,1)-elem_cord(1,1)) (elem_cord(2,2)-elem_cord(1,2))];
len_elem=sqrt((E1(1)^2)+((E1(2))^2));
E1=E1/len_elem;
E2=[-E1(2) E1(1)];
sm_bend1=[(((12*E(i)*I(i)))/(len_elem^3)) (((6*E(i)*I(i)))/(len_elem^2)) -
(((12*E(i)*I(i)))/(len_elem^3))
(((6*E(i)*I(i)))/(len_elem^2));(((6*E(i)*I(i)))/(len_elem^2))
(((4*E(i)*I(i))/len_elem)) -(((6*E(i)*I(i)))/(len_elem^2))
(((2*E(i)*I(i))/len_elem));-(((12*E(i)*I(i)))/(len_elem^3))
-(((6*E(i)*I(i)))/(len_elem^2)) (((12*E(i)*I(i)))/(len_elem^3))
-(((6*E(i)*I(i)))/(len_elem^2));(((6*E(i)*I(i)))/(len_elem^2))
(((2*E(i)*I(i))/len_elem)) -(((6*E(i)*I(i)))/(len_elem^2))
(((4*E(i)*I(i))/len_elem))];
sm1=(((E(i)*A(i))/len_elem))*[1 -1;-1 1];
sm_loc1([1,4],[1,4])=sm1;
sm_loc1([2,3,5,6],[2,3,5,6])=sm_bend1;
trans=[E1; E2];
trans(3,3)=1 ;
n1=3*(elem_nodes(1)-1)+1:3*(elem_nodes(1));
n2=3*(elem_nodes(2)-1)+1:3*(elem_nodes(2));
mem_f([1,2,3],i)=((sm_loc1([1,2,3],[1,2,3])*trans*d(n1))+(sm_loc1([1,2,3],
[4,5,6])*trans*d(n2)))%member forces ij
mem_f([4,5,6],i)=((sm_loc1([4,5,6],[1,2,3])*trans*d(n1))+(sm_loc1([4,5,6],
[4,5,6])*trans*d(n2)))%member forces ji
end
f_temp([1,2,3],1)=mem_f([1,2,3],1);
f_temp([1,2,3],2)=mem_f([4,5,6],19);
f_temp([1,2,3],3)=mem_f([1,2,3],21);
f_temp([1,2,3],4)=mem_f([4,5,6],39);
f_temp([1,2,3],5)=mem_f([1,2,3],41);
f_temp([1,2,3],6)=mem_f([4,5,6],60);
f_temp([1,2,3],7)=mem_f([1,2,3],61);
p=[0,0,0,0,0,0,0;-14.6261,-14.6261,-14.6261,-14.6261,-14.6261,-14.6261,-
14.6261,;0,0,0,0,0,0,0];
tij=[0,1,0;-1,0,0;0,0,1];
reac(3,7)=zeros;
for i=1:7
reac(:,i)=(tij*f_temp([1,2,3],i))-p(:,i);
end
reac %Reaction at the supports
mulim1=447178;
vulim1=556.659;
pulim1=4993.832;
dup=1:no_elem;
dup([3,6,8,9,10,11,13,15,17,19,23,25,27,29,30,32,34,36,38,40,42,44,46,48,50,51,53,5
5,57,59])=[];
for i=1:size(dup',1)
r(dup(i))=(abs(mem_f(1,dup(i)))/vulim1)+(abs(mem_f(2,dup(i)))/pulim1)+
(abs(mem_f(3,dup(i)))/mulim1);
end
dup1=[3,6,8,9,10,11,13,15,17,19,23,25,27,29,30,32,34,36,38,40,42,44,46,48,50,51,53,
55,57,59];
mulim=670769;
vulim=694.29;
pulim=4993.832;
for i=1:size(dup1',1)
r(dup1(i))=(abs(mem_f(1,dup1(i)))/pulim)+(abs(mem_f(2,dup1(i)))/vulim)+
(abs(mem_f(3,dup1(i)))/mulim);
end
true=mem_f'
true1=r'
