Finite elements model using STABIL
We are working on a finite elements model using the plug-in STABIL. We get no errors but our displacements are all giving NaN as solution. Can someone explain the meaning of these NaNs? Our K and P matrix look logical so we are out of options…
%initializing
clc
close all
clear all
% units kN,m
unzip(‘stabil-3.1.zip’);
% CALCULATING THE NODES
%STEP N1: the deck
% Nodes=[NodID X Y Z]
Nodes= [1 21 0 4;
2 21+5.25*1 0 4;
3 21+5.25*2 0 4;
4 21+5.25*3 0 4;
5 21+5.25*4 0 4;
6 21+5.25*5 0 4;
7 21+5.25*6 0 4;
8 21+5.25*7 0 4;
9 21+5.25*8 0 4];
Nodes = reprow(Nodes, 1:9, 1, [9 0 4.2 0]);
Nodes = reprow(Nodes, 10:18, 1, [9 0 4.5 0]);
Nodes = reprow(Nodes, 19:27, 1, [9 0 4.2 0]);
Nodes = reprow(Nodes, 28:36, 1, [9 0 3.7 0]);
%STEP N2: the 3 top points of triangle, then copy to second triangle
Nodes= [Nodes;
46 31.5 0 14.5;
47 42 0 25;
48 52.5 0 14.5];
Nodes = reprow(Nodes, 46:48, 1, [3 0 12.9 0]);
%STEP N3: reprow the full triangle + deck two times
%no overlapping nodes
%second triangle + deck
Nodes = reprow(Nodes, 1:8,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 10:17,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 19:26,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 28:35,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 37:44,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 46:51,1,[51 42 0 0]);
%third triangle + deck
Nodes = reprow(Nodes, 1:8,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 10:17,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 19:26,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 28:35,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 37:44,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 46:51,1,[102 84 0 0]);
%STEP N4: add the other notes of the deck parts not directly supported by
%the triangles to the left and right
%The left
Nodes = [Nodes;
154 0 0 4;
155 5.25 0 4;
156 10.5 0 4;
157 15.75 0 4;
158 0 4.2 4;
159 5.25 4.2 4;
160 10.5 4.2 4;
161 15.75 4.2 4;
162 0 8.7 4;
163 5.25 8.7 4;
164 10.5 8.7 4;
165 15.75 8.7 4;
166 0 12.9 4;
167 5.25 12.9 4;
168 10.5 12.9 4;
169 15.75 12.9 4;
170 0 16.6 4;
171 5.25 16.6 4;
172 10.5 16.6 4;
173 15.75 16.6 4];
%The right
Nodes = [Nodes;
174 152.25 0 4;
175 157.5 0 4;
176 162.75 0 4;
177 168 0 4;
178 152.25 4.2 4;
179 157.5 4.2 4;
180 162.75 4.2 4;
181 168 4.2 4;
182 152.25 8.7 4;
183 157.5 8.7 4;
184 162.75 8.7 4;
185 168 8.7 4;
186 152.25 12.9 4;
187 157.5 12.9 4;
188 162.75 12.9 4;
189 168 12.9 4;
190 152.25 16.6 4;
191 157.5 16.6 4;
192 162.75 16.6 4;
193 168 16.6 4];
%STEP N5: reference nodes
%Longitudinal beams and diagonal beams triangle 1 to 47 and copies
Nodes= [Nodes;
500 0 0 100;
501 0 4.2 100;
502 0 8.7 100;
503 0 12.9 100;
504 0 16.6 100];
%Beam element from 47 to 9 and all the copies
%also the ones going from 46 to 5 and the copies
Nodes= [Nodes;
505 168 0 100;
506 168 12.9 100];
%Transversal beams
Nodes= [Nodes;
507 31.5 0 100;
508 42 0 100;
509 52.5 0 100;
510 31.5+42 0 100;
511 42+42 0 100;
512 52.5+42 0 100;
513 31.5+42*2 0 100;
514 42+42*2 0 100;
515 52.5+42*2 0 100];
%STEP E1: initializing
% Element types -> {EltTypID EltName}
Types= {1 ‘beam’;
2 ‘truss’};
% Sections=[SecID A ky kz Ixx Iyy Izz] !!!all in m
Sections = [1 0.265*0.990*2 Inf Inf 2*(1/12*0.265*0.990^3) 2*(1/12*0.265*0.990^3) 2*(1/12*0.265*0.990^3); % plywood king post beam
2 0.265*0.540*2 Inf Inf 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3); % plywood diagonal beam
3 0.265*0.540*2 Inf Inf 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3); % horizontal beam
4 ((0.314/2)^2*pi-(0.314/2-0.019)^2*pi) NaN NaN NaN NaN NaN; % steel vertical
5 ((0.314/2)^2*pi-(0.314/2-0.019)^2*pi) NaN NaN NaN NaN NaN; % diagonal steel truss (only A)
6 0.265*0.1350*2 Inf Inf 2*(1/12*0.265*0.1350^3) 2*(1/12*0.265*0.1350^3) 2*(1/12*0.265*0.1350^3); % outer deck beam
7 0.240*0.1350*2 Inf Inf 2*(1/12*0.240*0.1350^3) 2*(1/12*0.240*0.1350^3) 2*(1/12*0.240*0.1350^3); % inner deck beam
8 0.200*0.650 Inf Inf 1/12*0.200*0.650^3 1/12*0.200*0.650^3 1/12*0.200*0.650^3; % outer footbridge
9 0.100*0.100 NaN NaN NaN NaN NaN]; % deck truss steel
% Materials=[MatID E nu rho];
Materials= [1 7.75e6 0.2 2000; % plywood
2 210e6 0.3 7500]; % steel KN/m2
% Elements=[EltID TypID SecID MatID n1 n2 n3]
%STEP E2: deck longitudonal plywood beams
%First triangle
Elements= [1 1 6 1 1 2 500;
2 1 7 1 10 11 501;
3 1 7 1 19 20 502;
4 1 6 1 28 29 503;
5 1 8 1 37 38 504];
Elements=reprow(Elements,1:5,7,[5 0 0 0 1 1 0]);
%Second triangle
Elements=[Elements;
41 1 6 1 9 52 500;
42 1 7 1 18 61 501;
43 1 7 1 27 70 502;
44 1 6 1 36 79 503;
45 1 8 1 45 88 504];
Elements= [Elements;
46 1 6 1 52 53 500;
47 1 7 1 61 62 501;
48 1 7 1 70 71 502;
49 1 6 1 79 80 503;
50 1 8 1 88 89 504];
Elements=reprow(Elements,46:50,6,[5 0 0 0 1 1 0]);
%Third triangle
Elements=[Elements;
81 1 6 1 59 103 500;
82 1 7 1 68 112 501;
83 1 7 1 77 121 502;
84 1 6 1 86 130 503;
85 1 8 1 95 139 504];
Elements= [Elements;
86 1 6 1 103 104 500;
87 1 7 1 112 113 501;
88 1 7 1 121 122 502;
89 1 6 1 130 131 503;
90 1 8 1 139 140 504];
Elements=reprow(Elements,86:90,6,[5 0 0 0 1 1 0]);
%Deck on the right of the triangles
Elements= [Elements;
121 1 6 1 110 174 500;
122 1 7 1 119 178 501;
123 1 7 1 128 182 502;
124 1 6 1 137 186 503;
125 1 8 1 146 190 504];
Elements= [Elements;
126 1 6 1 174 175 500;
127 1 7 1 178 179 501;
128 1 7 1 182 183 502;
129 1 6 1 186 187 503;
130 1 8 1 190 191 504];
Elements=reprow(Elements,126:130,2,[5 0 0 0 1 1 0]);
%Deck on the left of the triangles
Elements= [Elements;
141 1 6 1 154 155 500;
142 1 7 1 158 159 501;
143 1 7 1 162 163 502;
144 1 6 1 166 167 503;
145 1 8 1 170 171 504];
Elements=reprow(Elements,141:145,2,[5 0 0 0 1 1 0]);
Elements= [Elements;
156 1 6 1 157 1 500;
157 1 7 1 161 10 501;
158 1 7 1 165 19 502;
159 1 6 1 169 28 503;
160 1 8 1 173 37 504];
%STEP E2: deck transversal steel trusses
%First triangle
Elements= [Elements;
161 2 9 1 1 10 NaN;
162 2 9 1 10 19 NaN;
163 2 9 1 19 28 NaN;
164 2 9 1 28 37 NaN];
Elements=reprow(Elements,161:164,8,[4 0 0 0 1 1 0]);
%Second triangle
Elements= [Elements;
197 2 9 1 52 61 NaN;
198 2 9 1 61 70 NaN;
199 2 9 1 70 79 NaN;
200 2 9 1 79 88 NaN];
Elements=reprow(Elements,197:200,7,[4 0 0 0 1 1 0]);
%Third triangle
Elements= [Elements;
229 2 9 1 103 112 NaN;
230 2 9 1 112 121 NaN;
231 2 9 1 121 130 NaN;
232 2 9 1 130 139 NaN];
Elements=reprow(Elements,229:232,7,[4 0 0 0 1 1 0]);
%Right of the triangles
Elements= [Elements;
261 2 9 1 174 178 NaN;
262 2 9 1 178 182 NaN;
263 2 9 1 182 186 NaN;
264 2 9 1 186 190 NaN];
Elements=reprow(Elements,261:264,3,[4 0 0 0 1 1 0]);
%Left of the triangles
Elements= [Elements;
277 2 9 1 154 158 NaN;
278 2 9 1 158 162 NaN;
279 2 9 1 162 166 NaN;
280 2 9 1 166 170 NaN];
Elements=reprow(Elements,277:280,3,[4 0 0 0 1 1 0]);
%STEP E3: Beams elements triangles
%Beam element from 1 to 47 and all the copies (king post truss up)
Elements= [Elements;
293 1 1 1 1 47 500;
294 1 1 1 9 98 500;
295 1 1 1 59 149 500;
296 1 1 1 28 50 503;
297 1 1 1 36 101 503;
298 1 1 1 86 152 503];
%Beam element from 47 to 9 and all the copies (king post truss down)
Elements= [Elements;
299 1 1 1 47 9 505;
300 1 1 1 98 59 505;
301 1 1 1 149 110 505;
302 1 1 1 50 36 506;
303 1 1 1 101 86 506;
304 1 1 1 152 137 506];
%Beam element from 5 to 48 and all the copies (side diagonal)
Elements= [Elements;
305 1 2 1 5 48 500;
306 1 2 1 55 99 500;
307 1 2 1 106 150 500;
308 1 2 1 32 51 503;
309 1 2 1 82 102 503;
310 1 2 1 133 153 503];
%Beam element from 46 to 5 and all the copies (side diagonal)
Elements= [Elements;
311 1 2 1 46 5 505;
312 1 2 1 97 55 505;
313 1 2 1 148 106 505;
314 1 2 1 49 32 506;
315 1 2 1 100 82 506;
316 1 2 1 151 133 506];
%Transversal beams between the triangles
Elements= [Elements;
317 1 3 1 46 49 507;
318 1 3 1 47 50 508;
319 1 3 1 48 51 509;
320 1 3 1 97 100 510;
321 1 3 1 98 101 511;
322 1 3 1 99 102 512;
323 1 3 1 148 151 513;
324 1 3 1 149 152 514;
325 1 3 1 150 153 515];
%STEP E4: Truss elements triangles
%Vertical steel trusses
%First triangle
Elements= [Elements;
326 2 4 2 46 3 NaN;
327 2 4 2 47 5 NaN;
328 2 4 2 48 7 NaN;
329 2 4 2 49 30 NaN;
330 2 4 2 50 32 NaN;
331 2 4 2 51 34 NaN];
%Second triangle
Elements= [Elements;
332 2 4 2 97 53 NaN;
333 2 4 2 98 55 NaN;
334 2 4 2 99 57 NaN;
335 2 4 2 100 80 NaN;
336 2 4 2 101 82 NaN;
337 2 4 2 102 84 NaN];
%Third triangle
Elements= [Elements;
338 2 4 2 148 104 NaN;
339 2 4 2 149 106 NaN;
340 2 4 2 150 108 NaN;
341 2 4 2 151 131 NaN;
342 2 4 2 152 133 NaN;
343 2 4 2 153 135 NaN];
%Diagonal trusses
%First triangle
Elements= [Elements;
344 2 5 2 46 50 NaN;
345 2 5 2 47 49 NaN;
346 2 5 2 47 51 NaN;
347 2 5 2 48 50 NaN];
%Second triangle
Elements= [Elements;
348 2 5 2 97 101 NaN;
349 2 5 2 101 99 NaN;
350 2 5 2 98 100 NaN;
351 2 5 2 98 102 NaN];
%Third triangle
Elements= [Elements;
352 2 5 2 148 152 NaN;
353 2 5 2 150 152 NaN;
354 2 5 2 149 151 NaN;
355 2 5 2 149 153 NaN];
%STEP E5: plot
figure
plotelem(Nodes,Elements(find(Elements(:,3)==1),:),Types,’bl’);
hold(‘on’);
plotelem(Nodes,Elements(find(Elements(:,3)==2),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==3),:),Types,’m’);
plotelem(Nodes,Elements(find(Elements(:,3)==4),:),Types,’r’);
plotelem(Nodes,Elements(find(Elements(:,3)==5),:),Types,’r’);
plotelem(Nodes,Elements(find(Elements(:,3)==6),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==7),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==8),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==9),:),Types,’r’);
title(‘Elements’)
figure
plotnodes(Nodes);
title(‘Nodes’)
%STEP D1: generate all possible DOF
DOFall = getdof(Elements,Types);
%STEP D2: remove the DOF limited by the boundary conditions
selnodes = [1;9;28;36;59;86;110;137;154;158;162;166;170;177;181;185;189;193];
dofpattern = [0.01 0.02 0.03]; % the three dof’s
seldof = selnodes + dofpattern; % for the selected nodes, the pattern of dof
seldof = seldof(:);
DOF = removedof(DOFall,seldof);
%STEP K1: assemble the K matrix
K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF);
%STEP K2: check the K matrix
%check symmetry
disp(K-K’)
det(K)
issymmetric(K)
%positive definite
find(eig(K)<0)
%to check wheter invertible
%Now we assign the correct loads to the structure
%STEP L1: compute the self weight of the structure
% Own weight
DLoadsOwn=accel([0 0 9.81e-3],Elements,Types,Sections,Materials);
%STEP L2: compute the distributed load working on the deck
%Combine with own weight and turn into elemloads
%DLoads=[EltID n1globalX n1globalY n1globalZ …]
DLoadsTraffic = zeros([160 7]);
DLoadsTraffic =[DLoadsTraffic;
161 0 0 -2.5 0 0 -2.5;
162 0 0 -2.5 0 0 -2.5;
163 0 0 -2.5 0 0 -2.5;
164 0 0 -9 0 0 -9];
DLoadsTraffic=reprow(DLoadsTraffic,161:164,32,[4 0 0 0 0 0 0]);
DLoadsTraffic = [DLoadsTraffic;
zeros([355-292 7])];
DLoadsOwn(161:292,1)=0;
DLoads = DLoadsOwn + DLoadsTraffic
P=elemloads(DLoads,Nodes,Elements,Types,DOF);
%P = P1+P2
%STEP L3: compute the point loads working on the deck
% Nodal loads: 5 kN horizontally on node 4.
%seldof=[4.01];
%PLoad= [5];
% Assembly of the load vectors:
%P=nodalvalues(DOF,seldof,PLoad);
%PLoad=[0]
%seldof=[]
%Pp=nodalvalues(DOF,seldof,PLoad);
%STEP L4: Compute P
%STEP U1: Compute the displacements
U=KP;
% Plot displacements
figure
plotdisp(Nodes,Elements,Types,DOF,U,’DispScal’,5)
title(‘Displacements’)
exportgraphics(gcf,’FE_DYB_Displ.jpg’,’Resolution’,300)
printdisp(Nodes,DOF,U);
% Compute element forces
Forces=elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U,DLoads);
printforc(Elements,Forces);
% Plot element forces
figure
plotforc(‘norm’,Nodes,Elements,Types,Forces,DLoads)
title(‘Normal forces’)
exportgraphics(gcf,’FE_DYB_Normal.jpg’,’Resolution’,300)
figure
plotforc(‘sheary’,Nodes,Elements,Types,Forces,DLoads)
title(‘Shear forces’)
exportgraphics(gcf,’FE_DYB_Shear.jpg’,’Resolution’,300)
figure
plotforc(‘momz’,Nodes,Elements,Types,Forces,DLoads)
title(‘Bending moments’)
exportgraphics(gcf,’FE_DYB_Moment.jpg’,’Resolution’,300)We are working on a finite elements model using the plug-in STABIL. We get no errors but our displacements are all giving NaN as solution. Can someone explain the meaning of these NaNs? Our K and P matrix look logical so we are out of options…
%initializing
clc
close all
clear all
% units kN,m
unzip(‘stabil-3.1.zip’);
% CALCULATING THE NODES
%STEP N1: the deck
% Nodes=[NodID X Y Z]
Nodes= [1 21 0 4;
2 21+5.25*1 0 4;
3 21+5.25*2 0 4;
4 21+5.25*3 0 4;
5 21+5.25*4 0 4;
6 21+5.25*5 0 4;
7 21+5.25*6 0 4;
8 21+5.25*7 0 4;
9 21+5.25*8 0 4];
Nodes = reprow(Nodes, 1:9, 1, [9 0 4.2 0]);
Nodes = reprow(Nodes, 10:18, 1, [9 0 4.5 0]);
Nodes = reprow(Nodes, 19:27, 1, [9 0 4.2 0]);
Nodes = reprow(Nodes, 28:36, 1, [9 0 3.7 0]);
%STEP N2: the 3 top points of triangle, then copy to second triangle
Nodes= [Nodes;
46 31.5 0 14.5;
47 42 0 25;
48 52.5 0 14.5];
Nodes = reprow(Nodes, 46:48, 1, [3 0 12.9 0]);
%STEP N3: reprow the full triangle + deck two times
%no overlapping nodes
%second triangle + deck
Nodes = reprow(Nodes, 1:8,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 10:17,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 19:26,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 28:35,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 37:44,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 46:51,1,[51 42 0 0]);
%third triangle + deck
Nodes = reprow(Nodes, 1:8,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 10:17,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 19:26,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 28:35,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 37:44,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 46:51,1,[102 84 0 0]);
%STEP N4: add the other notes of the deck parts not directly supported by
%the triangles to the left and right
%The left
Nodes = [Nodes;
154 0 0 4;
155 5.25 0 4;
156 10.5 0 4;
157 15.75 0 4;
158 0 4.2 4;
159 5.25 4.2 4;
160 10.5 4.2 4;
161 15.75 4.2 4;
162 0 8.7 4;
163 5.25 8.7 4;
164 10.5 8.7 4;
165 15.75 8.7 4;
166 0 12.9 4;
167 5.25 12.9 4;
168 10.5 12.9 4;
169 15.75 12.9 4;
170 0 16.6 4;
171 5.25 16.6 4;
172 10.5 16.6 4;
173 15.75 16.6 4];
%The right
Nodes = [Nodes;
174 152.25 0 4;
175 157.5 0 4;
176 162.75 0 4;
177 168 0 4;
178 152.25 4.2 4;
179 157.5 4.2 4;
180 162.75 4.2 4;
181 168 4.2 4;
182 152.25 8.7 4;
183 157.5 8.7 4;
184 162.75 8.7 4;
185 168 8.7 4;
186 152.25 12.9 4;
187 157.5 12.9 4;
188 162.75 12.9 4;
189 168 12.9 4;
190 152.25 16.6 4;
191 157.5 16.6 4;
192 162.75 16.6 4;
193 168 16.6 4];
%STEP N5: reference nodes
%Longitudinal beams and diagonal beams triangle 1 to 47 and copies
Nodes= [Nodes;
500 0 0 100;
501 0 4.2 100;
502 0 8.7 100;
503 0 12.9 100;
504 0 16.6 100];
%Beam element from 47 to 9 and all the copies
%also the ones going from 46 to 5 and the copies
Nodes= [Nodes;
505 168 0 100;
506 168 12.9 100];
%Transversal beams
Nodes= [Nodes;
507 31.5 0 100;
508 42 0 100;
509 52.5 0 100;
510 31.5+42 0 100;
511 42+42 0 100;
512 52.5+42 0 100;
513 31.5+42*2 0 100;
514 42+42*2 0 100;
515 52.5+42*2 0 100];
%STEP E1: initializing
% Element types -> {EltTypID EltName}
Types= {1 ‘beam’;
2 ‘truss’};
% Sections=[SecID A ky kz Ixx Iyy Izz] !!!all in m
Sections = [1 0.265*0.990*2 Inf Inf 2*(1/12*0.265*0.990^3) 2*(1/12*0.265*0.990^3) 2*(1/12*0.265*0.990^3); % plywood king post beam
2 0.265*0.540*2 Inf Inf 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3); % plywood diagonal beam
3 0.265*0.540*2 Inf Inf 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3); % horizontal beam
4 ((0.314/2)^2*pi-(0.314/2-0.019)^2*pi) NaN NaN NaN NaN NaN; % steel vertical
5 ((0.314/2)^2*pi-(0.314/2-0.019)^2*pi) NaN NaN NaN NaN NaN; % diagonal steel truss (only A)
6 0.265*0.1350*2 Inf Inf 2*(1/12*0.265*0.1350^3) 2*(1/12*0.265*0.1350^3) 2*(1/12*0.265*0.1350^3); % outer deck beam
7 0.240*0.1350*2 Inf Inf 2*(1/12*0.240*0.1350^3) 2*(1/12*0.240*0.1350^3) 2*(1/12*0.240*0.1350^3); % inner deck beam
8 0.200*0.650 Inf Inf 1/12*0.200*0.650^3 1/12*0.200*0.650^3 1/12*0.200*0.650^3; % outer footbridge
9 0.100*0.100 NaN NaN NaN NaN NaN]; % deck truss steel
% Materials=[MatID E nu rho];
Materials= [1 7.75e6 0.2 2000; % plywood
2 210e6 0.3 7500]; % steel KN/m2
% Elements=[EltID TypID SecID MatID n1 n2 n3]
%STEP E2: deck longitudonal plywood beams
%First triangle
Elements= [1 1 6 1 1 2 500;
2 1 7 1 10 11 501;
3 1 7 1 19 20 502;
4 1 6 1 28 29 503;
5 1 8 1 37 38 504];
Elements=reprow(Elements,1:5,7,[5 0 0 0 1 1 0]);
%Second triangle
Elements=[Elements;
41 1 6 1 9 52 500;
42 1 7 1 18 61 501;
43 1 7 1 27 70 502;
44 1 6 1 36 79 503;
45 1 8 1 45 88 504];
Elements= [Elements;
46 1 6 1 52 53 500;
47 1 7 1 61 62 501;
48 1 7 1 70 71 502;
49 1 6 1 79 80 503;
50 1 8 1 88 89 504];
Elements=reprow(Elements,46:50,6,[5 0 0 0 1 1 0]);
%Third triangle
Elements=[Elements;
81 1 6 1 59 103 500;
82 1 7 1 68 112 501;
83 1 7 1 77 121 502;
84 1 6 1 86 130 503;
85 1 8 1 95 139 504];
Elements= [Elements;
86 1 6 1 103 104 500;
87 1 7 1 112 113 501;
88 1 7 1 121 122 502;
89 1 6 1 130 131 503;
90 1 8 1 139 140 504];
Elements=reprow(Elements,86:90,6,[5 0 0 0 1 1 0]);
%Deck on the right of the triangles
Elements= [Elements;
121 1 6 1 110 174 500;
122 1 7 1 119 178 501;
123 1 7 1 128 182 502;
124 1 6 1 137 186 503;
125 1 8 1 146 190 504];
Elements= [Elements;
126 1 6 1 174 175 500;
127 1 7 1 178 179 501;
128 1 7 1 182 183 502;
129 1 6 1 186 187 503;
130 1 8 1 190 191 504];
Elements=reprow(Elements,126:130,2,[5 0 0 0 1 1 0]);
%Deck on the left of the triangles
Elements= [Elements;
141 1 6 1 154 155 500;
142 1 7 1 158 159 501;
143 1 7 1 162 163 502;
144 1 6 1 166 167 503;
145 1 8 1 170 171 504];
Elements=reprow(Elements,141:145,2,[5 0 0 0 1 1 0]);
Elements= [Elements;
156 1 6 1 157 1 500;
157 1 7 1 161 10 501;
158 1 7 1 165 19 502;
159 1 6 1 169 28 503;
160 1 8 1 173 37 504];
%STEP E2: deck transversal steel trusses
%First triangle
Elements= [Elements;
161 2 9 1 1 10 NaN;
162 2 9 1 10 19 NaN;
163 2 9 1 19 28 NaN;
164 2 9 1 28 37 NaN];
Elements=reprow(Elements,161:164,8,[4 0 0 0 1 1 0]);
%Second triangle
Elements= [Elements;
197 2 9 1 52 61 NaN;
198 2 9 1 61 70 NaN;
199 2 9 1 70 79 NaN;
200 2 9 1 79 88 NaN];
Elements=reprow(Elements,197:200,7,[4 0 0 0 1 1 0]);
%Third triangle
Elements= [Elements;
229 2 9 1 103 112 NaN;
230 2 9 1 112 121 NaN;
231 2 9 1 121 130 NaN;
232 2 9 1 130 139 NaN];
Elements=reprow(Elements,229:232,7,[4 0 0 0 1 1 0]);
%Right of the triangles
Elements= [Elements;
261 2 9 1 174 178 NaN;
262 2 9 1 178 182 NaN;
263 2 9 1 182 186 NaN;
264 2 9 1 186 190 NaN];
Elements=reprow(Elements,261:264,3,[4 0 0 0 1 1 0]);
%Left of the triangles
Elements= [Elements;
277 2 9 1 154 158 NaN;
278 2 9 1 158 162 NaN;
279 2 9 1 162 166 NaN;
280 2 9 1 166 170 NaN];
Elements=reprow(Elements,277:280,3,[4 0 0 0 1 1 0]);
%STEP E3: Beams elements triangles
%Beam element from 1 to 47 and all the copies (king post truss up)
Elements= [Elements;
293 1 1 1 1 47 500;
294 1 1 1 9 98 500;
295 1 1 1 59 149 500;
296 1 1 1 28 50 503;
297 1 1 1 36 101 503;
298 1 1 1 86 152 503];
%Beam element from 47 to 9 and all the copies (king post truss down)
Elements= [Elements;
299 1 1 1 47 9 505;
300 1 1 1 98 59 505;
301 1 1 1 149 110 505;
302 1 1 1 50 36 506;
303 1 1 1 101 86 506;
304 1 1 1 152 137 506];
%Beam element from 5 to 48 and all the copies (side diagonal)
Elements= [Elements;
305 1 2 1 5 48 500;
306 1 2 1 55 99 500;
307 1 2 1 106 150 500;
308 1 2 1 32 51 503;
309 1 2 1 82 102 503;
310 1 2 1 133 153 503];
%Beam element from 46 to 5 and all the copies (side diagonal)
Elements= [Elements;
311 1 2 1 46 5 505;
312 1 2 1 97 55 505;
313 1 2 1 148 106 505;
314 1 2 1 49 32 506;
315 1 2 1 100 82 506;
316 1 2 1 151 133 506];
%Transversal beams between the triangles
Elements= [Elements;
317 1 3 1 46 49 507;
318 1 3 1 47 50 508;
319 1 3 1 48 51 509;
320 1 3 1 97 100 510;
321 1 3 1 98 101 511;
322 1 3 1 99 102 512;
323 1 3 1 148 151 513;
324 1 3 1 149 152 514;
325 1 3 1 150 153 515];
%STEP E4: Truss elements triangles
%Vertical steel trusses
%First triangle
Elements= [Elements;
326 2 4 2 46 3 NaN;
327 2 4 2 47 5 NaN;
328 2 4 2 48 7 NaN;
329 2 4 2 49 30 NaN;
330 2 4 2 50 32 NaN;
331 2 4 2 51 34 NaN];
%Second triangle
Elements= [Elements;
332 2 4 2 97 53 NaN;
333 2 4 2 98 55 NaN;
334 2 4 2 99 57 NaN;
335 2 4 2 100 80 NaN;
336 2 4 2 101 82 NaN;
337 2 4 2 102 84 NaN];
%Third triangle
Elements= [Elements;
338 2 4 2 148 104 NaN;
339 2 4 2 149 106 NaN;
340 2 4 2 150 108 NaN;
341 2 4 2 151 131 NaN;
342 2 4 2 152 133 NaN;
343 2 4 2 153 135 NaN];
%Diagonal trusses
%First triangle
Elements= [Elements;
344 2 5 2 46 50 NaN;
345 2 5 2 47 49 NaN;
346 2 5 2 47 51 NaN;
347 2 5 2 48 50 NaN];
%Second triangle
Elements= [Elements;
348 2 5 2 97 101 NaN;
349 2 5 2 101 99 NaN;
350 2 5 2 98 100 NaN;
351 2 5 2 98 102 NaN];
%Third triangle
Elements= [Elements;
352 2 5 2 148 152 NaN;
353 2 5 2 150 152 NaN;
354 2 5 2 149 151 NaN;
355 2 5 2 149 153 NaN];
%STEP E5: plot
figure
plotelem(Nodes,Elements(find(Elements(:,3)==1),:),Types,’bl’);
hold(‘on’);
plotelem(Nodes,Elements(find(Elements(:,3)==2),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==3),:),Types,’m’);
plotelem(Nodes,Elements(find(Elements(:,3)==4),:),Types,’r’);
plotelem(Nodes,Elements(find(Elements(:,3)==5),:),Types,’r’);
plotelem(Nodes,Elements(find(Elements(:,3)==6),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==7),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==8),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==9),:),Types,’r’);
title(‘Elements’)
figure
plotnodes(Nodes);
title(‘Nodes’)
%STEP D1: generate all possible DOF
DOFall = getdof(Elements,Types);
%STEP D2: remove the DOF limited by the boundary conditions
selnodes = [1;9;28;36;59;86;110;137;154;158;162;166;170;177;181;185;189;193];
dofpattern = [0.01 0.02 0.03]; % the three dof’s
seldof = selnodes + dofpattern; % for the selected nodes, the pattern of dof
seldof = seldof(:);
DOF = removedof(DOFall,seldof);
%STEP K1: assemble the K matrix
K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF);
%STEP K2: check the K matrix
%check symmetry
disp(K-K’)
det(K)
issymmetric(K)
%positive definite
find(eig(K)<0)
%to check wheter invertible
%Now we assign the correct loads to the structure
%STEP L1: compute the self weight of the structure
% Own weight
DLoadsOwn=accel([0 0 9.81e-3],Elements,Types,Sections,Materials);
%STEP L2: compute the distributed load working on the deck
%Combine with own weight and turn into elemloads
%DLoads=[EltID n1globalX n1globalY n1globalZ …]
DLoadsTraffic = zeros([160 7]);
DLoadsTraffic =[DLoadsTraffic;
161 0 0 -2.5 0 0 -2.5;
162 0 0 -2.5 0 0 -2.5;
163 0 0 -2.5 0 0 -2.5;
164 0 0 -9 0 0 -9];
DLoadsTraffic=reprow(DLoadsTraffic,161:164,32,[4 0 0 0 0 0 0]);
DLoadsTraffic = [DLoadsTraffic;
zeros([355-292 7])];
DLoadsOwn(161:292,1)=0;
DLoads = DLoadsOwn + DLoadsTraffic
P=elemloads(DLoads,Nodes,Elements,Types,DOF);
%P = P1+P2
%STEP L3: compute the point loads working on the deck
% Nodal loads: 5 kN horizontally on node 4.
%seldof=[4.01];
%PLoad= [5];
% Assembly of the load vectors:
%P=nodalvalues(DOF,seldof,PLoad);
%PLoad=[0]
%seldof=[]
%Pp=nodalvalues(DOF,seldof,PLoad);
%STEP L4: Compute P
%STEP U1: Compute the displacements
U=KP;
% Plot displacements
figure
plotdisp(Nodes,Elements,Types,DOF,U,’DispScal’,5)
title(‘Displacements’)
exportgraphics(gcf,’FE_DYB_Displ.jpg’,’Resolution’,300)
printdisp(Nodes,DOF,U);
% Compute element forces
Forces=elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U,DLoads);
printforc(Elements,Forces);
% Plot element forces
figure
plotforc(‘norm’,Nodes,Elements,Types,Forces,DLoads)
title(‘Normal forces’)
exportgraphics(gcf,’FE_DYB_Normal.jpg’,’Resolution’,300)
figure
plotforc(‘sheary’,Nodes,Elements,Types,Forces,DLoads)
title(‘Shear forces’)
exportgraphics(gcf,’FE_DYB_Shear.jpg’,’Resolution’,300)
figure
plotforc(‘momz’,Nodes,Elements,Types,Forces,DLoads)
title(‘Bending moments’)
exportgraphics(gcf,’FE_DYB_Moment.jpg’,’Resolution’,300) We are working on a finite elements model using the plug-in STABIL. We get no errors but our displacements are all giving NaN as solution. Can someone explain the meaning of these NaNs? Our K and P matrix look logical so we are out of options…
%initializing
clc
close all
clear all
% units kN,m
unzip(‘stabil-3.1.zip’);
% CALCULATING THE NODES
%STEP N1: the deck
% Nodes=[NodID X Y Z]
Nodes= [1 21 0 4;
2 21+5.25*1 0 4;
3 21+5.25*2 0 4;
4 21+5.25*3 0 4;
5 21+5.25*4 0 4;
6 21+5.25*5 0 4;
7 21+5.25*6 0 4;
8 21+5.25*7 0 4;
9 21+5.25*8 0 4];
Nodes = reprow(Nodes, 1:9, 1, [9 0 4.2 0]);
Nodes = reprow(Nodes, 10:18, 1, [9 0 4.5 0]);
Nodes = reprow(Nodes, 19:27, 1, [9 0 4.2 0]);
Nodes = reprow(Nodes, 28:36, 1, [9 0 3.7 0]);
%STEP N2: the 3 top points of triangle, then copy to second triangle
Nodes= [Nodes;
46 31.5 0 14.5;
47 42 0 25;
48 52.5 0 14.5];
Nodes = reprow(Nodes, 46:48, 1, [3 0 12.9 0]);
%STEP N3: reprow the full triangle + deck two times
%no overlapping nodes
%second triangle + deck
Nodes = reprow(Nodes, 1:8,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 10:17,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 19:26,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 28:35,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 37:44,1,[51 47.25 0 0]);
Nodes = reprow(Nodes, 46:51,1,[51 42 0 0]);
%third triangle + deck
Nodes = reprow(Nodes, 1:8,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 10:17,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 19:26,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 28:35,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 37:44,1,[102 89.25 0 0]);
Nodes = reprow(Nodes, 46:51,1,[102 84 0 0]);
%STEP N4: add the other notes of the deck parts not directly supported by
%the triangles to the left and right
%The left
Nodes = [Nodes;
154 0 0 4;
155 5.25 0 4;
156 10.5 0 4;
157 15.75 0 4;
158 0 4.2 4;
159 5.25 4.2 4;
160 10.5 4.2 4;
161 15.75 4.2 4;
162 0 8.7 4;
163 5.25 8.7 4;
164 10.5 8.7 4;
165 15.75 8.7 4;
166 0 12.9 4;
167 5.25 12.9 4;
168 10.5 12.9 4;
169 15.75 12.9 4;
170 0 16.6 4;
171 5.25 16.6 4;
172 10.5 16.6 4;
173 15.75 16.6 4];
%The right
Nodes = [Nodes;
174 152.25 0 4;
175 157.5 0 4;
176 162.75 0 4;
177 168 0 4;
178 152.25 4.2 4;
179 157.5 4.2 4;
180 162.75 4.2 4;
181 168 4.2 4;
182 152.25 8.7 4;
183 157.5 8.7 4;
184 162.75 8.7 4;
185 168 8.7 4;
186 152.25 12.9 4;
187 157.5 12.9 4;
188 162.75 12.9 4;
189 168 12.9 4;
190 152.25 16.6 4;
191 157.5 16.6 4;
192 162.75 16.6 4;
193 168 16.6 4];
%STEP N5: reference nodes
%Longitudinal beams and diagonal beams triangle 1 to 47 and copies
Nodes= [Nodes;
500 0 0 100;
501 0 4.2 100;
502 0 8.7 100;
503 0 12.9 100;
504 0 16.6 100];
%Beam element from 47 to 9 and all the copies
%also the ones going from 46 to 5 and the copies
Nodes= [Nodes;
505 168 0 100;
506 168 12.9 100];
%Transversal beams
Nodes= [Nodes;
507 31.5 0 100;
508 42 0 100;
509 52.5 0 100;
510 31.5+42 0 100;
511 42+42 0 100;
512 52.5+42 0 100;
513 31.5+42*2 0 100;
514 42+42*2 0 100;
515 52.5+42*2 0 100];
%STEP E1: initializing
% Element types -> {EltTypID EltName}
Types= {1 ‘beam’;
2 ‘truss’};
% Sections=[SecID A ky kz Ixx Iyy Izz] !!!all in m
Sections = [1 0.265*0.990*2 Inf Inf 2*(1/12*0.265*0.990^3) 2*(1/12*0.265*0.990^3) 2*(1/12*0.265*0.990^3); % plywood king post beam
2 0.265*0.540*2 Inf Inf 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3); % plywood diagonal beam
3 0.265*0.540*2 Inf Inf 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3) 2*(1/12*0.265*0.540^3); % horizontal beam
4 ((0.314/2)^2*pi-(0.314/2-0.019)^2*pi) NaN NaN NaN NaN NaN; % steel vertical
5 ((0.314/2)^2*pi-(0.314/2-0.019)^2*pi) NaN NaN NaN NaN NaN; % diagonal steel truss (only A)
6 0.265*0.1350*2 Inf Inf 2*(1/12*0.265*0.1350^3) 2*(1/12*0.265*0.1350^3) 2*(1/12*0.265*0.1350^3); % outer deck beam
7 0.240*0.1350*2 Inf Inf 2*(1/12*0.240*0.1350^3) 2*(1/12*0.240*0.1350^3) 2*(1/12*0.240*0.1350^3); % inner deck beam
8 0.200*0.650 Inf Inf 1/12*0.200*0.650^3 1/12*0.200*0.650^3 1/12*0.200*0.650^3; % outer footbridge
9 0.100*0.100 NaN NaN NaN NaN NaN]; % deck truss steel
% Materials=[MatID E nu rho];
Materials= [1 7.75e6 0.2 2000; % plywood
2 210e6 0.3 7500]; % steel KN/m2
% Elements=[EltID TypID SecID MatID n1 n2 n3]
%STEP E2: deck longitudonal plywood beams
%First triangle
Elements= [1 1 6 1 1 2 500;
2 1 7 1 10 11 501;
3 1 7 1 19 20 502;
4 1 6 1 28 29 503;
5 1 8 1 37 38 504];
Elements=reprow(Elements,1:5,7,[5 0 0 0 1 1 0]);
%Second triangle
Elements=[Elements;
41 1 6 1 9 52 500;
42 1 7 1 18 61 501;
43 1 7 1 27 70 502;
44 1 6 1 36 79 503;
45 1 8 1 45 88 504];
Elements= [Elements;
46 1 6 1 52 53 500;
47 1 7 1 61 62 501;
48 1 7 1 70 71 502;
49 1 6 1 79 80 503;
50 1 8 1 88 89 504];
Elements=reprow(Elements,46:50,6,[5 0 0 0 1 1 0]);
%Third triangle
Elements=[Elements;
81 1 6 1 59 103 500;
82 1 7 1 68 112 501;
83 1 7 1 77 121 502;
84 1 6 1 86 130 503;
85 1 8 1 95 139 504];
Elements= [Elements;
86 1 6 1 103 104 500;
87 1 7 1 112 113 501;
88 1 7 1 121 122 502;
89 1 6 1 130 131 503;
90 1 8 1 139 140 504];
Elements=reprow(Elements,86:90,6,[5 0 0 0 1 1 0]);
%Deck on the right of the triangles
Elements= [Elements;
121 1 6 1 110 174 500;
122 1 7 1 119 178 501;
123 1 7 1 128 182 502;
124 1 6 1 137 186 503;
125 1 8 1 146 190 504];
Elements= [Elements;
126 1 6 1 174 175 500;
127 1 7 1 178 179 501;
128 1 7 1 182 183 502;
129 1 6 1 186 187 503;
130 1 8 1 190 191 504];
Elements=reprow(Elements,126:130,2,[5 0 0 0 1 1 0]);
%Deck on the left of the triangles
Elements= [Elements;
141 1 6 1 154 155 500;
142 1 7 1 158 159 501;
143 1 7 1 162 163 502;
144 1 6 1 166 167 503;
145 1 8 1 170 171 504];
Elements=reprow(Elements,141:145,2,[5 0 0 0 1 1 0]);
Elements= [Elements;
156 1 6 1 157 1 500;
157 1 7 1 161 10 501;
158 1 7 1 165 19 502;
159 1 6 1 169 28 503;
160 1 8 1 173 37 504];
%STEP E2: deck transversal steel trusses
%First triangle
Elements= [Elements;
161 2 9 1 1 10 NaN;
162 2 9 1 10 19 NaN;
163 2 9 1 19 28 NaN;
164 2 9 1 28 37 NaN];
Elements=reprow(Elements,161:164,8,[4 0 0 0 1 1 0]);
%Second triangle
Elements= [Elements;
197 2 9 1 52 61 NaN;
198 2 9 1 61 70 NaN;
199 2 9 1 70 79 NaN;
200 2 9 1 79 88 NaN];
Elements=reprow(Elements,197:200,7,[4 0 0 0 1 1 0]);
%Third triangle
Elements= [Elements;
229 2 9 1 103 112 NaN;
230 2 9 1 112 121 NaN;
231 2 9 1 121 130 NaN;
232 2 9 1 130 139 NaN];
Elements=reprow(Elements,229:232,7,[4 0 0 0 1 1 0]);
%Right of the triangles
Elements= [Elements;
261 2 9 1 174 178 NaN;
262 2 9 1 178 182 NaN;
263 2 9 1 182 186 NaN;
264 2 9 1 186 190 NaN];
Elements=reprow(Elements,261:264,3,[4 0 0 0 1 1 0]);
%Left of the triangles
Elements= [Elements;
277 2 9 1 154 158 NaN;
278 2 9 1 158 162 NaN;
279 2 9 1 162 166 NaN;
280 2 9 1 166 170 NaN];
Elements=reprow(Elements,277:280,3,[4 0 0 0 1 1 0]);
%STEP E3: Beams elements triangles
%Beam element from 1 to 47 and all the copies (king post truss up)
Elements= [Elements;
293 1 1 1 1 47 500;
294 1 1 1 9 98 500;
295 1 1 1 59 149 500;
296 1 1 1 28 50 503;
297 1 1 1 36 101 503;
298 1 1 1 86 152 503];
%Beam element from 47 to 9 and all the copies (king post truss down)
Elements= [Elements;
299 1 1 1 47 9 505;
300 1 1 1 98 59 505;
301 1 1 1 149 110 505;
302 1 1 1 50 36 506;
303 1 1 1 101 86 506;
304 1 1 1 152 137 506];
%Beam element from 5 to 48 and all the copies (side diagonal)
Elements= [Elements;
305 1 2 1 5 48 500;
306 1 2 1 55 99 500;
307 1 2 1 106 150 500;
308 1 2 1 32 51 503;
309 1 2 1 82 102 503;
310 1 2 1 133 153 503];
%Beam element from 46 to 5 and all the copies (side diagonal)
Elements= [Elements;
311 1 2 1 46 5 505;
312 1 2 1 97 55 505;
313 1 2 1 148 106 505;
314 1 2 1 49 32 506;
315 1 2 1 100 82 506;
316 1 2 1 151 133 506];
%Transversal beams between the triangles
Elements= [Elements;
317 1 3 1 46 49 507;
318 1 3 1 47 50 508;
319 1 3 1 48 51 509;
320 1 3 1 97 100 510;
321 1 3 1 98 101 511;
322 1 3 1 99 102 512;
323 1 3 1 148 151 513;
324 1 3 1 149 152 514;
325 1 3 1 150 153 515];
%STEP E4: Truss elements triangles
%Vertical steel trusses
%First triangle
Elements= [Elements;
326 2 4 2 46 3 NaN;
327 2 4 2 47 5 NaN;
328 2 4 2 48 7 NaN;
329 2 4 2 49 30 NaN;
330 2 4 2 50 32 NaN;
331 2 4 2 51 34 NaN];
%Second triangle
Elements= [Elements;
332 2 4 2 97 53 NaN;
333 2 4 2 98 55 NaN;
334 2 4 2 99 57 NaN;
335 2 4 2 100 80 NaN;
336 2 4 2 101 82 NaN;
337 2 4 2 102 84 NaN];
%Third triangle
Elements= [Elements;
338 2 4 2 148 104 NaN;
339 2 4 2 149 106 NaN;
340 2 4 2 150 108 NaN;
341 2 4 2 151 131 NaN;
342 2 4 2 152 133 NaN;
343 2 4 2 153 135 NaN];
%Diagonal trusses
%First triangle
Elements= [Elements;
344 2 5 2 46 50 NaN;
345 2 5 2 47 49 NaN;
346 2 5 2 47 51 NaN;
347 2 5 2 48 50 NaN];
%Second triangle
Elements= [Elements;
348 2 5 2 97 101 NaN;
349 2 5 2 101 99 NaN;
350 2 5 2 98 100 NaN;
351 2 5 2 98 102 NaN];
%Third triangle
Elements= [Elements;
352 2 5 2 148 152 NaN;
353 2 5 2 150 152 NaN;
354 2 5 2 149 151 NaN;
355 2 5 2 149 153 NaN];
%STEP E5: plot
figure
plotelem(Nodes,Elements(find(Elements(:,3)==1),:),Types,’bl’);
hold(‘on’);
plotelem(Nodes,Elements(find(Elements(:,3)==2),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==3),:),Types,’m’);
plotelem(Nodes,Elements(find(Elements(:,3)==4),:),Types,’r’);
plotelem(Nodes,Elements(find(Elements(:,3)==5),:),Types,’r’);
plotelem(Nodes,Elements(find(Elements(:,3)==6),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==7),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==8),:),Types,’bl’);
plotelem(Nodes,Elements(find(Elements(:,3)==9),:),Types,’r’);
title(‘Elements’)
figure
plotnodes(Nodes);
title(‘Nodes’)
%STEP D1: generate all possible DOF
DOFall = getdof(Elements,Types);
%STEP D2: remove the DOF limited by the boundary conditions
selnodes = [1;9;28;36;59;86;110;137;154;158;162;166;170;177;181;185;189;193];
dofpattern = [0.01 0.02 0.03]; % the three dof’s
seldof = selnodes + dofpattern; % for the selected nodes, the pattern of dof
seldof = seldof(:);
DOF = removedof(DOFall,seldof);
%STEP K1: assemble the K matrix
K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF);
%STEP K2: check the K matrix
%check symmetry
disp(K-K’)
det(K)
issymmetric(K)
%positive definite
find(eig(K)<0)
%to check wheter invertible
%Now we assign the correct loads to the structure
%STEP L1: compute the self weight of the structure
% Own weight
DLoadsOwn=accel([0 0 9.81e-3],Elements,Types,Sections,Materials);
%STEP L2: compute the distributed load working on the deck
%Combine with own weight and turn into elemloads
%DLoads=[EltID n1globalX n1globalY n1globalZ …]
DLoadsTraffic = zeros([160 7]);
DLoadsTraffic =[DLoadsTraffic;
161 0 0 -2.5 0 0 -2.5;
162 0 0 -2.5 0 0 -2.5;
163 0 0 -2.5 0 0 -2.5;
164 0 0 -9 0 0 -9];
DLoadsTraffic=reprow(DLoadsTraffic,161:164,32,[4 0 0 0 0 0 0]);
DLoadsTraffic = [DLoadsTraffic;
zeros([355-292 7])];
DLoadsOwn(161:292,1)=0;
DLoads = DLoadsOwn + DLoadsTraffic
P=elemloads(DLoads,Nodes,Elements,Types,DOF);
%P = P1+P2
%STEP L3: compute the point loads working on the deck
% Nodal loads: 5 kN horizontally on node 4.
%seldof=[4.01];
%PLoad= [5];
% Assembly of the load vectors:
%P=nodalvalues(DOF,seldof,PLoad);
%PLoad=[0]
%seldof=[]
%Pp=nodalvalues(DOF,seldof,PLoad);
%STEP L4: Compute P
%STEP U1: Compute the displacements
U=KP;
% Plot displacements
figure
plotdisp(Nodes,Elements,Types,DOF,U,’DispScal’,5)
title(‘Displacements’)
exportgraphics(gcf,’FE_DYB_Displ.jpg’,’Resolution’,300)
printdisp(Nodes,DOF,U);
% Compute element forces
Forces=elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U,DLoads);
printforc(Elements,Forces);
% Plot element forces
figure
plotforc(‘norm’,Nodes,Elements,Types,Forces,DLoads)
title(‘Normal forces’)
exportgraphics(gcf,’FE_DYB_Normal.jpg’,’Resolution’,300)
figure
plotforc(‘sheary’,Nodes,Elements,Types,Forces,DLoads)
title(‘Shear forces’)
exportgraphics(gcf,’FE_DYB_Shear.jpg’,’Resolution’,300)
figure
plotforc(‘momz’,Nodes,Elements,Types,Forces,DLoads)
title(‘Bending moments’)
exportgraphics(gcf,’FE_DYB_Moment.jpg’,’Resolution’,300) stabil, finite elements, structures, nan MATLAB Answers — New Questions
