Info

此问题已关闭。 请重新打开它进行编辑或回答。

Why is this happening?

3 次查看（过去 30 天）
Scott Banks
Scott Banks 2023-11-18
关闭： Walter Roberson 2023-11-18
Hi Guys,
I firstly want to apologise, because I am going to post a lot of code, but hopefully you won't need to go through all of it.
I am a post grad Civil Engineer and I am using a method called The Direct Stiffness Method to solve for a portal gable frame. I don't know if it will help, but here is what it looks like:
The thing that makes this tricky to solve is that of the non-prismatic eaves haunch (seen in the pink). To model this I have had to divide the section into a number slices that in reality would give a staggered apperance. I have divided it into six equal segements, each with a smaller cross-section in the directon towards the centre of the frame.
Anyway, back to the programming side of things, I am getting some strange results when I do my computations. Please see the following code, but in particular go to where I have lablled 'Relevant Code'.
clear, clc, close all
% Haunch length
L = 3.0165;
% Length of far end of haunch and sharp end
h1 = 754/1000;
h5 = 453.4/1000;
% Haunch angle at column intersection
theta = atand(L/(h1-h5));
% Cut haunch into six segments
d = L/6;
% increment that the haunch decreases
hs = d/tand(theta);
% for loop to calculate height of sections
for i = 1:6
H(i) = h1 - hs;
h1 = H(i);
endfor
% properties of rafter sections
b = 189.8/1000;
tf = 12.7/1000;
tw = 8.5/1000;
% calculating of second moment of area
for i = 1:6
Iy(i) = 1/12.*b.*H(i).^3 - 1/12.*(b-tw).*(H(i) - 2*tf).^3;
endfor
%calculating cross-sectional area
for i = 1:6
Ac(i) = H(i).*b -(b - tw)*(H(i) - 2*tf);
endfor
Iy_start = 9.49961E-04;
Ac_start = 0.01112;
Iy = [Iy_start Iy(1) Iy(2) Iy(3) Iy(4) Iy(5) Iy(6)];
Ach = [Ac_start Ac(1) Ac(2) Ac(3) Ac(4) Ac(5) Ac(6)];
for i = 1:6
Ix(i) = (Iy(i) + Iy(i+1))/2
Ach(i) = (Ach(i) + Ach(i+1))/2
end
E = 2.1E+08;
Ic = 8.74999E-04;
Acr = 0.01442;
Lc = 8;
n = 51;
alpha = 6;
matrix = zeros(6,n);
Lr = 15/(cosd(alpha)) - 3.0165;
Ir = 2.94331E-04;
Ar = 0.00856;
% Local stiffness matrix for columns
Kc = [1 0 0 -1 0 0
0 12 6*Lc 0 -12 6*Lc
0 6*Lc 4*Lc^2 0 -6*Lc 2*Lc^2
-1 0 0 1 0 0
0 -12 -6*Lc 0 12 -6*Lc
0 6*Lc 2*Lc^2 0 -6*Lc 4*Lc^2];
% Local stiffness matrix for rafters
Kr = [1 0 0 -1 0 0
0 12 6*Lr 0 -12 6*Lr
0 6*Lr 4*Lr^2 0 -6*Lr 2*Lr^2
-1 0 0 1 0 0
0 -12 -6*Lr 0 12 -6*Lr
0 6*Lr 2*Lr^2 0 -6*Lr 4*Lr^2];
% Local stiffness matrix for eaves haunch
Kh = [1 0 0 -1 0 0
0 12 6*d 0 -12 6*d
0 6*d 4*d^2 0 -6*d 2*d^2
-1 0 0 1 0 0
0 -12 -6*d 0 12 -6*d
0 6*d 2*d^2 0 -6*d 4*d^2];
for i = 1:6
if i == 1 || i == 4
Kc(i,:) = Kc(i,:)*E*Acr/Lc;
Kr(i,:) = Kr(i,:)*E*Ar/Lr;
Kh1(i,:) = Kh(i,:)*E*Ach(1)/d;
Kh2(i,:) = Kh(i,:)*E*Ach(2)/d;
Kh3(i,:) = Kh(i,:)*E*Ach(3)/d;
Kh4(i,:) = Kh(i,:)*E*Ach(4)/d;
Kh5(i,:) = Kh(i,:)*E*Ach(5)/d;
Kh6(i,:) = Kh(i,:)*E*Ach(6)/d;
Kh7(i,:) = Kh(i,:)*E*Ach(6)/Lr;
else
Kc(i,:) = Kc(i,:)*E*Ic/Lc^3;
Kr(i,:) = Kr(i,:)*E*Ir/Lr^3;
Kh1(i,:) = Kh(i,:)*E*Ix(1)/d^3;
Kh2(i,:) = Kh(i,:)*E*Ix(2)/d^3;
Kh3(i,:) = Kh(i,:)*E*Ix(3)/d^3;
Kh4(i,:) = Kh(i,:)*E*Ix(4)/d^3;
Kh5(i,:) = Kh(i,:)*E*Ix(5)/d^3;
Kh6(i,:) = Kh(i,:)*E*Iy(6)/d^3;
Kh7(i,:) = Kh(i,:)*E*Iy(6)/Lr^3;
end
end
Kh = cat(6, Kh1, Kh2, Kh3, Kh4, Kh5, Kh6)
% Transformation matries
% Left-hand column
T1 = matrix;
T1(2,1) = -1;
T1(1,2) = 1;
T1(3,3) = 1;
T1(5,4) = -1;
T1(4,5) = 1;
T1(6,6) = 1;
% Start of haunch left side
T2 = matrix;
T2(1,4) = cosd(alpha);
T2(1,5) = sind(alpha);
T2(2,4) = -sind(alpha);
T2(2,5) = cosd(alpha);
T2(3,6) = 1;
T2(4,7) = cosd(alpha);
T2(4,8) = sind(alpha);
T2(5,7) = -sind(alpha);
T2(5,8) = cosd(alpha);
T2(6,9) = 1;
% Remaing Sections of eaves haunch
for i = 1:5
ThL(:,:,i) = circshift(T2,[0,i*3])
endfor
% Rest of Left Rafter
TrL = matrix;
TrL(1,22) = cosd(alpha);
TrL(1,23) = sind(alpha);
TrL(2,22) = -sind(alpha);
TrL(2,23) = cosd(alpha);
TrL(3,24) = 1;
TrL(4,25) = cosd(alpha);
TrL(4,26) = sind(alpha);
TrL(5,25) = -sind(alpha);
TrL(5,26) = cosd(alpha);
TRL(6,27) = 1;
% Start of Right Rafter
TrR = matrix;
TrR(1,25) = cosd(-alpha);
TrR(1,26) = sind(-alpha);
TrR(2,25) = -sind(-alpha);
TrR(2,26) = cosd(-alpha);
TrR(3,27) = 1;
TrR(4,28) = cosd(-alpha);
TrR(4,29) = sind(-alpha);
TrR(5,28) = -sind(-alpha);
TrR(5,29) = cosd(-alpha);
TrR(6,30) = 1;
% Start of Eaves Haunch
Th3 = matrix;
Th3(1,28) = cosd(-alpha);
Th3(1,29) = sind(-alpha);
Th3(2,28) = -sind(-alpha);
Th3(2,29) = cosd(-alpha);
Th3(3,30) = 1;
Th3(4,31) = cosd(-alpha);
Th3(4,32) = sind(-alpha);
Th3(5,31) = -sind(-alpha);
Th3(5,32) = cosd(-alpha);
Th3(6,33) = 1;
% Remaining sections of eaves haunch
for i = 1:5
ThR(:,:,i) = circshift(Th3,[0,i*3])
endfor
% Right-hand column
T4 = matrix;
T4(2,49) = -1;
T4(1,50) = 1;
T4(3,51) = 1;
T4(5,46) = -1;
T4(4,47) = 1;
T4(6,48) = 1;
% Transformation from local to global
Km1 = T1'*Kc*T1;
Km2 = T2'*Kh(1)*T2;
Km3 = TrL'*Kh(7)*TrL;
Km4 = TrR'*Kh(7)*TrR;
Km5 = Th3'*Kh(6)*Th3;
Km6 = T4'*Kc*T4;
(Relevant Code) for i = 1:5
KmhL(:,:,i) = ThL(:,:,i)'*Kh(i+1)*ThL(:,:,i);
KmhR(:,:,i) = ThR(:,:,i)'*Kh(6-i)*ThR(:,:,i);
endfor
% Assembly
Ks = KmhL(:,:,1) + KmhL(:,:,2) + KmhL(:,:,3)+ KmhL(:,:,4)+ KmhL(:,:,5)+ KmhR(:,:,1)+ KmhR(:,:,2)+KmhR(:,:,3)+KmhR(:,:,4)+KmhR(:,:,5) + Km1 + Km2 + Km3 + Km4 + Km5 + Km6
Ks([1,2,49,50],:)=[];
Ks(:,[1,2,49,50])= [];
F = zeros(n,1);
F(4) = 0.484;
F(46) = 0.565;
F([1,2,49,50]) = [];
% Solving for displacements
U = inv(Ks)*F
Ux = [0 0 U(1:46)' 0 0 U(47)]'
For some reason, the KmhL matrix is producing a matrix with all values as zero - like this:
There should be some numbers in this matrix!!
So I really hope you guys don't have to go through the whol;e code, but if you could help, that would be great.
Many Thanks,
Scott
  1 个评论
Walter Roberson
Walter Roberson 2023-11-18
This is not a matlab question. matlab does not use endfor.
The code appears to be Octave. The purpose of Octave is to force Mathworks out of business (or at the very least to release literally all Mathworks proprietary code as Open Source.)

此问题已关闭。

回答（0 个）

此问题已关闭。

标签

Community Treasure Hunt

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

Start Hunting!

Translated by