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

THE OBJECTIVE OF THE PROBLEM :

To solve the given 2D heat conduction problem for the steady and unsteady stade using different iterative solvers like Jacobi, Gauss-Seidel And Successive Over Relaxations iterative solvers and with the help of the given Boundary condition.

Given Data and Assumptions to be considered :

1. Area of computation: SQUARE

2. Number of grid points along the x-direction = Number of grid points along the y-direction

3. Initial Temperature of all nodes on the grid: 300k

4. Given Boundary Condition

• Temperature on Left boundary = 400k
• Temperature on Right boundary = 600k
• Temperature on Top boundary = 800k
• Temperature on Bottom boundary = 900k

Methodology :

The given 2D heat conduction equation is solved in two different states

1. Steady-state Conduction
2. Un-Steady State conduction

and also By using Two different Methods

1. Implicit
2. Explicit

Note: The above methods are used in Un-Steady State conduction

and also by using Three Iterative schemes

1. Jacobi Iterative Schemes
2. Gauss-Seidel Iterative Schemes
3. Successive Over Relaxations Iterative Schemes

Governing Equations and Derivations :

Steady-State Heat Conduction

Conduction is a type of heat Conduction that occurs in materials as a result of temperature differences caused by free electron and lattice vibration.

The heat equation is a partial differential equation that describes how the distribution of a quantity (such as heat) in a solid medium varies over time as it flows spontaneously from higher to lower points.

The General Heat Conduction Equation is given by

`frac(del^2T)(delx^2)+frac(del^2T)(dely^2)+frac(del^2T)(delz^2) + frac(q)(k) = frac(1)(alpha).frac(delT)(del tau)`

Where ,

T = Temperature

q = Heat Generation

`alpha` = Thermal diffusivity

Thermal diffusivity is a thermophysical property defined as the ratio of a material's thermal conductivity to its thermal capacity. Because the problem is 2D and no heat generation is taken into account, those variables are ignored. As a result, the equation depicts unstable 2D heat conduction (parabolic).

`(frac(del^2T)(delx^2)+frac(del^2T)(dely^2))alpha = frac(delT)(del tau)`

For the given Instinct, steady-state conduction temperature does not change with time, Therefore The time derivation of the above equation can be equated to zero

we get,

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

Since we are using the Central Difference method to solve the differential equation

`frac(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))(dx^2) + frac(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))(dy^2) = 0`

From the Given Data and Assumptions considered

Jacobi Iterative Method :

`T_(i,j) = 0.25*(Temp-old_(i+1,j)+Temp-old_(i-1,j)+Temp-old_(i,j+1)+Temp-old_(i,j-1))`

Gauss Seidel Iterative method :

`T_(i,j) = 0.25*(Temp-old_(i+1,j)+Temp_(i-1,j)+Temp-old_(i,j+1)+Temp_(i,j-1))`

Sucess Over Relaxation Iterative method :

`T_(i,j) = Omega*(0.25*(Temp-old_(i+1,j)+Temp_(i-1,j)+Temp-old_(i,j+1)+Temp_(i,j-1)) + (1-Omega)*(Temp-Old_(i,j))`

`Omega` = Relaxation Factor

MatLab Code: Solving Steady-State Heat Conduction

Main Code: From the below code we can call given three Iterative schemes for solving Steady-State Heat Conduction

clear all
close all
clc

%Number of Nodes
nx=20; %Number Of Nodes in x-Direction5
ny=20; %Number of Nodes In y- direction

x=linspace(0,1,nx);
y=linspace(0,1,ny);

dx=x(2)-x(1);
dy=y(2)-y(1);

error = 9e9; %Initial Error Value 
tolerance=1e-4; %Absolute Error Value

%Boundary Condition
Temp_Left= 400;  %Initial Condition Of Temperature at Left Boundary
Temp_Right=800; %Initial Condition Of Temperature at Right Boundary
Temp_Top=600; %Initial Condition Of Temperature at Top Boundary
Temp_Bottom=900; %Initial Condition Of Temperature at Bottom Boundary

Temp = 300*ones(nx,ny); %Initial Temperature of all the node points 

Temp(2:ny-1,1) =  Temp_Left; %Initial Temperature at Left Boundary Nodes 
Temp(2:ny-1,nx)= Temp_Right; %Initial Temperature at Right Boundary Nodes
Temp(1,2:nx-1) = Temp_Top; %Initial Temperature at Top Boundary Nodes
Temp(ny,2:nx-1) = Temp_Bottom; %Initial Temperature at Bottom Boundary Nodes

Temp(1,1) = (Temp_Top+Temp_Left)/2; %Calculating The temperature of the top-left node of the grid
Temp(1,ny) = (Temp_Top+Temp_Right)/2;%Calculating The temperature of the top-Right node of the grid
Temp(nx,1) = (Temp_Bottom+Temp_Left)/2; %Calculating The temperature of the Bottom-left node of the grid
Temp(nx,ny) = (Temp_Bottom+Temp_Right)/2; %Calculating The temperature of the Bottom-Left node of the grid

Temp_Initial = Temp;
Temp_Old = Temp;

%Enter the input for required Output ( 1 -jacobi_iteration ,
%2-gauss_seidal_iteration , 3 - Successive_over_relaxation_iteration)

iteration_solver = input('Enter The number : ')

if iteration_solver == 1 
    tic;
    Temp_jacobi = jacobi_iteration(error,tolerance,Temp,Temp_Old,x,y,nx,ny);
    Time_jacobi_iteration = toc;
end

if iteration_solver == 2 
    tic;
    Temp_gauss_seidel = gauss_seidal_iteration(error,tolerance,Temp,Temp_Old,x,y,nx,ny);
    Time_gauss_seidal_iteration = toc;
end

if iteration_solver == 3
    tic;
    Temp_Successive_over_relaxation = Successive_over_relaxation_iteration(error,tolerance,Temp,Temp_Old,x,y,nx,ny);
    Time_Successive_over_relaxation_iteration = toc;
end

Function codes :

Jacobi Iterative Method :

function Temp = jacobi_iteration(error,tolerance,Temp,Temp_Old,x,y,nx,ny)
tic % timer On
jacobi_itration=1;
while (error>tolerance)
    for i=2:nx-1
        for j =2:ny-1
            Temp(i,j) = 0.25*( Temp_Old(i-1,j) + Temp_Old(i+1,j) + Temp_Old(i,j-1) + Temp_Old(i,j+1));
        end
    end
    error = max(max(abs(Temp_Old-Temp)));
    Temp_Old = Temp;
    jacobi_itration = jacobi_itration+1;
end

time=toc %timer Off
figure(1)
contourf(x,y,Temp)
clabel(contourf(x,y,Temp))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
grid on
title(sprintf('Steady State Solution JACOBI ITERATION METHOD\n No.Of Iterations : %d Computational Time : %f',jacobi_itration,time))
end

OUTPUT :

Gauss Seidel Iterative method :

function Temp = gauss_seidal_iteration(error,tolerance,Temp,Temp_Old,x,y,nx,ny)
tic
gauss_seidal_iteration =1;
while (error>tolerance)
    for i=2:nx-1
        for j =2:ny-1
            Temp(i,j) = 0.25*(Temp(i-1,j) + Temp_Old(i+1,j) + Temp(i,j-1) + Temp_Old(i,j+1));
        end
    end
    error = max(max(abs(Temp_Old-Temp)));
    Temp_Old = Temp;
    gauss_seidal_iteration = gauss_seidal_iteration+1;
end

time=toc
figure(1)
contourf(x,y,Temp)
clabel(contourf(x,y,Temp))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X-Axis')
ylabel('Y-Axis')
grid on
title(sprintf('Steady State Solution GAUSS SEIDEL ITERATION METHOD\n No.Of Iterations : %d Computational Time : %f',gauss_seidal_iteration,time))
end

OUTPUT:

Sucess Over Relaxation Iterative method :

function Temp = Successive_over_relaxation_iteration(error,tolerance,Temp,Temp_Old,x,y,nx,ny)
tic
Successive_over_relaxatio_itration=1;
omega=1.2;
while (error>tolerance)
    for i=2:nx-1
        for j =2:ny-1
            Temp(i,j) = omega*(0.25*( Temp(i-1,j) + Temp_Old(i+1,j) + Temp(i,j-1) + Temp_Old(i,j+1))) + (1-omega)*Temp_Old(i,j);
        end
    end
    error = max(max(abs(Temp_Old-Temp)));
    Temp_Old = Temp;
    Successive_over_relaxatio_itration = Successive_over_relaxatio_itration+1;
end

time=toc
figure(1)
contourf(x,y,Temp)
clabel(contourf(x,y,Temp))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X-Axis')
ylabel('Y-Axis')
grid on
title(sprintf('Steady State Solution SOR ITERATION METHOD\n No.Of Iterations : %d Computational Time : %f',Successive_over_relaxatio_itration,time))
end

OUTPUT :

UNSTEADY /TRANSIENT HEAT CONDUCTION :

unsteady-state Heat transfer, as we all know, is the transfer of thermal energy from one body to another at a lower temperature. This thermal energy transfer might take place in either a steady or an unstable state. The temperature within the system does not change over time when it is in a steady state. Under unstable state conditions, on the other hand, the temperature within the system varies over time.
Stable state conditions are preceded by unsteady-state conditions. Under steady-state conditions, no system exists throughout time. After heat transmission begins, the system must establish a steady-state after some time has passed.

During this moment of transition, the system is in an unstable state.
Clearly, no system can exist in an unsettled state indefinitely. The system's temperature will eventually reach the heat source's temperature, and the system will be in a steady-state at that point. Even if the quantity of heat supplied into the system is increased, the system eventually reaches its critical temperature, at which point the energy transferred into the system begins to cause phase shifts within the system rather than temperature increases.
As a result, the 2D steady-state heat conduction equation is as follows:

`(frac(del^2T)(delx^2)+frac(del^2T)(dely^2))alpha = frac(delT)(del tau)`

When discretizing the above equation we get :

`frac(T_p^(n+1)-T_p^n)(Deltat) = alphafrac(T_L-2.T_P+T_R)(delx^2)+frac(T_T-2.T_P+T_B)(dely^2)`

`T_P^(n+1) = frac(1)(1+k1+k2).(T_P^n+k1(T_L+T_R)^(n+1)+k2(T_T+T_B)^(n+1))`

Where,

`k1= alpha.frac(Deltat)(Deltax^2)`

`k2 = alpha.frac(Deltat)(Deltay^2)`

Converting the above equation in grid formate we get

`T_(i,j) = frac(1)(1+k1+k2)(T_(i,j)+k1.(T_(i-1,j)+T_(i+1,j))+k2.(T_(i,j-1)+T_(i,j+1)))`

MatLab Code: Solving UNSTEADY /TRANSIENT HEAT CONDUCTION

Main Code: From the below code we can call given three Iterative schemes for solving UNSTEADY /TRANSIENT HEAT CONDUCTION

clear all
close all
clc

%Number of Nodes
nx=50;%Number Of Nodes in x-Direction5
ny=50;%Number of Nodes In y- direction
nt=2000;

x=linspace(0,1,nx);
y=linspace(0,1,ny);

dx=x(2)-x(1);
dy=y(2)-y(1);

Temp = 300*ones(nx,ny);%Initial Temperature of all the node points 
dt=0.1;
alpha =1.65634e-4;  %Thermal Diffusivity of Copper 
k1 = alpha*(dt/(dx^2));
k2 = alpha*(dt/(dy^2));
cfl = k1+k2;

error = 9e9;%Initial Error Value 
tolerance=1e-4;%Absolute Error Value


%Boundary Condition
Temp_Left= 400;%Initial Condition Of Temperature at Left Boundary
Temp_Right=800;%Initial Condition Of Temperature at Right Boundary
Temp_Top=600;%Initial Condition Of Temperature at Top Boundary
Temp_Bottom=900;%Initial Condition Of Temperature at Bottom Boundary



Temp(2:ny-1,1) =  Temp_Left;%Initial Temperature at Left Boundary Nodes 
Temp(2:ny-1,nx)= Temp_Right;%Initial Temperature at Right Boundary Nodes
Temp(1,2:nx-1) = Temp_Top;%Initial Temperature at Top Boundary Nodes
Temp(ny,2:nx-1) = Temp_Bottom;%Initial Temperature at Bottom Boundary Nodes

Temp(1,1) = (Temp_Top+Temp_Left)/2;%Calculating The temperature of the top-left node of the grid
Temp(1,ny) = (Temp_Top+Temp_Right)/2;%Calculating The temperature of the top-Right node of the grid
Temp(nx,1) = (Temp_Bottom+Temp_Left)/2;%Calculating The temperature of the Bottom-left node of the grid
Temp(nx,ny) = (Temp_Bottom+Temp_Right)/2;%Calculating The temperature of the Bottom-Left node of the grid

Temp_Old = Temp;
Temp_Initial = Temp_Old;


%Enter the input for required Output ( 1 -jacobi_iteration ,
%2-gauss_seidal_iteration , 3 - Successive_over_relaxation_iteration,4- Explicit Iteration)
iteration_solver = input('Enter The number : ')

if iteration_solver ==1
    tic;
    Temp_jacobi = jacobi_iteration_Unsteady_Implicit(tolerance,Temp,Temp_Old,x,y,nx,ny,k1,k2,Temp_Initial);
    Time_jacobi_iteration = toc;
end

if iteration_solver ==2 
    tic;
    Temp_gauss_seidel = Gauss_sidel_Unsteady_Implicit(tolerance,Temp,Temp_Old,x,y,nx,ny,k1,k2,Temp_Initial);
    Time_gauss_seidal_iteration = toc;
end

if iteration_solver == 3
    tic;
    Temp_Successive_over_relaxation = SOR_Unsteady_Implicit(tolerance,Temp,Temp_Old,x,y,nx,ny,k1,k2,Temp_Initial);
    Time_Successive_over_relaxation_iteration = toc;
end

if iteration_solver == 4
    tic
    Temp_Unsteady_explicit_Iteration = Unsteady_Explicit(Temp,Temp_Old,x,y,nx,ny,nt,k1,tolerance);
    Temp_Unsteady_explicit = toc;
end

Function codes :

Jacobi Iterative Method :

function Temp = jacobi_iteration_Unsteady_Implicit(tolerance,Temp,Temp_Old,x,y,nx,ny,k1,k2,Temp_Initial)
tic
jacobi_itration=1;
conv = 9e9;
time_itration = 1;
while(conv>tolerance)
    error = 9e9;
    while (error>tolerance)
        for i=2:nx-1
            for j =2:ny-1
                Temp_Term_1 = 1/(1 + (2*k1) + (2*k2));
                Temp_Term_2 = k1*Temp_Term_1;
                Temp_Term_3 = k2*Temp_Term_1;
                
                a = (Temp_Old(i-1,j))+(Temp_Old(i+1,j));
                b = (Temp_Old(i,j-1))+(Temp_Old(i,j+1));
                
                Temp(i,j) = (Temp_Initial(i,j)*Temp_Term_1) + (a*Temp_Term_2) + (b*Temp_Term_3);
            end
        end
        error = max(max(abs(Temp_Old-Temp)));
        Temp_Old = Temp;
        jacobi_itration = jacobi_itration+1;
    end
      conv = max(max(abs(Temp_Initial-Temp)));
      Temp_Initial = Temp;
      time_itration = time_itration+1
end

time=toc
figure(1)
contourf(x,y,Temp)
clabel(contourf(x,y,Temp))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
grid on
title(sprintf('Un-Steady State Solution JACOBI ITERATION METHOD\n No.Of Iterations : %d Computational Time : %f',jacobi_itration,time))
end

OUTPUT :

Gauss Seidel Iterative method :

function Temp = Gauss_sidel_Unsteady_Implicit(tolerance,Temp,Temp_Old,x,y,nx,ny,k1,k2,Temp_Initial)
tic
Gauss_sidel_Unsteady_itration=1;
conv = 9e9;
time_itration = 1;
while(conv>tolerance)
    error = 9e9;
    while (error>tolerance)
        for i=2:nx-1
            for j =2:ny-1
                Temp_Term_1 = 1/(1 + (2*k1) + (2*k2));
                Temp_Term_2 = k1*Temp_Term_1;
                Temp_Term_3 = k2*Temp_Term_1;
                
                a = (Temp(i-1,j))+(Temp_Old(i+1,j));
                b = (Temp(i,j-1))+(Temp_Old(i,j+1));
                
                Temp(i,j) = (Temp_Initial(i,j)*Temp_Term_1) + (a*Temp_Term_2) + (b*Temp_Term_3);
            end
        end
        error = max(max(abs(Temp_Old-Temp)));
        Temp_Old = Temp;
        Gauss_sidel_Unsteady_itration = Gauss_sidel_Unsteady_itration+1;
    end
      conv = max(max(abs(Temp_Initial-Temp)));
      Temp_Initial = Temp;
      time_itration = time_itration+1;
end

time=toc
figure(1)
contourf(x,y,Temp)
clabel(contourf(x,y,Temp))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
grid on
title(sprintf('Un-Steady State Solution GAUSS SIDEL ITERATION METHOD\n No.Of Iterations : %d Computational Time : %f',Gauss_sidel_Unsteady_itration,time))
end

OUTPUT :

Sucess Over Relaxation Iterative method :

function Temp = SOR_Unsteady_Implicit(tolerance,Temp,Temp_Old,x,y,nx,ny,k1,k2,Temp_Initial)
tic
SOR_Unsteady=1;
omega=1.1;
conv = 9e9;
time_itration = 1;
while(conv>tolerance)
    error = 9e9;
    while (error>tolerance)
        for i=2:nx-1
            for j =2:ny-1
                Temp_Term_1 = 1/(1 + (2*k1) + (2*k2));
                Temp_Term_2 = k1*Temp_Term_1;
                Temp_Term_3 = k2*Temp_Term_1;
                
                a = (Temp(i-1,j))+(Temp_Old(i+1,j));
                b = (Temp(i,j-1))+(Temp_Old(i,j+1));
                
                Temp(i,j) = (1-omega)*(Temp_Old(i,j)) + omega*((Temp_Initial(i,j)*Temp_Term_1)+(a*Temp_Term_2)+(b*Temp_Term_3));
            end
        end
        error = max(max(abs(Temp_Old-Temp)));
        Temp_Old = Temp;
        SOR_Unsteady = SOR_Unsteady+1;
    end
      conv = max(max(abs(Temp_Initial-Temp)));
      Temp_Initial = Temp;
      time_itration = time_itration+1;
end

time=toc
figure(1)
contourf(x,y,Temp)
clabel(contourf(x,y,Temp))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
grid on
title(sprintf('Un-Steady State Solution SOR ITERATION METHOD\n No.Of Iterations : %d Computational Time : %f',SOR_Unsteady,time))
end

OUTPUT :

Unsteady state Explicit Scheme :

function Temp = Unsteady_Explicit(Temp,Temp_Old,x,y,nx,ny,nt,k1,tolerance)
tic
Unsteady_Explicit=1;
for k = 1:nt
    error = 9e9
    while (error>tolerance)
        for i=2:nx-1
            for j =2:ny-1
                Temp(i,j) = Temp_Old(i,j)*(1-4*k1)+k1*(Temp_Old(i-1,j)+Temp_Old(i+1,j)+Temp_Old(i,j-1)+Temp_Old(i,j+1));
            end
        end
        error = max(max(abs(Temp_Old-Temp)));
        Temp_Old = Temp;
        Unsteady_Explicit = Unsteady_Explicit+1;
    end
end

time=toc
figure(1)
contourf(x,y,Temp)
clabel(contourf(x,y,Temp))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
grid on
title(sprintf('Un-Steady State Solution EXPLICIT ITERATION METHOD\n No.Of Iterations : %d Computational Time : %f',Unsteady_Explicit,time))
end

Results And Discussion

For Solving Steady State heat conduction equation :

Iterative Scheme	No Of Iteration	Computational Time
Jacobi Iterative Method	826	0.023334
Gauss Seidel Iterative method	442	0.009939
Sucess Over Relaxation Iterative method	304	0.007959

For Solving Unsteady State HEat conduction equation :

Iterative Scheme	No Of iterations	Computational Time
Jacobi Iterative Method	78343	5.651383
Gauss Seidel Iterative method	68917	3.639329
Sucess Over Relaxation Iterative method	60640	3.004065
Explicit Method	25228	1.079713

We can Conclud from the above outputs :

1. Jacobi > Gauss-seidel > Successive over-relaxation, based on the No. of iterations computed for all methods and Schemes.
2. When compared to Jacobi and Gauss-Siedel, SOR. The SOR iteration method is the most successful because it took the least number of iterations
3. When compared to steady-state equations, transient equations require more time to solve. As a result, solving the 2d heat conduction equation under steady-state circumstances is advised.
4. When solving a transient state equation, explicit methods are faster than implicit methods, although implicit methods are suggested when CFL is less than 0.5
5. When solving the transient state equation implicitly, Jacobi requires less time than the other two approaches, although SOR is more stable than gauss seidel.