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

Solving the 2-D Heat Conduction Problem

PPritam Pritam19 Jul 2019

Introduction:

The 2-D Heat Conduction to be solved:

`alpha*((del^2T)/(delx^2)+(del^2T)/(dely^2))=(delT)/(delt)`

This is the unsteady form of 2D heat conduction equation. When the steady-state analysis is to be done:

`(delT)/(delt) = 0`

Steady state equation:

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

For Jacobi Solver, old values will be used in the scheme while for the Gauss-Seidel updated values are to be used. In SOR method, a factor of 1.7 is used to speed up the time required for reaching convergence.

Case 1(Steady State Implicit Solver):

The Central Difference Scheme is used as the discretization scheme here for the spatial terms:

`(del^2T)/(delx^2) = (T(i+1,j)-2*T(i,j)+T(i-1,j))/(delx^2)`

`(del^2T)/(dely^2) = (T(i,j+1)-2*T(i,j)+T(i,j-1))/(delx^2)`

Equating the sum of the spatial terms to zero,

`T(i,j) = (T(i,j-1)+T(i-1,j)+T(i+1,j)+T(i,j+1))/4`

Mesh used here is a 20x20 elements mesh.

Code for the main program that calls the appropriate functions for solving:

clc
clear all
close all
% number of elements in x and y direction
nx = 20;
ny = 20;
% tempreture distribution on the square plate
t = ones(20,20);
t(1,:)=600;
t(:,1)=400;
t(:,end)=800;
t(end,:)=900;
t(1,end)=(600+800)/2;
t(1,1)=(600+400)/2;
t(end,1)=(900+400)/2;
t(end,end)=(900+800)/2;
% calling the methods for performing the specific approach
my_jacobi(nx,ny,t)
my_siedel (nx,ny,t)
my_sor (nx,ny,t)

Steady State Jacobi Solver:

function [] = my_jacobi(n_x,n_y,temp)

told = temp;
t_allow = 1e-4;
error = 10;
iter=1;

while (error > t_allow)
    for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = (told(i+1,j)+told(i-1,j)+told(i,j+1)+told(i,j-1))/4;
        end
    end
    error = max(max(abs(temp-told)));
    figure(1)
    C=contourf(temp);
    clabel(C)
    colormap(jet)
    colorbar
    pause(0.03)
    title([\'Iteration number =\',num2str(iter)]);
    told = temp;
    iter = iter + 1;
end

Steady State Gauss-Seidel Solver:

function [] = my_siedel(n_x,n_y,temp)

told = temp;
t_allow = 1e-4;
error = 10;
iter=1;

while (error > t_allow)
    for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = (told(i+1,j)+temp(i-1,j)+told(i,j+1)+temp(i,j-1))/4;
        end
    end
    error = max(max(abs(temp-told)));
    figure(2)
    C=contourf(temp);
    clabel(C)
    colormap(jet)
    colorbar
    title([\'Iteration number =\',num2str(iter)]);
    told = temp;
    iter = iter + 1;
end

Steady State SOR Solver:

function [] = my_sor(n_x,n_y,temp)

told = temp;
t_allow = 1e-4;
error = 10;
iter=1;
fou=1.7;

while (error > t_allow)
    for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = (1-fou)*told(i,j)+fou*((told(i+1,j)+temp(i-1,j)+told(i,j+1)+temp(i,j-1))/4);
        end
    end
    error = max(max(abs(temp-told)));
    figure(3)
    C=contourf(temp);
    clabel(C)
    colormap(jet)
    colorbar
    title([\'Iteration number =\',num2str(iter)]);
    pause (0.03)
    told = temp;
    iter = iter + 1;
end

Plots Obtained:

Jacobi plot:-

Gauss-Seidel Plot:-

SOR Plot:

Case 2 (Unsteady Implicit Solver):

The Central difference Scheme is used here as the above case.

Discretizing the unsteady 2D conduction,

`T^(n+1)(i,j)=T^(n)(i,j)+alpha*Deltat((T(i-1,j)-2*T(i,j)+T(i+1,j))/(Deltax^2)+(T(i,j-1)-2*T(i,j)+T(i,j+1))/(Deltay^2))^n`

which gives,

`T^(n+1)(i,j)=T^(n)(i,j)+k1(T(i-1,j)-2*T(i,j)+T(i+1,j))+k2(T(i,j-1)-2*T(i,j)+T(i,j+1))`

where `k1 = (alpha * Deltat)/(Deltax^2)` and `k2 = (alpha *Deltat)/(Deltay^2)`

`=> T^(n+1)(i,j) = T^(n)(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)`

Code for the main program that calls the appropriate functions for solving:

clc
clear all
close all
% meshing
nx = 20;
msh = linspace(0,1,nx);
dx = msh(2)-msh(1);
dy=dx;
ny = 20;
% time steps 
dt = 1e-4;
time=10;
n=time/dt;
% tempreture initialization
t = 300*ones(20,20);
t(1,:)=600;
t(:,1)=400;
t(:,end)=800;
t(end,:)=900;
t(1,end)=(600+800)/2;
t(1,1)=(600+400)/2;
t(end,1)=(900+400)/2;
t(end,end)=(900+800)/2;
% constants definition
alpha = 0.000097; 
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
% function calling
my_implicit_jacobi(nx,ny,term1,term2,term3,t,n)
my_implicit_seidel(nx,ny,term1,term2,term3,t,n)
my_implicit_sor(nx,ny,term1,term2,term3,t,n)

Jacobi Solver:

% unsteady  jacobi implicit
function [] = my_implicit_jacobi(n_x,n_y,term_1,term_2,term_3,temp,tstep)

told = temp;
tprev =temp;
t_allow = 1e-4;
error = 1e5;
iter =1;
iter_timestep=1;

for k = 1:tstep
  
  error = 1e5;
  while (error > t_allow)
    
      for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = term_1*tprev(i,j) + (told(i-1,j)+told(i+1,j))*term_2 +(told(i,j-1)+told(i,j+1))*term_3;
        end
      end
      
    error = max(max(abs(temp-told)));
    told = temp;
    iter = iter + 1;
  end
  
  iter_timestep = iter_timestep+1;
  
  if(temp-tprev < 1e-5)
    break
  end
  tprev = temp;
end
figure(1)
C=contourf(temp);
clabel(C)
colormap(jet)
colorbar
title([\'Iteration number =\',num2str(iter)]);

Gauss Seidel Solver:

function [] = my_implicit_seidel(n_x,n_y,term_1,term_2,term_3,temp,tstep)

told = temp;
tprev =temp;
t_allow = 1e-4;
error = 1e5;
iter =1;
iter_timestep=1;

for k = 1:tstep
  error = 1e5;
  
  while (error > t_allow)
    
      for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = term_1*tprev(i,j) + (temp(i-1,j)+told(i+1,j))*term_2 +(temp(i,j-1)+told(i,j+1))*term_3;
        end
      end
      error = max(max(abs(temp-told)));
      told = temp;
      iter = iter + 1;
  end
  
  iter_timestep = iter_timestep+1;
  
  if(temp-tprev < 1e-5)
    break
  end
  tprev = temp;
end
figure(2)
C=contourf(temp);
clabel(C)
colormap(jet)
colorbar
title([\'Iteration number =\',num2str(iter)]);

SOR Solver:

function [] = my_implicit_sor(n_x,n_y,term_1,term_2,term_3,temp,tstep)

told = temp;
tprev =temp;
t_allow = 1e-4;
error = 1e5;
iter =1;
iter_timestep=1;
sl = 1.4;

for k = 1:tstep
  
  error = 1e5;
  while (error > t_allow)
    
      for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = (1-sl)*told(i,j) + sl*(term_1*tprev(i,j) + (temp(i-1,j)+told(i+1,j))*term_2 +(temp(i,j-1)+told(i,j+1))*term_3);
        end
      end
      
    error = max(max(abs(temp-told)));
    told = temp;
    iter = iter + 1;
  end
  
  iter_timestep = iter_timestep+1;
  
  if(temp-tprev < 1e-5)
    break
  end
  tprev = temp;
end
figure(3)
C=contourf(temp);
clabel(C)
colormap(jet)
colorbar
title([\'Iteration number =\',num2str(iter)];

Plots Obtained:

Jacobi:

Gauss-Seidel:

SOR:

Case 3(Unsteady Explicit Solver):

The Central difference Scheme is used here as the above case.

Discretizing the unsteady 2D conduction,

`(T^(n+1)(i,j) - T^(n)(i,j`T^(n+1)(i,j) = T^(n)(i,j)+alpha*Deltat((T(i-1,j)-2*T(i,j)+T(i+1,j))/(Deltax^2)+(T(i,j-1)-2*T(i,j)+T(i,j+1))/(Deltay^2))^n`))/(Deltat) = alpha*((T(i-1,j)-2*T(i,j)+T(i+1,j))/(Deltax^2)+(T(i,j-1)-2*T(i,j)+T(i,j+1))/(Deltay^2))^n`

Code for the main program that calls the appropriate functions for solving:

clc
clear all
close all

nx = 30;
ny = 30;
msh = linspace(0,1,nx);
dx = msh(2)-msh(1);
alpha = 0.000097;
dt = 1e-3;

t = ones(30,30);
t(1,:)=600;
t(:,1)=400;
t(:,end)=800;
t(end,:)=900;
t(1,end)=(600+800)/2;
t(1,1)=(600+400)/2;
t(end,1)=(900+400)/2;
t(end,end)=(900+800)/2;

term = (alpha*dt)/(dx*dx);

my_jacobi_time(nx,ny,term,t);
my_seidel_time(nx,ny,term,t);
my_sor_time(nx,ny,term,t);

Jacobi Solver:

function [] = my_jacobi_time(n_x,n_y,k,temp)

alpha = 0.000097;
told = temp;
t_allow = 1e-4;
error = 10;
iter=1;

while (error > t_allow)
    for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = told(i,j)+ k*(told(i-1,j)-4*told(i,j)+told(i+1,j)+told(i,j-1)+told(i,j+1));
        end
    end
    error = max(max(abs(temp-told)));
    
    told = temp;
    iter = iter + 1;
end
figure(1)
C=contourf(temp);
clabel(C)
colormap(jet)
colorbar
title([\'Iteration number =\',num2str(iter)]);

Gauss-Seidel Solver:

function [] = my_seidel_time(n_x,n_y,k,temp)

alpha = 0.000097;
told = temp;
t_allow = 1e-4;
error = 10;
iter=1;

while (error > t_allow)
    for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = told(i,j)+ k*(temp(i-1,j)-4*told(i,j)+told(i+1,j)+temp(i,j-1)+told(i,j+1));
        end
    end
    error = max(max(abs(temp-told)));
    
    told = temp;
    iter = iter + 1;
end
figure(2)
C=contourf(temp);
clabel(C)
colormap(jet)
colorbar
title([\'Iteration number =\',num2str(iter)]);

SOR Solver:

function [] = my_sor_time(n_x,n_y,k,temp)

alpha = 0.000097;
told = temp;
t_allow = 1e-4;
error = 10;
iter=1;
fou = 1.35;

while (error > t_allow)
    for i=2 : n_x-1
        for j= 2 : n_y-1
           temp(i,j) = (1-fou)*told(i,j)+fou*(told(i,j)+ k*(temp(i-1,j)-4*told(i,j)+told(i+1,j)+temp(i,j-1)+told(i,j+1)));
        end
    end
    error = max(max(abs(temp-told)));
    
    told = temp;
    iter = iter + 1;
end
figure(3)
C=contourf(temp);
clabel(C)
colormap(jet)
colorbar
title([\'Iteration number =\',num2str(iter)]);

Plots Obtained:

Jacobi:-

Gauss-Seidel:

SOR:

Observation:

1. Steady-state solvers take fewer iteration steps to converge than unsteady solvers. This happens because the solution is independent of time and it does not change its properties with time, as a result, we get the solution easily and at a much faster rate.

2. For an explicit solver, we need to take care of the Courant Number which should be less than one for providing a better convergence rate while in implicit solver the Courant Number does not hold a major factor in determining the convergence.

3. The SOR provides faster convergence than Jacobi or Gauss-Seidel as it helps in changing guess values more rapidly than latter solvers.

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.