Analysis of a plane truss using MATLAB (FINAL INDEPENDENT PROJECT)
In this project, we will write a MATLAB program in order to determine the internal member forces of a plane truss.
Now let us consider a prismatic bar whose length is L, Area of the cross-section is A, made up of a material of Young's modulus E, subjected to external axial load say P, Now considering the bar to be elastic,
According to Hook's Law
`delta = (P*L)/(A*E)`
`implies P=delta*(AE)/L`
This `(AE)/L` is known as the stiffness of the bar.
In other words, Stiffness is the force required to produce a unit amount of deflection in the bar.
If we know the stiffness of a member we can get the deflection.
In order to analyze the truss, we will follow the stiffness approach.
In a plane truss, the truss element can have any orientation, unlike the one-dimensional bar, the local coordinates of a plane Truss element is generally not be aligned in the same direction as the Global Axis for the truss although they are in the same plane. In order to define the displaced configuration of the plane Truss element wrt it's original configuration, it is necessary to define two degrees of freedom at each end of the element does there are four degrees of freedom to be considered.
The local coordinates are denoted as 1 and 2 in the X and Y directions respectively at the starting node that is left node and as 3 and 4 in the X and Y direction respectively at the ending node that is the right node. Elastic deformations are possible only along the longitudinal axis direction that is x-axis in a truss element, whereby based on Simple compatibility requirements the expression for member elongation e in terms of the member end Axial translations `D_1` and `D_3` is easily obtained as
e = `D_3` - `D_1` _______ eq 1
Axial tension developed in the element is given by
`N =k*e=(AE)/L*e`
Now from the figure
`F_1 = -(AE)/L*e` and
`F_3 = (AE)/L*e`
The truss element has zero stiffness in y-direction shows rigid body displacement in the y-direction.
`impliesF_1=-(AE)/L*(D_3-D_1)`
`impliesF_1=(AE)/L*(D_1-D_3)`
`andF_2=0`
`andF_3=(AE)/L*(D_3-D_1)`
`andF_4=0`
so in matrix from,
`[[F_1],[F_2],[F_3],[F_4]]=(AE)/L*[[1,0,-1,0],[0,0,0,0],[-1,0,1,0],[0,0,0,0]]*[[D_1],[D_2],[D_3],[D_4]]`
this is the stiffness matrix in terms of local axis, to convert it into interms of Global axis of the entire structure we have to do some manupulation
Now consider that the local axis is making an angle of `theta` with global axis
So, `D_1=D_(1G)costheta+D_(2G)sin theta`
`D_2=-D_(1G)sintheta+D_(2G)cos theta`
`D_3=D_(3G)costheta+D_(4G)sin theta`
`D_4=-D_(3G)sintheta+D_(4G)cos theta`
and similarly
`F_1=F_(1G)costheta+F_(2G)sin theta`
`F_2=-F_(1G)sintheta+F_(2G)cos theta`
`F_3=F_(3G)costheta+F_(4G)sin theta`
`F_4=-F_(3G)sintheta+F_(4G)cos theta`
`[[F_1],[F_2],[F_3],[F_4]]=(AE)/L*[[1,0,-1,0],[0,0,0,0],[-1,0,1,0],[0,0,0,0]]*[[D_1],[D_2],[D_3],[D_4]]`
`implies[[F_(1G)],[F_(2G)],[F_(3G)],[F_(4G)]]*[[costheta,sintheta,0,0],[-sintheta,costheta,0,0],[0,0,costheta,sintheta],[0,0,-sintheta,costheta]]=(AE)/L*[[1,0,-1,0],[0,0,0,0],[-1,0,1,0],[0,0,0,0]]*[[D_(1G)],[D_(2G)],[D_(3G)],[D_(4G)]]*[[costheta,sintheta,0,0],[-sintheta,costheta,0,0],[0,0,costheta,sintheta],[0,0,-sintheta,costheta]]`
`implies[[F_(1G)],[F_(2G)],[F_(3G)],[F_(4G)]]=(AE)/L*[[cos^2theta,costhetasintheta,-cos^2theta,-costhetasintheta],[costhetasintheta,sin^2theta,-costhetasintheta,-sin^2theta],[-cos^2theta,-costhetasintheta,cos^2theta,costhetasintheta],[-costhetasintheta,-sin^2theta,costhetasintheta,sin^2theta]]*[[D_(1G)],[D_(2G)],[D_(3G)],[D_(4G)]]`
Now we have our stiffness matrix and our transformation matrix which will convert our local stiffness matrix into global stiffness matrix.
Now let us consider an actual truss problem as shown in figure
So this is a truss having 5 members
now there are 4 joints so each of them will have two degrees of freedom so there will be total of 8 degrees of freedom but as the supports are restrained at 6,7 and 8 so these are not free to move so these are restrained degrees of freedom.
So Global labels of members and there lengths are
| Member | Labels | Length |
| AB | [ 6 4 5 3 ] | 3m |
| BC | [ 4 8 3 7 ] | 3m |
| BD | [ 4 2 3 1 ] | 3m |
| AD | [ 6 2 5 1 ] | 4.243m |
| DC | [ 2 8 1 7 ] | 4.243m |
Let us do the first phase of the program that is defining the input parameter
%% stiffness matrix method for planar truss
% Input
clc
clear
n = 5; % number of members
L = [3 3 3 4.243 4.243]; % length in m
A = [1 1 1 1 1]; % Area
theta= [0 0 90 45 -45]; % angle in degrees
uu = 5; % Number of unrestrained degrees of freedom
ur = 3; % Number of restrained degrees of freedom
uul = [1 2 3 4 5]; % global labels of unrestained dof
url = [6 7 8]; % global labels of restained dof
l1 = [6 4 5 3]; % Global labels for member 1
l2 = [4 8 3 7]; % Global labels for member 2
l3 = [4 2 3 1]; % Global labels for member 3
l4 = [6 2 5 1]; % Global labels for member 4
l5 = [2 8 1 7]; % Global labels for member 5
l= [l1; l2; l3; l4; l5];
dof = uu + ur; % Degrees of freedom
Ktotal = zeros (dof);
Tt1 = zeros (4); % Transformation matrix for member 1
Tt2 = zeros (4); % Transformation matrix for member 2
Tt3 = zeros (4); % Transformation matrix for member 3
Tt4 = zeros (4); % Transformation matrix for member 4
Tt5 = zeros (4); % Transformation matrix for member 5
fem1= [0; 0; 0; 0]; % Local Fixed end forces of member 1
fem2= [0; 0; 0; 0]; % Local Fixed end forces of member 2
fem3= [0; 0; 0; 0]; % Local Fixed end forces of member 3
fem4= [0; 0; 0; 0]; % Local Fixed end forces of member 4
fem5= [0; 0; 0; 0]; % Local Fixed end forces of member 5now let us go for the second phase that is the formulation of stiffness matrix
%% stiffness matrix 4 by 4
for i = 1:n
Knew = zeros (dof);
k1 = [0; 0; 0; 0];
k2 = [0; 0; 0; 0];
k3 = [0; 0; rc(i); -rc(i)];
k4 = -k3;
K = [k1 k2 k3 k4];
fprintf ('Member Number =');
disp (i);
fprintf ('Local Stiffness matrix of member, [K] = \n');
disp (K);
T1 = [cx(i); 0; cy(i); 0];
T2 = [0; cx(i); 0; cy(i)];
T3 = [-cy(i); 0; cx(i); 0];
T4 = [0; -cy(i); 0; cx(i)];
T = [T1 T2 T3 T4];
fprintf ('Tranformation matrix of member, [T] = \n');
disp (T);
Ttr = T';
fprintf ('Tranformation matrix Transpose, [T] = \n');
disp (Ttr);
Kg = Ttr*K*T;
fprintf ('Global Matrix, [K global] = \n');
disp (Kg);
for p = 1:4
for q = 1:4
Knew((l(i,p)),(l(i,q))) =Kg(p,q);
end
end
Ktotal = Ktotal + Knew;
if i == 1
Tt1= T;
Kg1=Kg;
fembar1= Tt1'*fem1;
elseif i == 2
Tt2 = T;
Kg2 = Kg;
fembar2= Tt2*fem2;
elseif i == 3
Tt3 = T;
Kg3 = Kg;
fembar3= Tt3'*fem3;
elseif i == 4
Tt4 = T;
Kg4 = Kg;
fembar4= Tt4'*fem4;
else
Tt5 = T;
Kg5=Kg;
fembar5= Tt5*fem5;
end
end
fprintf ('Stiffness Matrix of complete structure, [Ktotal] = \n');
disp (Ktotal);
Kunr = zeros(ur);
for x=1:uu
for y=1:uu
Kunr(x,y)= Ktotal(x,y);
end
end
fprintf ('Unrestrained Stiffness sub-matrix, [Kuu] = \n');
disp (Kunr);
KuuInv= inv(Kunr);
fprintf ('Inverse of Unrestrained Stiffness sub-matrix, [KuuInverse] =\n');
disp (KuuInv);Now we have obtained the local stiffness matrix for the individual members and then from that, we have deduced the global stiffness matrix, now as we already know that
load = stiffness X deflection
`implies` deflection = load X stiffnes`s^-1`
now as we know that labels 6,7 and 8 are restrained so there will be no deflection so what we will do we will from a sub-matrix of global stiffness matrix only for the unrestrained labels and we will do it's inverse and multiplying it with global load vector we can obtain the deflection at respective nodes.
so let define the load vector and obtain the displacement at respective nodes.
once we have the displacement at the respective nodes it is quite easy to get the member force in respective member
as Load = stiffness X deflection
here the obtained load will be member force acting on that member due to the external load, so let's do this
%% Creation of joint load vector
jl= [40; -20; 0; 0; 0; 0; 0; 0]; % values given in kN
jlu = [40; -20; 0; 0; 0]; % load vector in unrestrained dof
delu = KuuInv*jlu;
fprintf ('Joint Load vector, [Jl] = \n');
disp (jl');
fprintf ('Unrestrained displacements, [DelU] = \n');
disp (delu');
delr = zeros(ur,1);
del = zeros (dof,1);
del = [delu; delr];
deli= zeros (4,1);
for i = 1:n
for p = 1:4
deli(p,1) = del((l(i,p)),1) ;
end
if i == 1
delbar1 = deli;
mbar1= (Kg1 * delbar1)+fembar1;
fprintf ('Member Number =');
disp (i);
fprintf ('Global displacement matrix [DeltaBar] = \n');
disp (delbar1');
fprintf ('Global End force matrix [MBar] = \n');
disp (mbar1');
elseif i == 2
delbar2 = deli;
mbar2= (Kg2 * delbar2)+fembar2;
fprintf ('Member Number =');
disp (i);
fprintf ('Global displacement matrix [DeltaBar] = \n');
disp (delbar2');
fprintf ('Global End force matrix [MBar] = \n');
disp (mbar2');
elseif i ==3
delbar3 = deli;
mbar3= (Kg3 * delbar3)+fembar3;
fprintf ('Member Number =');
disp (i);
fprintf ('Global displacement matrix [DeltaBar] = \n');
disp (delbar3');
fprintf ('Global End force matrix [MBar] = \n');
disp (mbar3');
elseif i == 4
delbar4 = deli;
mbar4= (Kg4 * delbar4)+fembar4;
fprintf ('Member Number =');
disp (i);
fprintf ('Global displacement matrix [DeltaBar] = \n');
disp (delbar4');
fprintf ('Global End force matrix [MBar] = \n');
disp (mbar4');
else
delbar5 = deli;
mbar5= (Kg5 * delbar5)+fembar5;
fprintf ('Member Number =');
disp (i);
fprintf ('Global displacement matrix [DeltaBar] = \n');
disp (delbar5');
fprintf ('Global End force matrix [MBar] = \n');
disp (mbar5');
end
endso now let see the outputs
Member Number =
1
Local Stiffness matrix of member, [K] =
0 0 0 0
0 0 0 0
0 0 0.3333 -0.3333
0 0 -0.3333 0.3333
Tranformation matrix of member, [T] =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Tranformation matrix Transpose, [T] =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Global Matrix, [K global] =
0 0 0 0
0 0 0 0
0 0 0.3333 -0.3333
0 0 -0.3333 0.3333
Member Number =
2
Local Stiffness matrix of member, [K] =
0 0 0 0
0 0 0 0
0 0 0.3333 -0.3333
0 0 -0.3333 0.3333
Tranformation matrix of member, [T] =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Tranformation matrix Transpose, [T] =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Global Matrix, [K global] =
0 0 0 0
0 0 0 0
0 0 0.3333 -0.3333
0 0 -0.3333 0.3333
Member Number =
3
Local Stiffness matrix of member, [K] =
0 0 0 0
0 0 0 0
0 0 0.3333 -0.3333
0 0 -0.3333 0.3333
Tranformation matrix of member, [T] =
0 0 -1 0
0 0 0 -1
1 0 0 0
0 1 0 0
Tranformation matrix Transpose, [T] =
0 0 1 0
0 0 0 1
-1 0 0 0
0 -1 0 0
Global Matrix, [K global] =
0.3333 -0.3333 0 0
-0.3333 0.3333 0 0
0 0 0 0
0 0 0 0
Member Number =
4
Local Stiffness matrix of member, [K] =
0 0 0 0
0 0 0 0
0 0 0.2357 -0.2357
0 0 -0.2357 0.2357
Tranformation matrix of member, [T] =
0.7071 0 -0.7071 0
0 0.7071 0 -0.7071
0.7071 0 0.7071 0
0 0.7071 0 0.7071
Tranformation matrix Transpose, [T] =
0.7071 0 0.7071 0
0 0.7071 0 0.7071
-0.7071 0 0.7071 0
0 -0.7071 0 0.7071
Global Matrix, [K global] =
0.1178 -0.1178 0.1178 -0.1178
-0.1178 0.1178 -0.1178 0.1178
0.1178 -0.1178 0.1178 -0.1178
-0.1178 0.1178 -0.1178 0.1178
Member Number =
5
Local Stiffness matrix of member, [K] =
0 0 0 0
0 0 0 0
0 0 0.2357 -0.2357
0 0 -0.2357 0.2357
Tranformation matrix of member, [T] =
0.7071 0 0.7071 0
0 0.7071 0 0.7071
-0.7071 0 0.7071 0
0 -0.7071 0 0.7071
Tranformation matrix Transpose, [T] =
0.7071 0 -0.7071 0
0 0.7071 0 -0.7071
0.7071 0 0.7071 0
0 0.7071 0 0.7071
Global Matrix, [K global] =
0.1178 -0.1178 -0.1178 0.1178
-0.1178 0.1178 0.1178 -0.1178
-0.1178 0.1178 0.1178 -0.1178
0.1178 -0.1178 -0.1178 0.1178
Stiffness Matrix of complete structure, [Ktotal] =
0.2357 0 0 0 -0.1178 -0.1178 -0.1178 0.1178
0 0.5690 0 -0.3333 -0.1178 -0.1178 0.1178 -0.1178
0 0 0.6667 0 -0.3333 0 -0.3333 0
0 -0.3333 0 0.3333 0 0 0 0
-0.1178 -0.1178 -0.3333 0 0.4512 0.1178 0 0
-0.1178 -0.1178 0 0 0.1178 0.1178 0 0
-0.1178 0.1178 -0.3333 0 0 0 0.4512 -0.1178
0.1178 -0.1178 0 0 0 0 -0.1178 0.1178
Unrestrained Stiffness sub-matrix, [Kuu] =
0.2357 0 0 0 -0.1178
0 0.5690 0 -0.3333 -0.1178
0 0 0.6667 0 -0.3333
0 -0.3333 0 0.3333 0
-0.1178 -0.1178 -0.3333 0 0.4512
Inverse of Unrestrained Stiffness sub-matrix, [KuuInverse] =
5.7430 1.5000 1.5000 1.5000 3.0000
1.5000 5.7430 1.5000 5.7430 3.0000
1.5000 1.5000 3.0000 1.5000 3.0000
1.5000 5.7430 1.5000 8.7430 3.0000
3.0000 3.0000 3.0000 3.0000 6.0000
Joint Load vector, [Jl] =
40 -20 0 0 0 0 0 0
Unrestrained displacements, [DelU] =
199.7200 -54.8600 30.0000 -54.8600 60.0000
Member Number =
1
Global displacement matrix [DeltaBar] =
0 -54.8600 60.0000 30.0000
Global End force matrix [MBar] =
0 0 10.0000 -10.0000
Member Number =
2
Global displacement matrix [DeltaBar] =
-54.8600 0 30.0000 0
Global End force matrix [MBar] =
0 0 10.0000 -10.0000
Member Number =
3
Global displacement matrix [DeltaBar] =
-54.8600 -54.8600 30.0000 199.7200
Global End force matrix [MBar] =
1.0e-15 *
0.1691 -0.1691 0 0
Member Number =
4
Global displacement matrix [DeltaBar] =
0 -54.8600 60.0000 199.7200
Global End force matrix [MBar] =
-10.0000 10.0000 -10.0000 10.0000
Member Number =
5
Global displacement matrix [DeltaBar] =
-54.8600 0 199.7200 0
Global End force matrix [MBar] =
-30.0000 30.0000 30.0000 -30.0000now got the member forces, so get the picture more clear the results are shown in the figure as a free body
Now to cross-check the solution we can perform two operation
`1.sumf_X=0`
`2.sumf_Y=0`
And clearly we can see that our solution is justified.
So now conclusion, in order to design the truss member it is very much essential to determine the member forces as we can clearly see that the Member AD is subjected to the most critical loading, so we have to design the structure according to that. By this method using MATLAB, we can get the member forces for critical truss arrangements also.