MATLAB Answers

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

조회 수: 4(최근 30일)
Shimon Katzman
Shimon Katzman 4 Dec 2019
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

Shimon Katzman
Shimon Katzman 4 Dec 2019
Can anyone help please please?
Star Strider
Star Strider 4 Dec 2019
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.

로그인 to comment.

답변 수 (2)

Star Strider
Star Strider 4 Dec 2019
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

표시 이전 댓글 수: 3
Star Strider
Star Strider 5 Dec 2019
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.
Shimon Katzman
Shimon Katzman 5 Dec 2019
Thank you very much!!!!

로그인 to comment.


Walter Roberson
Walter Roberson 5 Dec 2019
Repaired code is attached. Some optimizations were made as well.

  댓글 수: 1

Shimon Katzman
Shimon Katzman 5 Dec 2019
Wow Thank you very much.
Thank you.

로그인 to comment.

이 질문에 답변하려면 로그인을(를) 수행하십시오.


Translated by