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

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

[Mohammad Farhadi]

PLEASE HELP ME :'(

[darova]

Can you attach original equations?

[Mohammad Farhadi]

Hi darova

Thanks for replay.

https://www.youtube.com/watch?v=z5tBlHeSpSI

After numerical discretization:

After rearranging:

I try to solve this matrix but my result is worng.

[darova]

It's too complicated for this forum. What about method of lines?

[Mohammad Farhadi]

sorry i am new in this forum.

This method use for solving Partial differential equations like heat equation. we consider a domain like this.

we open the equation in time step. < n > is time step.

This matrix solve iteratively over time.

Finally we have temperature in each time.

For solving this equation we need boundary and initial condition. I use "diag" function for creating left hand side. Then i change temperature matrix to a vector and solve matrix with Gauss Seidel Method for any time step.

Sorry for my bad english.

[Vineet Sengar]

Bro, could you send me your code? As I'm having the project related to 2-D transient heat conduction. Please do reply asap. Thanking u in the advance😊.

[Leroy Coelho]

I have written the code for this problem

I am more then happy to send the code your way for your refrence

[Ragul Kumar]

Hello experts,

I am trying to solve the finite difference methof for crank nicolson scheme to 2d heat equation. please let me know if you have any MATLAB CODE for this

Boundary condition are

If you can kindly send me the matlab code, it will be very useful for my research work . thank you very much.

Answer 1

[Leroy Coelho]

%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'])
        
        

[Mohammad Farhadi]

[Leroy Coelho]

HTH

Keep coding,keep learnning ,keep sharing

[Héctor Antonio González Pomares]

¿Cómo podría simular esto? Cuando el flujo es bidimensional de calor por conducción en condiciones transitorias para esta placa plana por medio de la aplicacióndel Método de Diferencias Finitas.

temperaturas en °C donde t1=250, t2=270, t3=230, t4=150, t5=210, t6=210, t7=270; t8=110; t9=300; t10=100; t11=290; t12=130 Cada cuadricula con dimensiones a x a donde a=0.25 m

[Soukaina Fekkar]

¿Puedes darme el código para este problema?@Héctor Antonio González Pomares