Problem Using Finite Difference Method to Simulate 1D Heat Conduction (UPDATED)

I am writing a script to perform a 1D heat transfer simulation on a system of two materials (of different k) with convection from a flame on one side and free convection (assumed room temperature) at the other. I am trying to employ central finite difference method to solve the general equation for conduction through the material. I am confident in my boundary conditions, though my constants still need to be tweaked (not the problem at hand).

The basic overview of the system is an aluminum disk with a coating or barrier designed to reduce heat transfer from hot gases (such as combusting gases in an ICE). The script must solve the convection to the barrier, conduction through the barrier, and conduction through the disk.

Here is the script:

dt = 0.001; %time step is 1 millisecond
dx_B = 1E-6; %depth step is one micron
dx_W = 1E-4;
t = [0:dt:15]; %total time is 15 seconds
Nt = 15/dt; %number of time steps
Tgas = 550 + 200*sin(2*pi*t); %gas temperature fluxuates sinusoidally
Nx_B = 1000; %depth steps; thickness of barrier (1mm = 1000microns)
Nx_W = 1000; %thickness of wall is 10 cm

x_B = [0:dx_B:0.001];
x_W = [0:dx_W:0.1];

T_B = 300*ones(Nx_B+1,Nt+1);
T_W = 300*ones(Nx_W+1,Nt+1);

h = 50;
cp = 50000;
k_B = 10;
k_W = 270;

maxiter = 500;
for iter = 1:maxiter
    %%T_B(:,1) = T_W(:,end); %Initialize the temp at t=0 to the last temp

      for i=2:Nt+1
      T_B(1,i) = T_B(1,i-1) + (h/cp)*(Tgas(i-1) - T_B(1,i-1)); %Convection boundary layer between surface and gas
      depth_2D = (T_B(1:end-2,i-1) - 2*T_B(2:end-1,i-1) + T_B(3:end,i-1))/dx_B^2; %use central FDM to solve d^2(T)/dx^2
          time_1D = k_B*depth_2D; %dT/dt = k*d^2(T)/dx^2
          T_B(2:end-1,i) = time_1D*dt + T_B(2:end-1,i-1);

          %T_B(end,i) = T_B(end-1,i) - (k_W/k_B)*(T_W(2,i) - T_W(1,i));
          T_W(1,i) = T_W(2,i) - (k_B/k_W)*(T_B(end,i) - T_B(end-1,i)); %Conduction into boundary for barrier is equal to conduction out of boundary for wall

          %depth_2D = (T_W(1:end-2,i)-2*T_W(2:end-1,i)+T_W(3:end,i))/dx_W^2;
          %time_1D = k_W*depth_2D;
          %T_W(2:end-1,i) = time_1D*dt + T_W(2:end-1,i-1);
          %T_W(end,i) = 300;
      end
  err = T_W(end-1,end) - T_W(end,end);
  if err<1 break;
  end
  if iter == maxiter;
      warning('convergence not reached')
  end
  end

figure(1)
plot(t,Tgas,'b',t,T_B(1,:),'k')

figure(2)
plot(t,T_B(end-1,:))

figure(3)
depth = [0:dx_B:Nx_B*dx_B];
contourf(t,depth,T_B); title('Temperature plot (contourf)')
colorbar

You can see that I have commented out the second part of the loop where FDM is used for the second material. The problem I am having is in the first implementation of FDM for conduction through the barrier. The surface temperature (first row of T_B matrix) fluctuates appropriately according to the convection equation, however all the other values in the matrix remain constant. For some reason the script is not modeling any heat transfer past the surface through the solid. I am certain the issue is in either my syntax or indexing, as other users employ FDM for conduction problems quite frequently.

If anyone has any ideas or knows how to fix this please help. Or if you see a general improvement to be made in my script (any part of it) feedback is always appreciated. Thanks in advance for all help.

Cheers, David

[David]

UPDATE: I tried adding another for loop inside the time loop to solve for the temp change over x. It did not give me better results, but maybe I'm moving in the right direction?

 for i=2:Nt+1
    T_B(1,i) = T_B(1,i-1) + (h/cp)*(Tgas(i-1) - T_B(1,i-1)); %Convection boundary layer between surface and gas
    for j = 2:Nx_B
    depth_2D = (T_B(j-1,i-1) - 2*T_B(j,i-1) + T_B(j+1,i-1))/(dx_B^2) %use central FDM to solve d^2(T)/dx^2
        time_1D = k_B*depth_2D; %dT/dt = k*d^2(T)/dx^2
        T_B(j,i) = time_1D*dt + T_B(j-1,i-1);
    end

        T_W(1,i) = T_W(2,i) - (k_B/k_W)*(T_B(end,i) - T_B(end-1,i)); %Conduction into boundary for barrier is equal to conduction out of boundary for wall

        %depth_2D = (T_W(1:end-2,i)-2*T_W(2:end-1,i)+T_W(3:end,i))/dx_W^2;
        %time_1D = k_W*depth_2D;
        %T_W(2:end-1,i) = time_1D*dt + T_W(2:end-1,i-1);
        %T_W(end,i) = 300;
    end