Solving the steady and unsteady 2D heat conduction problem

RESULT:

steady-state analysis 

The 2D heat conduction equation by using the point iterative techniques are implemented for  the following methods

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

where absolute error criteria are 1e-4

Derivation of the different Numerical scheme is done  before starting the program

Here results for iterative methods is the same but the no of iterations for all methods and the time is taken for convergence differs and is explained below. The result of 2D heat conduction is shown below in a contour plot below.

-----------------------------------------------------------------------------------------
Gaussian Siedel Method
-----------------------------------------------------------------------------------------
function [out] = gs();
Lx = 1;
Ly = Lx;
nx = 20;
ny = nx;

x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
[XX YY]=meshgrid(x,y);

% Initialize the temperature field

T = zeros(nx,ny);

% Boundary Conditions

T(1,:) = 900;       % Temperature at the top face
T(end,:) = 600;     % Temperature at the bottom face
T(:,1) = 400;       % Temperature at the left face
T(:,end) = 800;     % Temperature at the right face

error = 9e9;
toler = 1e-5;
Told = T;
T_initial = T;
k = 2*(dx^2 + dy^2)/(dx*dy)^2;

    gs_iterations =1 ;
    while (error > toler)
        for i=2:nx-1
            for j=2:ny-1
                T(i,j) = ((T(i,j-1) + T(i,j+1))*(1/dx^2) + (T(i-1,j) + T(i+1,j))*(1/dy^2))*(1/k);
            end
        end
        error_t = max(abs(Told - T));
        error = max(error_t);
        Told = T;
  gs_iterations = gs_iterations + 1;
  figure(1)
contourf(XX,YY,T)
title(\'Gaussian Siedel Method\')
colorbar

pause(0.003)
    end

out = gs_iterations;

end



----------------------------------------------------------------------------------------

 jacobi method
-----------------------------------------------------------------------------------------



function[out]= gs0;
 %Defining a function with time stepas argument 
clear all;
close all;
clc
Lx=1;
Ly=1;
nx = 10;
ny = 10;

x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
[XX YY]=meshgrid(x,y);

% Initialize the temperature field

T = zeros(nx,ny);

% Boundary Conditions
% Temperature at the top face
T(1,:) = 900;
 % Temperature at the bottom face
T(end,:) = 600;    
% Temperature at the left face
T(:,1) = 400;  
% Temperature at the right face
T(:,end) = 800;     

iterative_solver = 2;       
error = 9e9;
toler = 1e-4;
Told = T;
T_initial = T;
k = 2*(dx^2 + dy^2)/(dx*dy)^2;




    jacobi_iterations =1 ;
    while (error > toler)
        for i=2:nx-1
            for j=2:ny-1
               
                T(i,j) = ((Told(i,j-1) + Told(i,j+1))*(1/dx^2) + (Told(i-1,j) + Told(i+1,j))*(1/dy^2))*(1/k);
               
                end
        end
        error_t = max(abs(Told - T));
        error = max(error_t);
        Told = T;
        jacobi_iterations = jacobi_iterations + 1;
        figure(1)
contourf(XX,YY,T)

title(\'Jacobi Method\')
colorbar
 title(sprintf(\'Temprature distribution on a %d x %d  grid\\n using Jacobi iteration no %d\', nx, ny, jacobi_iterations),\'FontSize\',12);
 pause(0.003)
    end

out = jacobi_iterations;    

end
 
----------------------------------------------------------------------------------------






function [out] = sor()
Lx = 1;
Ly = 1;
nx = 20;
ny = nx;

x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
[XX YY]=meshgrid(x,y);

% Initialize the temperature field

T = zeros(nx,ny);

% Boundary Conditions

T(1,:) = 900;       % Temperature at the top face
T(end,:) = 600;     % Temperature at the bottom face
T(:,1) = 400;       % Temperature at the left face
T(:,end) = 800;     % Temperature at the right face

error = 9e9;
toler = 1e-4;
Told = T;
T_initial = T;
k = 2*(dx^2 + dy^2)/(dx*dy)^2;


    sor_iterations =1 ;
    while (error > toler)
        for i=2:nx-1
            for j=2:ny-1
              
                T(i,j) = ((Told(i,j-1) + Told(i,j+1))*(1/dx^2) + (Told(i-1,j) + Told(i+1,j))*(1/dy^2))*(1/k);
               
            end
 
end
        error_t = max(abs(Told - T));
        error = max(error_t);
        Tinitial = Told;
        Told = T;
        sor_iterations = sor_iterations + 1;
        figure(1)
contourf(XX,YY,T)
title(\'Successive Over Relaxation Method (SOR) Method\')
colorbar
pause(0.003)
    end

out = sor_iterations;

end
----------------------------------------------------------------------------------------

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ -------- The explicit method makes use of same time level and hence takes less no of iterations so solve, but takes more time. 

clear all
clc
%Initializing V23ariables:
Nx = 21;
Ny = Nx;
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 1e-4;
t_start = 0;
t_end = 1/alpha;
dt = 0.01;
error = 6e9;
tol = 1e-4;

% initilize
T = ones(Nx,Ny);
T(1,:) = 900; % Bottom Temperature
T(end,:) = 600; % Top Temperature
T(:,1) = 400; % Left Temperature
T(:,end) = 800; % Right Temperature


 %CFL Number calculation and adaptive Time-step control:
CFL_Number = alpha*dt/(dx^2);


%%% Calculation of Temperature Distribution explicitly using iterative solvers:

k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;
tic;
counter = 1;
for z = 1:500000
    %jacobi
    for i = 2:Nx-1
        for j = 2:Ny-1
            term1 = Told(i,j);
            term2 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
            term3 = k2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
            T(i,j) = term1+term2+term3;
        end
    end
    Told = T;
    
end
time_counter = toc;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf(\'2D steady state heat conduction (Explicit Method) \\n Number of iterations: %d (Time: %0.2fs) \\n Computation Time: %0.4f \\n CFL Number: %0.4f\',counter-1,z*dt,time_counter,CFL_Number));
xlabel(\'xx\');
ylabel(\'yy\');

implicit or explicit method schemes for the Transient Analysis are given above 

Here a unit Square plate is assumed for heat conduction under the given boundary conditions.we have also assumed here that nx=ny, which is grid point in x and y-direction 

The boundary conditions are given below for the plate

1.Top wall = 600k  

2.Right wall = 800k

3. left wall    =400k

4. Bottom wall = 900k

 At initially, the boundary condition of the plate is subjected at time t=0 And temperature of the plate is considered to be ambient temperature i.e T= 298k

from the above figures,

1. It can be seen that Gauss-Siedel Iteration takes less time when compared to Jacobi Iteration.

2. It can also be observed that the Implicit method takes lesser time when compared to the explicit method for solving the heat conduction. Also, in the implicit method, we require only the stable condition at the end of the output rather than the process happening in the intermediate stage. In the explicit method, it takes care of each and every step and hence the stability is a major concern in the case of the explicit method.