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

Aim

1. Steady state analysis & Transient State Analysis

Solve the 2D heat conduction equation by using the point iterative techniques . The Boundary conditions for the problem are as follows;

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

You will implement the following methods for solving implicit equations.

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Your absolute error criteria is 1e-4

2. write code for both implicit and explicit schemes (This is for the Transient Analysis only)

2D heat conduction equation

(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )=1/α (∂T)/(∂t)

The heat conduction equation can be solved in two ways 

1. Steady state

2. Transient state

A)Solving 2D steady state heat conduction equation 

 Here the time term gets eliminated ie,    $\frac{\partial \mathrm{T/}}{\partial t}=0$

so the eqution becomes,(∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )=0

1. Jacobi method

In numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges.

2. Gauss siedel method

Gauss–Seidel method is an iterative method to solve a set of linear equations and very much similar to Jacobi’s method. This method is also known as Liebmann method or the method of successive displacement. The name successive displacement is because the second unknown is determined from the first unknown in the current iteration, the third unknown is determined from the first and second unknowns. 

3. Successive over relaxed method

The method of successive over relaxation (SOR) is a variant of the Gauss-siedel method for solving a linear system of equations, resulting in faster convergence. A similar method can be used for any slowly converging iterative process

 for converting into matlab code' I've changed temprature terms as follows:

$T_p$ = T(i,j)

$T_L$ = T(i-1,j)

$T_R$ = T(i+1,j)

$T_B$ = T(i,j-1)

$T_T$ = T(i,j+)

• For jacobi method all T on the Right hand side will be the old calculated values of T , Hence I'll be using Told
• For Gauss-siedel method only updated values to be substituted for T for the previous nodes but can use either T or Told for successive nodes 
• And for successive over relaxed method it will be of form

         T(i,j) = (1-w)*Told + w*(Tgs)

          Tgs is the value found through gauss siedel method

           here when omega, w >1 then it is known as over relaxed 

            and whe w<1 it will be under relaxed.

So the main code for solving steady state will be;

clear all
close all
clc

%inputs
L=1;
nx =20;
ny=20;
error = 8e8;
tol = 1e-4;

iteration_Method =2 ;

if iteration_Method == 1;
    steady_state_jacobian(L,nx,ny,error,tol);
end

if iteration_Method == 2;
    steady_state_GS(L,nx,ny,error,tol)
end

if iteration_Method ==3 ;
    steady_state_SOR(L,nx,ny,error,tol)
end
 

Code for steady state Jacobian 

%function for jacobian method
function steady_state_jacobian(L,nx,ny,error,tol)

tic

x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = L/(nx-1);
dy = L/(ny-1);
%Temperature matrix
T = 300*ones(nx,ny);

Tleft = 400;
Tright = 800;
Ttop = 600;
Tbottom = 900;


T(:,1)= Tleft;
T(1,:) = Ttop;
T(:,end) = Tright ;
T(end,:) = Tbottom;
%defining the edge values
T(1,1) = (Ttop + Tleft)/2;
T(end,end) = (Tright+Tbottom)/2;
T(1,end) = (Ttop + Tright)/2;
T(end,1) = (Tbottom+Tleft)/2;

k = 2*(dx^2+dy^2)/((dx^2)*(dy^2));
Told = T;
j_iter = 1


while (error>tol)
    for i = 2:(nx-1);
    for j = 2:(ny-1) ;
        T(i,j)=(1/k)*(((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j-1)+Told(i,j+1))/dy^2))
    end
    end
    
    error = max(max(abs(Told-T)));
    Told = T;  %replacing the old value to new
    j_iter = j_iter + 1; %for continueing the itration process
    %plotting
    figure(1)
    contourf(x,y,T) 
    
    colormap(jet);
    colorbar %for showing color scale
    toc;
    pause (0.03)
    titletext = sprintf('jacobi method with iteration number=%d',j_iter);
    title(titletext);
    xlabel('x')
    ylabel('y')
   
end

end

Code for steady state Gauss-siedel

%function for Gauss-siedel method
function steady_state_GS(L,nx,ny,error,tol)

tic

x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = L/(nx-1);
dy = L/(ny-1);
%Temperature matrix
T = 300*ones(nx,ny);

Tleft = 400;
Tright = 800;
Ttop = 600;
Tbottom = 900;


T(:,1)= Tleft;
T(1,:) = Ttop;
T(:,end) = Tright ;
T(end,:) = Tbottom;
%defining the edge values
T(1,1) = (Ttop + Tleft)/2;
T(end,end) = (Tright+Tbottom)/2;
T(1,end) = (Ttop + Tright)/2;
T(end,1) = (Tbottom+Tleft)/2;

k = 2*(dx^2+dy^2)/((dx^2)*(dy^2));
Told = T;
GS_iter = 1


while (error>tol)
    for i = 2:nx-1;
    for j = 2:ny-1 ;
        T(i,j)=(1/k)*(((T(i-1,j)+T(i+1,j))/dx^2)+((T(i,j-1)+T(i,j+1))/dy^2))
    end
    end
    
    error = max(max(abs(Told-T)));
    Told = T;  %replacing the old value to new
    GS_iter = GS_iter + 1; %for continueing the itration process
    %plotting
    
    contourf(x,y,T) 
    
    colormap(jet);
    colorbar %for showing color scale
    toc;
    pause (0.03)
    titletext = sprintf('Gauss siedel method with iteration number=%d',GS_iter);
    title(titletext);
    xlabel('x')
    ylabel('y')
end
end
 

Code for steady state SOR method

%function for sor method
function steady_state_SOR(L,nx,ny,error,tol)

tic

x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = L/(nx-1);
dy = L/(ny-1);
%Temperature matrix
T = 300*ones(nx,ny);

Tleft = 400;
Tright = 800;
Ttop = 600;
Tbottom = 900;


T(:,1)= Tleft;
T(1,:) = Ttop;
T(:,end) = Tright ;
T(end,:) = Tbottom;
%defining the edge values
T(1,1) = (Ttop + Tleft)/2;
T(end,end) = (Tright+Tbottom)/2;
T(1,end) = (Ttop + Tright)/2;
T(end,1) = (Tbottom+Tleft)/2;

k = 2*(dx^2+dy^2)/((dx^2)*(dy^2));
Told = T;
w = 1.7;
SOR_iter = 1


while (error>tol)
    for i = 2:nx-1;
    for j = 2:ny-1 ;
        Tgs(i,j)=(1/k)*(((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2));
        T(i,j) = (1-w)*(Told(i,j))+w*(Tgs(i,j));
    end
    end
    
    error = max(max(abs(Told-T)));
    Told = T;  %replacing the old value to new
    SOR_iter = SOR_iter + 1; %for continueing the itration process
    %plotting
    figure(3)
    contourf(x,y,T) 
    
    colormap(jet);
    colorbar %for showing color scale
    toc;
    pause(0.03)
    titletext = sprintf('SOR method with iteration number=%d',SOR_iter );
    
    title(titletext)
    xlabel('x')
    ylabel('y')
end
end

Results obtained

1. Jacobian method

Elapsed time is 490.967885 seconds.

2. Gauss-siedel method

Elapsed time is 295.191002 seconds.

3. Successive over relaxed method

Elapsed time is 14.475782 seconds.

Explanation

As you can see from the above results Jacobi method has the highest number of iteration with an elapsed time of 490.96 seconds. Gauss-siedel method has less number of iterations than jacobi method .Also it can be seen that SOR method is most efficient since it hace less number of iteration and the elapsed time is very less. The SOR method really depend on the relaxation factor which I used here as omega (w) . When we use the values more or less of 1.7 for omega calculated value will be away from actual value, so iteration gets incresed . so the preffered values are 1.7 and values around it .

B) Transient state explicit method

%function for solving transient state heat equation explicitly
function transient_state(L,nx,ny,error,tol)

tic

nx= 20;
L=1;
ny =nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = L/(nx-1);
dy = L/(ny-1);
%Temperature matrix
T = 300*ones(nx,ny);
t =0.14;
dt = 1e-4;
nt =t/dt;

Tleft = 400;
Tright = 800;
Ttop = 600;
Tbottom = 900;


T(:,1)= Tleft;
T(1,:) = Ttop;
T(:,end) = Tright ;
T(end,:) = Tbottom;
%defining the edge values
T(1,1) = (Ttop + Tleft)/2;
T(end,end) = (Tright+Tbottom)/2;
T(1,end) = (Ttop + Tright)/2;
T(end,1) = (Tbottom+Tleft)/2;
alpha =1;
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);

Tpre =T;



for k = 1:nt
    for i = 2:(nx-1);
    for j = 2:(ny-1) ;
        T(i,j)=Tpre(i,j)*(1-2*k1-2*k2)+k1*(Tpre(i-1,j)+Tpre(i+1,j))+k2*(Tpre(i,j-1)+Tpre(i,j+1));
    end
    end
    
   
   
    Tpre = T
    %plotting
    figure(1)
    contourf(x,y,T) 
    
    colormap(jet);
    colorbar %for showing color scale
    toc;
    pause (0.03)
    
    
    xlabel('x')
    ylabel('y')
   
end

end
 

Result obtained

Elapsed time is 374.642637 seconds.