1D transient heat conduction

Hello everybody, i am currently working on a simple modeling of a transient 1D heat conduction in a plate. I use the following script:

clc 
clear all; 
%-----------------------------Input---------------------------------------- 
% Inputs 
L=0.05;                          % wall thickness(m)
lambda=5.0;                           % conductivity(W/m-K)
rho=2000;                        % density(kg/m^3)                   
cp=200;                          % specific heat capacity(J/kg-K)

T_ini=293.2;                      % initial temperature(K) 
T_inf=273.2;                      % external temperature(K) 
alpha=500;                        % heat transfer coefficient(K) 

%-----------------------------Time stepping-------------------------------- 
% Setup time steps 
M=100;                            % number of time steps 
t=40; 
DELTA_t=t/(M);                 % time step duration(s) 
for j=1:M 
    time(j)=(j-1)*DELTA_t; 
end 

%-----------------------------Grid----------------------------------------- 
%Setup grid 
N=10;                              % number of nodes 1.Layer 

DELTA_x=L/(N);                   % distance between adjacent nodes 1.Layer(m) 

x=0:DELTA_x:L;                   % position of each node 1.Layer(m) 

%--------------------------Stability criterion----------------------------- 

DELTA_t_crit_N = DELTA_x*rho*cp/(2*(lambda/DELTA_x+alpha)); 
if DELTA_t > DELTA_t_crit_N 
    disp(' ') 
    disp('Time step exeeds the limit'); 
    return 
end 

%----------------------Initial wall temperatures--------------------------- 
%Initial wall temperatures T(i,1) 
for i=1:N+1   
   T(i,1)=T_ini;       
end 

% Step trough time 
for j=1:(M-1)     
    % Heat flux condition(q=n*(-k*dT/dx))[W/m^2] 
    T(1,j+1)=T(1,j)+2*lambda*(T(2,j)-T(1,j))*DELTA_t/(rho*cp*DELTA_x^2); 

      for i=2:(N)    
           T(i,j+1)=T(i,j)+lambda*(T(i-1,j)+T(i+1,j)-2*T(i,j))*DELTA_t/(rho*cp*DELTA_x^2); 
      end 

      % Heat flux condition(q=n*(-k*dT/dx))[W/m^2] + heat transfer coefficient(hout*(Tinf-T))[W/(m^2*K] 
      T(N+1,j+1)=T(N,j)+(2*lambda*(T(N-1,j)-T(N,j))/(rho*cp*DELTA_x^2)+2*alpha*(T_inf-T(N,j))/(rho*cp*DELTA_x))*DELTA_t; 

end 

%-----------------------------Plot----------------------------------------- 
plot(time,T) 
figure 
plot(x,T) 
figure 
g=T(:,M);   
plot(x,g)

Now the results look reasonible at first sight. However when i increase the number of time steps, the temperature difference between left and right side of the plate are getting lower and lower. With 10000 time steps i only get a difference of 0.2 K, while getting 9 K with 100 time steps. Where is my mistake?

I would be really gratful for your help. Thanks in advance,

Malte

[Juan Federico Herrera Ruiz]

The problem is that decreasing the time steps affects the Fourier number, making the solution less accurate

[Tyler Else]

Does anyone have an answer?

[三丰]

The heat transfer coefficient alpha=lambda/(rho*cp),why is it 500?

[Torsten]

@三丰

lambda/(rho*cp) is the thermal diffusivity, not the heat transfer coefficient.

The heat transfer coefficient influences the heat flux at the boundary via

-lambda*dT/dn = alpha*(T-Tinf)

where n is the outer normal.