Truss Analysis script - matrix close to singular, incorrect results

4 views (last 30 days)
Marcin Morawski
Marcin Morawski on 21 Feb 2016
Edited: Marcin Morawski on 21 Feb 2016
I am trying to model a simple bridge in 2D based on a catenary and an A – truss. I used the Truss Analysis file from the file exchange ( http://uk.mathworks.com/matlabcentral/fileexchange/14313-truss-analysis ), however the results I get do not make any sense.
The input structure is:
Coord =
-4.0000 0.0000 0
-2.7386 0.5988 0
-1.3910 0.9641 0
-0.0000 1.0862 0
1.3910 0.9641 0
2.7386 0.5988 0
4.0000 0.0000 0
-6.6667 0 0
-4.0000 0.9567 0
-2.7386 1.4093 0
-1.3910 1.8928 0
-0.0000 2.3918 0
1.3910 1.8928 0
2.7386 1.4093 0
4.0000 0.9567 0
6.6667 0 0
Con =
1 2
2 3
3 4
4 5
5 6
6 7
8 9
9 10
10 11
11 12
12 13
13 14
14 15
15 16
1 9
2 10
3 11
4 12
5 13
6 14
7 15
Re =
1 1 1
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
1 1 1
1 1 1
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
0 0 1
1 1 1
E =
1.0e+11 *
Columns 1 through 8
2.1000 2.1000 2.1000 2.1000 2.1000 2.1000 0.0900 0.0900
Columns 9 through 16
0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.1100 0.1100
Columns 17 through 21
0.1100 0.1100 0.1100 0.1100 0.1100
A =
Columns 1 through 8
0.0100 0.0100 0.0100 0.0100 0.0100 0.0100 0.0400 0.0400
Columns 9 through 16
0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0100 0.0100
Columns 17 through 21
0.0100 0.0100 0.0100 0.0100 0.0100
I think that the connections and coordinates are correct, because the truss analysis script plots the correct structure:
The problem is, that while for a load (10kN) placed at a single node, the script gives results which look correct:
F =
1.0e+04 *
Columns 1 through 8
-0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -1.4806 -1.4806
Columns 9 through 16
-1.4806 -1.4806 -1.4806 -1.4806 -1.4806 -1.4806 0.0000 -0.0000
Columns 17 through 21
-0.0000 1.0000 -0.0000 -0.0000 -0.0000
U =
1.0e-03 *
Columns 1 through 8
0 -0.0357 -0.0870 0.0224 0.1769 -0.0791 0 0
0 0.0751 0.2646 -0.9813 0.7780 -0.1666 0 0
0 0 0 0 0 0 0 0
Columns 9 through 16
-0.1238 -0.2093 -0.3398 0.0000 0.5240 0.1226 0.1238 0
0.0000 0.0751 0.2646 -0.8626 0.7780 -0.1666 -0.0000 0
0 0 0 0 0 0 0 0
R =
1.0e+04 *
Columns 1 through 8
0.0000 0 0 0 0 0 -0.0000 1.3936
0.0000 0 0 0 0 0 0.0000 0.5000
0 0 0 0 0 0 0 0
Columns 9 through 16
0 0 0 0 0 0 0 -1.3936
0 0 0 0 0 0 0 0.5000
0 0 0 0 0 0 0 0
But if I use two much smaller loads (100N) placed on two nodes, it looks like the whole thing collapses (the deflections are of the order 10^9 m), which, physically, just doesn't make sense:
F =
Columns 1 through 8
-256.0000 -128.0000 -96.0000 -64.0000 -32.0000 -24.0000 -145.1047 -96.0000
Columns 9 through 16
-192.0000 -160.0000 -136.0000 -128.0000 -135.5000 -135.2881 12.3376 0
Columns 17 through 21
0 89.7262 -7.0598 -4.0451 0.1431
U =
1.0e+09 *
Columns 1 through 8
0 0.6130 -0.0095 0.0788 0.1199 -0.0161 0 0
0 -1.2912 1.0058 -0.0000 0.4681 -0.0340 0 0
0 0 0 0 0 0 0 0
Columns 9 through 16
-0.0000 0.4632 -0.3608 0.0000 0.1680 -0.0122 0.0000 0
0.0000 -1.2912 1.0058 -0.0000 0.4681 -0.0340 0.0000 0
0 0 0 0 0 0 0 0
R =
Columns 1 through 8
128.0000 0 0 0 0 0 -24.0000 136.5807
115.6624 0 0 0 0 0 9.8569 49.0011
0 0 0 0 0 0 0 0
Columns 9 through 16
0 0 0 0 0 0 0 -127.3407
0 0 0 0 0 0 0 45.6861
0 0 0 0 0 0 0 0
It also gives the following error with every calculation I make:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =
1.881839e-18.
> In ST at 8
I tried checking the whole structure a number of times, changing the dimensions/stiffness of members to very high values etc. but the result is always the same - 1 force, it holds, 2 forces - it breaks.
The code I used to generate the data is:
% calculating the points on the catenary
x1 = linspace (-4, 4, 1000);
y1 = -2.0*(cosh(x1/4.0))+3.0862;
% calculating the location of equally spaced joints on the catenary
pt1 = interparc (7, x1, y1, 'spline');
% calculating the data for plotting the triangular section
theta = atan(1.5/4); % calculating the angle
% calculating the joints on the triangular section
pt2 = [-6.666667, 0; pt1(1:3,1), theta.*(pt1(1:3,1)+6.666667); ...
pt1(4:7,1), -theta.*(pt1(4:7,1)-6.666667); 6.666667, 0];
% matrix containing the coordinates of all the joints in 3D
Coord = [pt1, zeros(size(pt1,1),1); pt2, zeros(size(pt2,1),1)];
% matrix of connections
Con = [];
for n=1:15
if n==7
continue
else
Con = [Con; n, n+1];
end
end
for n=1:7
Con = [Con; n, n+8];
end
% matrix of degrees of freedom
n = 2;
Re = NaN(16, 3);
Re (1,:) = [1,1,1];
while n>1 & n<7
Re (n,:) = [0,0,1];
n=n+1;
end
while n==7 | n==8
Re (n,:) = [1,1,1];
n=n+1;
end
while n>8 & n<16
Re (n,:) = [0,0,1];
n=n+1;
end
Re (16,:) = [1,1,1];
% Matrix containing the nodal loads
Load = zeros(size(Coord));
Load (4, :) = [0, -100, 0];
Load (3, :) = [0, -100, 0];
% Matrix containing the modulus of elasticity
E = NaN(1, 21);
for n=1:6
E(n)= 210E9;
end
for n=7:14
E(n)= 9E9;
end
for n=15:21
E(n)= 11E9;
end
% Definition of cross sectional areas
A = NaN(1,21);
for n=1:6
A(n)= (10e-2)^2;
end
for n=7:14
A(n)= (20e-2)^2;
end
for n=15:21
A(n)= (10e-2)^2;
end
% Convert to structure array
D=struct('Coord',Coord','Con',Con','Re',Re','Load',Load','E',E','A',A');
Does anyone know what might be causing the issue?

Sign in to answer this question.

Answers (0)

Categories

Find more on Structural Analysis in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!