Infinite loop error in moment-curvature diagram program.
Show older comments
Hello,
I have a task with submission in today. The source code is added. Program first take a concrete strain (EpsCtop) in for loop. For the while loop, initial variables are taken which are cmin,cmax,c,Hsb and NumSlice. I aspect from code in while loop, after taking the initial c value, do the calculation and check the absolute value of error is less than error limit or not, if it is not, define another c with interpolation method and do the same until abs(Error)<ErrorLimit. For the easier understanding, a simple sketch is added. 
My problem is that c is always equal to the 2500 and loop does not end. Probably, it is a simple mistake in structure but I cannot see.
Thank u all for your answers, good healthy days.
%--------Moment_Curvature_Relationship_Program--------
% Basic_Assumptions
% 1_Plane sections remain plane.
% 2_Perfect bond btw. concrete and steel.
% 3_Material constituve models==>Hognasted curve for concrete and
% elastic-perfectly plastic model for concrete
% 4_Effect of confinement on core concrete, strain hardening in steel
% and tensile strength of concrete are neglected.
% Note: This program runs with error less than 1.
%---------------------------------------------------------------------
%Inputs
% Load(N)__N
% Width of cross-section(b)__mm
% Height of cross-section(h)__mm
% Characteristic strength of concrete(fck)__MPa
% Characteristic strength of steel(fyk)__MPa
% Total number of steel layers in cross section(NumSteel)
% Diameter of steel bars in ith layer(d)__mm
% Number of steel bars in ith layer(l)
% Steel layer distances from top of cross-section(y(i))__mm
%---------------------------------------------------------------------
%Some Parameters
%Es:______Modulus of elasticity for steel__MPa
%EpsCo:___Concrete strain corresponding the maximum stress
%EpsCu:___Concrete ultimate strain
%EpsSy:___Yield strain of steel
%EpsSu:___Ultimate strain of steel
%EpsCtop:_Extreme fiber concrete strain
%c=_______Neutral axis depth__mm
%i=_______Steel layer number from top to bottom
%As(i)=___Area of steel bar__mm^2
%EpsS(i)=_Strain of ith steel layer
%SgmS(i)=_Stress of ith steel layer__MPa
%Fs(i)=___Force of ith steel layer__kN
%TFs=_____Total force of steel__kN
%Hsb=_____Height of concrete stress block__mm
%j=_______Number of slices
%NumSlice=Total slice number
%Y(j)=____Distance from top to concrete slice's midpoint__mm
%EpsC(j)=_Strain of jth concrete slice
%SgmC(j)=_Stress of jth steel layer__MPa
%Fc(j)=___Force of concrete slice__kN
%TFc=_____Total concrete block force__kN
%Ms(i)=___Moment of steel layer i about centroid__kN/m^2
%Mc(j)=___Moment of concrete slice j about centroid__kN/m^2
%TMc=_____Total concrete moment__kN/m^2
%TMs=_____Total steel moment__kN/m^2
%M=_______Moment__kN/m^2
%K=_______Curvature__rad/m
%--------------------------------------------------------------------
clear
clc
%Entering Inputs---------------------------------------------------
N=input('Enter load on section:');
b=input('Enter cross-section width:');
h=input('Enter the cross-section height:');
fck=input('Enter the characteristic strength of concrete:');
fyk=input('Enter the characteristic strength of steel:');
NumSteel=input('Enter total number of steel layers:');
d=zeros(NumSteel,1);
l=zeros(NumSteel,1);
As=zeros(NumSteel,1);
y=zeros(NumSteel,1);
for i=1:NumSteel
i
d(i,1)=input('Enter steel diameter of layer i:');
l(i,1)=input('Enter how many steel bar in layer i:');
y(i,1)=input('Please enter the y(i):');
As(i,1)=(pi)*l(i,1)*d(i,1)^2/4;
end
%Constants-------------------------------------------------
EpsCo=0.002;
EpsCu=0.0038;
Es=200000;
EpsSu=0.1;
%Algorithm-------------------------------------------------
EpsSy=fyk/Es;
ErrorLimit=5;
EpsS=zeros(NumSteel,1);
SgmS=zeros(NumSteel,1);
Fs=zeros(NumSteel,1);
Ms=zeros(NumSteel,1);
cmin=0;
cmax=10.*h
c=0.5*(cmin+cmax)
Hsb=min(c,h);
NumSlice=Hsb/10;
Y=zeros(NumSlice,1);
EpsC=zeros(NumSlice,1);
SgmC=zeros(NumSlice,1);
Fc=zeros(NumSlice,1);
Mc=zeros(NumSlice,1);
M=zeros(99,1);
K=zeros(99,1);
for EpsCtop=0.0002:0.0002:0.02
Error=7;
while abs(Error)>ErrorLimit
for i=1:NumSteel
EpsS(i,1)=EpsCtop*(c-y(i,1))/c;
x=EpsS(i,1);
if abs(x)<=EpsSy
SgmS(i,1)=EpsS(i,1)*Es;
else if (EpsSy<abs(x))&&(abs(x)<=EpsSu)
SgmS(i,1)=fyk*EpsS(i,1)/(abs(EpsS(i,1)));
else
SgmS(i,1)=0;
end
end
Fs(i,1)=SgmS(i,1)*As(i,1)/10^3;
end
TFs=sum(Fs);
for j=1:NumSlice
Y(j,1)=Hsb*(j-0.5)/NumSlice;
EpsC(j,1)=EpsCtop*(c-Y(j,1))/c;
if (0<EpsC(j,1))&&(EpsC(j,1)<=EpsCo)
SgmC(j,1)=0.85*fck*(2*(EpsC(j,1)/EpsCo)-(EpsC(j,1)/EpsCo)^2);
else if (EpsCo<EpsC(j,1))&&(EpsC(j,1))<=EpsCu
SgmC(j,1)=0.85*fck-(0.15*(0.85*fck))*(EpsC(j,1)-EpsCo)/(EpsCu-EpsCo);
else
SgmC(j,1)=0;
end
end
Fc(j,1)=(SgmC(j,1)*b*Hsb/NumSlice)/10^3;
end
TFc=sum(Fc);
if Error<0
cmax=c
c=0.5*(cmin+c)
else
cmin=c
c=0.5*(c+cmax)
end
Error=N-TFc-TFs
Hsb=min(c,h);
NumSlice=Hsb/10;
end
for i=1:NumSteel
Ms(i,1)=Fs(i,1)*(h/2-y(i,1))/10^3;
end
TMs=sum(Ms);
for j=1:NumSlice
Mc(j,1)=Fc(j,1)*(h/2-Y(j,1))/10^3;
end
TMc=sum(Mc);
for i=1:99
M(i,1)=TMc+TMs;
K(i,1)=10^3*EpsCtop/c;
end
end
%Outputs----------------------------------------------
M
K
%Plotting---------------------------------------------
plot(M,K)
xlabel('K(rad)')
ylabel('M(kN/m^2)')
title('Moment vs. Curvature Diagram')
ln.LineWidth = 3;
ln.Color = [47 79 79];
Answers (1)
I have no idea what your code is intented to do. The missing comments are a serious problem. That all inputs must by typed in instead of providing them as input arguments to a function, makes the debugging very hard.
But one guess:
if Error<0
cmax=c
c=0.5*(cmin+c)
else
cmin=c
c=0.5*(c+cmax)
end
Error=N-TFc-TFs % Are you sure, that you evaluate Error at first and
% want to calculate it afterwards?
% What about moving this line before the if() block?
If this is not the problem, set a breakpoint in the code and use the debugger to step through your function line by line.
A hint: 10^3 is a very expensive power operation, while 1e3 is a cheap constant.
3 Comments
IHSAN GENCTURK
on 8 Apr 2021
Edited: IHSAN GENCTURK
on 8 Apr 2021
Jan
on 8 Apr 2021
"Code calculate moment and curvature values for different strain levels of a reinforced concrete cross-section, related to the civil engineering, to the understand strength for flexural buckling of section."
:-) Thanks, this is pure poetry for me. I like it. I have no idea what this means. For Matlab this means a bunch of numbers.
I admit, that code blocks like this are very frightening for me:
Y(j,1)=Hsb*(j-0.5)/NumSlice;
EpsC(j,1)=EpsCtop*(c-Y(j,1))/c;
if (0<EpsC(j,1))&&(EpsC(j,1)<=EpsCo)
SgmC(j,1)=0.85*fck*(2*(EpsC(j,1)/EpsCo)-(EpsC(j,1)/EpsCo)^2);
else if (EpsCo<EpsC(j,1))&&(EpsC(j,1))<=EpsCu
SgmC(j,1)=0.85*fck-(0.15*(0.85*fck))*(EpsC(j,1)-EpsCo)/(EpsCu-EpsCo);
else
SgmC(j,1)=0;
end
end
Fc(j,1)=(SgmC(j,1)*b*Hsb/NumSlice)/10^3;
It looks, like someone has rolledn an angry armadillo over the keyboard. I'm not able to debug this.
You could simplify some code, e.g.:
for i=1:NumSteel
Ms(i,1)=Fs(i,1)*(h/2-y(i,1))/10^3;
end
TMs=sum(Ms);
to
TMs = sum(Fs .* (0.5 * h - y)) / 1e3;
and
for i=1:NumSteel
EpsS(i,1)=EpsCtop*(c-y(i,1))/c;
x=EpsS(i,1);
if abs(x)<=EpsSy
SgmS(i,1)=EpsS(i,1)*Es;
else if (EpsSy<abs(x))&&(abs(x)<=EpsSu)
SgmS(i,1)=fyk*EpsS(i,1)/(abs(EpsS(i,1)));
else
SgmS(i,1)=0;
end
end
Fs(i,1)=SgmS(i,1)*As(i,1)/10^3;
end
to
for i = 1:NumSteel
x = EpsCtop * (1 - y(i) / c);
if abs(x) <= EpsSy
SgmS(i) = x * Es;
elseif abs(x) <= EpsSu % EpsSy < abs(x) excluded before already
SgmS(i) = fyk * sign(x); % Or is NaN wanted for x==0? Then: x/abs(x)
else
SgmS(i) = 0;
end
EpsS(i) = x;
Fs(i) = SgmS(i) * As(i) * 1e-3;
end
For me the spaces improve the readability and this help to debug the code.
IHSAN GENCTURK
on 8 Apr 2021
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
on 8 Apr 2021
on 8 Apr 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
