Numerical Simulation of Steady and Transient 2D heat conduction problem using different iterative solvers.

Abstract

The given problem demonstrates the flow of heat along a unit square domain at a defined length. Using both explicit and implicit scheme, the steady state and transient state conduction is computed. Applying Jacobi, gauss seidel and SOR iterative solver with a defined CFL value, a comparison study is done to find out the optimal iterative solver in terms of no of iterations performed and computational time.

Governing equations:

The general formula for 2D heat conduction is given by

The first part of (1) represents transient or unsteady state and the second part represents steady state of heat conduction.

Now considering the steady state equation let us discretize using numerical schemes

From (1) we can represent the steady state as

Discretizing the above equation, we get

Denoting

Rearranging the equation, we can write as

or

Finally, the equation can be represented as

Using implicit solvers, we can write the equation with respect to the number of iterations m.

1.Jacobi solver

2.Gauss-Siedel solver

3.Succesive over relaxation solver

Where w is the relaxation factor taken as 1.1

Now considering the unsteady state equation let us discretize using two numerical schemes:

1.Explicit schemes

Discretizing the equation

we get

Where n denote the time step

Rearranging it we can write as

Here k1 and k2 represents CFL no or diffusion probability

If CFL no>0.5, the solution becomes unstable.

If CFL<0.5, the solution is stable.

2.Implicit schemes

Discretizing the equation for (n+1) time steps

Using iterative solvers on (3) with respect to iteration no m.

1.Jacobi solver

2.Gauss-Siedel solver

3.Succesive over relaxation solver

Objective:

To observe the computational time, no of iterations and thermal gradient effect by different iterative solvers both for steady and transient case.

Initial conditions and boundary conditions:

1. Assume the domain to be a unit square.
2. Assume nx=ny [ no of grid points in x axis is equal to y axis].
3. Boundary conditions for steady and transient case

Left 400k

Right 800k

Top 600k

Bottom 900k

Workflow of the problem:

Steady state conduction

Transient state conduction

Algorithm

Algorithm for steady state conduction

1. Define the no of grid points in x and y direction.
2. Define the temperature domain matrix by applying boundary conditions.
3. Perform implicit iterations by Jacobi, gauss seidel and SOR solvers by function call method.
4. In each function call it takes temperature and nodal points as input
5. Inside the iterative solver function, length of the domain, Told, mesh grid and value of k is defined.
6. Starting with convergence criteria, two nodal loops in x and y axis is executed.
7. With each space marching the temperature value at the next nodal point is computed.
8. After each iteration, the error is computed and the old values of temperature are updated.
9. When the convergence criteria fail, the loop terminates out and the function sends back the new updated values of temperature, computational time and the no of iterations performed to the main program.
10. In the main program the temperature plots are plotted in a contour fashion.
11. Finally, a bar plot of simulation time and no of iterations performed for different iterative solvers are plotted to compare the results.

Algorithm for transient state conduction explicit solvers

1. Define the no of grid points in x and y direction.
2. Define the temperature domain matrix by applying boundary conditions.
3. Perform explicit iterations by function call method.
4. In each function call it takes temperature and nodal points as input
5. Inside the explicit solver function, length of the domain, Told, mesh grid and value of CFL no is defined.
6. Starting with time marching step, two for loops for inner nodal points at x and y axis are defined.
7. After nodal iteration the old temperature value are updated.
8. After execution of the time loop the final updated temperature, no of iteration performed and computational time are returned back to the main program.
9. Temperature plots in contour fashion are plotted.
10. Finally, a bar plot of simulation time and no of are plotted to compare the results.

Algorithm for transient state conduction implicit solvers

1. Define the no of grid points in x and y direction.
2. Define the temperature domain matrix by applying boundary conditions.
3. Perform implicit iterations by function call method.
4. In each function call it takes temperature and nodal points as input
5. Inside the implicit solver function, length of the domain, Told, mesh grid and value of CFL no is defined.
6. Starting with time marching step, the convergence criteria is applied. Following the convergence loop, two for loops for the inner nodes in x and y direction are executed.
7. After nodal iteration, the new temperature value gets updated.
8. After the convergence criteria fails the temperature value at the previous timestep gets updated.
9. With every increment of time marching step, the convergence loop is executed.
10. After the execution of the time steps the final updated temperature, computational time and the no of iteration performed are returned back to the main program.
11. Finally, a bar plot of simulation time and no of are plotted to compare the results.

Algorithm for computing the optimal value of omega for SOR solvers

1. Define the range of omega.
2. Perform a for loop iteration for the entire length of omega
3. Pass the value of omega at each index position to the SOR solver function. The function returns the no of iterations performed for each value of omega.
4. Initialize a for loop for the entire length of iteration variable. Extract the minimum iteration value using min function.
5. Use find function to extract the nodal position at which the minimum no of iteration takes place.
6. At the optimal position, the optimal omega value is extracted.

MATLAB CODE

Steady state solver code

clear all
clc
close all

%%% define the number of grid points required in x and y direction
nx=10;%no of grids in x axis
ny=nx;%no of grids in y axis

%define the temperature domain matrix
T=ones(nx,ny)*300;
%boundary conditions of temperature domain
TL=400;
TR=800;
TT=600;
TB=900;
%defining the temperatures at the boundaries
T(2:ny-1,1)=TL;% temperature at the left boundary
T(2:ny-1,nx)=TR;% temperature at the right boundary
T(1,2:nx-1)=TT;% temperature at the top boundary
T(nx,2:ny-1)=TB;% temperature at the bottom boundary
%define the temperature at the corner points
T(1,1)=(TT+TL)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% inputs for the three function call are [T,nx,ny]
%%%%%%%%%%% function call for Jacobi iteration solver%%%%%%%%%%%
[jacobi_computation_time, TOLD, X, Y,jacobi_iteration]=steady_state_jacobi_solver(T,nx,ny);
figure(1)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
titletext1=sprintf('Thermal gradient plot of steady state heat conduction by Jacobi iterative solver\n');
titletext2=sprintf('Jacobi iteration no=%d\n',jacobi_iteration);
titletext3=sprintf('Jacobi computation time=%d',jacobi_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%% function call for Gauss-Siedel iteration solver%%%%%%%%%%%
[gs_computation_time, TOLD, X, Y,gs_iteration]=steady_state_gauss_siedel_solver(T,nx,ny);
figure(2)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
titletext1=sprintf('Thermal gradient plot of steady state heat conduction by Gauss-Siedel iterative solver\n');
titletext2=sprintf('Gauss-Siedel iteration no=%d\n',gs_iteration);
titletext3=sprintf('Gauss-Siedel computation time=%d',gs_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%% function call for Successive Over Relaxation iteration solver%%%%%%%%%%%
[sor_computation_time, TOLD, X, Y,sor_iteration]=successive_over_relaxation(T,nx,ny);
figure(3)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
titletext1=sprintf('Thermal gradient plot of steady state heat conduction by successive over relaxation iterative solver\n');
titletext2=sprintf('successive over relaxation iteration no=%d\n',sor_iteration);
titletext3=sprintf('successive over relaxation computation time=%d',sor_computation_time);
title([titletext1 titletext2 titletext3]);

function call for steady state implicit time integration jacobi solver

function [time1, TOLD, X, Y,jacobi_iteration]=steady_state_jacobi_solver(T,nx,ny)
 l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
error=9e9;
tolerance=1e-4;
%define the constant k
k=2*(dx^2+dy^2)/(dx^2*dy^2);
[X,Y]= meshgrid(x,y);% generating mesh grid
TOLD=T;
jacobi_iteration=1;
   tic
   while(error>tolerance)
        for i=2:nx-1%space marching step for inner loop in x direction
            for j=2:ny-1%space marching step for inner loop in y direction
 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;
       jacobi_iteration=jacobi_iteration+1;% computes the no of iteration performed
   end
time1=toc;% computes the time to run the iteration
end

function call for steady state implicit time integration gauss siedel solver

function [time2, TOLD, X, Y,gs_iteration]=steady_state_gauss_siedel_solver(T,nx,ny)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
%define the constant k
k=2*(dx^2+dy^2)/(dx^2*dy^2);
[X,Y]= meshgrid(x,y);% generating mesh grid
TOLD=T;
error=9e9;
gs_iteration=1;
   tic
   while(error>tolerance)
        for i=2:nx-1%space marching step for inner loop in x direction
            for j=2:ny-1%space marching step for inner loop in y direction
 T(i,j)=(1/k)*[((T(i-1,j)+TOLD(i+1,j))/(dx^2))+((TOLD(i,j+1)+T(i,j-1))/(dy^2))];
            end
        end
        error=max(max(abs(TOLD-T)));
        TOLD=T;
       gs_iteration=gs_iteration+1;% computes the no of iteration performed
   end
    time2=toc;% computes the time to run the iteration
end

function call for steady state implicit time integration SOR solver

function [time4, TOLD, X, Y,sor_iteration]=successive_over_relaxation(T,nx,ny)
l=1;%length of the domain
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction

tolerance=1e-4;
%define the constant k
k=2*(dx^2+dy^2)/(dx^2*dy^2);
[X,Y]= meshgrid(x,y);% mesh grid generation
TOLD=T;
error=9e9;
w=1.1;% relaxation factor
    sor_iteration=1;
   tic
   while(error>tolerance)
        for i=2:nx-1%space marching step for inner loop in x direction
            for j=2:ny-1%space marching step for inner loop in y direction
 T(i,j)=(1-w)*TOLD(i,j)+w*[(1/k)*[((T(i-1,j)+TOLD(i+1,j))/(dx^2))+((TOLD(i,j+1)+T(i,j-1))/(dy^2))]];
            end
        end
        error=max(max(abs(TOLD-T)));
        TOLD=T;
       sor_iteration=sor_iteration+1;% computes the no of iteration performed
   end
    time4=toc;% computes the time to run the iteration
end

Transient state solver code

clear all
clc
close all
%%% define the number of grid points required in x and y direction
nx=10;%no of grids in x axis
ny=nx;%no of grids in y axis
alpha=0.1;%thermal diffusivity factor
dt=1e-2;
nt=1400;%no of timesteps
%define the temperature domain matrix
T=ones(nx,ny)*300;
%boundary conditions of temperature domain
TL=400;
TR=800;
TT=600;
TB=900;
%defining the temperatures at the boundaries
T(2:ny-1,1)=TL;% temperature at the left boundary
T(2:ny-1,nx)=TR;% temperature at the right boundary
T(1,2:nx-1)=TT;% temperature at the top boundary
T(nx,2:ny-1)=TB;% temperature at the bottom boundary
%define the temperature at the corner points
T(1,1)=(TT+TL)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% inputs for the function call are [T,nx,ny,alpha]

%%%%%%%%%%%%%%%%%%%%%%%%%EXPLICIT SOLVER%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% function call for explicit iteration solver%%%%%%%%%%%
[computation_time, TOLD, X, Y,explicit_iteration,CFL_explicit]=transient_conduction_explicit_solver(T,nx,ny,alpha,dt,nt);
figure(1)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Explicit iterative solver\n');
titletext2=sprintf('Explicit iteration no=%d\n',explicit_iteration);
titletext3=sprintf('Explicit solver computation time=%d',computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%IMPLICIT SOLVER%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% function call for implicit Jacobi iteration solver%%%%%%%%%%%
[jacobi_computation_time, TOLD, X, Y,jacobi_iteration,CFL_jacobi]=transient_conduction_implicit_jacobi_solver(T,nx,ny,alpha,dt,nt);

figure(2)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Jacobi solver\n');
titletext2=sprintf('Jacobi iteration no=%d\n',jacobi_iteration);
titletext3=sprintf('Jacobi solver computation time=%d',jacobi_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%% function call for implicit Gauss-Siedel iteration solver%%%%%%%%%%%
[gs_computation_time, TOLD, X, Y,gs_iteration,CFL_gs]=transient_conduction_implicit_gauss_siedel_solver(T,nx,ny,alpha,dt,nt);
figure(3)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Gauss-Siedel solver\n');
titletext2=sprintf('Gauss-Siedel iteration no=%d\n',gs_iteration);
titletext3=sprintf('Gauss-Siedelsolver computation time=%d',gs_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%SUCCESSIVE OVER RELAXATION%%%%%%%%%%%%%
%%%%%%%%%%% function call for SOR Jacobi iteration solver%%%%%%%%%%%
[sor_jacobi_computation_time, TOLD, X, Y,sor_jacobi_iteration,CFL_jacobi_sor]=transient_conduction_implicit_sor_jacobi_solver(T,nx,ny,alpha,dt,nt);
figure(4)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Jacobi SOR solver\n');
titletext2=sprintf('SOR jacobi iteration no=%d\n',sor_jacobi_iteration);
titletext3=sprintf('SOR jacobi solver computation time=%d',sor_jacobi_computation_time);
title([titletext1 titletext2 titletext3]);

%%%%%%%%%%% function call for SOR Gauss-Siedel iteration solver%%%%%%%%%%%
[sor_gs_computation_time, TOLD, X, Y,sor_gs_iteration,CFL_gs_sor]=transient_conduction_implicit_sor_gauss_siedel_solver(T,nx,ny,alpha,dt,nt);
figure(5)
[M,N]=contourf(X,Y,TOLD);
set(gca,'YDir','reverse')
clabel(M,N)
colormap(jet)
c=colorbar;
c.Label.String='temperature range';
c.Label.FontSize=10;
c.Label.FontWeight='bold';
xlabel('length of domain in x axis')
ylabel('length of domain in y axis')
titletext1=sprintf('Thermal gradient plot of transient heat conduction by Implicit Gauss Siedel SOR solver\n');
titletext2=sprintf('SOR gauss siedel iteration no=%d\n',sor_gs_iteration);
titletext3=sprintf('SOR gauss siedel solver computation time=%d',sor_gs_computation_time);
title([titletext1 titletext2 titletext3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%plot of simulation time and no of iterations%%%%%%%%%%%%%%
iterations=[explicit_iteration jacobi_iteration gs_iteration sor_jacobi_iteration sor_gs_iteration];
simulationtime=[computation_time jacobi_computation_time gs_computation_time sor_jacobi_computation_time sor_gs_computation_time];
figure(6)
c=categorical({'explicit','jacobi','Gauss Siedel','SOR jacobi','SOR Gauss Siedel'});
bar(c,iterations,'facecolor',[1 0.5 0.5])
ylabel('No of iterations')
figure(7)
c=categorical({'explicit','jacobi','Gauss Siedel','SOR jacobi','SOR Gauss Siedel'});
bar(c,simulationtime,'facecolor',[1 0.1 0.8])
ylabel('simulation time')

function call for transient state explicit time integration solver

function [time, TOLD, X, Y,iteration,k]=transient_conduction_explicit_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain

x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);% generating meshgrid
TOLD=T;
error=9e9;

%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
iteration=0;
   tic
for n=1:nt  % time marching loop
  
for i=2:nx-1%space marching loop at inner nodes of x direction
    for j=2:ny-1%space marching loop at inner nodes of y direction
        T(i,j)=TOLD(i,j)+k1*(TOLD(i+1,j)-2*TOLD(i,j)+TOLD(i-1,j))+k2*(TOLD(i,j-1)-2*TOLD(i,j)+TOLD(i,j+1));
    end
end
TOLD=T;% updates the old values of temperature
iteration=iteration+1;% computes the no of iteration performed
end
  
time=toc;% computes the time to run the simulation
end

function call for transient state implicit time integration jacobi solver

function [time, TOLD, X, Y,jacobi_iteration,k]=transient_conduction_implicit_jacobi_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain

x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction 
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating meshgrid
TOLD=T;


%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;

jacobi_iteration=1;
tic
 for n=1:nt % time marching loop
     error=9e9;
   while(error>tolerance)
        for i=2:nx-1%space marching loop at inner nodes of x direction
            for j=2:ny-1%space marching loop at inner nodes of y direction
 T(i,j)=term1*TPREV(i,j)+term2*(TOLD(i+1,j)+TOLD(i-1,j))+term3*(TOLD(i,j-1)+TOLD(i,j+1));
            end
        end
        error=max(max(abs(TOLD-T)));% computes the error at each iteration
        TOLD=T;% updates the value of temperature
       jacobi_iteration=jacobi_iteration+1;% computes the no of iteration performed
   end
       TPREV=T;% stores the value of temperature at the previous time steps
 end
time=toc;% computes the time to run the simulation
end

function call for transient state implicit time integration gauss siedel solver

function [time, TOLD, X, Y,gs_iteration,k]=transient_conduction_implicit_gauss_siedel_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain

x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);% generating meshgrid
TOLD=T;


%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%define the constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;

gs_iteration=1;
tic
for n=1:nt  % time marching loop
     error=9e9;
   while(error>tolerance)%convergence loop
      
        for i=2:nx-1%space marching loop at inner nodes of x direction
            for j=2:ny-1%space marching loop at inner nodes of y direction
T(i,j)=term1*TPREV(i,j)+term2*(TOLD(i+1,j)+T(i-1,j))+term3*(T(i,j-1)+TOLD(i,j+1));
            end
        end
        error=max(max(abs(TOLD-T)));% computes the error at each iteration
        TOLD=T;% updates the old values of temperature
       gs_iteration=gs_iteration+1;% computes the no of iteration performed
   end
     TPREV=T;% stores the value of temperature at the previous time steps  
   end
  time=toc; % computes the time to run the simulation
end                                              
                
                       

function call for transient state implicit time integration SOR jacobi solver

function [time, TOLD, X, Y,sor_iteration,k]=transient_conduction_implicit_sor_jacobi_solver(T,nx,ny,alpha,dt,nt)
l=1;%length of the domain

x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;


%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;

sor_iteration=1;
w=1.1;% relaxation factor
tic
for n=1:nt  % time marching loop
    error=9e9;
   while(error>tolerance)% convergence loop
       
        for i=2:nx-1%space marching loop at inner nodes of x direction
            for j=2:ny-1%space marching loop at inner nodes of y direction
                
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+TOLD(i-1,j))+term3*(TOLD(i,j-1)+TOLD(i,j+1)));
            end
        end
        error=max(max(abs(TOLD-T)));% computes the error at each iteration
        TOLD=T;% updates the value of temperature
       sor_iteration=sor_iteration+1;% computes the no of iteration performed
   end
       TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end

function call for transient state implicit time integration SOR gauss siedel solver

function [time, TOLD, X, Y,sor_iteration,k]=transient_conduction_implicit_sor_gauss_siedel_solver(T,nx,ny,alpha,dt,nt)

l=1;%length of the domain

x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction 
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;


%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;

sor_iteration=1;
w =1.1;% relaxation factor
tic
for n=1:nt  % time marching loop
    error=9e9;
   while(error>tolerance)% convergence loop
       
        for i=2:nx-1%space marching loop at inner nodes of x direction
            for j=2:ny-1%space marching loop at inner nodes of y direction
                
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+T(i-1,j))+term3*(T(i,j-1)+TOLD(i,j+1)));
            end
        end
        error=max(max(abs(TOLD-T)));% computes the error at each iteration
        TOLD=T;% updates the value of temperature
       sor_iteration=sor_iteration+1;% computes the no of iteration performed
   end
       TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end

Code for computing optimal omega

%%%%%%%%%%%% Code to compute the optimal value of omega for SOR iterative solver%%%%%
clear all
clc
close all
%%% define the number of grid points required in x and y direction
nx=10;%no of grids in x axis
ny=nx;%no of grids in y axis
alpha=0.1;%thermal diffusivity factor
dt=1e-2;
nt=1400;%no of timesteps
%define the temperature domain matrix
T=ones(nx,ny)*300;
%boundary conditions of temperature domain
TL=400;
TR=800;
TT=600;
TB=900;
%defining the temperatures at the boundaries
T(2:ny-1,1)=TL;% temperature at the left boundary
T(2:ny-1,nx)=TR;% temperature at the right boundary
T(1,2:nx-1)=TT;% temperature at the top boundary
T(nx,2:ny-1)=TB;% temperature at the bottom boundary
%define the temperature at the corner points
T(1,1)=(TT+TL)/2;
T(nx,ny)=(TR+TB)/2;
T(1,ny)=(TT+TR)/2;
T(nx,1)=(TL+TB)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%IMPLICIT SOLVER%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%SUCCESSIVE OVER RELAXATION%%%%%%%%%%%%%
omega_range = linspace(0.5, 1.6, 100);%define the range of omega
%%%%%%%%%%% function call for SOR Gauss-Siedel and SOR Jacobi solver%%%%%%%%%%%
for i = 1:length(omega_range)%marching loop for the entire range of omega
    omega = omega_range(i);
    %stores the no of iteration value for each omega value for jacobi and GS solver
    
iter_gs(i)=transient_implicit_sor_gauss_siedel_optimal_omega(T,nx,ny,alpha,dt,nt, omega);%
iter_jacobi(i)= transient_implicit_sor_jacobi_optimal_omega(T,nx,ny,alpha,dt,nt, omega);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%Computing the optimal node and omega value for SOR Gauss Siedel and SOR Jacobi %%%%%%%%%%%%%%%%%%%%%%%%
for j=1:length(iter_gs)% marching loop for the entire length of iter
 minimum_iteration_gs=min(iter_gs);%stores the minimum value of iteration along the entire range of omega 
 minimum_iteration_jacobi=min(iter_jacobi);%stores the minimum value of iteration along the entire range of omega   
end
position_of_optimal_omega_gs=find(iter_gs==minimum_iteration_gs);% stores the nodal position of the minimum iteration value
optimal_value_of_omega_gs=omega_range(position_of_optimal_omega_gs);%stores the optimum omega value for optimum nodal position of omega

position_of_optimal_omega_jacobi=find(iter_jacobi==minimum_iteration_jacobi);% stores the nodal position of the minimum iteration value
optimal_value_of_omega_jacobi=omega_range(position_of_optimal_omega_jacobi);%stores the optimum omega value for optimum nodal position of omega

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%plotting function%%%%%%%%%%%%%%%%%%
figure(1)%for SOR Gauss Siedel
plot(omega_range, iter_gs);
xlabel ('omega values');
ylabel ('No of iterations taken for convergence');
titletext1=sprintf('Computing optimal omega for SOR Gauss Siedel iterative solver\n');
titletext2=sprintf('Optimal value of omega=%d\n',optimal_value_of_omega_gs);
titletext3=sprintf('Optimal nodal position of omega=%d\n',position_of_optimal_omega_gs);
titletext4=sprintf('Optimal no of iteration for optimal omega value=%d\n',minimum_iteration_gs);
title([titletext1 titletext2 titletext3 titletext4]);
%%%%%%%%%%%%%%%%
figure(2)% for SOR Jacobi
plot(omega_range, iter_jacobi);
xlabel ('omega values');
ylabel ('No of iterations taken for convergence');
titletext1=sprintf('Computing optimal omega for SOR Jacobi iterative solver\n');
titletext2=sprintf('Optimal value of omega=%d\n',optimal_value_of_omega_jacobi);
titletext3=sprintf('Optimal nodal position of omega=%d\n',position_of_optimal_omega_jacobi);
titletext4=sprintf('Optimal no of iteration for optimal omega value=%d\n',minimum_iteration_jacobi);
title([titletext1 titletext2 titletext3 titletext4]);

function call for transient state implicit time integration SOR gauss siedel solver for optimal omega

function sor_iteration = transient_implicit_sor_gauss_siedel_optimal_omega(T,nx,ny,alpha,dt,nt, omega)

l=1;%length of the domain

x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction 
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;


%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;

sor_iteration=1;
w = omega;% relaxation factor
tic
for n=1:nt  % time marching loop
    error=9e9;
   while(error>tolerance)% convergence loop
       
        for i=2:nx-1%space marching loop at inner nodes of x direction
            for j=2:ny-1%space marching loop at inner nodes of y direction
                
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+T(i-1,j))+term3*(T(i,j-1)+TOLD(i,j+1)));
            end
        end
        error=max(max(abs(TOLD-T)));% computes the error at each iteration
        TOLD=T;% updates the value of temperature
       sor_iteration=sor_iteration+1;% computes the no of iteration performed
   end
       TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end

function call for transient state implicit time integration SOR jacobi solver for optimal omega

function sor_iteration = transient_implicit_sor_jacobi_optimal_omega(T,nx,ny,alpha,dt,nt, omega)

l=1;%length of the domain

x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);%grid spacing in x direction
dy=y(2)-y(1);%grid spacing in y direction
tolerance=1e-4;
[X,Y]= meshgrid(x,y);%generating mesh grid
TOLD=T;


%define the constant k1 and k2
k1=alpha*dt/(dx^2);%CFL no at x direction
k2=alpha*dt/(dy^2);%CFL no at y direction
k=k1+k2;%%%total CFL number
TPREV=T;
%defining constant terms
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;

sor_iteration=1;
w=omega;% relaxation factor
tic
for n=1:nt  % time marching loop
    error=9e9;
   while(error>tolerance)% convergence loop
       
        for i=2:nx-1%space marching loop at inner nodes of x direction
            for j=2:ny-1%space marching loop at inner nodes of y direction
                
T(i,j)=(1-w)*TOLD(i,j)+w*(term1*TPREV(i,j)+term2*(TOLD(i+1,j)+TOLD(i-1,j))+term3*(TOLD(i,j-1)+TOLD(i,j+1)));
            end
        end
        error=max(max(abs(TOLD-T)));% computes the error at each iteration
        TOLD=T;% updates the value of temperature
       sor_iteration=sor_iteration+1;% computes the no of iteration performed
   end
       TPREV=T;% stores the value of temperature at the previous time steps
end
time=toc; % computes the time to run the simulation
end

Results and analysis

1.Steady state solution

Thermal plot of Jacobi solver

Thermal plot of Gauss Seidel solver

Thermal plot of SOR solver

Bar plot with respect to no of iterations

Bar plot with respect to simulation time

Iteration and computational table

2.Transient state solution

Thermal plot of Jacobi solver

Thermal plot of Gauss Seidel solver

Thermal plot of SOR Jacobi solver

Thermal plot of SOR Gauss Seidel solver

Bar plot with respect to no of iterations

Bar plot with respect to simulation time

Iteration and computational table

3.Optimal omega plot for SOR Jacobi solver

4.Optimal omega plot for SOR Gauss siedel solver

5.Optimal omega computational table

Conclusion

1. In steady state conduction from implicit solver time integration scheme it is seen that SOR solver converges to steady state solution at a much faster rate than Jacobi and gauss seidel solver.
2. The computational time required for convergence is minimum for SOR and maximum for Jacobi solver.
3. In case of transient state implicit solver, SOR Jacobi and SOR gauss seidel solver converges to steady state solution faster with less no of iterations.
4. The value of omega for SOR solver in transient case is optimized within the range (0.5-1.6), where the optimal omega reaches 1.3 for Jacobi and 1.08 for gauss seidel with minimum no of iterations.
5. At omega equal to 1, no of iteration for jacobi solver is equal to gauss siedel solver.
6. Taking too high omega (>1.9) makes the solution divergent instead of converging to a particular solution. Such divergence gives rise to instability in a problem and must be avoided.