%% Calculating backstress and the fraction of plastically deforming UFG in harmonic-structure material
% This is an open-source code for free use.
% If you will publish the results of calculations using this code, please
% (i)  acknowledge the use of this code, and
% (ii) cite the following article where respective mathematical model was developed:
% D.Orlov, R.Kulagin, Y.Beygelzimer, Strain partitioning and back-stress evaluation in
% harmonic-structure materials, Materials Letters (2020) 275, 128126; DOI: 10.1016/j.matlet.2020.128126

%% If you have any questions, comments, suggestions or requests, feel free to contact:
% Dmytro Orlov, E-Mail: dmytro.orlov@material.lth.se

%% Please enter INPUT data in the form of coefficients for Sherby equation for strain hardening curves
% s=ss-(ss-sy)*exp(-n*e^p)
% ss - steady-state flow stress
% sy - yield strength
% n, p - material parameters
% ro - initial fraction of CG phase in the HS material
global ro ss1 sy1 n1 p1 ss2 sy2 n2 p2 ss3 sy3 n3 p3  % DO NOT CHANGE ANYTHING IN THIS LINE
de=0.001;
e=(0:de:0.5);
ro=0.6;
ss1=1319; sy1=194; n1=1.022; p1=0.715;  %stress-strain curve of 100% CG material
ss2=1731; sy2=1502; n2=128.6; p2=0.941;%stress-strain curve of 100% UFG material
ss3=1450; sy3=253; n3=1.838; p3=0.857;  %stress-strain curve of HS material
%
% DO NOT CHANGE ANYTHING BELOW THIS LINE
%
%% calculation of coefficients c1 and c2
x=fminsearch(@sfcc,[0.5,0.5]);
c1=round(x(1),3);
c2=round(x(2),3);
e1=e./ro;
e2=e.^2./(ro+e).*((8-3.1416)/2/3.1416/ro);
sig1=ss1-(ss1-sy1).*exp(-n1.*e1.^p1);
sigb=c1.*sig1;
sig2=ss2-(ss2-sy2).*exp(-n2.*e2.^p2);
s_hs=sig1.*(1+c1).*ro.*(1+e)./(ro+e)+...
    sig2.*c2.*(1-ro).*(8-3.1416).*(2*ro+e).*(1+e).*e./((ro+e).*(2*3.1416*ro.*(ro+e)+e.^2.*(8-3.1416)));

%% preparing data for the plot
sig1=ss1-(ss1-sy1).*exp(-n1.*e.^p1);
dsig1(1)=10*sig1(1);
i=1;
while sig1(i)<dsig1(i)
    e1p(i)=e(i);
    sig1p(i)=sig1(i);
    dsig1(i+1)=(sig1(i+1)-sig1(i))/de;
    i=i+1;
end

sig2=ss2-(ss2-sy2).*exp(-n2.*e.^p2);
dsig2(1)=10*sig2(1);
i=1;
while sig2(i)<dsig2(i)
    e2p(i)=e(i);
    sig2p(i)=sig2(i);
    dsig2(i+1)=(sig2(i+1)-sig2(i))/de;
    i=i+1;
end

sig3=ss3-(ss3-sy3).*exp(-n3.*e.^p3); sig4=s_hs;
dsig3(1)=10*sig3(1);
i=1;
while sig3(i)<dsig3(i)
    e3p(i)=e(i);
    sig3p(i)=sig3(i);
    sig4p(i)=sig4(i);
    sigbp(i)=sigb(i);
    dsig3(i+1)=(sig3(i+1)-sig3(i))/de;
    i=i+1;
end

%% Plot
clf
axes1 = axes('FontSize',16,'PlotBoxAspectRatio',[15 15 15]);
%axis(axes1,[0 0.5 0 2000]);
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');
xlabel(axes1,'Engineering Strain'); 
ylabel(axes1,'Engineering Stress, MPa');
plot(e1p,sig1p,'m-','LineWidth',2);
plot(e2p,sig2p,'g-','LineWidth',2);
plot(e3p,sig3p,'k-','LineWidth',2);
plot(e3p,sig4p,'b-','LineWidth',2);
plot(e3p,sigbp,'r-','LineWidth',2);
text(0.1,300,"C1 = "+cellstr(num2str(c1))+";",'FontSize',16);
text(0.1,150,"C2 = "+cellstr(num2str(c2))+";",'FontSize',16);
legend({'CG','UFG','HS','HSmodel','backstress'},'FontSize',16);

clear axes1 c1 c2 e ro de dsig1 dsig2 dsig3 e1 e1p e2 e2p e3p i n1 n2 n3 p1 p2 p3
clear s_hs sig1 sig1p sig2 sig2p sig3 sig3p sig4 sig4p sigb sigbp ss1 ss2 ss3
clear sy1 sy2 sy3 x

function rss = sfcc(c)
global ro ss1 sy1 n1 p1 ss2 sy2 n2 p2 ss3 sy3 n3 p3
rss=0;
for i=1:100
    ei=0.2*(i-0.2)/(100-1);
    e1=ei/ro;
    e2=ei^2/(ro+ei)*((8-3.1416)/2/3.1416/ro);
    sig1=ss1-(ss1-sy1)*exp(-n1*e1^p1);
    sig2=ss2-(ss2-sy2)*exp(-n2*e2^p2);
    sig3=ss3-(ss3-sy3)*exp(-n3*ei^p3);
    sig4=sig1*(1+c(1))*ro*(1+ei)/(ro+ei)+...
        sig2*c(2)*(1-ro)*(8-3.1416)*(2*ro+ei)*(1+ei)*ei/((ro+ei)*(2*3.1416*ro*(ro+ei)+ei^2*(8-3.1416)));
    rss=rss+(sig3-sig4)^2;
end
end

