Finding a value and marking it in a graph

1 view (last 30 days)
Hi everyone,
I am trying to find the value of M when epss=epssh=0.009 for every graph (there are 4 of them) and to mark it on the graph unsucssesfully.
Can anyone help me please with this please?
b=400; %mm
d=500; %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30 30 90 90]; %Mpa
num_fckv = numel(fckv);
epscmv = linspace(0.1, 100, 5000)*1E-3;
num_epscvm = numel(epscmv);
fsolve_opts = optimoptions('fsolve','Display','off');
M = zeros(num_epscvm, num_fckv);
phi = zeros(num_epscvm, num_fckv);
idx = zeros(1, num_fckv);
Mmax = zeros(1, num_fckv);
phiAtMmax = zeros(1, num_fckv);
epsAtMmax = zeros(1, num_fckv);
tic
for k = 1:num_fckv
fck = fckv(k);
Ecshah=57000/145*(fck*145)^0.5; %Mpa
Es=200000; %Mpa
Esh=8500; %Mpa
fy=500; %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck*145+0.00195;
kshah=0.025*fck*10^3;
A=Ecshah*eps0/fck;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
As = Asv(k);
c0 = 1000;
c = zeros(1, num_epscvm);
for i=1:num_epscvm
epscm = epscmv(i);
epss=@(c) (d-c)/c*epscm;
funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
compression=@(c) b*fck*c/epscm*integral(funCshah,0,epscm)/1000;
sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
tension=@(c) sigmaSteel(c).*As/1000;
c(i)=fsolve(@(c) compression(c)-tension(c), c0, fsolve_opts);
c0 = c(i);
funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
phi(i,k)=epscm/c(i);
epss1(k,i)= (d-c(i))/c(i)*epscm;
sigmaSteel1(k,i)= Es*epss1(k,i) .* (epss1(k,i)<=epsy) + fy .* (epss1(k,i)>epsy & epss1(k,i)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss1(k,i))./(epssu-epssh)).^(1/P)) .* (epss1(k,i)>epssh & epss1(k,i)<=epssu) + 0 .* (epss1(k,i)>epssu);
c_mtx(i,k) = c(i);
end
idx(k) = find(diff(M(:,k)) < 0, 1, 'first')
[Mmax(k),idx(k)]=max(M(:,k)) %[kNm]
phiAtMmax(k)=phi(idx(k),k) %[1/mm]
epsAtMmax(k)=epscmv(idx(k))*1000 %[promil]
end
toc
% epss_c = fsolve(@(c) epss(c)-epssh, 100, fsolve_opts)
% c_idx = find(c >= epss_c, 1, 'first')
% epscm_exact = interp1(c([-5 5]+c_idx), epscmv([-5 5]+c_idx), epss_c)
figure
hold all
hp = zeros(1, num_fckv);
lgdstr = cell(1, num_fckv);
for k = 1:num_fckv
hp(k) = plot(phi(1:idx(k)+50,k), M(1:idx(k)+50,k));
% plot([1 1]*epss_c, ylim, ':k')
lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
plot(phi(idx(k),k),M(idx(k),k),'r*') % marking the Mmax
end
hold off
grid on
xlabel('phi [1/mm]')
ylabel('Moment [kNm]')
hl = legend(hp,lgdstr);
hl.FontSize = 7;
% text(epss_c, 400, sprintf('\\leftarrow epss_c = %.1f', epss_c), 'HorizontalAlignment','left')

Accepted Answer

Star Strider
Star Strider on 5 Dec 2019
The only way I can see to calculate ‘M’ where ‘epss’ is equal to 0.009 is to duplicate the first nested loop, defining:
epss = @(c) epssh;
If you do that:
for k = 1:numel(fckv)
fck = fckv(k);
Ecshah=57000/145*(fck*145)^0.5; %Mpa
Es=200000; %Mpa
Esh=8500; %Mpa
fy=500; %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck*145+0.00195;
kshah=0.025*fck*10^3;
A=Ecshah*eps0/fck;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
epscmv = linspace(0.1, 100, 5000)*1E-3;
As = Asv(k);
for i=1:numel(epscmv);
epscm = epscmv(i);
epss = @(c) epssh;
% epss=@(c) (d-c)/c*epscm;
funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
compression=@(c) b*fck*c/epscm*integral(funCshah,0,epscm)/1000;
sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
tension=@(c) sigmaSteel(c).*As/1000;
c(i)=fsolve(@(c) compression(c)-tension(c),1000);
funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
M009(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
phi(i,k)=epscm/c(i);
c_mtx(i,k) = c(i);
end
end
toc
figure
mesh(M009)
grid on
and also do the surface plot, you can see the result.
I have no idea what to do with this result, so I leave that to you.
  11 Comments
Shimon Katzman
Shimon Katzman on 8 Dec 2019
I know but I still cant find what's wrong :(
VBBV
VBBV on 25 Dec 2022
b=400; %mm
d=500; %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30 30 90 90]; %Mpa
num_fckv = numel(fckv);
fsolve_opts = optimoptions('fsolve','Display','off');
c0 = 100;
tic
for k = 1:num_fckv
fck = fckv(k);
Ecshah=57000/145*(fck*145)^0.5; %Mpa
Es=200000; %Mpa
Esh=8500; %Mpa
fy=500; %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck*145+0.00195;
kshah=0.025*fck*10^3;
A=Ecshah*eps0/fck;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
As = Asv(k);
%Values at yield
epssAtMy= epsy;
epscmAtMy=@(cAtMy) epssAtMy*cAtMy./(d-cAtMy);
funCshahAtMy=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
compressionAtMy=@(cAtMy) b*fck*cAtMy/epscmAtMy(cAtMy)*integral(funCshahAtMy,0,epscmAtMy(cAtMy))/1000;
%compressionAtMy=@(cAtMy) b*fck*cAtMy/epscmAtMy*integral(funCshahAtMy,0,epscmAtMy)/1000;
sigmaSteelAtMy= fy;
tensionAtMy= sigmaSteelAtMy*As/1000;
cAtMy(k)=fsolve(@(cAtMy) compressionAtMy(cAtMy)-tensionAtMy, c0, fsolve_opts);
c0=cAtMy(k);
funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c0+(c0./epscmAtMy(c0)).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c0+(c0./epscmAtMy(c0)).*epsc) .* (epsc>eps0);
K = integral(funM,0,epscmAtMy(c0));
M(k)=b*fck*cAtMy(k)/epscmAtMy(c0)*K/1000000;
end
M
M = 1×4
1.0e+03 * 0.6333 0.9588 0.6696 1.0822
toc
Elapsed time is 0.099555 seconds.
you need to provide input value to variable to the function epscmAtMy defined using function handle when calling it in the lines
funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c0+(c0./epscmAtMy(c0)).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c0+(c0./epscmAtMy(c0)).*epsc) .* (epsc>eps0);
M(k)=b*fck*cAtMy(k)/epscmAtMy(c0)*K/1000000;

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!