Solving the steady state 2D heat conduction problem

The purpose of this project is to implement point iteratives techniques such as Jacobi,Gauss Seidel and Successive Over-Relaxation to the 2D heat conduction equation respectively.

STEADY STATE SOLUTION

Derivation

When in the case of a steady state solution, convection, internal heat conduction and dissipation is neglected.

Hence, the energy equation in 2D becomes

`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

Discretizing this equation by using the central differencing scheme,we get

`(del^2T)/(delx^2)=(T(i-1,j) -2T(i,j)+T(i+1,j))/ (Deltax^2)`

`(del^2T)/(dely^2)=(T(i,j-1)-2T(i,j)+T(i,j+1))/(Deltay^2)`

As a result we get

`(T(i-1,j)-2T(i,j)+T(i+1,j))/(Deltax^2)+ (T(i,j-1)-2T(i,j)+T(i,j+1))/(Deltay^2)=0`

simplifying , we substitute the following values

T(i-1,j)=T_l

T(i+1,j)=T_r

T(i,j-1)=T_b

T(i,j+1)=T_t

T(i,j)=T_p=k

which leaves us with the following equation

`Tp=1/k((Tl+Tr)/(Deltax^2)+(Tb+Tt)/(Deltay^2))`

where

`k=2(1/(Deltax^2)+1/(Deltay^2))`

clear all
close all
clc
%initial conditions
L=1;
nx=10;
ny=nx;
% grid spacing along x and y axes
x=linspace(0,L,nx);
dx=L/(nx-1);
y=linspace(0,L,ny);
dy=L/(ny-1);
error=7e7;
tolerance=1e-4;
k=2*(1/(dx^2)+1/(dy^2));
% BCs
T=ones(nx,ny);

T(:,1)=400;
T(end,:)=600;
T(:,end)=800;
T(1,:)=900;
% temperature at corners by refining edges
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(600+400)/2;
T(end,end)=(600+800)/2;

for method=[1 2 3]

if method==1
    jacobi_steady(nx,ny,dx,dy,x,y,T,error,tolerance,k);
else if method==2
        gs_steady(nx,ny,dx,dy,x,y,T,error,tolerance,k);
    else if method==3
            sor_steady(nx,ny,dx,dy,x,y,T,error,tolerance,k);
        end
    end 
end

end
        
FUNCTION CODE FOR JACOBI METHOD
function out= jacobi_steady(nx,ny,dx,dy,x,y,T,error,tolerance,k)
[X,Y]=meshgrid(x,y);
T_old=T;
tic
jacobi_iter=1;

while(error>tolerance)
    for i=2:nx-1
        for j=2:ny-1
            term_1=(T_old(i-1,j)+T_old(i+1,j))/(k*(dx^2));
            term_2=(T_old(i,j-1)+T_old(i,j+1))/(k*(dy^2));
            
            T(i,j)= term_1+term_2;
            
        end
    end
    error=max(max(abs(T-T_old)));
    T_old=T;
    jacobi_iter=jacobi_iter+1;
end
iteration_time=toc;
figure(1)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel('x-axis')
ylabel('U')
colorbar
title_name=sprintf('the number of iterations for jacobi method is %d', jacobi_iter);
title(title_name)
file_text=sprintf('jacobi_steady.jpg')
saveas(gca,file_text);
fprintf('time taken for solving jacobi method is %f \n',iteration_time);
out=figure(1)
end

FUNCTION CODE FOR GAUSS SEIDEL METHOD

function out= gs_steady(nx,ny,dx,dy,x,y,T,error,tolerance,k)
[X,Y]=meshgrid(x,y);
T_old=T;
tic
gs_iter=1;

while(error>tolerance)
    for i=2:nx-1
        for j=2:ny-1
            term_1=(T(i-1,j)+T(i+1,j))/(k*(dx^2));
            term_2=(T(i,j-1)+T(i,j+1))/(k*(dy^2));
            T(i,j)=term_1+term_2;
        end
    end
    error=max(max(abs(T-T_old)));
    T_old=T;
    gs_iter=gs_iter+1;
end
iteration_time=toc;
figure(2)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel('x-axis')
ylabel('U')
colorbar
title_name=sprintf('the number of iterations for gauss seidel method is %d', gs_iter);
title(title_name)
file_text=sprintf('gs_steady.jpg')
saveas(gca,file_text);
fprintf('time taken for solving gauss seidel method is %f \n',iteration_time);
out=figure(2)
end

FUNCTION CODE FOR SOR METHOD

function out= sor_steady(nx,ny,dx,dy,x,y,T,error,tolerance,k)
lambda=1.5;
T_old=T;
[X,Y]=meshgrid(x,y);
tic
sor_iter=1;

while(error>tolerance)
    for i=2:nx-1
        for j=2:ny-1
         term_1=(T(i-1,j)+T_old(i+1,j))/(k*(dx^2));
         term_2=(T(i,j-1)+T_old(i,j+1))/(k*(dy^2));
         T(i,j)=(term_1+term_2)*lambda+(1-lambda)*T_old(i,j);
        end
    end
    error=max(max(abs(T-T_old)));
    T_old=T;
    sor_iter=sor_iter+1;
end
iteration_time=toc;
figure(3)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel('x-axis')
ylabel('U')
colorbar
title_name=sprintf('the number of iterations for sor method is %d', sor_iter);
title(title_name)
file_text=sprintf('sor_steady.jpg')
saveas(gca,file_text);
fprintf('time taken for solving sor method is %f \n',iteration_time);
out=figure(3);
end

The following are the plots generated

The time taken for solving the problem using Jacobi method is 0.003922 seconds

The time taken to solve the problem using Gauss Seidel method is 0.003802 seconds and

The time taken to solve the problem using SOR method is 0.003278 seconds.

• It is seen that by implementing the Gauss Seidel method, the number of iterations required to reach convergence is almost halved when compared to the Jacobi method(217).

By setting the relaxation factor to 1.5, convergence is reached much faster at 30 iterations,however, on decreasing the value to 1, the iterations are of the same number as in the Gauss Seidel method(116).