clc
clear
close all
%% LOADING CONDITION %%
t = 30;

%% MATERIAL PARAMETERS %%
A = 3*10^-10;
sigma_cr = 47*10^6;
R = 8.314;
E_fiber = 224*10^9;

i = 0;
j = 0;
k = 0;

for stress = 10^4:10^5:6*10^7;
    i = i+1;
    for Temp = 400:0.1:600;
        j = j+1;
        
% power-law exponent n        
n = 3.1943019-0.002283*Temp;
        
% fit activition energy using bilinear interpolation
if Temp<=548
    Q = 123*10^3;
end

if Temp>548 && Temp<=558
    if stress<=0.45*sigma_cr
        Q = 123*10^3;
    end
    if stress>0.45*sigma_cr && stress<=0.6*sigma_cr
        a = 123-(stress/sigma_cr-0.45)/0.15*3;
        Q = (123-(Temp-548)/10*(123-a))*10^3;
    end
    if stress>0.6*sigma_cr
        a = 120-(stress/sigma_cr-0.6)/0.1*4;
        Q = (123-(Temp-548)/10*(123-a))*10^3;
    end
end

if Temp>558
    if stress<=0.3*sigma_cr
        Q = 123*10^3;
    end
    if stress>0.3*sigma_cr && stress<=0.45*sigma_cr
        a = 123-(stress/sigma_cr-0.3)/0.15*4;
        Q = (123-(Temp-558)/15*(123-a))*10^3;
    end
    if stress>0.45*sigma_cr && stress<=0.6*sigma_cr
        a = 119-(stress/sigma_cr-0.45)/0.15*5;
        b = 123-(stress/sigma_cr-0.45)/0.15*3;
        Q = (b-(Temp-558)/15*(b-a))*10^3;
    end
    if stress>0.6*sigma_cr
        a = 114-(stress/sigma_cr-0.6)/0.1*8;
        b = 120-(stress/sigma_cr-0.6)/0.1*4;
        Q = (b-(Temp-558)/15*(b-a))*10^3;
    end
end

x(i,j) = Temp;
y(i,j) = stress;
eps_ss = A*stress^n*exp(-Q/(R*Temp));
Q1(i,j) = Q;

if Temp<=523
slope = 0.00687.*523.^2-7.2317.*523+1877.184;
end
if Temp>523
slope = 0.00687.*Temp.^2-7.2317.*Temp+1877.184;
end

% elastic modulus
E_ave = -1.2148*10^7.*Temp+8.5038*10^9;
E = slope.*(stress-0.51*sigma_cr)+E_ave;
eps_elastic = stress/E;

% transient constant
if Temp<=548
    C = 10^(6.6184)*(stress/10^6)^(-2.1575);
end

if Temp>548
    T1 = 548;
    a = -3.8947+(573-Temp)/(573-T1)*1.7372;
    b = 8.3528-(573-Temp)/(573-T1)*1.7344;
    C = 10^(b)*(stress/10^6)^(a);
end

% transient strain
eps_t = ((1.022*10^-9*exp(0.0227*Temp))*exp((0.00689*exp(0.00487*Temp))*(stress/10^6)));

% plastic strain calcualte by power-law model
eps_plastic = eps_t*(1-exp(-C*eps_ss*t))+eps_ss*t;

% total creep strain on polymer matrix
eps_m = eps_elastic+eps_plastic;

% total loading stress
stress_total(i,j) = 0.6*E_fiber*eps_m+0.4*stress;

trial = eps_plastic; %total plastic strain

eps(i,j) = eps_m;   % total strain (equals to strain on polymer matrix)

if trial>0.001999 && trial<0.002
    k = k+1;
    x1(k) = Temp;
    y1(k) = stress; 
    stress_design(k) = (0.6*E_fiber*eps_m+0.4*stress)/1.25;   %design stress
end

    end
    
    j = 0;
    
end

%% TRANSIENT MAP AND DESIGN SPACE %%
temp = [400 600];
ultimate = [2070 2070];   %ultimate tensile strength of composite

figure
v = [6*10^-3 8*10^-3 10^-2 0.012 0.014 0.016 0.018 0.02 0.022];
c = contour(x,stress_total/10^6,eps,v);
clabel(c)
hold on
plot(x1,stress_design/10^6,'k','LineWidth',2)
hold on
plot(temp,ultimate,'--k','LineWidth',2)
axis([400 580 8*10^2 32*10^2])
xlabel('Temperature (K)')
ylabel('Stress \sigma (MPa)')
text(460,2170,'ultimate tensile strength in axial fiber direction')
title('T650-35/PMR-15 in axial fiber direction, for a time of 30 seconds')

%% DESIGN STRESS %%
figure
plot((x1-273.15)*9/5+32,stress_design/10^6,'LineWidth',2)
hold on
xlabel('Temperature (F)')
ylabel('Design Stress \sigma_D (MPa)')
