Info
This question is closed. Reopen it to edit or answer.
Tryin to perform a iteration, no sucess
1 view (last 30 days)
Show older comments
Hello Guys,
I'm tryin to perform a iteration but no sucess, The error should feed back a vector called epsilon, can anyone help me?
%Análise não-linear de seção de concreto armado%
%--------------------------------------------------------------------%
% DADOS DE ENTRADA %
%--------------------------------------------------------------------%
Lx = 20; Mxe = 5000; Ne = 12000;
Ly = 20; Mye = 1500;
n = 10; % Discretização da seção transversal
xci= zeros(10,1); ycj = zeros(10,1); acij = zeros(10,10);
for i = 1:n
for j = 1:n
xci(i,1) = (Lx/n)*(i-(1/2)) - Lx/2;
ycj(j,1) = (Ly/n)*(j-(1/2)) - Ly/2;
acij(i,j) = (Lx/n)*(Ly/n);
end
end
xsk = [-7 -7 7 7];
ysk = [-7 7 7 -7];
ask = [6; 6; 6; 6];
%--------------------------------------------------------------------%
% CHUTE INICIAL DA SEÇÃO DEFORMADA %
%--------------------------------------------------------------------%
epsilon0 = -0.00001;
kx = 0;
ky = 0;
erro = 0;
maxTimes = 0;
while max(erro) <= 10e-3 || maxTimes < 100
epsilon = [epsilon0; kx; ky];
for i = 1:n
for j = 1:n
ecij(i,j) = -ky*xci(i,1) + kx*ycj(j,1) + epsilon0;
end
end
esk = -ky*xsk + kx*ysk + epsilon0;
%--------------------------------------------------------------------%
% RELAÇÕES CONSTITUTIVAS %
%--------------------------------------------------------------------%
%Concreto
fc = 200;
fcd = fc/1.4;
Ec = fc/0.002;
if ecij >= 0
sigmacij = 0;
elseif 0 > ecij > -0.002
sigmacij = 0.85*fcd*(1 - (1 - ((ecij/0.002).^2)));
else
sigmacij = -0.85*fcd;
end
if ecij >= 0
dcij = 0;
elseif 0 > ecij > -0.002
dcij = 850*fcd*(1 - ecij);
else
dcij = 0;
end
%Aço
fy = 4200;
Es = 2.05*10^6;
epsilony = fy/Es;
if esk >= epsilony
sigmask = fy;
elseif -epsilony < esk < epsilony
sigmask = Es*esk;
else
sigmask = -fy;
end
if abs(esk) >= epsilony
dsk = 0;
else
dsk = Es;
end
%--------------------------------------------------------------------%
% MATRIZ DE RIGIDEZ TANGENTE %
%--------------------------------------------------------------------%
k11 = (sum(dcij*acij,'all')) + sum((dsk*ask),'all');
k22 = (sum(dcij*acij*(ycj.^2),'all'))+(sum(dsk*ask*(ysk.^2),'all'));
k33 = (sum(dcij*acij*(xci.^2),'all'))+(sum(dsk*ask*(xsk.^2),'all'));
k12 = (sum(dcij*acij*ycj,'all'))+(sum(dsk*ask*ysk,'all'));
k13 = (sum(dcij*acij*xci,'all'))+(sum(dsk*ask*xsk,'all'));
k23 = (sum(dcij*acij*xci.*ycj,'all'))+(sum(dsk*ask*xsk.*ysk,'all'));
k21 = k12;
k31 = k13;
k32 = k23;
K = [k11 k21 k31; k21 k22 k32; k31 k32 k33];
%--------------------------------------------------------------------%
% ESFORÇOS INTERNOS %
%--------------------------------------------------------------------%
N = (sum(sigmacij*acij,'all')+(sum(sigmask*ask,'all')));
Mx = (sum(sigmacij*acij*ycj,'all')+(sum(sigmask*ask*ysk,'all')));
My = (sum(sigmacij*acij*xci,'all')+(sum(sigmask*ask*xsk,'all')));
%--------------------------------------------------------------------%
% VETOR DE FORÇAS DESBALANCEADAS %
%--------------------------------------------------------------------%
F = [Ne;Mxe;Mye] - abs([N;Mx;My]);
%--------------------------------------------------------------------%
% APLICAÇÃO DO MÉTODO DE NEWTON %
%--------------------------------------------------------------------%
delta1 = N/k11 + Mx/k12 + My/k13;
delta2 = N/k21 + Mx/k22 + My/k23;
delta3 = N/k31 + Mx/k32 + My/k33;
delta = [delta1; delta2; delta3];
erro = epsilon - delta;
epsilon = epsilon + erro;
maxTimes = maxTimes + 1;
end
disp(erro)
5 Comments
Answers (1)
Jan
on 9 Apr 2019
Edited: Jan
on 9 Apr 2019
These lines will no do, what you expect:
elseif 0 > ecij > -0.002
Matlab evaluates them from left to right. In the first step it compares 0 > ecij. This is either true or false, which is treated as 1 or 0. Both are greater than -0.002, such that the condition is true in every case.
You want:
elseif 0 > ecij && ecij > -0.002
and the same in some other lines.
By the way, I'd reduce the number of parentheses. Compare
k11 = (sum(dcij*acij,'all')) + sum((dsk*ask),'all');
with
k11 = sum(dcij * acij, 'all') + sum(dsk * ask, 'all');
2 Comments
Jan
on 10 Apr 2019
Edited: Jan
on 10 Apr 2019
@Anax: Oh, I see. ecij is a matrix, but if requires a scalar condition. Then Matlab converts e.g.
if ecij >= 0
automatically to
if all(ecij >= 0, 'all') && ~isempty(ecij)
Most likely you want to apply the operation for each element matching this condition. Such details should be mentioned in the comments. Then replace e.g.:
if ecij >= 0
sigmacij = 0;
elseif 0 > ecij > -0.002
sigmacij = 0.85*fcd*(1 - (1 - ((ecij/0.002).^2)));
else
sigmacij = -0.85*fcd;
end
by
index1 = (ecij >= 0);
sigmacij(index1) = 0;
index2 = (0 > ecij) & (ecij > -0.002);
sigmacij(index2) = 0.85 * fcd * (ecij(index)/0.002) .^ 2;
index3 = ~index1 & ~index2;
sigmacij(index3) = -0.85 * fcd;
I've used the simplification: 1 - (1 - x) == x.
By the way, the block of code to create K looks like a simple matrix equation. Then the matrix notation is much easier and clearer. Compare e.g.
delta1 = N/k11 + Mx/k12 + My/k13;
delta2 = N/k21 + Mx/k22 + My/k23;
delta3 = N/k31 + Mx/k32 + My/k33;
delta = [delta1; delta2; delta3];
with
delta = [N, Mx, My] * (1 ./ K); % Maybe, just as demonstration
I guess, that all the sum(x, 'all') can be replaced by matrix operations such that the code is much easier to read, to understand and to debug.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!