% Transient temperature distribution in aninfinite plate with the thickness 2L.
% Initially temperature is Tinit when it is exposed to a fluid with
% the temperature Tfluid.
% Explicit method
%
clear all;
t=[];ap=[];aw=[];ae=[];d=[];
n_cv = 30; % number of control volumes
ni = n_cv+2; % number of node points
L = 0.05; % half thickness
dx = L/n_cv % increment
dt =0.1 % time step
conductivity = 40; % thermal conductivity
density = 8000; % density
cp = 500 % Heat capacity
alfa = 1600; % heat transfer coefficient
maxtime = 100; % total time
Tinit = 1; % initial temperature (dimensionless)
Tfluid = 0; % fluid temperature (dimensionless)
x(1) = 0;
x(2) = dx/2;
for k = 3:ni-1
   x(k) = (k-1)*dx;
end;
x(ni) = L;
% Explicit formulation
for k = 1:ni % coefficients
   t(k) = Tinit;
   told(k) = Tinit;
   ap(k) = density*cp*dx/dt;
   ae(k) = conductivity/dx;
   aw(k) = conductivity/dx;
   if (k==2);aw(k)=0;end;
   if (k==ni-1);ae(k)=0;end;  
   d(k) = ap(k)- aw(k)-ae(k);
end;
time = 0;
while (time < (maxtime+dt/2))
   told = t;
   for k = 2:ni-1
       dd = d(k)*told(k);
       if k == 2;
         dd = d(k)*told(k)+(Tfluid-told(k))/(1/alfa+2*dx/conductivity);  
       end;  
       t(k) = (ae(k)*told(k+1)+aw(k)*told(k-1)+dd)/ap(k);
   end;
   time = time +dt;
end; %while
% Boundaries
t(ni)=t(ni-1)
t(1)=(alfa*Tfluid+2*conductivity/dx*t(2))/(alfa+2*conductivity/dx);
%
plot(x,t,'*');ylabel('Temperature');xlabel('x-pos');
title('Transient temperature distribution');