Why Does MATLAB Show NaN for Truss Deflections?

  • Thread starter Thread starter nate9519
  • Start date Start date
  • Tags Tags
    Matlab Matrices
Click For Summary
SUMMARY

The discussion centers on the issue of receiving "NaN" (not a number) outputs for truss deflections in MATLAB due to a singular matrix problem. The warnings indicate that the system may be partially constrained and that the matrix is singular to working precision. The root cause is identified as the absence of specified initial displacements for at least one node, which is crucial for solving the stiffness matrix equations. Users are advised to ensure that initial conditions are set, particularly at support points, to avoid singular results.

PREREQUISITES
  • Understanding of MATLAB programming, specifically version R2023a.
  • Familiarity with structural analysis concepts, particularly truss design.
  • Knowledge of matrix operations and linear algebra.
  • Experience with defining boundary conditions in finite element analysis.
NEXT STEPS
  • Review MATLAB's matrix manipulation functions, particularly inv() and det().
  • Learn about specifying boundary conditions in finite element models.
  • Study the implications of singular matrices in structural analysis.
  • Explore resources on truss analysis and stiffness matrix formulation.
USEFUL FOR

Structural engineers, MATLAB users involved in finite element analysis, and students studying mechanics of materials will benefit from this discussion.

nate9519
Messages
47
Reaction score
0

Homework Statement


I'm using MATLAB to design a truss and the output is suppose to be the deflections caused by the applied loads. But I keep getting "not a number" for my deflections . here is the warnings I get.

Warning: System may be partially constrained.
> In truss3 (line 222)
In truss3ex2 (line 26)
Warning: Matrix is singular to working precision.
> In truss3 (line 224)
In truss3ex2 (line 26)

Homework Equations


relevant lines from truss3 starting at line 219

Matlab:
% Solve the force-displacement equations.AA = [a,zero1;e,d]; BB = [b,zero2]'; detAA = det(AA);

if abs(detAA) < 1e-16

warning('System may be partially constrained.')

end

AAinv = inv(AA); sol = AAinv*BB;

% Prepare the output quantities.

pp = 0;

for i = 1:nJforce(i,1) = i; Jdispl(i,1) = i;

isol = 2*m+p+3*(i-1);

Jdispl(i,2) = sol(isol+1); Jdispl(i,3) = sol(isol+2);

Jdispl(i,4) = sol(isol+3);

[\code][h2]The Attempt at a Solution[/h2] 
[/B]
here is the code.

[code=matlab]
% class design project example%

% all of the members are quenched steel

% k = 2000 k-lb, Pmax = 1500 k-lb (A = 100 sq-inch)

%

clear

n = 24; m = 60; LOADZ = -20000; LOADY = 1000; A = 100;

joint = [0,0,0;0,0,-10;0,0,-20;12,0,0;12,0,-10;12,0,-20;...

0,10,0;0,10,-10;12,10,0;12,10,-10;0,20,0;0,20,-10;...

12,20,0;12,20,-10;0,30,0;0,30,-10;12,30,0;12,30,-10;...

0,40,0;0,40,-10;0,40,-20;12,40,0;12,40,-10;12,40,-20];

assembly = [1,2;1,4;1,7;2,3;2,4;2,5;2,6;2,7;2,8;3,6;3,8;...

4,5;4,9;5,6;5,9;5,10;6,10;7,8;7,9;7,11;8,9;8,10;8,11;...

8,12;9,10;9,13;10,13;10,14;11,12;11,13;11,14;11,15;11,16;...

12,14;12,16;13,14;13,17;13,18;14,18;15,16;15,17;15,19;15,20;...

16,17;16,18;16,20;16,21;17,18;17,22;17,23;18,23;18,24;19,20;...

19,22;20,22;20,23;20,24;21,24;22,23;23,24];

forceJ = [3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;...

3,1,1,1;-1,0,0,0;-1,0,0,0;...

-1,0,0,0; -1,0,0,0;-1,0,LOADY,LOADZ;-1,0,0,0;...

-1,0,LOADY,LOADZ;-1,0,0,0;-1,0,0,0;-1,0,0,0;-1,0,0,0;-1,0,0,0;...

3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1];

for i = 1:m; stretch(i) = 2000*1000; end; stretch(1);

%

index = 1;

[Jforce,Mforce,Jdispl,Mdispl] = ...                                                ([B]LINE 26)[/B]

truss3(n,m,joint,assembly,forceJ,stretch,index);

%

peak_klb = 18*A

maxMforce_klb = max(abs(Mforce/1000))

maxJdispl = max(abs(Jdispl*12));

maxDX_in = maxJdispl(2),maxDY_in = maxJdispl(3),maxDZ_in = maxJdispl(4)

[\code]for forceJ, (3,1,1,1) indicates a external joint with 3 reaction forces. (-1,0,0,0) indicates an internal joint with no loads.
Can anyone help me figure out this problem with the matrices? I have no idea what to do
 
Last edited:
Physics news on Phys.org
Hi, you want to use ##[##code=matlab##]##
##[##code=matlab##]##

(code goes here)​

##[##\code##]##
to keep it legible:

Matlab:
% Solve the force-displacement equations.

AA = [a,zero1;e,d]; BB = [b,zero2]'; detAA = det(AA);

if abs(detAA) < 1e-16
    warning('System may be partially constrained.')
end

AAinv = inv(AA); sol = AAinv*BB;

% Prepare the output quantities.

pp = 0;

for i = 1:n

   Jforce(i,1) = i; Jdispl(i,1) = i;

   isol = 2*m+p+3*(i-1);

   Jdispl(i,2) = sol(isol+1); Jdispl(i,3) = sol(isol+2);

   Jdispl(i,4) = sol(isol+3);

end

% 3. The Attempt at a Solution 
% here is the code.

% class design project example
%
% all of the members are quenched steel
% k = 2000 k-lb, Pmax = 1500 k-lb (A = 100 sq-inch)
%

clear

n = 24; m = 60; LOADZ = -20000; LOADY = 1000; A = 100;

joint = [0,0,0;0,0,-10;0,0,-20;12,0,0;12,0,-10;12,0,-20;...

0,10,0;0,10,-10;12,10,0;12,10,-10;0,20,0;0,20,-10;...

12,20,0;12,20,-10;0,30,0;0,30,-10;12,30,0;12,30,-10;...

0,40,0;0,40,-10;0,40,-20;12,40,0;12,40,-10;12,40,-20];

assembly = [1,2;1,4;1,7;2,3;2,4;2,5;2,6;2,7;2,8;3,6;3,8;...

4,5;4,9;5,6;5,9;5,10;6,10;7,8;7,9;7,11;8,9;8,10;8,11;...

8,12;9,10;9,13;10,13;10,14;11,12;11,13;11,14;11,15;11,16;...

12,14;12,16;13,14;13,17;13,18;14,18;15,16;15,17;15,19;15,20;...

16,17;16,18;16,20;16,21;17,18;17,22;17,23;18,23;18,24;19,20;...

19,22;20,22;20,23;20,24;21,24;22,23;23,24];

forceJ = [3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;...

3,1,1,1;-1,0,0,0;-1,0,0,0;...

-1,0,0,0; -1,0,0,0;-1,0,LOADY,LOADZ;-1,0,0,0;...

-1,0,LOADY,LOADZ;-1,0,0,0;-1,0,0,0;-1,0,0,0;-1,0,0,0;-1,0,0,0;...

3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1;3,1,1,1];

for i = 1:m; stretch(i) = 2000*1000; end; stretch(1);

%

index = 1;

[Jforce,Mforce,Jdispl,Mdispl] = ...                                     %           ([B]LINE 26)[/B]

truss3(n,m,joint,assembly,forceJ,stretch,index);

%

peak_klb = 18*A

maxMforce_klb = max(abs(Mforce/1000))

maxJdispl = max(abs(Jdispl*12));

maxDX_in = maxJdispl(2),maxDY_in = maxJdispl(3),maxDZ_in = maxJdispl(4)
 
Seems like this is old news.

In order for the stiffness matrix equations to have a solution, at least one or more nodes need to have an initial displacement specified. Thus usually is accomplished by specifying a zero displacement at the support points of the structure.

It's not clear from your code where such initial conditions are specified by the user. Without being able to specify an initial displacement for at least one node, the stiffness matrix equations are probably giving a singular result when the deflections are being solved for.

See p. 11 of the following article:

http://www.engr.sjsu.edu/ragarwal/ME273/pdf/Chapter 4 - Beam Element.pdf

More details can be found in the attached file below.
 

Attachments

Last edited:

Similar threads

  • · Replies 1 ·
Replies
1
Views
2K
Replies
1
Views
2K
Replies
2
Views
3K