I used this code for another truss and it worked and now I am getting an error
12 次查看(过去 30 天)
显示 更早的评论
This is the error I keep getting:
Index in position 1 exceeds array bounds (must not exceed 4).
Error in CEE210_CP4_Franke_IDK (line 289)
SupportReaction(j,:))
Here is my code: (Any help would be greatly appreciated!)
%% **********************************************************************
% Start of Input Section
% Input data for Benchmark Truss........................
% Units used:
LUnit = 'ft'; % Length = ft or m
FUnit = 'lb'; % Force = lb or kip or N or kN
%.. Input problem data
nNodes = 5; % Number of nodes (pin joints) in the truss
nMembers = 9; % Number of members
nReactions = 4; % Number of force vector reactions
nLoads = 2; % Number of externally applied loads
nDim = 3; % Dimensions 2=2D, 3=3D
% *** Once your code is modified and tested, nDim must always equal 3 and
% the variable nDim cannot be used in your code.
nDOF = nDim * nNodes; % Number of degrees of freedom
tol = 1.e-8; % Remainder tolerance on norm of residual force
% Input coordinates of each node
Node(1,:) = [0, 0, 12]; % A
Node(2,:) = [5, -5, 0]; % B
Node(3,:) = [5, 5, 0]; % C
Node(4,:) = [-5, 5, 0]; % D
Node(5,:) = [-5, -5, 0]; % E
% Input member connections
%.. First input is the start node; Second input is the end node
Member(1,:) = [1,2]; % AB
Member(2,:) = [1,3]; % AC
Member(3,:) = [1,4]; % AD
Member(4,:) = [1,5]; % AE
Member(5,:) = [2,3]; % BC
Member(6,:) = [3,4]; % CD
Member(7,:) = [4,5]; % DE
Member(8,:) = [5,2]; % DG
Member(9,:) = [2,4]; % EB
% Input support reactions
%.. First input is the node number
%.. Second input is resistance in the e1 direction
%.. Third input is resistance in the e2 direction
%.. Fourth input is resistance in the e3 direction (zero for 2-D trusses)
%*** 1 for resistance in that direction, 0 for no resistance ***
Reaction(1,:) = [2, 1,1,1]; % Pin at Node 1
Reaction(2,:) = [3, 0,0,1]; % Roller at Node 5
Reaction(3,:) = [4, 1,0,0];
Reaction(4,:) = [5, 0,0,0];
% Input externallly applied loads
%.. First input is the node number, remaining inputs are the magnitude of
%.. the force in the e1, e2 and e3 directions, respectively.
Load(1,:) = [1, 0, 0, -2400]; % 1000 lb load at Node 2
Load(2,:) = [4, 0, -1200, 0]; % 1000 lb load at Node 4
% End of Input Section for Benchmark Truss......
% ************************************************************************
%% Compute the unit vector (direction) for every member
PosVec = zeros(nMembers, 3); %Initialize Position Vector matrix
UnVec = zeros(nMembers, 3); %Initialize Unit Vector matrix
OppUnVec = zeros(nMembers, 3); %Initialize opposite direction matrix
Length = zeros(nMembers, 3); %Initialize Length matrix
for i = 1 : nMembers
SN = Member(i,1);
EN = Member(i,2);
PosVec(i,:) = Node(EN,:)- Node(SN,:);
Length(i,:) = norm(PosVec(i,:));
UnitVec(i,:) = PosVec(i,:)/Length(i,:);
negUnitVec(i,:) = - UnitVec(i,:);
OppUnVec(i,:) = -PosVec(i,:) / norm(PosVec(i,:));
end
%% The coefficient matrix, C, contains the unit vectors (direction vectors)
% for each member force associated with each node. The unit vectors are
% then separated into e1, e2, e3 components and placed in separate rows.
% Therefore, each node has 2 rows (for 2-D) or 3 rows (for 3-D) associated
% with it; one for each of the e1, e2, e3 components.
% This corresponds with the Degrees of Freedom, nDOF, which is 2*nNodes
% for 2-D or 3*nNodes for 3-D.
% So the coefficient matrix has 2 or 3 rows for each node and one column
% for each member force.
C = zeros(3*nNodes, nMembers); %Initialize coefficient matrix
% Loop through all nodes, create the vector force equation for that node
% and store it in the proper row of the coefficient matrix, C.
for i = 1 : nNodes
for j = 1 : nMembers
SN = Member(j,1);
EN = Member(j,2);
if (SN == i) % If member j starts at node i
C(3*i-2, j) = UnVec(j,1);
C(3*i-1, j) = UnVec(j,2);
C(3*i, j) = UnVec(j,3);
elseif (EN == i) % If member j ends at node i
C(3*i-2, j) = OppUnVec(j,1);
C(3*i-1, j) = OppUnVec(j,2);
C(3*i, j) = OppUnVec(j,3);
end % if SN
end % for nMembers
end % for nNodes
%% The External Force matrix, F, contains the magnitude of all externally
% applied loads, (input as e1, e2, e3 components), stored in the
% proper rows to corespond with the node it is applied at.
% Therefore, the F matrix has 2 or 3 rows for each node and one column.
F = zeros(3*nNodes, 1); %Initialize External Force matrix
% Loop through all externally applied loads and place each component of the
% load (e1, e2, e3) in the proper rows in the external load matrix, F.
for i = 1 : nLoads
j = Load(i,1); % Which node, j is the load i on
F(3*j-2) = Load(i,2);
F(3*j-1) = Load(i,3);
F(3*j) = Load(i,4);
end % for i,nLoads
%% The Restraint matrix, Crest, contains the unit vectors (directions)
% for each member force associated with those nodes that are restrained by
% external supports.
% The Crest matrix has one row for each reaction component and one column
% for each member.
%.... Loop through all nodes and determine if there is a reaction at that
%.... node and seperate the reaction node equations from the force load
%.... equations.
nReaction = 1; %Used to count the number of reaction equations
%.... The nReactionEquation vector determines if a row should be put in the
%.... Crest matrix or the Cfree matrix
nReactionEquation = zeros(1, 3*nNodes);
for (j = 1 : nReactions)
for (i = 1 : nNodes)
if (Reaction(j,1) == i)
if (Reaction(j,2) == 1)
Crest(nReaction,:) = C(3*i-2,:);
nReactionEquation(3*i-2) = 1;
% If nReactionEquation is 1 it should be in the Crest
% matrix and not the Cfree matrix.
nReaction = nReaction+1;
end
if (Reaction(j,3) == 1)
Crest(nReaction,:) = C(3*i-1,:);
nReactionEquation(3*i-1) = 1;
% If nReactionEquation is 1 it should be in the Crest
% matrix and not the Cfree matrix.
nReaction = nReaction+1;
end
if (Reaction(j,4) == 1)
Crest(nReaction,:) = C(3*i,:);
nReactionEquation(3*i) = 1;
% If nReactionEquation is 1 it should be in the Crest
% matrix and not the Cfree matrix.
nReaction = nReaction+1;
end
end
end
end
nForce = 1; %Used to count the number of Cfree equations
for (i = 1 : 3*nNodes)
if (nReactionEquation(i) == 1)
else
Cfree(nForce,:) = C(i,:);
Ffree(nForce) = F(i);
nForce = nForce + 1;
end
end
%% Solve first for the member forces and then solve for the support
% reactions.
% All of the values in the Cfree and Ffree matrices are known so you can
% now solve for the T matrix, which is the tension force in each member.
T = -(Cfree)\(Ffree'); % finish this equation
% Since all values in the Crest and T matrices are now known, you can solve
% for the support reactions.
Reactions = -Crest*T; % finish this equation
% Calculate the residual to find the amount of error and ensure that
% equilibrium was satisfied.
Residual1 = (Cfree * T + Ffree');
Residual2 = (Crest * T + Reactions);
Res = norm(Residual1) + norm(Residual2);
% Create the support reaction matrix.
nReaction = 1;
for (i = 1 : nNodes)
if (nReactionEquation(3*i-2) == 1)
SupportReaction(i,1) = i;
SupportReaction(i,2) = Reactions(nReaction);
nReaction = nReaction+1;
end
if (nReactionEquation(3*i-1) == 1)
SupportReaction(i,1) = i;
SupportReaction(i,3) = Reactions(nReaction);
nReaction = nReaction+1;
end
if (nReactionEquation(3*i) == 1)
SupportReaction(i,1) = i;
SupportReaction(i,4) = Reactions(nReaction);
nReaction = nReaction+1;
end
end
%% Create output for command window
fprintf('%s\n' , '----------------------------------------')
fprintf('%s\n' , '------------- Truss --------------')
fprintf('%s\n' , '----------------------------------------')
fprintf('%s',' Node Coordinates (', LUnit, ')')
fprintf ('\n')
for i = 1 : nNodes
fprintf('%s%4i%8.3f%8.3f%8.3f\n',' Node: ', i, Node(i,:)')
end % for i, nNodes
fprintf ('\n')
fprintf('%s',' Member Forces (', FUnit, ')')
fprintf ('\n')
for i = 1 : nMembers
fprintf('%s%4i%12.3f\n',' Member Force: ',i, T(i))
end % for i, nMembers
fprintf ('\n')
for i = 1 : nReactions
for j = 1 : nNodes
if (Reaction(i,1) == j)
fprintf('%s%s%s%4i%12.3f%12.3f%12.3f\n',...
' Support Reaction (',FUnit,') at Node ',...
SupportReaction(j,:))
fprintf ('\n')
end % if Reaction
end % for j, nNodes
end % for i, nReactions
fprintf('\n%s%8.3f\n',' Residual Error: ', Res)
fprintf('\n\n\n')
%% Create the plot............................................
% Once your code is modified the plot will be 3-D and so you may need the
% "Rotate 3D" command to view the truss figure properly.
%.. Plot members with color indicating tension or compression
Marker = ceil(25 / sqrt(nNodes));
fig1 = figure(1); clf; grid on; axis equal;
hold on;
xlabel(cat(2, 'X (',LUnit,')' ));
ylabel(cat(2, 'Y (',LUnit,')' ));
zlabel(cat(2, 'Z (',LUnit,')' ));
title('Truss Analysis');
small = max(T)*tol;
for m = 1 : nMembers
if (T(m) < -small)
MColor = 'b'; % Color compression members blue
elseif (T(m) > small)
MColor = 'r'; % Color tension members red
else
MColor = 'k'; % Color "zero force members" black
end % if T
SN = Member(m,1);
EN = Member(m,2);
X = [Node(SN,1); Node(EN,1)];
Y = [Node(SN,2); Node(EN,2)];
Z = [Node(SN,3); Node(EN,3)];
p = plot3(X,Y,Z);
set(p, 'Color', MColor, 'LineWidth', 6/nDim);
end % for m, nMembers
view(3);
%.. Plot external reaction forces
Rlength = 0.1 * max(max(Node)); % establish size of line
e = [ 1, 0, 0 ; 0, 1, 0 ; 0, 0, 1];
for i = 1 : nReactions
for j = 1 : nNodes
if (Reaction(i,1) == j)
for k = 1 : 2
if(Reaction(i, k+1) == 1)
xR(1,:) = Node(j,:);
xR(2,:) = Node(j,:) - Rlength * e(k, 1:nDim);
xR(3,:) = Node(j,:) - Rlength * e(k, 1:nDim);
p = plot3(xR(:,1),xR(:,2),xR(:,3));
set(p,'Color','g','LineWidth',8/nDim);
end % if Reaction == 1
end % for k
end % if Reaction == j
end % for j, nNodes
end % for i, nReactions
%.. Plot loads
Llength = 0.1 * max(max(Node)); % establish size of line
Fmax = max(max(abs(Load))); % establish maximum force level
e = [ 1, 0, 0 ; 0, 1, 0 ; 0, 0, 1];
for i = 1 : nLoads
for j = 1 : nNodes
if (Load(i,1) == j)
for k = 1 : 2
F = Load(i, k+1) / Fmax;
xL(1,:) = Node(j,:);
xL(2,:) = Node(j,:) + Llength * F * e(k, 1:nDim);
xL(3,:) = Node(j,:) + Llength * F * e(k, 1:nDim);
p = plot3(xL(:,1), xL(:,2), xL(:,3));
set(p,'Color','c','LineWidth',8/nDim);
end % for k
end % if Load == j
end % for j, nNodes
end % for i, nLoads
%.. Plot nodes
plot3(Node(:,1), Node(:,2),Node(:,3), 'o', 'LineWidth', 2,...
'MarkerFaceColor','w',...
'MarkerEdgeColor','k',...
'MarkerSize',Marker);
view(3);
% End of Program
0 个评论
回答(1 个)
Walter Roberson
2020-11-11
for (i = 1 : nNodes)
if (nReactionEquation(3*i-2) == 1)
SupportReaction(i,1) = i;
SupportReaction(i,2) = Reactions(nReaction);
nReaction = nReaction+1;
end
if (nReactionEquation(3*i-1) == 1)
SupportReaction(i,1) = i;
SupportReaction(i,3) = Reactions(nReaction);
nReaction = nReaction+1;
end
if (nReactionEquation(3*i) == 1)
SupportReaction(i,1) = i;
SupportReaction(i,4) = Reactions(nReaction);
nReaction = nReaction+1;
end
end
That code only writes into SupportReaction(i,:) if at least one of those three conditions is true, and you are not pre-allocating SupportReaction. If there is an i for which none of the three tests is true, then SupportReaction(i,:) will not be written that time. However if there is a later i for which one of the tests is true, then writing into that later SupportReaction(i,:) would implicitly extend SupportReaction to have that many rows.
Therefore, SupportReaction will have as many rows as the last i for which one of the conditions hold. And that implies that SupportReaction will have fewer rows than nNodes if there are trailing i for which none of the conditions match.
K>> reshape(nReactionEquation,3,[])
ans =
0 1 0 1 0
0 1 0 0 0
0 1 1 0 0
We can see from this at a glance that none of the conditions match for the final i in your data.
You could get around this problem by simply pre-allocating
SupportReaction = zeros(3, nNodes);
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Structural Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!