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

Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

HHaji Mohamed29 Jan 2023

2D Steady Heat Conduction Equation

`(d^2T)/(dx^2)+(d^2T)/(dy^2)=0`

Objective:

Primary objective is to compare convergence rate i.e No. of iteration of steady and unsteady 2D heat conduction problem using Explicit and Implicit form of solutions for different iterative solvers: 1. Jacobi method, 2. Gauss Seidal method and 3. Successive over Relaxation (SOR) method.

The Matlab Code used for solving 2D Steady Heat Conduction Problem is given below:

 clear all;close all;clc;

L=1;
% Discritizing nodes
nx = 20;
ny = 20;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
k=2*((dx^2)+(dy^2))/((dx^2)*(dy^2));
% Assigning B.C's
Tx=ones(1,nx);Ty=ones(1,ny);
T=meshgrid(Tx,Ty);
T(:,1)=400;T(:,end)=800;T(1,:)=600;T(end,:)=900;
T(1,1)=500;T(1,end)=750;T(end,1)=650;T(end,end)=850;
T_old=T;
w=1.2; % Relaxation Factor
tol=1e-4; % Residual tolerance
error=9e9;

% for specifying Iterative solvers
% 1--> Jacobi || 2--> Gauss Seidal || 3--> SOR
iteration=3; 

% If the iterative solver is Jacobi
iter_no=1;

if iteration==1
while (error>tol)

    for i=2:nx-1
        
        for j=2:ny-1
            
            T(i,j) = ((T_old((i-1),j)+T_old((i+1),j))/(k*(dx^2)))+((T_old(i,(j-1))+T_old(i,(j+1)))/(k*(dy^2)));
                     
        end
    end
    
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;
    
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)
end
end

%  If the iterative solver is Gauss Seidal
if iteration==2
while (error>tol)

    for i=2:nx-1
        
        for j=2:ny-1
            
            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)));
                     
        end
    end
    
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;
    
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)
end
end

%  If the iterative solver is SOR
if iteration==3
while (error>tol)

    for i=2:nx-1
        
        for j=2:ny-1
            
            T(i,j) = ((1-w)*T_old(i,j))+w*(((T((i-1),j)+T((i+1),j))/(k*(dx^2)))+((T(i,(j-1))+T(i,(j+1)))/(k*(dy^2))));
                     
        end
    end
    
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;
    
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)
    
end
end

The Outputs obtained with this code is shown in following figures:

From the above figures we can infer that

jacobi iteration> Gauss seidal iteration> SOR Iteration

2D Unsteady Heat Conduction is given below:

`(dT)/dt=alpha*((d^2T)/dx^2+(d^2T)/dy^2)`

The Code used for solving 2D unsteady heat conduction problem in EXPLICIT Form is shown below:

clear all;close all;clc;

set(groot,'DefaultLineLineWidth',3);
set(0,'DefaultAxesFontSize',20);
set(0,'DefaultAxesGridAlpha',0.6);
set(0,'DefaultAxesMinorGridAlpha',0.3);
set(0,'DefaultAxesXminorTick','On');
set(0,'DefaultAxesYminorTick','On');
set(gcr,'YDir','rev');

% Given Space and Time info
L=1;Time=0.4;
% Discritizing nodes
nx = 20;ny = 20;nt=1200;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
t=linspace(0,Time,nt);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=t(2)-t(1);
a=1.4; % thermal diffusivity
k1=a*(dt/(dx^2));
k2=a*(dt/(dy^2));
% Assigning B.C's
Tx=ones(1,nx);Ty=ones(1,ny);
T=meshgrid(Tx,Ty);
T(:,1)=400;T(:,end)=800;T(1,:)=600;T(end,:)=900;
T(1,1)=500;T(1,end)=750;T(end,1)=650;T(end,end)=850;
T_old=T;
w=1.2; % Relaxation factor
tol=1e-4; % Residual tolerance
error=9e9;

% for specifying Iterative solvers
% 1--> Jacobi || 2--> Gauss Seidal || 3--> SOR
iteration=3;

% If the iterative solver is Jacobi
iter_no=1;

if iteration==1
    % Time loop
for it=1:nt    
    % Residual loop
while (error>tol)
    % Space loop
    for i=2:nx-1
        
        for j=2:ny-1
            
            T(i,j) = T_old(i,j) + (k1*(T_old((i-1),j)-2*(T_old(i,j))+T_old((i+1),j)))+(k2*(T_old(i,(j-1))-2*(T_old(i,j))+T_old(i,(j+1))));
                     
        end
    end
    
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;
    
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)
end
    T_old=T;
end
end 

%  If the iterative solver is Gauss Seidal
if iteration==2
    % Time loop
for it=1:nt
    % Residual loop
while (error>tol)
    % Space loop
    for i=2:nx-1
        
        for j=2:ny-1
            
            T(i,j) = T(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))));
                     
        end
    end
    
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;
    
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)
end
    T_old=T;
end
end

%  If the iterative solver is SOR
if iteration==3
    %  Time loop
for it=1:nt
    %  Residual loop
while (error>tol)
    % Space loop
    for i=2:nx-1
        
        for j=2:ny-1
            
            T(i,j) = ((1-w)*T_old(i,j))+w*(T(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)))));
                     
        end
    end
    
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;
    
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)
end
    T_old=T;
end
end

The Corresponding output figures obtained from the above is code are shown below:

Output Plot for Jacobi Iterative Solver in Explicit form:

Output Plot for Gauss Seidal Iterative Solver in Explicit form:

Output Plot for Successive Over Relaxation (SOR) Iterative Solver in Explicit form:

From the above figures we can infer that by solving 2D unsteady Heat conduction problem in explicit form, we have

jacobi iteration> Gauss seidal iteration> SOR Iteration

The Code used for solving 2D unsteady heat conduction problem in IMPLICIT Form is shown below:

clear all;close all;clc;

% Given Space and Time info
L=1;Time=0.4;
% Discritizing nodes
nx = 20;ny = 20;nt=1500;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
t=linspace(0,Time,nt);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=t(2)-t(1);
a=1.4;  % thermal diffusivity
% Assigning B.C's
Tx=ones(1,nx);Ty=ones(1,ny);
T=meshgrid(Tx,Ty);
T(:,1)=400;T(:,end)=800;T(1,:)=600;T(end,:)=900;
T(1,1)=500;T(1,end)=750;T(end,1)=650;T(end,end)=850;
T_old=T;
T_prev=T;
w=1.2;  % Relaxation factor
tol=1e-4;  % Residual tolerance
error=9e9;
k1=a*dt/dx^2;
k2=a*dt/dy^2;
term1=1/(1+(2*k1)+(2*k2));
term2=k1*term1;
term3=k2*term1;

% for specifying Iterative solvers
% 1--> Jacobi || 2--> Gauss Seidal || 3--> SOR
iteration=1;

% If the iterative solver is Jacobi
iter_no=1;

if iteration==1
    % Time loop
for it=1:nt
  error=9e9; % Specifying Error Randomly
  % Residual loop
while (error>tol)
    % Space loop
    for i=2:nx-1
        for j=2:ny-1
            T(i,j)=(term1*T_prev(i,j))+((T_old(i-1,j)+T_old(i+1,j))*term2)+((T_old(i,j-1)+T_old(i,j+1))*term3);
        end
    end
    E=abs(T_old-T);
    error=max(max(E)); %  Computing error after iteration
    T_old=T;           %  Updating Time for next iteration
    iter_no=iter_no+1; % Counting no. of iteration

end
T_prev=T; % Updating Time for next time step
end
    % Plotting Contour plot
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)

end

%  If the iterative solver is Gauss Seidal
if iteration==2
    % Time loop
for it=1:nt
  error=9e9;
      % Residual loop
while (error>tol)
    % Space loop
    for i=2:nx-1
        for j=2:ny-1
            T(i,j)=(term1*T_prev(i,j))+((T(i-1,j)+T(i+1,j))*term2)+((T(i,j-1)+T(i,j+1))*term3);
        end
    end
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;

end
T_prev=T;
end
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)

end

%  If the iterative solver is SOR
if iteration==3
    %  Time loop
for it=1:nt
  error=9e9;
      %  Residual loop
while (error>tol)
    % Space loop
    for i=2:nx-1
        for j=2:ny-1
            T(i,j)=(1-w)*T_old(i,j)+(w*((term1*T_prev(i,j))+((T(i-1,j)+T(i+1,j))*term2)+((T(i,j-1)+T(i,j+1))*term3)));
        end
    end
    E=abs(T_old-T);
    error=max(max(E));
    T_old=T;
    iter_no=iter_no+1;

end
T_prev=T;
end
    colormap(jet);
    A=contourf(x,y,T);
    xlabel('X');
    ylabel('Y');
    tit=sprintf('Iteration No= %d',iter_no);
    title(tit);
    clabel(A);
    set(gca,'YDir','reverse');
    pause(0.03)

end

The Corresponding output figures obtained from the above is code are shown below:

Output Plot for Jacobi Iterative Solver in Implicit form:

Output Plot for Gauss Seidal Iterative Solver in Implicit form:

Output Plot for Successive Over Relaxation (SOR) Iterative Solver in Implicit form:

From the above figures we can infer that by solving 2D unsteady Heat conduction problem in explicit form, we have

jacobi iteration> Gauss seidal iteration> SOR Iteration

Summary:

From the obtained output plots by solving 2D Steady and Unsteady Heat Conduction Equation by different iterative solvers, we can infer the following:

  • Convergence rate (No. of iteration) is higher for Jacobi iterative solver followed by Gauss Seidal and then SOR iterative solver. This is because in jacobi iterative solver, the value at (i,j) has not been using the updated values at next space loop whereas in Gauss Seidal and SOR the value at (i,j) has been using the updated values for the further iterations, this makes the solution to converge faster in Gauss Seidal and SOR than that of in Jacobi iterative solver.

    jacobi iteration> Gauss seidal iteration> SOR Iteration

  • Relaxation factor (`omega`), Thermal Diffusivity (`alpha`) and CFL term plays a vital role in convergence for the Unsteady solution where those parameters are defined in code as follow for stable convergence:Thermal Diffusivity (`alpha`) = 1.4Relaxation factor (`omega`) = 1.2CFL term (k1 + k2) = 0.26, which is `le`0.5 for stable convergence

For solving 2D Steady Heat Conduction equation, the number of iteration taken be different iterative solver is mentioned below:

  1. Jacobi Iterative Solver = 869
  2. Gauss Seidal Iterative Solver = 463
  3. Successive Over Relaxation (SOR) = 318

For solving 2D Unsteady Heat Conduction equation in Explicit form, the number of iteration taken be different iterative solver is mentioned below:

  1. Jacobi Iterative Solver = 1248
  2. Gauss Seidal Iterative Solver = 861
  3. Successive Over Relaxation (SOR) = 660

For solving 2D Unsteady Heat Conduction equation in Implicit form, the number of iteration taken be different iterative solver is mentioned below:

  1. Jacobi Iterative Solver = 9857
  2. Gauss Seidal Iterative Solver = 7623
  3. Successive Over Relaxation (SOR) = 5771
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.