MATLAB: Bad result in 2D Transient Heat Conduction Problem Using BTCS Finite Difference Method implicitly

2d transient heat conductionimplicit

Everything seem's ok, but my solution's is wrong.
(Thank's to youtuber Sam R. & Bhartendu for Gauss-Seidel Method)
%% 2D HEAT EQUATION WITH CONSTANT TEMP. AT BC'S
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////




%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%% INITIALIZING
clc
clear 
close all
%% PARAMETERIZATION
a  = 1e-4;       % THERMAL DIFFUSIVITY
t  = 200;        % TOTAL TIME
nt = 2;          % TOTAL NUMBER OF TIME STEPS
dt = t/nt;       % TIMESTEP
L  = 1;          % X = Y
nx = 4;          % TOTAL NUMBER OF SPATIAL GRIDS
dx = L/nx;       % dX = dY
%% BC'S
T1 = 100;        % BCS OF DOMAIN
T2 = 100;
T3 = 100;
T4 = 100;
T5 = 200;         % INITIAL TEMP.
Tmax = max([T1,T2,T3,T4,T5]);
%% DEFINITION
m = (L/dx)+1;    % NO. OF X STEP (X = Y)
r = (t/dt)+1;    % NO. OF TIME STEP
d = a*dt/(dx)^2; % DIFF. NO.
%% CREATING BC's & IC'S
B = zeros(m-2,m-2);
[nr,nc] = size(B);
nT = nr * nc;      
B = zeros(m,m);
B(:,end) = T1;
B(end,:) = T2;
B(1,:) = T3;
B(:,1) = T4;
B(2:end-1,2:end-1) = T5; %% EDGES MUST BE REFINED
B1 = B(2:end-1,2:end-1)';
BC_v = B1(:);             %IC's
%% CREATE TDMA LHS
AA2(1:nT)   = 1+4*d;
AA3(1:nT-1) = -d;
AA5(1:nT-3) = -d;
AA = diag (AA2,0) + diag (AA3, -1) + diag (AA3,1) + ...
    diag(AA5,-3)+diag(AA5,3);   %% LHS HAS BEEN CREATEAD
%% CREATE RHS     
T13(:,:,:) = zeros(size(B,1),size(B,1),r); % Define a cubic matrix
n = size(AA,1);
T13(:,:,1) = B;
BB(n,1) = 0;
for i = 1:n
        if mod(i,2) == 1
            BB(i,1) = T1+T2;
        end
        if mod(i,2) ~= 1
            BB(i,1) = T1;
        end
          if mod(n,2)~=1 
              BB((n+1)/2,1) = 0;
          elseif mod(n,2)==1 
              BB((n+1)/2,1) = 0;
          end
end  
RHS = BC_v + d*BB;
for k = 1:r
%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////



%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////
C = RHS;     % constants vector
n = length(C);
X = zeros(n,1);
Error_eval = ones(n,1);
% Start the Iterative method
iteration = 0;
while max(Error_eval) > 1e-7
    iteration = iteration + 1;
    Z = X;              % save current values to calculate error later
    for i = 1:n
        j = 1:n;        % define an array of the coefficients' elements
        j(i) = [];      % eliminate the unknow's coefficient from the remaining coefficients
        Xtemp = X;      % copy the unknows to a new variable
        Xtemp(i) = [];  % eliminate the unknown under question from the set of values
        X(i) = (C(i) - sum(AA(i,j) * Xtemp)) / AA(i,i);
    end
    Error_eval = sqrt((X - Z).^2);
end
%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////
%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////
    RHS = X + d*BB;
    T13(2:end-1,2:end-1,k+1) = (vec2mat(RHS',m-2)); % VEC. TO MAT.
    T13(:,end,k+1) = T1;
    T13(end,:,k+1) = T2;
    T13 (1,:, k + 1) = T3;
    T13 (:, 1, k + 1) = T4;
end
        
%% Animation
x = 1:m;
y = 1:m;
T21 = T13 (:,:, 1);
 for k = 1:r-1
     figure(1)
     h = imagesc(x,y,T21);
     shading interp
     axis([0 m 0 m 0 Tmax])
     title({['Transient Heat Conduction'];['time = ',...
         num2str((k)*dt),' s']})
     colorbar;
     drawnow;
     xlabel(['T = ',num2str(T2),'C'])
     yyaxis left
     ylabel(['T = ',num2str(T4),'C'])
     yyaxis right
     str = 'T0 = 20C' ;
     dim = [.6 0.55 .3 .3];
     pause(0.1);
     refreshdata(h);
     if k~=r
         T21 = T13(:,:,k+1);
     else
         break;
     end
 end

Best Answer

%2-D Transient Heat Conduction With No Heat Generation [FDM][BTCS]
clc;
clear all;
%% Variable Declaration
n = 21;                 %number of nodes
L = 1;                  %length of domain
W = 1;                  %width of domain
alpha = 1e-4;           %thermal diffusivity (m^2/s)
m =(n-2)*(n-2);         %variable used to construct penta diagonal matrix
sm = sqrt(m);
dx = L/(n-1);           %delta x domain 
dy = W/(n-1);           %delta y domain
x = linspace(0,L,n);    %linearly spaced vectors x direction
y = linspace(0,W,n);    %linearly spaced vectors y direction
[X,Y]=meshgrid(x,y);
Tin = 200;              %internal temperature
T  = ones(m,1)* Tin;    %initilizing Space
Ta = zeros(n,n);
A  = zeros(m,m);         
Ax = zeros(m,1);
B  = zeros(sm,sm);
dt = 0.1;                %time step
tmax = 500;              %total Time steps (s)
t = 0 : dt : tmax;
r = alpha * dt /(dx^2);  %for stability, must be 0.5 or less
Tt = 100;                %Top Wall

Tb = 100;                %Bottom Wall

Tl = 100;                %Left Wall

Tr = 100;                %Right Wall

%% Boundry Conditions
Ta(1,1:n) = Tt;         %Top Wall
Ta(n,1:n) = Tb;         %Bottom Wall
Ta(1:n,1) = Tl;         %Left Wall
Ta(1:n,n) = Tr;         %Right Wall
%% Setup Matrix
for ix = 1 : 1 : m
    
    for jx =1 : 1 : m
        
        if (ix == jx)
                A(ix,jx) = (1 + 4*r); 
                
        elseif ( (ix == jx + 1) && ( (ix - 1) ~= sm * round( (ix-1)/sm) ) ) %RHS
                A(ix,jx) = -r;
                
        elseif ( (ix == jx - 1) && ( ix ~= sm*round(ix/sm) ) )
                A(ix,jx) = -r;
                
        elseif (ix == jx + sm)
                A(ix,jx) = -r;
                
        elseif (jx == ix + sm)
                A(ix,jx) = -r;
        else
                A(ix,jx) = 0;
        end
    end
end
 
for iy = 1 : 1 : sm
    
    for jy = 1 : 1 : sm
        
        if (iy == 1) && (jy == 1)
            B(iy,jy) = r * (Tl + Tt);
            
        elseif (iy == 1) && (jy == sm)
            B(iy,jy) = r * (Tt + Tr);                           %LHS
            
        elseif (iy == sm) && (jy == sm)
            B(iy,jy) = r * (Tb + Tr);
            
        elseif (iy == sm) && (jy == 1)
            B(iy,jy) = r * (Tb + Tl);
            
        elseif (iy == 1) && (jy == sm)
            B(iy,jy) = r * (Tt + Tr);
            
        elseif (iy == 1)&&(jy > 1 || jy < sm)
            B(iy,jy) = r * Tt;
            
        elseif (jy == sm) && (iy > 1 || iy < sm)
            B(iy,jy) = r * Tr;
            
        elseif (iy == sm) && ( jy > 1 || jy < sm)
            B(iy,jy) = r * Tb;
            
        elseif (jy == 1) && ( iy > 1 || jy < sm)
            B(iy,jy) = r * Tl;
        else
            B(iy,jy) = 0;
        end
                    
    end
end
 Bx = reshape(B,[],1);  %Convert matrix to vector
 
%% Solution
for l = 2 : length(t)   %time steps
    
     Xx = ( T + Bx );
     
           Ax = A \ Xx;
           
                 T( 1 : m ) = Ax( 1 : m );
                       fprintf('Time t=%d\n',l-1);
    
end
Tx = reshape( Ax , sm , sm); %convert vector to matrix
for i = 2 : 1 : n-1
    
    for j = 2 : 1 :n-1
        
        Ta(i,j) = Tx ( i-1 , j-1 );
        
    end
    
end
%% Plot
   contourf(X,Y,Ta,50,'edgecolor','none');
        h = colorbar;
        ylabel(h, 'Temperature °C')
        colormap jet
        axis equal
            
        title(['Top (Tt)= ',num2str(Tt),'°C']);
            xlabel(['Bottom (Tb)= ',num2str(Tb),'°C'])
            
            yyaxis left
            ylabel(['Left (Tl)= ',num2str(Tl),'°C'])
            
            yyaxis right
            ylabel(['Right (Tr)= ',num2str(Tr),'°C'])
Related Question