clc
clear
close all

%% INPUT PARAMETERS %%
% steady-state parameters fitted from experimental data
A = 8*10^16;   
C = 4.5*10^6;
n = 3.7;
% parameters from literature
Q = 450000;
Tm = 1350;
mu_0 = 54000*10^6;
mu_1 = -17*10^6;
k = 8.314;
delta = 0.035;
f = 0.6;
h = 50*10^6;
l = 2.5;
poisson = 0.34;
alpha = 16.5*10^-6;
delta_T = 0;

% initial drag strength
D_0 = C/17;

% initial temperature & thermal parameters
T_0 = 773;
rho = 9.131*10^3;
c = 3860;

% temperature-dependent function
T = 773;
if  T>0.5*Tm
    v = exp(-Q/(k*T));
end 

if T<0.5*Tm
    v = exp(-Q/(k*0.5*Tm)*(log(0.5*Tm/T)+1));
end

% shear moudulus and elastic modulus
mu = mu_0+mu_1*T;
E = 2*mu*(1+poisson);

%% CALCULATE CREEP STRAIN %%
i = 0;
j = 0;
m = 0;

% initial guess value
sigma = 0; 
beta = 0;
epsdot_p_mag = 0;
epsdot_p = sign(sigma-beta)*epsdot_p_mag/sqrt(3);
D = D_0;
L = 0;
xi = 0;
lambda = l*xi;
r = 0;
H = 0.1*mu;

% input loading stress
sigma_creep = 72.9*10^6;

   for t = 0:0.01:8;
        delta_t = 0.01;
        i = i+1
         
if  T>0.5*Tm
    v = exp(-Q/(k*T));
end 

if T<0.5*Tm
    v = exp(-Q/(k*0.5*Tm)*(log(0.5*Tm/T)+1));
end

mu = mu_0+mu_1*T;
E = 2*mu*(1+poisson);

%% INPUT STRESS %%

if t<8
    sigma = sigma_creep/8*t;
end

if t>=8
    sigma = sigma_creep;
end
    
%% VALUES FROM LAST TIME STEP %%
L = f*(D-D_0)*(C-D)/(delta*C);

if L == 0 && sigma == 0
    xi = 0;
end

if L ~= 0
xi = 0.5*abs(sign(sigma-beta)-beta/(sqrt(3)*L));
end

H = (0.1+0.9*xi)*mu;

lambda = 0;
r = v*(sinh(sign(((D-D_0)/(delta*C)))*(abs((D-D_0)/(delta*C)))))^n;

if L == 0 && sigma == 0 
    beta_dot = 0;
end

if L ~= 0
beta_dot = 3*H*(epsdot_p-beta/(3*L)*epsdot_p_mag);
end

D_dot = h*((1-lambda)*epsdot_p_mag-r);
    
if i == 2
    S_magnitude = sqrt(1/3)*100*10^6;   %initial magnitude of deviatoric stress
    beta = f*(C-D_0-delta*S_magnitude)*S_magnitude/C/sqrt(1/3);
    D = D_0+delta*S_magnitude;
end

%% VALUES FOR CURRENT TIME STEP %%
beta = beta+delta_t*beta_dot;
D = D+D_dot*delta_t;
Y = (1-f)*(D-D_0)*(C-D)/(delta*C);

beta1(i) = beta;
D1(i) = D;

%% CALCULATE STRAIN %%
if (abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D) >= 0
Z = A*(sinh((abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D)))^n;

epsdot_p_mag = v*Z;
epsdot_p = epsdot_p_mag/sqrt(3)*sign(sigma-beta);

if i==1
    eps_p(i) = 0;
end
if i>1
eps_p(i) = delta_t*epsdot_p+eps_p(i-1);
end

eps = sigma/E+eps_p(i)+alpha*delta_T;

end

if (abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D) < 0
Z = 0;

epsdot_p_mag = 0;
epsdot_p = 0;

if i==1
    eps_p(i) = 0;
end
if i>1
eps_p(i) = eps_p(i-1);
end

eps = sigma/E+eps_p(i)+alpha*delta_T;

end


trial_3(i) = epsdot_p_mag; % don't change
trial_4(i) = delta_T;

stress(i) = sigma/10^6;
strain(i) = eps;
time(i) = t;


%% DISSIPATION RATE %%
sigma_dot(1) = stress(1)*10^6/delta_t;
if i>1
    sigma_dot(i) = (stress(i)-stress(i-1))*10^6/delta_t;
end

alpha_dot(1) = alpha*trial_4(1)/delta_t;
if i>1
    alpha_dot(i) = (alpha*trial_4(i)-alpha*trial_4(i-1))/delta_t;
end

eps_dot(1) = 0;
if i>1
    eps_dot = (strain(i)-strain(i-1))/delta_t;
end

W_dot(i) = sigma*(eps_dot-sigma_dot(i)/E-alpha_dot(i));

%% ENERGY DISSIPATION %%
W(1) = 0;
if i>=2
t_1(i) = t;   % time vector
W(i) = trapz(t_1,W_dot);
end

%% TEMPERATURE INCREASE %%
Tem(i) = T_0+W(i)/(rho*c);
T = Tem(i);
delta_T = T-T_0;
trial(i) = delta_T;

trial_5(i) = eps_dot;

    end
    
for t = 8:1:50000;
        i = i+1
        
delta_t = 1;
         
if  T>0.5*Tm
    v = exp(-Q/(k*T));
end 

if T<0.5*Tm
    v = exp(-Q/(k*0.5*Tm)*(log(0.5*Tm/T)+1));
end

mu = mu_0+mu_1*T;
E = 2*mu*(1+poisson);

%% INPUT STRESS %%
if t<8
    sigma = sigma_creep/8*t;
end

if t>=8
    sigma = sigma_creep;
end

%% VALUES FROM LAST TIME STEP %%
L = f*(D-D_0)*(C-D)/(delta*C);

if L == 0
    xi = 0;
end
if L ~= 0
xi = 0.5*abs(sign(sigma-beta)-beta/(sqrt(3)*L));
end

H = (0.1+0.9*xi)*mu;

lambda = 0;
r = v*(sinh(sign(((D-D_0)/(delta*C)))*(abs((D-D_0)/(delta*C)))))^n;

if L == 0 
    beta_dot = 0;
end
if L ~= 0
beta_dot = 3*H*(epsdot_p-beta/(3*L)*epsdot_p_mag);
end

D_dot = h*((1-lambda)*epsdot_p_mag-r);
    
if i == 2
    S_magnitude = sqrt(1/3)*sigma_creep;
    beta = f*(C-D_0-delta*S_magnitude)*S_magnitude/C/sqrt(1/3);
    D = D_0+delta*S_magnitude;
end

%% VALUES FOR CURRENT TIME STEP %%
beta = beta+delta_t*beta_dot;
D = D+D_dot*delta_t;
Y = (1-f)*(D-D_0)*(C-D)/(delta*C);

beta1(i) = beta;
D1(i) = D;


%% CALCULATE STRAIN %%
  
if (abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D) >= 0
Z = A*(sinh((abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D)))^n;

epsdot_p_mag = v*Z;
epsdot_p = epsdot_p_mag/sqrt(3)*sign(sigma-beta);

if i==1
    eps_p(i) = 0;
end
if i>1
eps_p(i) = delta_t*epsdot_p+eps_p(i-1);
end

eps = sigma/E+eps_p(i)+alpha*delta_T;

end

if (abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D) < 0
Z = 0;

epsdot_p_mag = 0;
epsdot_p = 0;

if i==1
    eps_p(i) = 0;
end
if i>1
eps_p(i) = eps_p(i-1);
end

eps = sigma/E+eps_p(i)+alpha*delta_T;

end


trial_3(i) = epsdot_p_mag; % don't change
trial_4(i) = delta_T;

stress(i) = sigma/10^6;
strain(i) = eps;
time(i) = t;


%% DISSIPATION RATE %%
sigma_dot(1) = stress(1)*10^6/delta_t;
if i>1
    sigma_dot(i) = (stress(i)-stress(i-1))*10^6/delta_t;
end

alpha_dot(1) = alpha*trial_4(1)/delta_t;
if i>1
    alpha_dot(i) = (alpha*trial_4(i)-alpha*trial_4(i-1))/delta_t;
end

eps_dot(1) = 0;
if i>1
    eps_dot = (strain(i)-strain(i-1))/delta_t;
end

W_dot(i) = sigma*(eps_dot-sigma_dot(i)/E-alpha_dot(i));

%% ENERGY DISSIPATION %%
W(1) = 0;
if i>=2
t_1(i) = t;   % time vector
W(i) = trapz(t_1,W_dot);
end

%% TEMPERATURE INCREASE %%
Tem(i) = T_0+W(i)/(rho*c);
T = Tem(i);
delta_T = T-T_0;
trial(i) = delta_T;


trial_5(i) = eps_dot;

    end
 

%% PLOT FIGURES %%    
figure
plot(strain,stress,'LineWidth',2)
xlabel('strain \epsilon')
ylabel('stress \sigma (MPa)')
grid on

figure
plot(time,strain,'LineWidth',2)
xlabel('time')
ylabel('strain \epsilon')
grid on

figure
plot(time,W_dot)
hold on
grid on

figure
plot(time,Tem,'LineWidth',2)
xlabel('time')
ylabel('temperature (K)')
grid on
