MATLAB: Problems with 1D heat diffusion with the Crank Nicholson method

1d heat diffusioncrank nicholson method

Hi everyone, I am trying to implement the code for 1D heat diffusion with the Crank Nicholson method. but it turned out to an error at the temperature matrix. could you please help me with the code? I would really appreciate it.
% Heat Diffusion in 1D within the Crank-Nicholson Method
%du/dt=cond*d2u/d2t
%dU/dt is the net change in the total energy of the system
%the flux of heat is proportional to the gradient of the temperature
clear all;
clc;
% Parameters to define space and time
L = 0.05; % Lenght/Thickness [m]
Tmax =900.; % Final time [s]
% Parameters needed to solve the equation
Nt = 1000; % Number of time steps
t = Tmax/Nt; %time differential
N = 25.; % Number of space steps
h = L/N; %space differential
cond = 0.4; % Conductivity
rho= 900; % density
ce= 800; %specific heat
b = cond*t/(rho*ce*h*h); % Parameter of the method 
%solution "u" is temperature , 
%that is why I use specific heat (ce) and density (rho)
% Initial temperature 
Tamb= 25; % constant ambient temperature
for i0 = 1:N+1
      dx(i0) =(i0-1)*h; %position vector (used in graphics)
      u(i0,1) =Tamb; % initial ambient temperature 
end
% BC: Temperature at the boundary
T1=25; %inital temperature for boundary condition
T2=35; %final temperature for boundary condition
T0= linspace(T1,T2,Nt+1); 
%linearly gradient between 2 known temperatures in x=0
TL= linspace(T1,T2,Nt+1);
%linearly gradient between 2 known temperatures in x=L
for j0=1:Nt+1
      u(1,j0) = T0(j0); %BC in the left extreme x=0
      u(N+1,j0) = TL(j0); %BC in the right extreme x=L
      dt(j0) = (j0-1)*t; %time vector (used in graphics)
end
% Defining the Matrices M_right and M_left in the CN method
aal(1:N-2)=-b;
bbl(1:N-1)=2.+2.*b;
ccl(1:N-2)=-b;
MMl=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);
aar(1:N-2)=b;
bbr(1:N-1)=2.-2.*b;
ccr(1:N-2)=b;
MMr=diag(bbr,0)+diag(aar,-1)+diag(ccr,1);
% Implementation of the Crank-Nicholson method
for j0=2:Nt+1             % Time Loop
      uu=u(2:N,j0-1);
      u(2:N,j0)=inv(MMl)*MMr*uu;
end
The problem is that the solution vector "u" seems to end with an error (illogical solution), or divided by 0.
Thank you in advanced for the help.

Best Answer

rho and ce are undefined. Thus b is undefined.
Best wishes
Torsten.
Related Question