BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Project

Steady and unsteady state analysis of 2D heat conduction equation using different iterative Methods

PPushkal Awasthi11 Nov 2019

Aim: Write a code for 2D heat conduction equation for steady and unsteady state.

Introduction: The 2D heat conduction is given below:

`(∂T)/(∂t)= α*[ (∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 ) ] …………………. (i)`

1) For steady state:

`(∂T)/(∂t)=0`

so `(∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 )=0 ……….. (ii)`

Now we descritized the equation using central differencing scheme…

`(T(i+1,j)– 2*T(i,j)+T(i-1,j))/(dx^2) + (T(i,j+1)– 2*T(i,j)+T(i,j-1))/(dy^2) = 0`

`T(i,j)=(T(i+1,j)+T(i-1,j))/(K*dx^2) + (T(i,j+1)+T(i,j-1))/(K*dy^2)`

Where

`K=2* ([dx^2+ dy^2])/(dx^2*dy^2 )`

2) For unsteady state:

`(∂T)/(∂t)= α*[ (∂^2 T)/(∂x^2 ) + (∂^2 T)/(∂y^2 ) ]`

Now we descritized the equation using central differencing scheme…

  1. Explicit method

`(T_(i,j)^(n+1)-T_(i,j)^n)/dt=α*[(T_(i+1,j)^n – 2*T_(i,j)^n +T_(i-1,j)^n)/(dx^2 ) + (T_(i,j+1)^n – 2*T_(i,j)^n +T_(i,j-1)^n)/(dy^2 ) ]`

`K1 = ([α*dt])/(dx^2 ) , K2 = ([α*dt])/(dy^2 )`

`T_(i,j)^(n+1)=T_(i,j)^n+ K1*[T_(i+1,j)^n – 2*T_(i,j)^n +T_(i-1,j)^n ]+ K2*[T_(i,j+1)^n – 2*T_(i,j)^n +T_(i,j-1)^n]`

  1. Implicit method

`(T_(i,j)^(n+1)-T_(i,j)^n)/dt=α*[(T_(i+1,j)^(n+1) – 2*T_(i,j)^(n+1) +T_(i-1,j)^(n+1))/(dx^2 ) + (T_(i,j+1)^(n+1) – 2*T_(i,j)^(n+1) +T_(i,j-1)^(n+1))/(dy^2 ) ]`

where `K3 = ([α*dt])/(dx^2 ) , K4 = ([α*dt])/(dy^2 )`

`T_(i,j)^(n+1)=[T_(i,j)^n+ K3*[T_(i+1,j)^(n+1) +T_(i-1,j)^(n+1) ]+ K4*[T_(i,j+1)^(n+1)+T_(i,j-1)^(n+1)]]/[1+2*K3+2*K3]`

There are three methods of iterative solver….

  1. Jaicobi method
  2. Gauss seidel method
  3. Sucessive over relaxation method (SOR)

Matlab code for Steady state 2D heat conduction case:

% write a code to destritized a steady state 2D heat conduction equation
clear all
close all
clc

% Steady state case
% First we create the grid space points
nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X Y] = meshgrid(x,y);
% dx and dy are the grid spacing.
dx = x(2) - x(1);
dy = y(2) - Y(1);

% error
error = 9e9;
tolerance = 1e-4;

% Now we define the boundary conditions
T = ones(nx, ny);
% Initial boundary condition given
% T_left face = 400; T_right face = 800;
% T_top face = 600;  T_bottom face = 900;
T_t = 600; T_b = 900; T_l = 400; T_r = 800;

% Now we create the temperature at boundary surface
T(:,1)= 400;
T(:,nx)= 800;
T(1,:)= 600;
T(ny,:)= 900;
T(1,1) = (T_t + T_l)/2;
T(1,nx) = (T_t + T_r)/2;
T(ny,1) = (T_l + T_b)/2;
T(ny,nx) = (T_r + T_b)/2;
Told = T;
% Descritization of 2D steady state heat conduction
% T(i,j) = (1/(k*dx^2))(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1))
% k = 2*(dx^2 + dy^2)/(dx^2*dy^2)
k = 2*(dx^2 + dy^2)/(dx^2*dy^2);
A = (1/(k*dx^2));
% Here above equation is solved by using iterative solver method.
% There are three iterative solver method: Jacobi, Gauss seidel, SOR
% Here first iterative solver is Jacobi method.
%%%%%%%%%%%%%%%%%%%%%%%%% Jacobi Method %%%%%%%%%%%%%%%%%%%%%%%%%
  Iterative_solver = 1
  tic;
   count = 1;
   while(error > tolerance)
    for i = 2:nx-1
      for j = 2:ny-1
          H = Told(i-1,j) + Told(i+1,j);
          V = Told(i,j-1) + Told(i,j+1);
          T(i,j) = A*( H + V);
      end
    end
    error = max(max(abs(Told -T)));
    Told = T;
    count = count +1;
    end
   
   toc;
   % plot
   figure(1)
   [a b] = contourf(X, Y, T);
   clabel(a , b);
   colormap(jet);
   colorbar
   time1 = toc;
   title_text1 = sprintf(\'steady state jacobi iterations = %d\',count);
   title_text2 = sprintf(\'Estimated time = %d\',time1);
   title({title_text1; title_text2})
   xlabel(\'x axis\')
   ylabel(\'y axis\')
   fprintf(\'Number of Jaicobi iterations : %d\', count);

   
%%%%%%%%%% Gauss seidel Method %%%%%%%%%%%%%%%%%%%%%%%%%%%
nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X Y] = meshgrid(x,y);
% dx and dy are the grid spacing.
dx = x(2) - x(1);
dy = y(2) - Y(1);
% error
error = 9e9;
tolerance = 1e-4;

% Now we define the boundary conditions
T = ones(nx, ny);
% Initial boundary condition given
% T_left face = 400; T_right face = 800;
% T_top face = 600;  T_bottom face = 900;
T_t = 600; T_b = 900; T_l = 400; T_r = 800;

% Now we create the temperature at boundary surface
T(:,1)= 400;
T(:,nx)= 800;
T(1,:)= 600;
T(ny,:)= 900;
T(1,1) = (T_t + T_l)/2;
T(1,nx) = (T_t + T_r)/2;
T(ny,1) = (T_l + T_b)/2;
T(ny,nx) = (T_r + T_b)/2;
Told = T;

% Descritization of 2D steady state heat conduction by gauss seidel method
% T(i,j) = (1/(k*dx^2))(T(i-1,j) + Told(i+1,j) + T(i,j-1) + Told(i,j+1))
% k = 2*(dx^2 + dy^2)/(dx^2*dy^2)
k = 2*(dx^2 + dy^2)/(dx^2*dy^2);
A = (1/(k*dx^2));
tic;
Iterative_solver = 2
   count1 = 1;
   while(error > tolerance)
    for i = 2:nx-1
      for j = 2:ny-1
        T(i,j) = A*( T(i-1,j) + Told(i+1,j) + T(i,j-1) + Told(i,j+1));
      end
    end
    error = max(max(abs(Told -T)));
    Told = T;
    count1 = count1 +1;
    end
   toc;
   time2 = toc;
   % plot
   figure(2)
   [c d] = contourf(X, Y, T);
   clabel(c, d);
   colormap(jet);
   colorbar
   title_text1 = sprintf(\'steady state Gauss seidel iterations = %d\',count1);
   title_text2 = sprintf(\'Estimated time = %d\',time2);
   title({title_text1; title_text2})
   xlabel(\'x axis\')
   ylabel(\'y axis\')
   fprintf(\'Number of G S iterations : %d\', count1);
   
   
%%%%%%%%%% SOR Method %%%%%%%%%%%%%
nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X Y] = meshgrid(x,y);
% dx and dy are the grid spacing.
dx = x(2) - x(1);
dy = y(2) - Y(1);
% error
error = 9e9;
tolerance = 1e-4;
alpha = 1.5;  % SOR coefficient
% Now we define the boundary conditions
T = ones(nx, ny);
% Initial boundary condition given
% T_left face = 400; T_right face = 800;
% T_top face = 600;  T_bottom face = 900;
T_t = 600; T_b = 900; T_l = 400; T_r = 800;

% Now we create the temperature at boundary surface
T(:,1)= 400;
T(:,nx)= 800;
T(1,:)= 600;
T(ny,:)= 900;
T(1,1) = (T_t + T_l)/2;
T(1,nx) = (T_t + T_r)/2;
T(ny,1) = (T_l + T_b)/2;
T(ny,nx) = (T_r + T_b)/2;
Told = T;

% Descritization of 2D steady state heat conduction USING SOR method
% T(i,j) = (1-alpha)*Told(i,j) + alpha*(A)(T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1))
% k = 2*(dx^2 + dy^2)/(dx^2*dy^2)
k = 2*(dx^2 + dy^2)/(dx^2*dy^2);
A = (1/(k*dx^2));
% calculations
tic;
Iterative_solver = 3
count2 = 1;
   while(error > tolerance)
    for i = 2:nx-1
      for j = 2:ny-1
        T(i,j) =(1-alpha)*Told(i,j) + alpha*A*( T(i-1,j) + Told(i+1,j) + T(i,j-1) + Told(i,j+1));
      end
    end
    error =max( max(abs(Told -T)));
    Told = T;
    count2 = count2 +1;
  
end
   toc;
   time3 = toc;
   % plot
   figure(3)
   [e f] = contourf(X, Y, T);
   clabel(e, f);
   colormap(jet);
   colorbar
   title_text1 = sprintf(\'steady state SOR iterations = %d\',count2);
   title_text2 = sprintf(\'Estimated time = %d\',time3);
   title({title_text1; title_text2})
   xlabel(\'x axis\')
   ylabel(\'y axis\')
   fprintf(\'Number of SOR iterations : %d\', count2);

OUTPUT:

Conclusions:

Iterative method

No. of iterations

Time taken (sec)

Jacobi method

1884

6.41

Gauss seidel method

1005

0.239

SOR method

363

0.116

From above table we can say that the least time taken by SOR method to converge the solution also have less number of iterations required. Whereas Jacobi method has maximum time taken to converge the solution amongst all iterative method.

Matlab code for unsteady state 2D heat conduction case:

a) Code for explicit method:

% write a code to destritized a 2D heat conduction second order equation by Explicit method
clear all
close all
clc

% Unsteady state case by using explicit method
% First we create the grid space points
nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X Y] = meshgrid(x,y);
% dx and dy are the grid spacing.
dx = x(2) - x(1);
dy = y(2) - Y(1);

% error
error = 9e9;
tolerance = 1e-4;
% time step 
dt = 10;
% Thermal diffusity for steel material
alpha = 1.88e-5; 

% Now we define the boundary conditions
T = ones(nx, ny);
% Initial boundary condition given
% T_left face = 400; T_right face = 800;
% T_top face = 600;  T_bottom face = 900;
T_t = 600; T_b = 900; T_l = 400; T_r = 800;

% Now we create the temperature at boundary surface
T(:,1)= 400;
T(:,nx)= 800;
T(1,:)= 600;
T(ny,:)= 900;
T(1,1) = (T_t + T_l)/2;
T(1,nx) = (T_t + T_r)/2;
T(ny,1) = (T_l + T_b)/2;
T(ny,nx) = (T_r + T_b)/2;
Told = T;
% Descritization of 2D unsteady state heat conduction
% T(i,j) = Told(i,j) + k1*(Told(i-1,j) + Told(i+1,j) - 2*Told(i,j)) + k2*(Told(i,j-1) + Told(i,j+1) - 2*Told(i,j))
% T(i,j) = A + B + C
% A = Told(i,j)
% B = k1*(Told(i-1,j) + Told(i+1,j) - 2*Told(i,j))
% C = k2*(Told(i,j-1) + Told(i,j+1) - 2*Told(i,j))
% k1 = alpha*dt/(dx^2),  k2 = alpha*dt/(dy^2)
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
% In Explicit case we do not need any iterative solver method.
 tic;
 count = 1
  
   
    error = 1e9;
    while(error>tolerance)
    for i = 2:nx-1
      for j = 2:ny-1
          A =  Told(i,j);
          B = k1*(Told(i-1,j) + Told(i+1,j) - 2*Told(i,j));
          C = k2*(Told(i,j-1) + Told(i,j+1) - 2*Told(i,j));
          T(i,j) = A + B + C;
      end
    end
    error = max(max(abs(Told -T)));
    Told = T;
    count = count +1;
   end
   
   time = toc;
   % plot
   figure(1)
   [a b] = contourf(X, Y, T);
   clabel(a , b);
   colormap(jet);
   colorbar
   title_text1 = sprintf(\'Unsteady state explcit iterations = %d\',count);
   title_text2 = sprintf(\'Estimated time = %d\',time);
   title({title_text1; title_text2})
   xlabel(\'x axis\')
   ylabel(\'y axis\')
   fprintf(\'Number of  iterations : %d\', count)

Output:

b) Code for implicit method:

% write a code to destritized a 2D heat conduction second order equation by Implicit method
clear all
close all
clc

% Unsteady state case by using explicit method
% First we create the grid space points
nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X Y] = meshgrid(x,y);
% dx and dy are the grid spacing.
dx = x(2) - x(1);
dy = y(2) - Y(1);

% error
error = 9e9;
tolerance = 1e-4;
% time step 
dt = 10;
nt = 2000;
% Thermal diffusity
alpha = 1.88e-5;

% Now we define the boundary conditions
T = ones(nx, ny);
% Initial boundary condition given
% T_left face = 400; T_right face = 800;
% T_top face = 600;  T_bottom face = 900;
T_t = 600; T_b = 900; T_l = 400; T_r = 800;

% Now we create the temperature at boundary surface
T(:,1)= 400;
T(:,nx)= 800;
T(1,:)= 600;
T(ny,:)= 900;
T(1,1) = (T_t + T_l)/2;
T(1,nx) = (T_t + T_r)/2;
T(ny,1) = (T_l + T_b)/2;
T(ny,nx) = (T_r + T_b)/2;
Told = T;
Told1 = T;
% Descritization of 2D unsteady state heat conduction
% T(i,j) = [ Told1(i,j) + k1*(T(i-1,j) + T(i+1,j)) + k2*(T(i,j-1) + T(i,j+1))]/(1 + 2*k1 + 2*k2)
% T(i,j) = (A + B + C)/D
% A = Told1(i,j)
% B = k1*(T(i-1,j) + T(i+1,j))
% C = k2*(T(i,j-1) + T(i,j+1)
% D = (1 + 2*k1 + 2*k2)
% k1 = alpha*dt/(dx^2),  k2 = alpha*dt/(dy^2)
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
% Here above equation is solved by using iterative solver method.
% There are three iterative solver method: Jacobi, Gauss seidel, SOR
% Here first iterative solver is Jacobi method.
% 1-Jacobi Method 
 Iterative_solver = 1
   tic;
   count = 1;
  for k = 1:nt
      error = 9e9;
   while(error > tolerance)
     for i = 2:nx-1
      for j = 2:ny-1
          A = Told1(i,j);
          B = k1*(Told(i-1,j) + Told(i+1,j));
          C = k2*(Told(i,j-1) + Told(i,j+1));
          D = (1 + 2*k1 + 2*k2);
          T(i,j) = (A + B + C)/D;
      end
    end
    error = max(max(abs(Told -T)));
    Told = T;
    count = count +1;
   end
  Told1 = T;
  end
   time = toc;
   % plot
   figure(1)
   [a b] = contourf(X, Y, T);
   clabel(a , b);
   colormap(jet);
   colorbar
   title_text1 = sprintf(\'Unsteady state implicit jacobi iterations = %d\',count);
   title_text2 = sprintf(\'Estimated time = %d\',time);
   title({title_text1; title_text2})
   xlabel(\'x axis\')
   ylabel(\'y axis\')
   fprintf(\'Number of Jaicobi iterations : %d\', count)
   
   
%%%%%%%%%%%%%%%%%%%%% Gauss seidel method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X Y] = meshgrid(x,y);
% dx and dy are the grid spacing.
dx = x(2) - x(1);
dy = y(2) - Y(1);

% error
error = 9e9;
tolerance = 1e-4;
% time step 
dt =10;
nt = 2000;
% Thermal diffusity
alpha = 1.88e-5;

% Now we define the boundary conditions
T = ones(nx, ny);
% Initial boundary condition given
% T_left face = 400; T_right face = 800;
% T_top face = 600;  T_bottom face = 900;
T_t = 600; T_b = 900; T_l = 400; T_r = 800;

% Now we create the temperature at boundary surface
T(:,1)= 400;
T(:,nx)= 800;
T(1,:)= 600;
T(ny,:)= 900;
T(1,1) = (T_t + T_l)/2;
T(1,nx) = (T_t + T_r)/2;
T(ny,1) = (T_l + T_b)/2;
T(ny,nx) = (T_r + T_b)/2;
Told = T;
Told1 = T;
% Descritization of 2D unsteady state heat conduction
% T(i,j) = [ Told1(i,j) + k1*(T(i-1,j) + T(i+1,j)) + k2*(T(i,j-1) + T(i,j+1))]/(1 + 2*k1 + 2*k2)
% T(i,j) = (A + B + C)/D
% A = Told1(i,j)
% B = k1*(T(i-1,j) + T(i+1,j))
% C = k2*(T(i,j-1) + T(i,j+1)
% D = (1 + 2*k1 + 2*k2)
% k1 = alpha*dt/(dx^2),  k2 = alpha*dt/(dy^2)
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

% Gauss seidel Method 
 Iterative_solver = 2
   tic;
   count1 = 1;
   for k = 1:nt
       error = 9e9;
   while(error > tolerance)
     for i = 2:nx-1
      for j = 2:ny-1
          A = Told1(i,j);
          B = k1*(T(i-1,j) + Told(i+1,j));
          C = k2*(T(i,j-1) + Told(i,j+1));
          D = (1 + 2*k1 + 2*k2);
          T(i,j) = (A + B + C)/D;
      end
    end
    error = max(max(abs(Told -T)));
    Told = T;
    count1 = count1 +1;
   end
  Told1 = T;
   end
   time1 = toc;
   % plot
   figure(2)
   [c d] = contourf(X, Y, T);
   clabel(c , d);
   colormap(jet);
   colorbar
   title_text1 = sprintf(\'Unsteady state implicit GS iterations = %d\',count1);
   title_text2 = sprintf(\'Estimated time = %d\',time1);
   title({title_text1; title_text2})
   xlabel(\'x axis\')
   ylabel(\'y axis\')
   fprintf(\'Number of GS iterations : %d\', count1)
   
   
   
%%%%%%%%%%%%%%%%%%%%%%%%%% SOR method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X Y] = meshgrid(x,y);
% dx and dy are the grid spacing.
dx = x(2) - x(1);
dy = y(2) - Y(1);

% error
error = 9e9;
tolerance = 1e-4;
% time step 
dt =10;
nt = 2000;
% Thermal diffusity
alpha = 1.88e-5;
beta = 2/(1 + sin(pi*dx));

% Now we define the boundary conditions
T = ones(nx, ny);
% Initial boundary condition given
% T_left face = 400; T_right face = 800;
% T_top face = 600;  T_bottom face = 900;
T_t = 600; T_b = 900; T_l = 400; T_r = 800;

% Now we create the temperature at boundary surface
T(:,1)= 400;
T(:,nx)= 800;
T(1,:)= 600;
T(ny,:)= 900;
T(1,1) = (T_t + T_l)/2;
T(1,nx) = (T_t + T_r)/2;
T(ny,1) = (T_l + T_b)/2;
T(ny,nx) = (T_r + T_b)/2;
Told = T;
Told1 = T;
   
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

% SOR Method 
 Iterative_solver = 3
   tic;
   count2 = 1;
   for k = 1:nt
       error = 9e9;
   while(error > tolerance)
     for i = 2:nx-1
      for j = 2:ny-1
         term1= (1 + 2*k1 + 2*k2)^-1;
         term2=  k1*term1;
         term3=  k2*term1;
          H = T(i-1,j) + Told(i+1,j);
          V = T(i,j-1) + Told(i,j+1);
          T(i,j) = Told1(i,j)*term1 + (H*term2) + (V*term3);
          T(i,j) = (Told1(i,j)*(1-beta)) + beta*T(i,j);
       end
    end
    error = max(max(abs(Told -T)));
    Told = T;
    count2 = count2 +1;
   end
  Told1 = T;
   end
   time2 = toc;
   % plot
   figure(3)
   [e f] = contourf(X, Y, T);
   clabel(e , f);
   colormap(jet);
   colorbar
   title_text1 = sprintf(\'Unsteady state implicit SOR iterations = %d\',count2);
   title_text2 = sprintf(\'Estimated time = %d\',time2);
   title({title_text1; title_text2})
   xlabel(\'x axis\')
   ylabel(\'y axis\')
   fprintf(\'Number of SOR iterations : %d\', count2)

Output:

Conclusions:

Iterative method

No. of iterations

Time taken (sec)

Jacobi method(i)

16485

 5.45

Gauss seidel method(i)

12329

 4.04

SOR method(i)

9800

4.25

Explicit

2858

 1.2
  • From above table we can say that explicit method have taken least time to converge the solution for transient 2D heat conduction equation.
  • In case of implicit method of discretization gauss seidel have taken less time to converge the solution among other iterative solver. If we compare the solution in terms of required iterations to converge the solutions, the successive over relaxation (SOR) had required minimum number of iterations.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.