MATLAB Answers

0

Marking a specific point in a graph made out of a loop

Asked by Shimon Katzman on 4 Dec 2019 at 4:59
Latest activity Commented on by Star Strider
on 5 Dec 2019 at 17:13
Hi everyone.
Trying to mark 2 points at the 4 curves in the graph:
1) Marking the peak in each graph (the Mmax point)
2) Marking the point where the value of M belongs to epss=epssh in each graph.
Thank You very much.
b=400; %mm
d=500; %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30 30 90 90]; %Mpa
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) (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);
M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
phi(i,k)=epscm/c(i);
end
end
[Mmax,idx]=max(M) %[kNm]
phiAtMmax=phi(idx) %[1/mm]
epsAtMmax=epscmv(idx)*1000 %[promil]
figure
hold all
for k = 1:numel(fckv)
plot(phi(1:idx+50,k), M(1:idx+50,k))
lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
plot(phi(idx,k),M(idx,k),'r*') % marking the Mmax
end
hold off
grid on
xlabel('phi [1/mm]')
ylabel('Moment [kNm]')
legend(lgdstr)

  2 Comments

Can anyone help please please?
Working on it. Windows 10 gave me a BSOD a few minutes ago for the second time in as many days (I have no kind works for Micro$oft just now), and it took a while to put the pieces of everything back together.

Sign in to comment.

2 Answers

Star Strider 님의 답변 4 Dec 2019 23:02

Try this:
b=400; %mm
d=500; %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30 30 90 90]; %Mpa
tic
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) (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);
M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
phi(i,k)=epscm/c(i);
end
idx(k) = find(diff(M(:,k)) < 0, 1, 'first')%[kNm]
Mmax(k) = M(idx(k),k) %[kNm]
phiAtMmax(k)=phi(idx(k)) %[1/mm]
epsAtMmax(k)=epscmv(idx(k))*1000 %[promil]
end
toc
% [Mmax,idx]=max(M) %[kNm]
% phiAtMmax=phi(idx) %[1/mm]
% epsAtMmax=epscmv(idx)*1000 %[promil]
epss_c = fsolve(@(c) epss(c)-epssh, 100)
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
for k = 1:numel(fckv)
hp(k) = plot(phi(1:idx+50,k), M(1:idx+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')
This correctly plots the maxima of every curve (although the green curve maximum is beyond the green curve plotted limit (30 Mpa, 3000 mm^3), and I corrected the legend problem as well, so that it only refers to the lines, not the markers.
I do not understand about ‘epss_c’ with respect to every curve. It appears to be beyond every curve. In any event, I have no idea how to calculate it with respect to plotting it as a function of ‘phi’.
The Windows 10 crash (that I did not immediately recognise) put me abouit two hours behind.

  6 Comments

As always, my pleasure!
In your code, ‘c’ is a vector, not a matrix. To create a matrrix for it, the easiest way would be to add:
c_mtx(i,k) = c(i);
just before the end of the ‘i’ loop, then change the plot call to:
hp(k) = plot(epscmv(1:idx(k)), c_mtx(1:idx(k),k));
EDIT (5 Dec 2019 at 14:03)
That appears to plot correctly. You will need to change the related calls to show the correct xlabel and ylabel strings, and remove or change the second plot call.
Thank you very much!!!!
My pleasure!

Sign in to comment.


Walter Roberson 님의 답변 5 Dec 2019 8:32

Repaired code is attached. Some optimizations were made as well.

  1 Comment

Wow Thank you very much.
Thank you.

Sign in to comment.