Having trouble modelling 2-dimensional transient heat transfer

I'm having trouble modelling 2-dimensional transient heat transfer in MATLAB. I've reduced the problem down to a very simple model, simply an insulated box split into a 2x2 matrix with a heat generation term inside, but the values I'm getting out are way off what should be correct. The code makes use of the ODE45 command in MATLAB. This is the code I'm using, any help is much appreciated! The version of MATLAB that I'm using is 7.10 (R2010a)

This is the code I'm using:

Main file

%%Basic 2 dimensional heat transfer code
clear all
close
clc

global l h k mass spec_heat T_amb m n z

% Box parameters
l = 1; % length
h = 1; % height
k = 0.2; % W/m.K
mass = 100; % kg
spec_heat = 4000; % J/kg

% Initial conditions
T_amb = 288; % Kelvin

% Number of elements
m = 2; % x direction
n = 2; % y direction
z = m*n; % total number of elements

T0 = zeros(m*n, 1);

for i = 1:(m*n)
  T0(i,1) = T_amb; % T0(i,1) is initial temp of the material in K
end

tdur = 3456000;                         % Duration of simulation in seconds (s).
tstep = 3600;                           % Time step (of 1 hour) in seconds.
tspan = [0:tstep:tdur];  

[t, T] = ode45(@twod_function, tspan, T0);

plot (t,T(:,1),'-r',t,T(:,2),'b',t,T(:,3),'-y',t,T(:,4),'--')

And for the function file

function dTdt = twod_function(t,T)

global l h k mass spec_heat T_amb m n z
R = (l/m)/(k*l*h); % resistance of material in box
Q_gen = 50;

 Q_cond = zeros(1,(m*n));
 for i = 1:(m*n)
    if i == 1
        Q_cond(i) = (T(i) - 288)/R;
    else
    Q_cond(i) = (T(i)-T(i-1))/R; % heat transfer by conduction
    end
end

 dTdt = zeros(4,1);

 for i = 1:(m*n)
    if i == 1
        dTdt(i,1) = (Q_gen - Q_cond(i+1) - Q_cond(i+2))/((1/m)*mass*spec_heat);
    elseif i == 2
        dTdt(i,1) = (Q_cond(i-1) - Q_cond(i+2))/((1/m)*mass*spec_heat);
    elseif i == 3
        dTdt(i,1) = (Q_cond(i-2) - Q_cond(i+1))/((1/m)*mass*spec_heat);
    else
        dTdt(i,1) = (Q_cond(i-2) + Q_cond(i-1))/((1/m)*mass*spec_heat);
    end
 end
 end

[Youssef Khmou]

Hi , i tried your code, but it gives the same oscillatory Temp in the end, i think you should work on a 2d matrix , can you ?

[Youssef Khmou]

try this as preliminary test :

 surf(10*log10(abs(T)))
 >> shading interp

Answer 1

[Youssef Khmou]

hi, please, change the duration as :

 tspan=[0:4:7200];

and after running the file, type in the command prmpt :

 surf(T);

Does that mesh make sens as Decreasing Gradient of Temperature ??

[Mark McNallly]

Hi Youssef, thanks for the reply. What do you mean by I should work on a 2nd matrix though?

The solution should end up with the temperature reaching a steady state eventually. One element in the matrix is reducing in temperature using

 tspan=[0:4:7200];

Which doesn't really make sense, although I guess neither dies the oscillating temperatures that I currently have

[Mark McNallly]

Oh wait, I read your reply as 2nd, not 2-d. I will try it with a 2-d matrix.

[Youssef Khmou]

hi, Mark,

yes you try per example initializing a 2D T matrix , and put the equation you want to solve as Laplace(T) = constant ( just as example, you put those convection and conduction equations in general form ) , next it will be easy to write a function that differentiate the 2D T .