Tryin to perform a iteration, no sucess

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)

[Jan]

Your comments are very nice, but unfortunately not in English and a little bit lean. Please explain what "no sucess" means with any details.

Matrix operations are easier than loops in Matlab. Compare:

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

with

xci = (Lx / n) * (0.5:n-0.5) - Lx / 2;
ycj = (Ly / n) * (0.5:n-0.5) - Ly / 2;

and

	for i = 1:n
	   for j = 1:n
		   ecij(i,j) = -ky*xci(i,1) + kx*ycj(j,1) + epsilon0;
	   end
	end

with

 ecij = -ky * xci.' + kx * ycj + epsilon0

[Anax Souza]

Ok, wiil try to explain all this script is to calculate the tension in a section of concrete.

This part of the code is my input data:

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]; 

In the other part of the code I have a guess of value of my deformation, they are in this variable

epsilon = [epsilon0; kx; ky]

Where: epsilon0 = Deformation of compression

kx = curvature in x

ky = curvature in y

Then, The code calculates all tensions with the first guess, at the end the code generates 3 deltas:

delta1 = N/k11 + Mx/k12 + My/k13;
delta2 = N/k21 + Mx/k22 + My/k23;
delta3 = N/k31 + Mx/k32 + My/k33;

I want that those deltas make a new epsilon doing this:

erro = epsilon - delta;
epsilon = epsilon + erro;

then the code performs the calculation with the new values unil the error becomes less than 10e-3.

[Jan]

Now mention, what is not working as expected.

[Anax Souza]

The code is not feeding back the variable epsilon with the new value, the script perform the analysis with the first value the stop.

[Jan]

The code stops with an error message. So how can it feed back anything?

Answer 1

[Jan]

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');

[Anax Souza]

I tried to fix as you suggest, but I got this error:

Operands to the || and && operators must be convertible to logical scalar values.

[Jan]

@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.