Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

Defining problem:

3D transient Heat conduction equation

`(deltaT)/(deltat) +alpha((delta^2T)/(deltax^2)+(delta^2T)/(deltay^2)+(delta^2T)/(deltaz^2))=0`

Where alpha is thermal diffusivity (m^2/s)

Given information

1. Top boundary = 600K

2. Bottom boundary = 900K

3. Left boundary = 400K

4. Right boundary = 800K

Following are the iterative techniques we are going to use -

1. Jacobi: The initial values at all nodes to calculate the new values, which will updated once we find the all the nodes values and this process will goes on untill we arrive with an converge solution. For different simulateneous equations with many unknowns first we have to make the equations into Diagonally dominant matrix form for converging results.

2. Gauss-seidel: In this iteration we assume the values obtained in Jacobi iteration for the next equation immediately without waiting for completion of first cycle.

3. Successive Over-relaxation (SOR): This method is an updated version of previous two. In addition we use the relaxation factor to ease the process. For this analysis we use the relaxation facor to be greater than 1 and less than 2.

2D transient state of equation

`(deltaT)/(deltat) +alpha((delta^2T)/(deltax^2)+(delta^2T)/(deltay^2))=0`

For steady state delta T is 0

`(delta^2T)/(deltax^2)+(delta^2T)/(deltay^2)=0`

For solving 2D steady state 2D above mentioned iterative techniques the steady state equation is discretised using 2nd order central differencing scheme

`(T_L-2T+T_R)/(Deltax^2)+(T_B-2T+T_T)/(Deltay^2)=0`

Assumption `k=2*(Deltax^2+Deltay^2)/(Deltax^2*Deltay^2)`

We can write

`T_P = 1/k((T_L+T_R)/(Deltax^2)+(T_B+T_T)/(DeltaY^2))`

all close 
all clear 
clc

%inputs
nx = 20; %number of grids in x direction 
ny = nx; %number of grids in y direction

tol = 1e-5; %tolerance 
error = 9e9; %intial guess 
T = 5; %time limit 
x = linspace(0,1,nx); %x axis mesh
y = linspace(0,1,ny); %y axis mesh

%grid spacing 
dx = x(2) -x(1);
dy = y(2) -y(1);

[xx,yy] = meshgrid(x,y);

%matrix definition 
t = 300*ones(nx,ny); %intial temp guess
t(end,:) = 900; %top boundary 
t(1,:) = 600; %bottom boundary
t(:,1) = 400; %left boundary
t(:,end) = 800; %right boundary

k = 2*(dx^2+dy^2)/(dx^2*dy^2);

t_old = t;
t_p = t;

%solving method
method = '\which method ?   \';
%\'1\' for jacobi
%\'2\' for gauss sidel
%\'3\' for SOR
solver = input(method); 

%Jacobi solver 

if solver ==1
    iter_1 = 1;
    tic;
    %convergence loop
    while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                TT = t_old(i,j+1);
                TB = t_old(i,j-1);
                
                H = (TL+TR)/(dx^2);
                V = (TB+TT)/(dy^2);
                
                t(i,j) = (H+V)/k;
            end 
        end 
        time = toc;
        error = max(max(abs(t_old-t)));
        t_old = t;
        iter_1 = iter_1 + 1;
        
        [x,y] = contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('\steady state using jacobi = %d\',iter_1);
        title_text2 = sprintf('\Simulation time +%f s\',time);
        title_text3 = sprintf('\Final error = %f\',error);
        title({title_text1;title_text2;title_text3})
        title(sprintf('\steady state using jacobi = %d\',iter_1))
    end 
end 

%GS solver  

if solver ==2
    iter_2 = 1;
    tic;
    %convergence loop
    while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                %TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                TT = t_old(i,j+1);
                %TB = t_old(i,j-1);
                
                H = (t(i-1,j)+TR)/(dx^2);
                V = (t(i,j-1)+TT)/(dy^2);
                
                t(i,j) = (H+V)/k;
            end 
        end 
        time2 = toc;
        error = max(max(abs(t_old-t)));
        t_old = t;
        iter_2 = iter_2 + 1;
        
        [x,y] = contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('\steady state using GS = %d\',iter_2);
        title_text2 = sprintf('\Simulation time +%f s\',time);
        title_text3 = sprintf('\Final error = %f\',error);
        title({title_text1;title_text2;title_text3})
    end 
end 

%SOR solver 
alpha = 0.3; 
omega = 1.2;
if solver == 3
    iter_3 = 1;
    tic;
    %convergence loop
    while(error>tol)
        for i = 2:(nx-1)
            for j = 2:(ny-1)
                %TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                TT = t_old(i,j+1);
                %TB = t_old(i,j-1);
                
                H = (t(i-1,j)+TR)/(dx^2);
                V = (t(i,j-1)+TT)/(dy^2);
                tgs(i,j) = (H+V)/k;
                
                t(i,j) = (1-omega)*(t_old(i,j))+omega*tgs(i,j);
            end 
        end 
        time3 = toc;
        error = max(max(abs(t_old-t)));
        t_old = t;
        iter_3 = iter_3 + 1;
        
        [x,y] = contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('\steady state using SOR = %d\',iter_1);
        title_text2 = sprintf('\Simulation time +%f s\',time);
        title_text3 = sprintf('\Final error = %f\',error);
        title({title_text1;title_text2;title_text3})
    end 
end 

The number of iteration required to converge the solution is more in Jacobi iteration scheme followed by Gauss Siedle iteraion and SOR scheme. Which means the computaional time required for a computer program to run the iteration is more in Jacobi as compare to GS and SOR scheme.