% Transient temperature distribution in an infinite plate with the thickness 2L.
% Initially the plate is at Tinit when it is exposed to a fluid with
% the temperature Tfluid.
% Implicit scheme
% Finite volume method
% Solver TDMA
t=[];ap=[];aw=[];ae=[];d=[];Q=[];P=[];
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 =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;
% Implicit formulation
for k = 1:ni % coefficients
   t(k) = Tinit;
   told(k) = Tinit;
   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) = density*cp*dx/dt;
   ap(k)=aw(k)+ae(k)+density*cp*dx/dt;
end;
time = 0;
maxres = 1.e-5;
maxit = 10;
while (time < (maxtime+dt/2))
   told = t;
   P(1) = 0;
   Q(1) = t(1);
   for k = 2:ni-1
      dd = d(k)*told(k);
      a = ap(k);
      if k == 2;
        a = ap(k)+1/(1/alfa+2*dx/conductivity);  
        dd = d(k)*told(k)+Tfluid/(1/alfa+2*dx/conductivity);        
      end;
      P(k) = ae(k)/(a-aw(k)*P(k-1));
      Q(k) = (aw(k)*Q(k-1)+dd)/(a-aw(k)*P(k-1));
   end;
   for k = ni-1:-1:2
     t(k) = P(k)*t(k+1)+Q(k);
   end;
   time = time +dt;
end; %while  time loop
% 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');