Finding a value and marking it in a graph

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')

Answer 1

[Star Strider]

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.

[Shimon Katzman]

Hi Star, it doesnt works..

The purpose of what i want to do is to add those points to the graph.

I need to know in what values of phi and M , epss=epssh=0.009.

example:

is it possible? Of course to calculate the values too.

Thank you very much.

[Star Strider]

My pleasure.

I gave you a reliable way of calculating ‘M when epss=epssh=0.009’, that calculates the ‘M009’ matrix. You have to take it from there. I have no idea how you want to plot it.

Try this:

plot(phi(1:idx(k)+50,k), M009(1:idx(k)+50,k));

or something similar. The ‘M009’ matrix is what you requested.

I am out of ideas.

[Shimon Katzman]

Hi Star,

I am trying this but it errors me. Can you please help me with the error?

Thank you very much.

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');
cAtMy = zeros(1, num_epscvm);
c0 = 1000;

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*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-cAtMy+(cAtMy./epscmAtMy).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-cAtMy+(cAtMy./epscmAtMy).*epsc) .* (epsc>eps0);
    %M(k)=b*fck*cAtMy/epscmAtMy*integral(funM,0,epscmAtMy)/1000000;
        
     
end 
toc

[Star Strider]

The error it throws when I run it is:

Unrecognized function or variable 'num_epscvm'.

cAtMy = zeros(1, num_epscvm);

I cannot correct that because I have no idea what ‘num_epscvm’ is.

[Shimon Katzman]

You can delete it.. still doesnt work.

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 = 1000;

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*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-cAtMy+(cAtMy./epscmAtMy).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-cAtMy+(cAtMy./epscmAtMy).*epsc) .* (epsc>eps0);
    %M(k)=b*fck*cAtMy/epscmAtMy*integral(funM,0,epscmAtMy)/1000000;
        
     
end 
toc

[Star Strider]

Using a function handle in an expression requires that the argument of the function be included.

Correct this function expression to:

compressionAtMy=@(cAtMy) b*fck*cAtMy/epscmAtMy(cAtMy)*integral(funCshahAtMy,0,epscmAtMy(cAtMy))/1000;

and the code runs without error, and produces:

cAtMy =
     1.2016e+003     1.2016e+003     1.0980e+003     1.0980e+003

[Shimon Katzman]

Hi Star.

1) It gives me back only one value of cAtMy (instead of 4) what is missing?

2) i added the M(k) and again it errors :(((

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-cAtMy+(cAtMy./epscmAtMy).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-cAtMy+(cAtMy./epscmAtMy).*epsc) .* (epsc>eps0);
    M(k)=b*fck*cAtMy(k)/epscmAtMy*integral(funM,0,epscmAtMy)/1000000
          
end 
toc

[Shimon Katzman]

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-cAtMy(k)+(cAtMy(k)./epscmAtMy).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-cAtMy(k)+(cAtMy(k)./epscmAtMy).*epsc) .* (epsc>eps0);
    M(k)=b*fck*cAtMy(k)/epscmAtMy*integral(funM,0,epscmAtMy)/1000000
          
end 
toc

Whatever i do it errors me :

Undefined operator '/' for input arguments of type 'function_handle'.

Error in Untitled (line 39)
    M(k)=b*fck*cAtMy(k)/epscmAtMy*integral(funM,0,epscmAtMy)/1000000

[Star Strider]

That is essentially the same problem that I described the correction for in my previous Comment.

[Shimon Katzman]

I know but I still cant find what's wrong :(

[VBBV]

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;