% 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;
n = 20; % number of increments
nk = n+1; % number of node points
L = 0.05; % half thickness
dx = L/n % increment
dt =0.1 % time step
conductivity = 40; % thermal conductivity
density = 8000; % density
cp = 500 % Heat capacity
alfa = 1600; % heat transfer coefficient
diffusivity = conductivity/(density*cp); % thermal diffusivity
maxtime = 100; % total time
Tinit = 1; % initial temperature (dimensionless)
Tfluid = 0; % fluid temperature (dimensionless)
for k = 1:nk
   x(k) = (k-1)*dx;
end;
% Explicit formulation
for k = 1:nk % coefficients
   t(k) = Tinit;
   told(k) = Tinit;
   a(k) = 1;
   b(k) = diffusivity*dt/dx^2;
   c(k) = diffusivity*dt/dx^2;
   d(k) = 1 - 2*diffusivity*dt/dx^2;
end;
% second order approximation at boundary
c(1) = 0;
b(1) = 4;
a(1) = 3 + 2*alfa*dx/conductivity;
a(nk) = 3;
b(nk) = 0;
c(nk) = 4;
time = 0
while (time < (maxtime+dt/2))
   told = t;
   d(1) = 2*Tfluid*alfa*dx/conductivity - t(3);
   t(1)=(b(1)*t(2)+d(1))/a(1);
   for k = 2:n
      t(k) = (b(k)*told(k+1)+c(k)*told(k-1)+d(k)*told(k))/a(k);
   end;
   d(nk)=-t(nk-2);
   t(nk) =(c(nk)*t(nk-1)+d(nk))/a(nk);
   time = time +dt;
end; %while
%
plot(x,t,'*');ylabel('Temperature');xlabel('x-pos');
title('Transient temperature distribution');