Simulation for steady-state Transient unsteady-state 2D heat conduction using various iterative techniques

Objective: To Solve the 2D heat conduction equation by using the point iterative techniques                  (Jacobi, Gauss-Seidel, SOR)

Theory:

A generalized form of the equation can be written as for 3D space:

where, dT/dt- time derivative term 

           α- constant for linear

PART-I (STEADY-STATE)

In this challenge, we are going to evaluate for 2D space, equation breakdown to

equation---1

We are going to use a 2nd order CD scheme

FORMULATION of Problem:

From equation 1 we will get a system of linear equations which will be solved by following iterative techniques

1. JACOBI METHOD
2. GAUSS-SEIDEL METHOD
3.Successive Over-Relaxation METHOD (SOR-GS method)

1.Jacobi Method:

In numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a diagonally dominant system of a linear equation. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges. This algorithm is a stripped-down version of the Jacobi transformation method of matrix diagonalization.

2. Gauss-seidel Method:

In numerical linear algebra, the Gauss-Seidel method, also known as the Liebmann method or the method of successive displacement, is an iterative method used to solve a linear system of equation and is similar to the Jacobi method. Though it can be applied to any matrix with non-zero elements on the diagonals, convergence is only guaranteed if the matrix is either diagonally dominant or symmetric and positive definite.

3.Successive Over-Relaxation method:

In numerical linear algebra, the method of successive over-relaxation (SOR) is a variant of the Gauss-seidel Method for solving a linear system of equation, resulting in faster convergence. A similar method can be used for any slowly converging iterative process.

Solving Equations for Steady-state: 

Assumptions:

1.domain is a unit square

2.nx = ny [Number of points along the x-direction is equal to the number of points along the     y direction]

3.Boundary conditions for steady and transient case

    1. Left:400K

    2. Top:600K

    3. Right:800K

    4. Bottom:900K

4. Absolute error criteria is 1e-4

MATLAB CODE:

% For Steady-state 2D-Heat conduction
%d^2T/dx^2+d^2T/dy^2=0

clear all
close all
clc

% inputs
nx=10; % for x-axis
ny=nx; % for y-axis


T=273.15*ones(nx,ny);

% Applying Boundary condition
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;

% creating old temp. value
%T_old=T;

x=linspace(0,2,nx);
dx=x(2)-x(1);

y=linspace(0,2,ny);
dy=y(2)-y(1);

% constant l=2*(1/dx^2+1/dy^2)
l=2*((1/dx^2)+(1/dy^2));

error=1e4;
tol=1e-4; % criteria for absolute error

alpha=1.5 % Over-relaxation factor for SOR-GS methode


[XX,YY]=meshgrid(x,y)
T_old=T;


Solving_method=3 % For selecting method

%Jacobi Method

if Solving_method ==1
k=1; % for Nos. of iterations

tic;
 while (error>tol)
   for i=2:(nx-1)
      for j=2:(ny-1)
        t_1 =(T_old(i-1,j)+T_old(i+1,j))/(l*dx^2);
        t_2 =(T_old(i,j-1)+T_old(i,j+1))/(l*dy^2);
        T(i,j)=t_1+t_2;
      end
   end 
   t(i,j)=toc;
   error = max(max((abs(T_old-T))))
   T_old=T;
   figure(1);
   [c,h]=contourf(XX,YY,T)  ;
   colorbar;
   clabel(c,h);
   xlabel(\'X-aixs\');
   ylabel(\'Y-axis\');
   k=k+1;
   time1=max(max(t));
   title_text= sprintf(\'Simulation for Steady-state 2D Heat conduction( by JACOBI solver)-Iteration number= %d  Simulation time= %0.3f s\',k,time1 );
   title(title_text);
   pause(0.003);
  
   
 end

end


% Gauss-seidal method

if Solving_method ==2
k=1;

tic;
  while (error>tol)
   for i=2:(nx-1)
      for j=2:(ny-1)
        t_1 =(T(i-1,j)+T_old(i+1,j))/(l*dx^2)
        t_2 =(T(i,j-1)+T_old(i,j+1))/(l*dy^2)
        T(i,j)=t_1+t_2;
      end
   end 
  t(i,j)=toc; 
   error = max(max((abs(T_old-T))))
   T_old=T;
   figure(1)
   [c,h]=contourf(XX,YY,T)  
   colorbar
   clabel(c,h)
   xlabel(\'X-aixs\')
   ylabel(\'Y-axis\')
   k=k+1;
   time2=max(max(t));
   title_text= sprintf(\'Simulation for Steady-state 2D Heat conduction( by GS solver)-Iteration number= %d  Simulation time= %0.3f s\',k,time2 )
   title(title_text)
   pause(0.003)
   
  end
end


% SOR method

if Solving_method ==3
k=1;

tic;
 while (error>tol)
   for i=2:(nx-1)
      for j=2:(ny-1)
        
        Tgs1=(T(i-1,j)+T_old(i+1,j))/(l*dx^2);
        Tgs2=(T(i,j-1)+T_old(i,j+1))/(l*dy^2);
        t_0 =((1-alpha)*T_old(i,j))+(alpha*(Tgs1+Tgs2));
        T(i,j)=t_0;                           
        
      end
   end 
  t(i,j)=toc; 
   error = max(max((abs(T_old-T))))
   T_old=T;
   figure(1)
   [c,h]=contourf(XX,YY,T)  
   colorbar
   clabel(c,h)
   xlabel(\'X-aixs\')
   ylabel(\'Y-axis\')
   k=k+1;
   time3=max(max(t));
   title_text= sprintf(\'Simulation for Steady-state 2D Heat conduction( by SOR-GS solver)-Iteration number= %d  Simulation time= %0.3f s\',k,time3 )
   title(title_text)
   pause(0.003)
  
 end
end


 

The results are as below:

Select
Solving_method=1 For JACOBI METHOD:

Select
Solving_method=2 For GAUSS-SEIDEL METHOD:

Select
Solving_method=3 For SOR-GS METHOD:

Above graphs are to visualize the difference of using different iterative techniques for solving  2D heat conduction equation for the STEADY-STATE condition.

Iteration numbers and simulation time is as below:

1.JACOBI METHOD:
   Iteration nos.=209, Simulation time=169.699 sec
2.GAUSS-SEIDEL METHOD:
   Iteration nos.=112, Simulation time=96.339 sec
2.SOR-GS METHOD:
   Iteration nos.=29, Simulation time=23.918 sec

Conclusion:

1. Iteration nos. and simulation time for SOR-GS method is quite less as compared to other iterative methods that mean SOR-GS has the ability to converge towards an exact solution at a faster rate than that of other iterative solving methods.

2. alpha(relaxation factor) plays an important role in the SOR method, here we use alpha=1.5, that\'s why solution change is a bit more aggressive towards an exact solution. hence iterative solver is overtuned. This condition we usually called \'Over-relaxed\'

3. If we use alpha= 0 to 1, then iteration time and simulation time considerably increases. As alpha value changes from 0 to 1, we are shifting new guess value from no updated value to updated guess value using GS or JACOBI (if we use JACOBI instead GS). This condition we called \'Under-relaxed\'

PART-II (TRANSIENT-STATE)

A generalized form of the equation can be written as for 3D space:

where, dT/dt- time derivative term 

           α- constant for linear

We will be considering for 2D space, equation breakdown to:

Time derivative will be solved by using FDS and for Spatial derivative 2nd order CDS.

Time derivative in 2D space:

Solving Equations for Transient-state: (IMPLICIT FORM)

Here in implicit form, RHS derivative will be computed by using n+1 terms which further will provide a set fo system of coupled linear equations.

This system of equations will be solved by following iterative solvers:

1. JACOBI SOLVER
2. GAUSS-SEIDEL SOLVER
3. SOR-GS SOLVER

MATLAB CODE:

% For Transient state 2D-Heat conduction- IMPLICIT form
%dT/dt+omega*(d^2T/dx^2+d^2T/dy^2)=0

clear all
close all
clc

% inputs
nx=10; % for x-axis
ny=nx; % for y-axis



% Applying Boundary condition
T=273.15*ones(nx,ny);

T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;



x=linspace(0,2,nx);
dx=x(2)-x(1);

y=linspace(0,2,ny);
dy=y(2)-y(1);


error=9e9;
tol=1e-4; % criteria for absolute error

alpha=1.5 % Over-relaxation factor for SOR-GS methode


[XX,YY]=meshgrid(x,y)

% creating old temp. value
T_old=T; %for space
T_prev=T; %for time

% inputs for Transient condition
time=4;
dt=1e-2;
nt=time/dt; % nos. of temporal nodes
omega=1.4; % constant for linear
I1= ((omega*dt)/dx^2);
I2= ((omega*dt)/dy^2);


Solving_method=1 % For selecting solver method

%Jacobi Method

if Solving_method ==1

 tic;
 k=1; % for Nos. of iterations
for l=1:nt
  
error = 9e9;
  while (error>tol)
   for i=2:(nx-1)
        for j=2:(ny-1)
         
          I=(1+(2*I1)+(2*I2))^(-1);
          t_1 =(T_old(i-1,j)+T_old(i+1,j));
          t_2 =(T_old(i,j-1)+T_old(i,j+1));
          T(i,j)=(T_prev(i,j)*I)+(I1*t_1*I)+(I2*t_2*I);
          
        endfor
       
    endfor
    
    error = max(max(abs(T_old-T)));
    T_old=T; 
    t(i,j)=toc;
    figure(1);
   [c,h]=contourf(XX,YY,T)  ;
   colorbar;
   clabel(c,h);
   xlabel(\'X-aixs\');
   ylabel(\'Y-axis\');
   
   k=k+1;
   time1=max(max(t));
   title_text= sprintf(\'Simulation for Transient 2D Heat conduction( by JACOBI solver)-Iteration number= %d  Simulation time= %0.3f s\',k,time1 );
   title(title_text);
   pause(0.003);

 endwhile
 T_prev=T; 

endfor


endif


% Gauss-seidal method

if Solving_method ==2

  tic;
  k=1; % for Nos. of iterations
 for l=1:nt
  
  error = 9e9;
  while (error>tol)
   for i=2:(nx-1)
      for j=2:(ny-1)
        I=(1+(2*I1)+(2*I2))^(-1);
          t_1 =(T(i-1,j)+T_old(i+1,j));
          t_2 =(T(i,j-1)+T_old(i,j+1));
          T(i,j)=(T_prev(i,j)*I)+(I1*t_1*I)+(I2*t_2*I);
      endfor
   endfor
   t(i,j)=toc; 
   error = max(max((abs(T_old-T))))
   T_old=T;
   figure(1)
   [c,h]=contourf(XX,YY,T);
   colorbar;
   clabel(c,h);
   xlabel(\'X-aixs\');
   ylabel(\'Y-axis\');
   k=k+1;
   time2=max(max(t));
   title_text= sprintf(\'Simulation for Transient state 2D Heat conduction( by GS solver)-Iteration number= %d  Simulation time= %0.3f s\',k,time2);
   title(title_text);
   pause(0.003);
  
   endwhile
  T_prev=T;
  
  endfor

endif


% SOR method

if Solving_method ==3
tic;
  k=1; % for Nos. of iterations
   #
  
 for l=1:nt
  
  error = 9e9;       
 while (error>tol)
   for i=2:(nx-1)
      for j=2:(ny-1)
        
        I=(1+(2*I1)+(2*I2))^(-1);
        t_1 =alpha*T(i-1,j)+(1-alpha)*T_old(i-1,j)+T_old(i+1,j);
        t_2 =alpha*T(i,j-1)+(1-alpha)*T_old(i,j-1)+T_old(i,j+1);
        T(i,j)=(T_prev(i,j)*I)+((I1*t_1*I)+(I2*t_2*I));
    
      endfor
   endfor 
  t(i,j)=toc; 
   error = max(max((abs(T_old-T))))
   T_old=T;
   figure(1)
   [c,h]=contourf(XX,YY,T);
   colorbar;
   clabel(c,h);
   xlabel(\'X-aixs\');
   ylabel(\'Y-axis\');
   k=k+1;
   time3=max(max(t));
   title_text= sprintf(\'Simulation for Transient-state 2D Heat conduction( by SOR-GS solver)-Iteration number= %d  Simulation time= %0.3f s\',k,time3 );
   title(title_text);
   pause(0.003);
  
 endwhile
T_prev=T;
endfor

endif

Results are as follows:

1. Select- solving method=1 for JACOBI solver

2. Select- solving method=2 for GAUSS-SEIDEL solver

2. Select- solving method=2 for GAUSS-SEIDEL solver

the results that we have is as

1. JACOBI solver- Simulation time= 3900.79 sec taken by this method to converge solution to error less than 1e-5. As this solver taking all old values to calculate next values for the next time step.

2. GAUSS-SEIDEL solver- Simulation time= 1419.487 sec took by this method to converge solution to error less than 1e-5. As this solver taking all old values as well as newly calculated values to calculate the next values for the next time step.

3.SOR- GS solver- Simulation time= 789.935 sec took by this method to converge solution to error less than 1e-5. As this solver taking all old values as well as newly calculated values with a relaxation factor of 1.5 to calculate the next values for the next time step.

Solving Equations for Transient-state: (EXPLICIT FORM)

Here in explicit form, RHS derivative will be computed by using n terms which further will provide a set fo system of uncoupled linear equations.

MATLAB CODE:

%For Transient state 2D-Heat conduction- EXPLICIT form
%dT/dt+omega*(d^2T/dx^2+d^2T/dy^2)=0

clear all
close all
clc

% inputs
nx=10; % for x-axis
ny=nx; % for y-axis



% Applying Boundary condition
T=273.15*ones(nx,ny);

T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;



x=linspace(0,2,nx);
dx=x(2)-x(1);

y=linspace(0,2,ny);
dy=y(2)-y(1);


error=9e9;
tol=1e-4; % criteria for absolute error


[XX,YY]=meshgrid(x,y)

% creating old temp. value
T_old=T; %for space
T_prev=T; %for time

% inputs for Transient condition
time=4;
dt=1e-2;
nt=time/dt; % nos. of temporal nodes
omega=0.4; % constant for linear
I1= ((omega*dt)/dx^2);
I2= ((omega*dt)/dy^2);


 tic;
 k=1; % for Nos. of iterations
for l=1:nt
  
error = 9e9;
  while (error>tol)
   for i=2:(nx-1)
        for j=2:(ny-1)
         
          t_1 =(T_old(i-1,j)-2*T_old(i,j)+T_old(i+1,j));
          t_2 =(T_old(i,j-1)-2*T_old(i,j)+T_old(i,j+1));
          T(i,j)=(T_prev(i,j))+(I1*t_1)+(I2*t_2);
          
        endfor
       
    endfor
    
    error = max(max(abs(T_old-T)))
    T_old=T; 
    t(i,j)=toc;
    figure(1);
   [c,h]=contourf(XX,YY,T)  ;
   colorbar;
   clabel(c,h);
   xlabel(\'X-axis\');
   ylabel(\'Y-axis\');
   
   k=k+1;
   time1=max(max(t));
   title_text= sprintf(\'Simulation for Transient 2D Heat conduction( by Explicit Method)-Iteration number= %d  Simulation time= %0.3f s\',k,time1 );
   title(title_text);
   pause(0.003);

 endwhile
T_prev=T;
endfor

The results are as below:

Conclusion:

From the above graphs we can easily understand that:

Simulation time in STEADY STATE= JACOBI  > GAUSS-SEIDEL > SOR-GS

Simulation time in TRANSIENT STATE (IMPLICIT)= JACOBI  > GAUSS-SEIDEL > SOR-GS

 But when we go for solving TRANSIENT STATE (EXPLICIT METHOD) solution is conditionally stable whereas Transient state (Implicit) solution is Unconditionally stable that means the method is stable for all values of dt, dx, and dy.