MATLAB: Having some issues with truss system code

direct stiffnessstructurestruss

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

Best Answer

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.
Related Question