MATLAB: Problem with graph MATLAB

errorgraphiterationMATLAB

I am trying to plot a creep x time graph using matlab. However, there are some mismatches between my graph and the reference paper.
This is my code:
% Parameters
c1= 6.48; %MPa



c2= 5.17; %GPa.d
a1= 1.29e-7; %MPa^-5.d^-1

a2= 1.6e-11; %MPa^-5.d^-1
b=66.28;
n=5;
r=2.6;
G=26.9; %GPa
u=0.4;
t(1)=0;
et(1)=0;
es(1)=0;
ed(1)=0;
etotal(1)=0;
sigma1=26; %MPa
sigma2=1; %MPa
sigma3=1; %MPa
% Loop
for i=1:355
t(i+1)=t(i)+0.01;
%Calculation
sigma= sigma1-sigma3;
sigmam=(sigma1+sigma2+sigma3)/3;
sigmadamage=sigma*((2/3)*(1+u+3*(1-2*u)*((sigmam/sigma)^2)))^0.5
% Damage Factor
D(i+1)= 1-(1-(r+1)*((sigmadamage/b)^r)*t(i+1))^(1/(r+1));
% Creep
%Primary Creep
et(i+1)=(sigma/c1)*(1-exp(-G*t(i+1)/c2));
%Secondary Creep
es(i+1)=a1*(sigma^n)*t(i+1);
%Tertiary Creep
ed(i+1)=a2*((sigma/(1-D(i+1)))^n)*t(i+1);
%Total Creep
etotal(i+1)= et(i+1)+es(i+1)+ed(i+1);
end
p = plot(t,etotal,'-b');
p(1).LineWidth = 1.5;
title('Creep x Time')
xlabel('Time(day)')
ylabel('Creep (%)')
grid on
legend([p(1)], 'Creep Evolution')
I really cannot see which type of error in my code is provoking this mismatch.
Thanks
Additional data:

Best Answer

Your code seems to be right, but the for loop is really not suitable in this case. However, there is a slight variation in the result when the for loop is omitted. You might want to check your parameters to be sure that they are entered correctly, for instance, I see that you entered both units of GPa and MPa in the same way. I think you should enter all pressure as either GPa or MPa to be consistent. Here the code without the for loop.
clear variables
close all
% Parameters
c1 = 6.48;    %MPa



c2 = 5.17;    %GPa.d
A1 = (1.29e-7); %MPa^-5.d^-1

A2 = (1.6e-11); %MPa^-5.d^-1
B = 66.280;
n = 5;
r = 2.6;
G = 26.9; %GPa
mue = 0.4;
sigma1 = 26; %MPa
sigma2 = 1; %MPa
sigma3 = 1; %MPa
sig_hat = sigma1-sigma3;
sigm =  (sigma1+sigma2+sigma3)/3;
a1 = 1+ mue;
a2 = 3*(1-2*mue);
a3 =(sigm/sig_hat)^2;
sigmastar = sig_hat*sqrt(2/3*(a1+a2*a3));
tf = 4;
N = 355;
t = 0:tf/(N-1):tf;
b1 = (1+r);
b2 = (sigmastar/B)^r;
D = 1-((1-b1*b2*t)).^(1/b1);
et = (sig_hat/c1)*(1-exp(-G*t/c2));
es =  A1*(sig_hat^n)*t;
ed = (A2*(sig_hat./(1-D)).^n).*t;
%Total Creep
ec = et+es+ed;
p = plot(t,real(ec),'-b');
hold on
p(1).LineWidth = 1.5;
title('Creep x Time')
xlabel('Time(day)')
ylabel('Creep (%)')
grid on
% ylim([0 15])
%%Your code without the for loop
% c1= 6.48; %MPa
% c2= 5.17; %GPa.d
% a1= 1.29e-7; %MPa^-5.d^-1
% a2= 1.6e-11; %MPa^-5.d^-1
% b=66.28;
% n=5;
% r=2.6;
% G=26.9; %GPa
% u=0.4;
% 


% sigma1=26; %MPa
% sigma2=1; %MPa
% sigma3=1; %MPa
% % Loop
% N = 350;
% t = 0:4/N:4;
% 
% %Calculation
% sigma= sigma1-sigma3;
% sigmam=(sigma1+sigma2+sigma3)/3;
% sigmadamage=sigma*((2/3)*(1+u+3*(1-2*u)*((sigmam/sigma)^2)))^0.5;
% % Damage Factor
% D= 1-(1-(r+1)*((sigmadamage/b)^r)*t).^(1/(r+1));
% % Creep
% %Primary Creep
% et=(sigma/c1)*(1-exp(-G*t/c2));
% %Secondary Creep
% es= a1*(sigma^n)*t;
% %Tertiary Creep
% ed=a2*((sigma./(1-D)).^n).*t;
% %Total Creep
% etotal= et+es+ed;
% 
% p = plot(t,real(etotal),'-b');
% p(1).LineWidth = 1.5;
% title('Creep x Time')
% xlabel('Time(day)')
% ylabel('Creep (%)')
% grid on
% legend([p(1)], 'Creep Evolution')
% % legend([p(1)], 'Creep Evolution')
Related Question