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

2D-Heat conduction problem using steady and unsteady state conditions

SShaik Moulali1 Dec 2019

2D-Heat conduction problem using steady and unsteady state conditions:

Here the objective of the problem is to solve the 2D-heat conduction equation for steady and unsteady state conditions using different iterative techniques like Jacobi Method, Gauss-Siedel Method and Successive Over Relaxation Method.

Given:

  • Assume that the domain is a unit square
  • Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
  • Boundary conditions (Temperatures) for steady and transient case
  1. Left:400K
  2. Top:600K
  3. Right:800K
  4. Bottom:900K
  • Given point iterative techniques are Jacobi Method, Gauss-Siedel Method and Successive Over Relaxation Method.
  • Error tolerance is 1e-4

2D-Heat conduction equation (Steady state):

  • The heat equation is a parabolic partial differential equation, which describes the transfer or the distribution of heat in a given unit area over a unit of time.
  • Since local heat conduction is proportional to the local temperature gradient, then the equation must be a function of temperature, location and time.
  • The 2D-Heat conduction equation is

x,y,z are the coordinate points

T is the temperature

t is the time

alpha is the diffusion coefficient (m2/sec)

  • But, the given condition in the problem is the domain is an unit square. So, the equation is going to be written as

This equation we called as 2D – transient heat conduction equation.

  • 2D – steady state heat conduction equation is

Here in this equation, we see that there is no convection and time parameters or terms. Which means the data we are taking is at a particular point of interest.

Detailed derivation for 2D – steady state heat conduction equation in discretization form:

Program:

main code:

% steady state

clear all
close all
clc

% number of grid points in X-direction
nx = 10

% number of grid points in y-direction
ny = nx

x = linspace(0,1,nx)

% grid spacing in x-direction
dx = x(2) - x(1)

y = linspace(0,1,ny)

% grid spacing in y-direction
dy = y(2) - y(1)

z =fliplr(y)
f = 1.5
iter = 1

T =200*ones(nx,ny)

% boundary condition
T(2:ny-1,1) = 400
T(2:ny-1,nx) = 800
T(1,1:nx) = 600
T(ny,1:nx) = 900

% updating temperatures
Told = T
Tinitial =T

 error = 1e-4
 tol = 1e-5

for iterative_solver = [ 1 2 3 ]
        while(error > tol)
            for i = 2:nx-1
            for j = 2:ny-1
                tic
                  if(iterative_solver == 1)
                      T(i,j) = jacobi(T,Told,i,j)
                      toc
                      
                      tic
                  elseif(iterative_solver == 2)
                      T(i,j) = gauss_seidel(T,Told,i,j)
                      toc
                      
                      tic
                  else
                      T(i,j) = SOR(T,Told,i,j,f)
                      toc
               end
end
end
    
error = max(max(T-Told))
% updating temperature
Told = T
[xx , yy] = meshgrid(x,z)

%plot

if(iterative_solver == 1)
      figure(1)
         [c,d] = contourf(xx , yy ,T) 
         colormap(jet)
         colorbar
         title_name = sprintf(\'jacobi iteration number (steady state) = %d\', iter)
         title(title_name)
         xlabel(\'x-direction\')
         ylabel(\'y-direction\')
         clabel(c,d)
         pause(0.03)
elseif(iterative_solver == 2)
    figure(2)
         [c,d] = contourf(xx , yy ,T) 
         colormap(jet)
         colorbar
         title_name = sprintf(\'Gauss seidel iteration number (steady state) = %d\', iter)
         title(title_name)
         xlabel(\'x-direction\')
         ylabel(\'y-direction\')
         clabel(c,d)
         pause(0.03)
else
        figure(3)
         [c,d] = contourf(xx , yy ,T) 
         colormap(jet)
         colorbar
         title_name = sprintf(\'SOR iteration number (steady state) = %d\', iter)
         title(title_name)
         xlabel(\'x-direction\')
         ylabel(\'y-direction\')
         clabel(c,d)
         pause(0.03)
end
iter = iter+1

              end
            
           error = 1
           iter = 1
           T = Tinitial
           Told = Tinitial
end

Function code for Jacobi method:

function output = jacobi(T,Told,i,j)
output = 0.25*((Told(i,j-1) + Told(i,j+1) + Told(i-1,j) + Told(i+1,j)))
end

Function code for Gauss-Seidel method:

function output = gauss_seidel(T,Told,i,j)
output = 0.25*(T(i,j-1)+ Told(i,j+1) + T(i-1,j) + Told(i+1,j))
end

Function code for SOR method:

function output = SOR(T,Told,i,j,f)
output = ((1-f)*Told(i,j) + (f*0.25*(T(i,j-1) + Told(i,j+1) + T(i-1,j) + Told(i+1,j))))
end

Plots or Results:

Jacobi method:

Gauss-Seidel method:

SOR method:

  • Here, we see the obtained plots or graphs we can say that all three solutions are same and have a similar trend in their temperature gradients at different locations in the domain.
  • The point we need to observe is in how many iterations, the solution is going to converge.
  • If we compare the three plots, we can observe the difference in their no. of iterations. As

Jacobi Method converges in = 247 Iterations

Gauss-Seidel Method converges in = 131 Iterations

Successive-Over-Relaxation Method converges in = 30 Iterations

  • The phenomena here we observe is, SOR Method converges quickly because a factor called scaling factor.
  • The concept behind the scaling factor is to accelerate our rate of converge. In this concept we need to maintain scaling factor always upon one (over relaxation), because relaxation factor under one will cause slowness in converge of solution.
  • We need to choose over relaxation factor carefully, cause it may lead your solution unstable also.

Detailed derivation for 2D – Transient state heat conduction equation (Explicit) in discretization form:

  • As required in computational simulations, Explicit approach used for obtaining numerical approximations to the solutions of time dependent ordinary and partial differential equations.
  • The current value depends only on previous time-step data. Hence, we can just evaluate the current value at any location.

Program:

% calculating the 2D unsteady or transint heat conduction equation (explicit) using three different schemes

clear all
close all
clc


% number of grid points in X-direction
nx = 10

% number of grid points in y-direction
ny = nx

x = linspace(0,1,nx)

% grid spacing in x-direction
dx = x(2) - x(1)

y = linspace(0,1,ny)

% grid spacing in y-direction
dy = y(2) - y(1)

% convection term
c = 1e-4

%no.of time steps
nt = 200

% time step to be calculated
dt = 33

% courant number
k = (c*dt)/dx^2

T = 500*ones(nx,ny)

% boundary condition
T(2:ny-1,1) = 400
T(2:ny-1,nx) = 800
T(1,1:nx) = 600
T(ny,1:nx) = 900

% updating temperatures
Told = T
[xx , yy] = meshgrid(x,y)

tic

iter = 1
for p = 1:nt
    for i  = 2:nx-1
        for j = 2:ny-1
            T(i,j) = k*(Told(i-1,j) + Told(i,j-1) + Told(i+1,j) + Told(i,j+1)) + (1-4*k)*Told(i,j)
        end
    end
    
    error = max(abs(Told-T))
    
    % updating temperatures
    Told = T
    iter = iter+1
end
time = toc
figure(1)
[c,d] = contourf(xx,yy,T)
colorbar
colormap(jet)
clabel(c,d)
title_text=(sprintf(\'2D transient condition (Explicit method) \\n computation time: %0.4f \\n CFL number: %0.4f \\n iteration=%d\', time,k,iter))
title(title_text)
xlabel(\'x-direction\')
ylabel(\'y-direction\')

Plots or Results:

  • Here the plot shows that how the temperature gradient in the unit square domain is varying with different temperatures at their boundary’s. Make sure that parameters you are giving are in the range or not. If not the solution gets unstable by them.
  • It takes more time to compute the solution compared to implicit method.

Detailed derivation for 2D – Transient state heat conduction equation (Implicit) in discretization form:

  • As required in computational simulations, Implicit approach used for obtaining numerical approximations to the solutions of time dependent ordinary and partial differential equations.
  • Here, you neglect the dependence to previous time-step. Hence, you rely on current time-step data which you don’t know yet! Hence, you have to simultaneously solve for all the places.

Program:

% calculating the 2D unsteady or transint heat conduction equation (Implicit) using three different schemes

%%%%%% Jacobi Method %%%%%%

close all
clear all
clc

% length of domain
l=1;

% number of grid points in X-direction
nx = 10;

% number of grid points in y-direction
ny = nx;

x = linspace(0,l,nx);

% grid spacing in x-direction
dx = x(2) - x(1);

y = linspace(0,l,ny);

% grid spacing in y-direction
dy = y(2) - y(1);

% time step to be calculated
dt = 2000;

error = 1;
tol = 1e-4;

% thermal diffusivity of copper
alpha = 1.11e-4 ;

% boundary condition
T_l = 400;
T_r = 800;
T_t = 600;
T_b = 900;

T = ones(nx,ny);

T(2:ny-1,1) = T_l;
T(2:ny-1,nx) = T_r;
T(1,2:nx-1) = T_t;
T(nx,2:ny-1) = T_b;

% average temperatures at corners
T(1,1) = (T_t + T_l)/2;
T(nx,ny) =(T_r + T_b)/2;
T(1,ny) = (T_t + T_r)/2;
T(nx,1) = (T_l + T_b)/2;

% courant number calculation in two directions
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

[x , y] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;


iterative_solver = 1:3;

term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;

for iterative_solver = 1
    tic;
    if iterative_solver == 1
        jacobi_iter = 1;
        while(error>tol)
            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) = term1*T_prev_dt(i,j)+term2*h+term3*v;
                end
            end
            error = max(max(abs(Told-T)));
            
            % update temperatures
            Told = T;
            
            jacobi_iter = jacobi_iter+1;
        end
    end
            F = toc;

figure(1)
[c,d] = contourf(x,y,T)
clabel(c,d)
colormap(jet)
colorbar
title_text = sprintf(\'jacobi method (Implicit) - Transient  iteration number =%d, total time F=%fs\' ,jacobi_iter,F)
title(title_text)
xlabel(\'x-direction\')
ylabel(\'y-direction\')
fprintf(\'No. of iterations = %d\\n\', jacobi_iter)
end

%%%%% Gauss-Siedel Method %%%%%

% length of domain
l=1;

% number of grid points in X-direction
nx = 10;

% number of grid points in y-direction
ny = nx;

x = linspace(0,l,nx);

% grid spacing in x-direction
dx = x(2) - x(1);

y = linspace(0,l,ny);

% grid spacing in y-direction
dy = y(2) - y(1);

% time step to be calculated
dt = 2000;

error = 1;
tol = 1e-4;

% thermal diffusivity of copper
alpha = 1.11e-4 ;

% boundary condition
T_l = 400;
T_r = 800;
T_t = 600;
T_b = 900;

T = ones(nx,ny);

T(2:ny-1,1) = T_l;
T(2:ny-1,nx) = T_r;
T(1,2:nx-1) = T_t;
T(nx,2:ny-1) = T_b;

% average temperatures at corners
T(1,1) = (T_t + T_l)/2;
T(nx,ny) =(T_r + T_b)/2;
T(1,ny) = (T_t + T_r)/2;
T(nx,1) = (T_l + T_b)/2;

% courant number calculation in two directions
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

[x , y] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;


iterative_solver = 1:3;

term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;

for iterative_solver = 2
    tic;
    if iterative_solver == 2
        gs_iter = 2;
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    h = (T(i-1,j)+Told(i+1,j));
                    v = (T(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_prev_dt(i,j)+term2*h+term3*v;
                end
            end
            error = max(max(abs(Told-T)));
            
            % update temperatures
            Told = T;
            
            gs_iter = gs_iter+1;
        end
    end
            F = toc;

figure(2)
[c,d] = contourf(x,y,T)
clabel(c,d)
colormap(jet)
colorbar
title_text = sprintf(\'Gauss-Seidel method (Implicit) - Transient  iteration number =%d, total time F=%fs\' ,gs_iter,F)
title(title_text)
xlabel(\'x-direction\')
ylabel(\'y-direction\')
fprintf(\'No. of iterations = %d\\n\',gs_iter)
end


%%%%% Successive Over Relaxation Method %%%%%

% length of domain
l=1;

% number of grid points in X-direction
nx = 10;

% number of grid points in y-direction
ny = nx;

x = linspace(0,l,nx);

% grid spacing in x-direction
dx = x(2) - x(1);

y = linspace(0,l,ny);

% grid spacing in y-direction
dy = y(2) - y(1);

% time step to be calculated
dt = 2000;

error = 1;
tol = 1e-4;

% thermal diffusivity of copper
alpha = 1.11e-4 ;

% Scaling factor in SOR
alpha1 = 1.8;

% boundary condition
T_l = 400;
T_r = 800;
T_t = 600;
T_b = 900;

T = ones(nx,ny);

T(2:ny-1,1) = T_l;
T(2:ny-1,nx) = T_r;
T(1,2:nx-1) = T_t;
T(nx,2:ny-1) = T_b;

% average temperatures at corners
T(1,1) = (T_t + T_l)/2;
T(nx,ny) =(T_r + T_b)/2;
T(1,ny) = (T_t + T_r)/2;
T(nx,1) = (T_l + T_b)/2;

% courant number calculation in two directions
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

[x , y] = meshgrid(x,y);
Told = T;
T_prev_dt = Told;


iterative_solver = 1:3;

term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;

for iterative_solver = 3
    tic;
    if iterative_solver == 3
        SOR_iter = 3;
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    h = (alpha1*T(i-1,j)+(1-alpha1)*Told(i-1,j)+Told(i+1,j));
                    v = (alpha1*T(i,j-1)+(1-alpha1)*Told(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_prev_dt(i,j)*term1+(term2*h+term3*v);
                end
            end
            error = max(max(abs(Told-T)));
            
            % update temperatures
            Told = T;
            
            SOR_iter = SOR_iter+1;
        end
    end
            F = toc;

figure(3)
[c,d] = contourf(x,y,T)
clabel(c,d)
colormap(jet)
colorbar
title_text = sprintf(\'Successive Over Relaxation method (Implicit) - Transient  iteration number =%d, total time F=%fs\' ,SOR_iter,F)
title(title_text)
xlabel(\'x-direction\')
ylabel(\'y-direction\')
fprintf(\'No. of iterations = %d\\n\',SOR_iter)
end

Plots or Results:

Jacobi method:

Gauss-Seidel method:

SOR method:

  • While selecting the over relaxation factor we need to be careful. Because solution may become unstable.
  • Temperature gradients changes quickly when compared to the explicit method.
  • In comparison of Jacobi, gauss-seidel, and SOR methods,
  • Jacobi Method converges in = 178 Iterations in 0.008507 sec
  • Gauss-Seidel Method converges in = 97 Iterations in 0.007226 sec
  • Successive-Over-Relaxation Method converges in = 82 Iterations in 0.007001 sec.
  • If we increase the no. of elements, we can clearly seen the time difference between the three techniques.

Conclusion:

  • By comparing explicit and implicit schemes we can say that implicit method is more accurate and taking less time in converging the solutions.
  • Implicit method uses all the values in the exisitig domain.
  • On comparing all three iterative methods, we can say that the SOR method gives us the faster convergence due to the scaling factor.
  • Gauss-Seidel is better option over jacobi because it always uses new guess values resulting for a previous iteration for the next iteration. This phenomena making GS more accurate.
  • Inverse matrix is an intresting technique for improving the accuracy implicit method in converging the solution.
  • So in conclusion, using SOR method in steady and transient conduction equation leads to a faster convergance rate over other techniques.
  • For the uniformity, no. of elements in x and y direction taken as 10 in all conditions.
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.