L3: Solution of Truss Using Stiffness Matrix Method + Solve Truss using MATLAB Code

% MATLAB Code: Analysis of Truss using Stiffness Matrix Method, generalised Truss
% Just change the co-ordinate and connectivity and properties of Materials
close all; clear all;
clc;
format short;
format compact;

%-----------------Data Input and Pre Processing

% Units: mm, Newtons, MPa

length_unit = '[m]'; % For printing the result
force_unit = '[N]';
stress_unit = '[Pa]';

E = 200E9; %(200GPa for steel)
A = 0.0005; % (b=50mm, d=10mm)

% nodes and connection
nodes = [0 0; -10 5.7735; -10 0];
conn = [1 3; 1 2; 2 3];

% The no of elements in the model
nel = length(conn);

% The no of nodes in the model
nnode = length(nodes);

% Total no of degrees of freedom in the model
ndof= 2*nnode;

%Allocate space for the vector of element length
lengths = zeros(1,nel);

% Allocate space for the vector of element inclination angle
angles = zeros(1,nel);

% Allocate memory for the stiffness matrix
K = zeros(ndof,ndof);

% vector and rotation matrix
kel_store = cell(nel,1);
index_store = cell(nel,1);
rot_store = cell(nel,1);

% Calculate the length and angle of each element
for i=1:nel
n1 = conn(i,1);
n2 = conn(i,2);
x1 = nodes(n1,1);
y1 = nodes(n1,2);
x2 = nodes(n2,1);
y2 = nodes(n2,2);
dx = x2-x1;
dy = y2-y1;
angles(i) = atan2(dy,dx);
lengths(i) = sqrt(dx^2+dy^2);
anglesDEg(i)= rad2deg(angles(i));
end
% To display length and angle of the element
for i = 1:nel
disp(['Length of the element', num2str(i), '=', num2str(lengths(i)), length_unit])
disp(['Angle of the element', num2str(i), '=', num2str(angles(i)*180/pi), '[degrees]'])
disp(' ')
end
f = logical([0 0 1 0 1 1]);
s = not(f); % specified (fixed) DOFs
D(s) = 0; % we specify the displacements of s dofs
% Index vector, used for populating local stiffness matrix into global stiffness matrix
% here in the following programme {1} means the no of element ande the no in the '[]' represent the DOFs of the end nodes of the respected element
index_store{1} = [5 6 1 2];
index_store{2} = [1 2 3 4];
index_store{3} = [3 4 5 6];

% Global stiffness matrix assembly

for i = 1:nel % Now for each element in turn
ind = index_store{i};
L = lengths(i);
theta = angles(i)
ElementNo(i)=i
anglesDEg(i)= rad2deg(theta)
S = sin(theta)
C = cos(theta)
% Generate the element rotation matrix
T = [C, S, 0, 0; -S, C, 0, 0; 0, 0, C, S; 0, 0, -S, C]
%TTT=transpose(T)*T
Tmat_store{i} = T;
% Generate the element stiffness matrix in the elements local coordinate system
% condition for using different sectional property
C1 = (A*E/L);
% Rotate the element's stiffness matrix to global coordinates\
kel = C1*[1, 0, -1, 0; 0, 0, 0, 0; -1, 0, 1, 0; 0, 0, 0, 0]
kelBtr_s{i} = kel; %Befroe Transformation
kel_trans = T'*kel*T; % Element stiffness matrix in the global coordinate system
% Store the stiffness matrix for later use:
kelAtr_s{i} = kel_trans
% Assemble system Stiffness matrix(add the element matrix in the right spot)
K(ind,ind) = K(ind,ind) + kel_trans;

end

fprintf('The stiffness matrix is as follows \n')
disp(K);
% Solution by Adding Stiffness in Diagonal term
K(2,2)=K(2,2)*10000;
K(5,5)=K(5,5)*10000;
K(6,6)=K(6,6)*10000;

force=zeros(6,1);
force(3,1)=10;
force(4,1)=-20;
U1M=K\force;

U=zeros(6,1)
U(1,1)=U1M(1,1);
U(3,1)=U1M(3,1);
U(4,1)=U1M(4,1);

RF1M=K*U
% by Row column elemination
KC(1,1)=K(1,1);
KC(2,1)=K(3,1);
KC(3,1)=K(4,1);

KC(1,2)=K(1,3);
KC(2,2)=K(3,3);
KC(3,2)=K(4,3);

KC(1,3)=K(1,4);
KC(2,3)=K(3,4);
KC(3,3)=K(4,4);

fr(1,1)=0;
fr(2,1)=10;
fr(3,1)=-20;
Kinv=inv(KC)

U2M=KC\fr;

RF2M=K*U