three_d_heatsim14.m

clear
%%
%This model attempts to combine the heat transfer model and kinetics model
%to examine how the mass loss of the sample changes with heat.  The
%differnce between this model and three_d_heatsim10.m is that this model
%attempts to account for the non-linear heating rates most nodes experience



%Data
dx = 0.002;%mesh resolution in meters
h = 0.05; %height
w = 0.05; %width
d = 0.01; %depth
Nt = 1700; %#of time steps
dt = 1;%time step 
wt_pct_input = 5; %weight percent (in %)
depth = 4; %depth at which average temperature and graphs of plane 
           %perpendicular to surface being heated are calculated

%constants that dictate change of thermal properties with temperature
c1=0; %controls change in K(default=0.00001)
c2=0; %controls change in cp(default=12.5)

%% Kinetic Properties
Ti=310;
%kinetic parameters obtained from "Synthesis and properties of 
    %monomer cast nylon-6-b-polyether amine copolymers with different
    %structures"
    
%Properties of polyethelyne, see "Pyrolysis of Low and High Density 
%Polyethylene. Part I: Non-isothermal Pyrolysis Kinetics"
Eai=284.37;  %activation energy (kJ/mol)
Eaf=483.37;
dEa=(Eaf-Eai)/25;
Ea=[Eai:dEa:Eaf];  %generate range of activation energy
Pi=1.9*10^(16);  %pre-exponential factor (1/s)
Pf=4.27*10^(22); % 8*10^(2);
dP=(Pf-Pi)/25;
P=[Pi:dP:Pf];  %generate range of activation energy

Bi=0.333;  %external heating rate (K/s or C/s)  NOTE THIS IS NOT THE HEATING
          %RATE EXPERIENCED BY EACH NODE
Bf=10.5;
dB=(Bf-Bi)/5;
B=[Bi:dB:Bf]/60;

R=0.00831;  %gas constant (J/mol*K)
n=1;


%% Material Properties
K_bulk = 0.33; %thermal conductivity (W/m*K)
rho_bulk = 920; %density (kg/m^3)
cp_bulk = 1900; %specific heat (J/kg*K)
K_add = 0.16; %thermal conductivity of additive
rho_add = 1084; %density of additive
cp_add = 7700; %specific heat of additive
H = h/dx + 1; %# of elements in each dimension
W = w/dx + 1; 
D= d/dx + 1;
%arrays of material properties
rho = rho_bulk*ones(H,W,D,Nt); 
cp = cp_bulk*ones(H,W,D,Nt);
K = K_bulk*ones(H,W,D,Nt);
%arrays of kinetic properties
Ea_array=Eai*ones(H,W,D,Nt);
P_array=Pi*ones(H,W,D,Nt);
B_array=Bi*ones(H,W,D,Nt);

%dx_ref=dx*ones(H,W,D,Nt);

Temp_prev = Ti*ones(H,W,D,Nt); %mesh has uniform temperature Ti initially
T = Ti*ones(H,W,D,Nt); %define Current temperature

wt_pct = wt_pct_input/100;
rands=round(((wt_pct)*(rho_bulk/rho_add)*H*W*D)*3); %generate the 
    %number of coordinates(each node has 3 coordinates so multiply by 3) 
    %that will have additive properties in accordance with wt percent
randsx=randi(H,1,rands);
randsy=randi(W,1,rands);
randsz=randi(D,1,rands);
%[rands]=randi([1 (h/dx+1)],1,length_rands);
X1=[];
X2=[];
X3=[];
for i=1:rands %filter random numbers into coordinates
    if mod(i,3)==0
        X1=[X1,randsx(i)];
    elseif mod(i,2)==0
        X2=[X2,randsy(i)];
    else
        X3=[X3,randsz(i)];
    end
end
%[size(X1),size(X2),size(X3)]
for i=1:(rands/3)   
    K(X1(i),X2(i),X3(i),:)=K_add;
    cp(X1(i),X2(i),X3(i),:)=cp_add;
    rho(X1(i),X2(i),X3(i),:)=rho_add;
    Ea_array(X1(i),X2(i),X3(i),:)=10^5; %assume very high Ea of nanoparticles
end
wt_percent_actual = ((rands/3)/(H*W*D))*(rho_bulk/rho_add)*100;

hc = 0.1;  %convection constant
heat_rate = 0; %in degrees per second
IC = 360; %initial temperature at the surface being heated

%%
%Boundary Conditions
T_boundary=410;
T_bottom = T_boundary;
T_left = T_boundary;
T_right = T_boundary;
T_back = T_boundary;
T_front = T_boundary;

%% 
%mesh refining around additives
dx_vectori=(0:dx:h);
dx_vectorj=(0:dx:w);
dx_vectork=(0:dx:d);
dx_array=dx*ones(H,W,D,Nt);
dx_ref=0.02;


%%
%calculations
td = K./(rho.*cp); %thermal diffusivity                                  
FO = (td*dt)./((dx_array).^2); %Fourier number
Biot = (h*dx_array)./K;
condition=FO.*(1+Biot);
T_inf=[];
for i=1:Nt
    if Bi==0;
        T_inf = IC;
    else
        T_inf = [T_inf,Bi*i+IC];
    end
end
if max(max(max(max(FO)))) > 1/6
    disp('Check the inputs')
end
if max(max(max(max(condition)))) > 0.5
    disp('Check the inputs')
end
if FO <= 1/6
    if FO.*(1+Biot) <= 0.5
        for l = 1:Nt-1
            %{
            %Assign bottom boundary condition
            T(:,:,end,:) = T_bottom;
            Temp_prev(:,:,end,:) = T_bottom;
            %Assign Left side BC
            T(:,1,:,:) = T_left;
            Temp_prev(:,1,:,:) = T_left;
            %Assign Right side BC
            T(:,end,:,:) = T_right;
            Temp_prev(:,end,:,:) = T_right;
            %Assign Back side BC
            T(1,:,:,:) = T_back;
            Temp_prev(1,:,:,:) = T_back;
            %Assign Front side BC
            T(end,:,:,:) = T_front;
            Temp_prev(end,:,:,:) = T_front;
            %}
            
            
           %Assign bottom boundary condition
            T(:,:,end,l+1) = T_inf(l);
            Temp_prev(:,:,end,l+1) = T_inf(l);
            %Assign Left side BC
            T(:,1,:,l+1) = T_inf(l);
            Temp_prev(:,1,:,l+1) = T_inf(l);
            %Assign Right side BC
            T(:,end,:,l+1) = T_inf(l);
            Temp_prev(:,end,:,l+1) = T_inf(l);
            %Assign Back side BC
            T(1,:,:,l+1) = T_inf(l);
            Temp_prev(1,:,:,l+1) = T_inf(l);
            %Assign Front side BC
            T(end,:,:,l+1) = T_inf(l);
            Temp_prev(end,:,:,l+1) = T_inf(l);
            %}
            
            T(:,:,1,l+1)=T_inf(l); %Top side is being heated
            Temp_prev(:,:,1,l+1)=T_inf(l);
            
            T(1,1,1,l)=(1/3)*(T_inf(l)+T_back+T_left); %Solve for nodes on vertices (average of three adjacent surfaces)
            T(end,1,1,l)=(1/3)*(T_inf(l)+T_front+T_left);
            T(1,end,1,l)=(1/3)*(T_inf(l)+T_back+T_right);
            T(end,end,1,l)=(1/3)*(T_inf(l)+T_front+T_right);
            T(1,1,end,l)=(1/3)*(T_back+T_left+T_bottom);
            T(1,end,end,l)=(1/3)*(T_back+T_right+T_bottom);
            T(end,1,end,l)=(1/3)*(T_front+T_left+T_bottom);
            T(end,end,end,l)=(1/3)*(T_front+T_right+T_bottom);
            for i = 2:H-1 %mesh calculation
                for j = 2:W-1
                    %for k = 1:2
                    %CurrentTemp(i,j,k,l+1) = %Finite difference on edge exposed to heat flux
                    %end
                    for k = 2:D-1
                       %{   
                        USE THIS WHEN ASSUMPTION THAT K,rho,cp=CONSTANT IS
                        NOT VALID
                        if K(i,j,k,l)~=K_add || rho(i,j,k,l)~=rho_add || cp(i,j,k,l)~=cp_add
                            K(i,j,k,l+1) = K_bulk+c1*T(i,j,k,l);
                            cp(i,j,k,l+1) = cp_bulk+c2*T(i,j,k,l);
                            %rate_const = A*exp(-E/(R.*T));
                            %rn = -rate_const.*rho;
                            %rho(i,j,k,l+1) = rho(i,j,k,l)+rn(i,j,k,l)*dt;
                        end
    
                        %MAKE SURE TO UNCOMMENT BELOW SO THAT PROPERTIES
                        %ARE UPDATED WITH TEMPERATURE WHEN NEEDED
                        
                        %td = K./(rho.*cp); %thermal diffusivity
                        %FO = (td*dt)./((dx_array).^2); %Fourier number
                        %Bi = (h*dx_array)./K;
                        %}
                        
                        %Explicit 3D finite difference heat equation
                        T(i,j,k,l+1) = FO(i,j,k,l)*(Temp_prev(i,j-1,k,l)+...
                        Temp_prev(i+1,j,k,l)+Temp_prev(i,j+1,k,l)+Temp_prev(i-1,j,k,l)+...
                        Temp_prev(i,j,k-1,l)+Temp_prev(i,j,k+1,l))+(1-6*FO(i,j,k,l))*Temp_prev(i,j,k,l);
                    end
                end
            end
            Temp_prev=T; %save last temperature
            %for making animation...
            %M(l)=getframe(gcf);
        end
    end
end
avg=T(:,:,depth,Nt-20);
Sum=sum(sum(avg));
average_T=(Sum)/(H*W)  %average T for a given plane at "depth"

%% kinetics calculations
dT=zeros(size(T));
dT=(T(:,:,:,end)-T(:,:,:,1))./Nt; 
dTt=[];

%Note that this model assumes constant Ea, P.  When
%assumption is invalid, update integrand matrix,
%integral matrix,etc to account for variability
Y=zeros(H,W,D,length(T)); %preallocated integrand matrix
A=zeros(H,W,D,length(T)); %preallocated integral matrix
Y=exp(-(Eai./(R.*T(:,:,:,:))));  %calculate values of integrand
for i=2:length(T)
    A(:,:,:,i)=A(:,:,:,i-1)+0.5*(Y(:,:,:,i)+Y(:,:,:,i-1)).*dT(:,:,:);
end
for i=1:H
    for j=1:W
        for k=1:D
            for l=2:Nt
                    dTt(i,j,k,l)=[(T(i,j,k,l)-T(i,j,k,l-1))];           
                %dTt is the variable heating rate at each temperature step at each node
            end
        end
    end
end

a=zeros(size(A));
for i=1:H
    for j=1:W
        for k=1:D
            for l=1:length(T)
                if dTt(i,j,k,l)<10^(-10)
                    a(i,j,k,l)=0;
                else
                    if n==1
                        a(i,j,k,l)= 1-exp(-(Pi/dTt(i,j,k,l))*(A(i,j,k,l)));
                    elseif n>1
                        a(i,j,k,l)= 1-real((1-(1-n)*((Pi/dTt(i,j,k,l))*(A(i,j,k,l))^(1/(1-n)))));
                    else
                        a(i,j,k,l)= 1-real((1-(1-n)*((Pi/dTt(i,j,k,l))*(A(i,j,k,l))^(1-n))));
                    end
                end
            end
        end
    end
end

for i=1:(rands/3)
    a(X1(i),X2(i),X3(i),:)=0;
end
avg_a=a(:,:,depth,Nt);
Sum_a=sum(sum(avg));
average_a=(Sum)/(H*W);

%{
for i=1:H
    for j=1:W
        for k=1:D
            for l=1:Nt
                if T(i,j,k,l)>550
                    a(i,j,k,l)=1;
                end
            end
        end
    end
end
%}

%% Weight Calcualtions
%initial_weight(:,:,:)=(dx^3)*rho_add*a(:,:,:,1));
%initial_weight(:,:,:)=(dx^3)*rho_add*(1-a(:,:,:,1));
%initial_weight_tot=sum(sum(sum(initial_weight)));
if wt_pct_input~=0
for i=1:(rands/3)
    initial_weight_add(i)=(dx^3)*rho_add*(1-a(X1(i),X2(i),X3(i),1));
end
initial_weight_tot_add=sum(initial_weight_add);
end
for i=1:H
    for j=1:W
        for k=1:D
            for l=1:Nt
                if K(i,j,k,l)~=K_add || rho(i,j,k,l)~=rho_add || cp(i,j,k,l)~=cp_add
                    weight_bulk(i,j,k,l)=(dx^3)*rho_bulk*(1-a(i,j,k,l));
                    xx(i,j,k,l)=1-a(i,j,k,l);
                end
            end
        end
    end
end
weight_tot_bulk=sum(sum(sum(weight_bulk)));
initial_weight_tot_bulk=max(weight_tot_bulk);

if wt_pct_input~=0
    %xx_tot=sum(sum(sum(xx)));
    weight_tot=sum(sum(sum(weight_bulk)))+initial_weight_tot_add;
    initial_weight_tot=initial_weight_tot_bulk+initial_weight_tot_add;
    pct_mass=squeeze(squeeze(weight_tot))./initial_weight_tot;
else
    xx_tot=sum(sum(sum(xx)));
    weight_tot_bulk=sum(sum(sum(weight_bulk)));
    initial_weight_tot_bulk=max(weight_tot_bulk);
    pct_mass=squeeze(squeeze(weight_tot_bulk))./initial_weight_tot_bulk;
end


for i=2:length(T)
%    dTt13(i)=(T(13,13,3,i)-T(13,13,3,i-1))/dT(13,13,3);
    %dTt is the variable heating rate at each temperature step at node
    %(13,13,3)
end
%{
plot(1:Nt,pct_loss)
plot(1:Nt,squeeze(squeeze(T(13,13,3,:))))
plot(1:Nt,dTt13)
plot(squeeze(squeeze(T(13,13,3,:))),squeeze(squeeze(1-a(13,13,3,:))))
%}
    
%% Graphs

x = 1:W
y = 1:H;
z = 1:D;
lx = length(x);
ly = length(y);
lz = length(z);
x1 = w*x/lx;
y1 = h*y/ly;
z1 = d*z/lz;

figure(1)
title('Temperature Profile At Constant Depth (dx=0.002)')
contourf(squeeze(T(:,:,3,Nt-20)))
cbar2=colorbar;
colormap(hot)
title(cbar2,'Temperature(K)')

figure(2)
surf(squeeze(T(2:H-1,2:W-1,2,Nt-5)))
title('Planar Temperature Parallel to Heated Side')
zlabel('Temperature (K)')

figure(3)
plot(squeeze(squeeze(T(13,13,3,:)-273.15)),squeeze(squeeze(1-a(13,13,3,:))))
title('TGA Curve for Polyethelyne Composite')
xlabel('Temperature (C)')
ylabel('Mass Fraction')

%{
figure(1)
subplot(2,2,1)
pcolor(squeeze(T(:,5,:,Nt-20)))
colorbar

figure(2)
subplot(2,2,2)
contourf(squeeze(T(:,5,:,Nt-20)))
colorbar

figure(3)
subplot(2,2,3)
pcolor(T(:,:,depth,Nt-20))
colorbar
%}

%{
figure(1)
contourf(T(:,:,depth,Nt-20))
colormap(hot)
colorbar
movie(M)
%movie(T(:,:,1:end,199))
figure(2)
title('Temperature Profile With Varying Depth(dx=0.002)')
contourf(squeeze(T(:,10,1:(depth+2),Nt-20)))
cbar1=colorbar;
title(cbar1,'Temperature(K)')



figure(3)
subplot(2,2,3)
pcolor(T(:,:,5,Nt-20))
%colorbar

%figure(4)
figure(3)
title('Temperature Profile At Constant Depth (dx=0.002)')
contourf(T(:,:,depth,Nt-5))
cbar2=colorbar;
%colormap(hot)
title(cbar2,'Temperature(K)')

figure(4)
title('Temperature Profile At Constant Depth (dx=0.002)')
contourf(T(round(0.1*H):round(0.9*H),round(0.1*W):round(0.9*W),depth,Nt-5))
cbar2=colorbar;
%colormap(hot)
title(cbar2,'Temperature(K)')

figure(5)
surf(squeeze(T(2:H-1,2:W-1,depth,Nt-5)))
title('Planar Temperature Parallel to Heated Side')
zlabel('Temperature (K)')

figure(6)
subplot(2,1,1)
surf(squeeze(a(round(0.1*H):round(0.9*H),round(0.1*W):round(0.9*W),depth+1,end)))
subplot(2,1,2)
surf(squeeze(T(round(0.1*H):round(0.9*H),round(0.1*W):round(0.9*W),depth+1,end)))

figure(7)
plot(squeeze(squeeze(T(13,13,3,:))),squeeze(squeeze(a(13,13,3,:))))
%}

%{
figure
plot(T(5,5,5,1:100))

subplot(2,2,1)
pcolor(x1,y1,CurrentTemp(:,5,:,t))
shading interp

figure(2)
subplot(2,2,2)
contourf(x1,y1,CurrentTemp(:,5,:,t),'ShowText','on')


length_x = length(1:W);
length_y = length(1:H);
subplot(2,2,3)
pcolor(length_x,length_y,CurrentTemp)
shading interp
subplot(2,2,3)
contour(length_x,length_y,CurrentTemp)
disp('Total Time:')
%}

%%
%{
ideas
CurrentTemp(size(1),size(2),d3,t+1) = Ti*ones(d1,d2,d3,t); %top, bottom, and right sides are isothermal
              CurrentTemp(size(1),:,t+1) = Ti*ones(1,size(2),t);
              TempNew(:,size(2)) = Ti*ones(size(1),size(2)) figure out how to
              make right side isothermal
              
              
 exiting the loop if constraints are not met    
 
 a = K./(rho.*cp); 
 FO = (a*dt)/(dx)^2; 
 Bi = (h*dx)/K;
 if FO > 1/6
    break
 end
 if FO.*(1+Bi) > 0.5
    break
  end
dx_vectori=[0:dx:h];
dx_vectorj=[0:dx:w]
dx_vectork=[0:dx:d];
for i=2:length(dx_vectori)
    for j=2:length(dx_vectorj)
        for k=2:length(dx_vectork)
  

  for i=2:length(dx_vectori)
     if (K(i,j,k,l)==K_add)
         for i=i-1:1+1
             dx_array(i,j,k,l)=[dx_array,dx_ref*ones(2,2,2,Nt)];
         end
     end
     for j=2:length(dx_vectorj)
         if (K(i,j,k,l)==K_add)
             for j=j-1:j+1
                 dx_array(i,j,k,l)=[dx_array,dx_ref*ones(2,2,2,Nt)];
             end
         end
         for k=2:length(dx_vectork)
             if (K(i,j,k,l)==K_add)
                 for k=k-1:k+1
                     dx_array(i,j,k,l)=[dx_array,dx_ref*ones(2,2,2,Nt)];
                 end
             end
         end
     end
end


for l=1:Nt
    for i=2:H-2
        for j=2:W-2
            for k=2:D-2
                if (K(i,j,k,l)==K_add)||(rho(i,j,k,l)==rho_add)||(cp(i,j,k,l)==cp_add);
                    for i=i-1:1+1
                        for j=j-1:j+1
                            for k=k-1:k+1
                                dx_array(i,j,k,l)=[dx_array(i,j,k,l),(dx_ref)*ones((dx_ref/dx),(dx_ref/dx),(dx_ref/dx),l)];
                            end
                        end
                    end
                end
            end
        end
    end
end
%}