Having some issues with truss system code

5 views (last 30 days)
I am trying to use http://www.mathworks.com/matlabcentral/fileexchange/38044-truss-analysis to verify my work on a project. Having a few issues with implementing. Not sure if they lie within the connectivity, TP plotting function or something else. Trying to solve a truss system with 12 nodes and 20 members. ST is returning=> Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.549154e-19. > In ST at 46
The TP function is incorrectly plotting 2/3 of the system. I tried messing with for i=1:n with no luck. Also unsure of the plotting of X(:,1)... Do i need to run that for the number of nodes? When I push to X(:,12) I still return an error of index excedes matrix dimension. Tried doing i=1:length(D.Coord(1,:)) with no luck also.
*Ends are pinned supports
Any help is appreciated.I may just be having a major brain fart somewhere and not catching it.
function D=Data
% Definition of Data
% Nodal Coordinates
Coord=[0 0 0;4 0 0;8 0 0;16 0 0;20 0 0;24 0 0;28 0 0;4 4 0;8 8 0;16 8 0;20 8 0; 24 4 0];
% Connectivity
Con=[10 11;12 11;11 4;4 3;3 2;1 2;1 8;8 10;4 12;5 12;5 9;6 5;9 6;7 6;9 7;12 9;3 8;10 3;2 8;4 5];
% Definition of Degree of freedom (free=0 & fixed=1); for 2-D trusses the last column is equal to 1
%Re=zeros(size(Coord));Re(7:10,:)=[1 1 1;1 1 1;1 1 1;1 1 1];
Re=[1 1 1;0 0 1;0 0 1;0 0 1;0 0 1;0 0 1;1 1 1;0 0 1;0 0 1;0 0 1; 0 0 1;0 0 1];
% Definition of Nodal loads
% Load=zeros(size(Coord));Load([1:3,6],:)=1e3*[1 -10 -10;0 -10 -10;0.5 0 0;0.6 0 0];
Load=1e3*[0 0 0;0 -150 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0; 0 0 0];
% Definition of Modulus of Elasticity
E=ones(1,size(Con,1))*45e9;
% or: E=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]*1e7;
% Definition of Area
A=ones(1,size(Con,1))*7.55;
% Convert to structure array
D=struct('Coord',Coord','Con',Con','Re',Re','Load',Load','E',E','A',A');
% D=Data31686814 ; [F,U,R]=ST(D)
% TP(D,U,1)
%function [F,U,R]=ST(D)
% Analyize a Truss using the direct stiffness method
%
% Input D defines the Truss structure as follows
% D.Coord -- N x 3 array of node coordinates
% D.Con -- N x 2 array of connector or member mapping
% D.Re -- N x 3 array of node freedom 1 = fixed 0 = free
% D.Load -- N x 3 array of load force vectors
% D.A -- M x 1 array of member cross section areas
% D.E -- M x 1 array of member Elasticity ( Youngs Modulous)
%
% Ouput F -- M x 1 array of force along members
% U -- N x 3 array of Node displacement vectors
% R -- N x 3 array of Reaction force vectors
History
Original code by Hossein Rahami 17 Mar 2007 (Updated 13 Apr 2007) Reformatted and comments added by Frank McHugh 06 Sep 2012
function [F,U,R]=ST(D)
w=size(D.Re); % 3 x number of nodes
S=zeros(3*w(2)); % stiffness matrix is 3*(number of nodes) square matrix
U=1-D.Re; % U is displacement matrix [
% column index by node
% x , y , z by rows
% initialize U to 1 for non fixed nodes 0 for fixed
f=find(U); % f index in U of free nodes
for i=1:size(D.Con,2) % Loop through Connectors (members)
H=D.Con(:,i);
C=D.Coord(:,H(2))-D.Coord(:,H(1)); % C is vector for connector i
Le=norm(C); % Le length of connector i
T=C/Le; % T is unit vector for connector i
s=T*T'; % Member Siffness matrix is of form
% k * | s -s |
% | -s s | in global truss coordinates
G=D.E(i)*D.A(i)/Le; % G aka k stiffness constant of member = E*A/L
Tj(:,i)=G*T; % Stiffness vector of this member
e=[3*H(1)-2:3*H(1),3*H(2)-2:3*H(2)];
% indexes into Global Stiffness matrix S for this member
S(e,e)=S(e,e)+G*[s -s;-s s];
% add this members stiffness to stiffness matrix
end
U(f)=S(f,f)\D.Load(f); % solve for displacements of free nodes
% ie solve F = S * U for U where S is stiffness
% matrix.
F=sum(Tj.*(U(:,D.Con(2,:))-U(:,D.Con(1,:))));
%project displacement of each node pair on to member
% between
% f = Tj dot ( U2j - U1j ). Then sum over all contributing
% node pairs.
R=reshape(S*U(:),w); % compute forces at all nodes = S*U
R(f)=0; % zero free nodes leaving only reaction.
function TP(D,U,Sc)
C=[D.Coord;D.Coord+Sc*U];
e=D.Con(1,:);
f=D.Con(2,:);
for i=1:6
M=[C(i,e);C(i,f); repmat(NaN,size(e))];
X(:,i)=M(:);
end
plot3(X(:,1),X(:,2),X(:,3),'k',X(:,4),X(:,5),X(:,6),'m');
axis('equal');
if D.Re(3,:)==1;
view(2);
end
  5 Comments
Andrew Stutz
Andrew Stutz on 1 Jul 2016
It was just the order of connectivity. I got it thanks. I'm using R2013b.
Ahmad Azn
Ahmad Azn on 25 Oct 2019
why can i not run this file?
how can i input the arguments of D? ST(D)

Sign in to comment.

Accepted Answer

KSSV
KSSV on 1 Jul 2016
Edited: KSSV on 1 Jul 2016
I have arranged Coord and Con, now you will not get the error. Use the following data.
Coord = [0 0 0
4 0 0
8 0 0
16 0 0
20 0 0
24 0 0
28 0 0
4 4 0
24 4 0
8 8 0
16 8 0
20 8 0 ] ;
Con = [1 2
2 3
3 4
4 5
5 6
6 7
7 9
9 12
12 11
11 10
10 8
8 1
8 2
8 3
10 3
10 4
11 4
12 4
12 5
9 5
9 6];
With the above data, you truss will look like the attached one.

More Answers (1)

KSSV
KSSV on 1 Jul 2016
I guess you have messed up the connectivity matrix. The given coordinates are not connecting properly...with the given data the truss looks like attached file. You have to give correct connections in order to get the required matrices (Mass, Stiffness etc) in FEM calculations.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!