Asked by Shimon Katzman
on 4 Dec 2019 at 4:59

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)

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.

Star Strider
on 5 Dec 2019 at 13:43

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
5 Dec 2019 16:49

Thank you very much!!!!

Star Strider
5 Dec 2019 17:13

My pleasure!

Sign in to comment.

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

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

Shimon Katzman
on 5 Dec 2019 at 9:06

Wow Thank you very much.

Thank you.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Shimon Katzman (view profile)

## 이 댓글에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/494721-marking-a-specific-point-in-a-graph-made-out-of-a-loop#comment_774527

## Star Strider (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494721-marking-a-specific-point-in-a-graph-made-out-of-a-loop#comment_774536

Sign in to comment.