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

Steady and Transient 2D Heat Conduction Equation

AAntim Gupta27 Jan 2021
MATLAB

Aim:

The major objective of this project was to solve the Steady and Transient 2D Heat Conduction Equation by using iterative solver such as Jacobi. gauss Seidel and Sucessive Over-Relaxation for both Implicit and Explicit Schemes using MATLAB.

Introduction:

The Heat Conduction Equation is a Partial Differential Equation that describes how the distribution of some quantity (Heat) evolves over time in a solid domain, as it spontaneous flows from places where it is higher towards place where its lower.

If the Temperature is a function of all three coordinates and the time, that is, T=T(x,y,z,t)

Where `alpha` is the thermal diffusivity of the material that tell us about the how fast the heat flow inside the medium. `alpha=k/(rhoc)(m^2/s)`

3D Steady State with Heat Generation

3D Steady State without Heat Generation

Laplace Equation

Steady State and Transient Equation

steady State means if the Properties changes with respect to time or the Time derivative term is envolved in the Equation to descrtize that Equation is generally called as Steady State heat transfer equation whereas if there is envolvement of time derivative or the property changes with respect to time those equation we call them as Transient equation or Unsteady State Heat Transfer Equation of Conduction.

Algorithm for solving this problem.

Problem Statement

Computation Area is Square

The Boundary Condition are:

  • Left=400k
  • Right=600k
  • Top=800k
  • Bottom=900k

Initial Temperature is 300k

Schemes used for solving this problem is

  • Jacobi Method
  • Gauss-Seidel
  • Successive Over Relaxation

The central difference scheme was chosen in order to calculate more accuarte result for space derivatives whereas for time derivative the forward scheme was choosen.

2D steady State Heat Conduction

3D Transient Heat Conduction Equation is represented by:-

The Pde's can be reduced to 2D Steady State Heat Conduction Equation.

Discritization using Central Difference Scheme for 2D Steady State Heat Conduction Equation.

The the Steady State Equation , the amount of heat entering any region of an object is equal to the amount of heat coming out.

If we further more discritize it we get,

`T_(i,j)=0.25(T_(i-1,j)+T_(i+1,j)+T_(i,j+1)+T_(i,j-1))`

Jacobi Method

%****************Antim Gupta************
%********2D Conduction Jacobbi Method**********
clc
clear all
close all

%Input Paramters

nx=100;% Number of Grid Points in x
ny=100;% Number of Grid Points in y

x=linspace(0,1,nx);% Method to Create The Grid Points in x
y=linspace(0,1,ny);% Method to Create The Grid Points in y

dx=x(2)-x(1);
dy=y(2)-y(1);

%Defining absolute error criteria 
error=9e9;
tolerance=1e-4;

%Defining the Boundary Condition 
T_L=400; %Left side of domain should be at 400K
T_R=800; %Right side of domain should be at 800K
T_T=600; %Top side of domain should be at 600K
T_B=900; %Botom side of domain should be at 900K

T=300*ones(nx,ny); %Defining the Inital Condition
T(2:ny-1,1)=T_L;
T(2:ny-1,nx)=T_R;
T(1,2:nx-1)=T_T;
T(ny,2:ny-1)=T_B;

% Now we have to assign the T value at corner
T(1,1)=(T_L+T_T)/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;

%keeping the Back Up of T 
T_ini=T;
T_old=T;

%Calculation of 2D Steady State Heat Conduction Equation by Jaccobi Method
is=1;

if is==1;
    ji=1;
    
    while(error>tolerance)
        
        for i=2:nx-1
            
            for j=2:ny-1
                T(i,j)=0.25*(T_old(i-1,j)+T_old(i+1,j)+T_old(i,j+1)+T_old(i,j-1));
            end
        end
        error=max(max(abs(T_old-T)));
        T_old=T;
        ji=ji+1;
           
    end
 
end

% Code for ploting the Values of T
figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')
title(sprintf('No. of Jacobi Iterations(Implicit)= %d',ji));
pause(0.03)

Fig 1.1 Steady State- Implicit Scheme using Jacobbi Method

Gauss-Seidel Method

%****************Antim Gupta************
%********2D Conduction Gauss Seidel Method**********
clc
clear all
close all

%Input Paramters

nx=100;% Number of Grid Points in x
ny=100;% Number of Grid Points in y

x=linspace(0,1,nx);% Method to Create The Grid Points in x
y=linspace(0,1,ny);% Method to Create The Grid Points in y

dx=x(2)-x(1);
dy=y(2)-y(1);

%Defining absolute error criteria 
error=9e9;
tolerance=1e-4;

%Defining the Boundary Condition 
T_L=400; %Left side of domain should be at 400K
T_R=800; %Right side of domain should be at 800K
T_T=600; %Top side of domain should be at 600K
T_B=900; %Botom side of domain should be at 900K

T=300*ones(nx,ny); %Defining the Inital Condition
T(2:ny-1,1)=T_L;
T(2:ny-1,nx)=T_R;
T(1,2:nx-1)=T_T;
T(ny,2:ny-1)=T_B;

% Now we have to assign the T value at corner
T(1,1)=(T_L+T_T)/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;

%keeping the Back Up of T 
T_ini=T;
T_old=T;

%Calculation of 2D Steady State Heat Conduction Equation by Gauss Seidel Method
is=2;

if is==2;
    gs=1;
    
    while(error>tolerance)
        
        for i=2:nx-1
            
            for j=2:ny-1
                T(i,j)=0.25*(T(i-1,j)+T_old(i+1,j)+T_old(i,j+1)+T(i,j-1));
            end
        end
        error=max(max(abs(T_old-T)));
        T_old=T;
        gs=gs+1;
           
    end
 
end

% Code for ploting the Values of T
figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')
title(sprintf('No. of Gauss-Seidel(Implicit)= %d',gs));
pause(0.03)

Fig 1.2 Steady State- Implicit Scheme using Gauss Seidel Method

Successive Over Relaxation

%****************Antim Gupta************
%********2D Conduction SOR Method**********
clc
clear all
close all

%Input Paramters

nx=100;% Number of Grid Points in x
ny=100;% Number of Grid Points in y

x=linspace(0,1,nx);% Method to Create The Grid Points in x
y=linspace(0,1,ny);% Method to Create The Grid Points in y

dx=x(2)-x(1);
dy=y(2)-y(1);

%Defining absolute error criteria 
error=9e9;
tolerance=1e-4;
omega=1.2;

%Defining the Boundary Condition 
T_L=400; %Left side of domain should be at 400K
T_R=800; %Right side of domain should be at 800K
T_T=600; %Top side of domain should be at 600K
T_B=900; %Botom side of domain should be at 900K

T=300*ones(nx,ny); %Defining the Inital Condition
T(2:ny-1,1)=T_L;
T(2:ny-1,nx)=T_R;
T(1,2:nx-1)=T_T;
T(ny,2:ny-1)=T_B;

% Now we have to assign the T value at corner
T(1,1)=(T_L+T_T)/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;

%keeping the Back Up of T 
T_ini=T;
T_old=T;

%Calculation of 2D Steady State Heat Conduction Equation by Successive Over Relaxation Method
is=3;

if is==3;
    sor=1;
    
    while(error>tolerance)
        
        for i=2:nx-1
            
            for j=2:ny-1
                T(i,j)=omega*(0.25*(T(i-1,j)+T_old(i+1,j)+T_old(i,j+1)+T(i,j-1)))+(1-omega)*T_old(i,j);
            end
        end
        error=max(max(abs(T_old-T)));
        T_old=T;
        sor=sor+1;
           
    end
 
end

% Code for ploting the Values of T
figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')
title(sprintf('No. of SOR(Implicit)= %d',sor));
pause(0.03)

Fig 1.3 Steady State- Implicit Scheme using SOR Method

Heat Conduction Transient Analysis

Heat Equation in 3D is given by:

Transient 2D Heat Conduction Equation is given by:

Discritization using Implicit Scheme

`(T_p^(n+1)-T_p^n)/dt=alpha[((T_l-2T_p+T_r)/dx^2)+(T_T-2T_p+T_B)/dy^2)]`

If we further more simplify it we get:

`Tp^(n+1)=1/(1+K1+K2)+(T_p^n+K1(T_L+T_R)^(n+1)+K2(T_T+T_B)^(n+1))`

`K1=alphadt/dx^2``K2=alpha*dt/dy^2`

Successive over Relaxation

%****************Antim Gupta************
%********Transient Heat Conduction SOR Method**********
clc
clear all
close all

%Input Paramters

nx=10;% Number of Grid Points in x
ny=10;% Number of Grid Points in y
nt=1000; 

x=linspace(0,1,nx);% Method to Create The Grid Points in x
y=linspace(0,1,ny);% Method to Create The Grid Points in y

dx=x(2)-x(1);
dy=y(2)-y(1);

%Defining absolute error criteria 
error=9e9;
tolerance=1e-4;
dt=1e-3; 
omega=1.1;

%Defining the Boundary Condition 
T_L=400; %Left side of domain should be at 400K
T_R=800; %Right side of domain should be at 800K
T_T=600; %Top side of domain should be at 600K
T_B=900; %Botom side of domain should be at 900K

T=300*ones(nx,ny); %Defining the Inital Condition
T(2:ny-1,1)=T_L;
T(2:ny-1,nx)=T_R;
T(1,2:nx-1)=T_T;
T(ny,2:ny-1)=T_B;

% Now we have to assign the T value at corner
T(1,1)=(T_L+T_T)/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;

%keeping the Back Up of T 
T_ini=T;
T_old=T;
T_end=T;

%Calculation of 2D Transient Heat Conduction Equation by GS
k1=1.1*(dt/(dx^2));
k2=1.1*(dt/(dy^2));

ct=1;   
gs=1;
for k=1:nt
    error=9e9;
    
    while(error>tolerance)
        
        for i=2:nx-1
            
            for j=2:ny-1
                term1=1/(1+2*k1+2*k2);
                term2=k1*term1;
                term3=k2*term1;
                
                T(i,j)=((T_ini(i,j)*term1))+(term2*(T(i-1,j)+T_old(i+1,j)))+(term3*(T(i,j+1)+T_old(i,j-1)));
            T_end(i,j)=((1-omega)*T_old(i,j))+(T(i,j)*omega);
            end
        end
        error=max(max(abs(T_old-T)));
        T_old=T_end;
        gs=gs+1;
           
    end
    T_ini=T_end;
    
    % Code for ploting the Values of T
 figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')
title(sprintf('No. of SOR(Implicit)= %d',gs));
pause(0.00000000000000003)
M(ct)=getframe(gcf);
    ct=ct+1;

end
figure(3)
movie(M)
videofile=VideoWriter('sor.avi','Uncompressed Avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

VID 2.1 UnSteady State- Implicit Scheme using SOR Method

Gauss Seidel Method

%****************Antim Gupta************
%********Transient Heat Conduction Gauss Sidel Method**********
clc
clear all
close all

%Input Paramters

nx=10;% Number of Grid Points in x
ny=10;% Number of Grid Points in y
nt=1000; 

x=linspace(0,1,nx);% Method to Create The Grid Points in x
y=linspace(0,1,ny);% Method to Create The Grid Points in y

dx=x(2)-x(1);
dy=y(2)-y(1);

%Defining absolute error criteria 
error=9e9;
tolerance=1e-4;
dt=1e-3; 
omega=1.1;

%Defining the Boundary Condition 
T_L=400; %Left side of domain should be at 400K
T_R=800; %Right side of domain should be at 800K
T_T=600; %Top side of domain should be at 600K
T_B=900; %Botom side of domain should be at 900K

T=300*ones(nx,ny); %Defining the Inital Condition
T(2:ny-1,1)=T_L;
T(2:ny-1,nx)=T_R;
T(1,2:nx-1)=T_T;
T(ny,2:ny-1)=T_B;

% Now we have to assign the T value at corner
T(1,1)=(T_L+T_T)/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;

%keeping the Back Up of T 
T_ini=T;
T_old=T;
T_end=T;

%Calculation of 2D Transient Heat Conduction Equation by GS
k1=1.1*(dt/(dx^2));
k2=1.1*(dt/(dy^2));

ct=1;   
ji=1;
for k=1:nt
    error=9e9;
    
    while(error>tolerance)
        
        for i=2:nx-1
            
            for j=2:ny-1
                term1=1/(1+2*k1+2*k2);
                term2=k1*term1;
                term3=k2*term1;
                
                T(i,j)=((T_ini(i,j)*term1))+(term2*(T(i-1,j)+T_old(i+1,j)))+(term3*(T(i,j+1)+T_old(i,j-1)));
            end
        end
        error=max(max(abs(T_old-T)));
        T_old=T;
        ji=ji+1;
           
    end
    T_ini=T;
    
    % Code for ploting the Values of T
 figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')
title(sprintf('No. of Gauss Seidel(Implicit)= %d',ji));
pause(0.00000000000000003)
M(ct)=getframe(gcf);
    ct=ct+1;

end
figure(3)
movie(M)
videofile=VideoWriter('GauusSeidel.avi','Uncompressed Avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

VID 2.2 UnSteady State- Implicit Scheme using Gauss Seidel Method

Jacobbi Method

%****************Antim Gupta************
%********Transient Heat Conduction Jacobbi Method**********
clc
clear all
close all

%Input Paramters

nx=10;% Number of Grid Points in x
ny=10;% Number of Grid Points in y
nt=1000; 

x=linspace(0,1,nx);% Method to Create The Grid Points in x
y=linspace(0,1,ny);% Method to Create The Grid Points in y

dx=x(2)-x(1);
dy=y(2)-y(1);

%Defining absolute error criteria 
error=9e9;
tolerance=1e-4;
dt=1e-3; 

%Defining the Boundary Condition 
T_L=400; %Left side of domain should be at 400K
T_R=800; %Right side of domain should be at 800K
T_T=600; %Top side of domain should be at 600K
T_B=900; %Botom side of domain should be at 900K

T=300*ones(nx,ny); %Defining the Inital Condition
T(2:ny-1,1)=T_L;
T(2:ny-1,nx)=T_R;
T(1,2:nx-1)=T_T;
T(ny,2:ny-1)=T_B;

% Now we have to assign the T value at corner
T(1,1)=(T_L+T_T)/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;

%keeping the Back Up of T 
T_ini=T;
T_old=T;

%Calculation of 2D Transient Heat Conduction Equation by Jaccobi Method
k1=1.1*(dt/(dx^2));
k2=1.1*(dt/(dy^2));

ct=1;   
ji=1;
for k=1:nt
    error=9e9;
    
    while(error>tolerance)
        
        for i=2:nx-1
            
            for j=2:ny-1
                term1=1/(1+2*k1+2*k2);
                term2=k1*term1;
                term3=k2*term1;
                
                T(i,j)=((T_ini(i,j)*term1))+(term2*(T_old(i-1,j)+T_old(i+1,j)))+(term3*(T_old(i,j+1)+T_old(i,j-1)));
            end
        end
        error=max(max(abs(T_old-T)));
        T_old=T;
        ji=ji+1;
           
    end
    T_ini=T;
    
    % Code for ploting the Values of T
 figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')
title(sprintf('No. of Jacobi Iterations(Implicit)= %d',ji));
pause(0.00000000000000003)
M(ct)=getframe(gcf);
    ct=ct+1;

end
figure(3)
movie(M)
videofile=VideoWriter('Jacobbi.avi','Uncompressed Avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

VID 2.3 UnSteady State- Implicit Scheme using Jacobbi Method

Explicit Scheme

Heat Conduction Equation in 3D

If we disscritize the equation for 2d

`(T_p^(n+1)-T_p)/dt=alpha[(T_l*2T_p+T_R)/dx^2+(T_T-2T_p+T_b)/dy^2]`

If we further more discritize it:

`T_(i,j)=T_(i,j)+K1(T_(i-1,j)-2T_(i,j)+K1(T_(i+1,j))+K2(T_(i,j+1)-2T_(i,j)+K1(T_(i,j-1))`

`K1=alphadt/dx^2``K2=alpha*dt/dy^2`

%****************Antim Gupta************
%********Transient Heat Conduction Explicit Scheme**********
clc
clear all
close all

%Input Paramters

nx=10;% Number of Grid Points in x
ny=10;% Number of Grid Points in y
nt=1000; 

x=linspace(0,1,nx);% Method to Create The Grid Points in x
y=linspace(0,1,ny);% Method to Create The Grid Points in y

dx=x(2)-x(1);
dy=y(2)-y(1);

%Defining absolute error criteria 
error=9e9;
tolerance=1e-4;
dt=1e-3; 
omega=1.1;

%Defining the Boundary Condition 
T_L=400; %Left side of domain should be at 400K
T_R=800; %Right side of domain should be at 800K
T_T=600; %Top side of domain should be at 600K
T_B=900; %Botom side of domain should be at 900K

T=300*ones(nx,ny); %Defining the Inital Condition
T(2:ny-1,1)=T_L;
T(2:ny-1,nx)=T_R;
T(1,2:nx-1)=T_T;
T(ny,2:ny-1)=T_B;

% Now we have to assign the T value at corner
T(1,1)=(T_L+T_T)/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;

%keeping the Back Up of T 
T_ini=T;
T_old=T;
T_end=T;

%Calculation of 2D Transient Heat Conduction Equation by Explicit
k1=1.1*(dt/(dx^2));
k2=1.1*(dt/(dy^2));
is=1;

fc=1;
if fc==1
    gs=1;
for k=1:nt
    error=9e9;
    
    while(error>tolerance)
        
        for i=2:nx-1
            
            for j=2:ny-1
                t1=T_old(i,j);
                t2=k1*(T_old(i+1,j)-2*T_old(i,j)+(T_old(i-1,j)));
                t3=k2*(T_old(i,j+1)-2*T_old(i,j)+(T_old(i,j-1)));
                T(i,j)=t1+t2+t3;
            end
        end
        error=max(max(abs(T_old-T)));
        T_old=T;
        gs=gs+1;
           
    end
end
end 
 figure(1)

contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')

xlabel('X -Axis')
ylabel('Y -Axis')
title(sprintf('No. of Transient(Explicit)= %d',gs));
pause(0.00000000000000003)

Img 3.1 UnSteady State- Explict Scheme using SOR Method

Conclusion

Types of Schemes & Analysis Methods Used Number of Iterations (for nx=ny=100)
Jacobbi 15940
IMPLICIT-Steady State Analysis Gauss Seidel 8672
Successive Over Relaxation 6053
Number of Iterations (for nx=ny=10)
Jacobbi 3471
IMPLICIT-unSteady State Analysis Gauss Seidel 3221
Successive Over Relaxation 2890
Explicit Analysis(Transient) 1941

From the above table it can be concluded that:

  • Jacobi>Gauss Seidel>Successive Over Relaxation according to the number of Iteration computed for all method and schemes
  • It can also be conclused that steady states is more sucessful as it takes less computational time because it dont have any time derivative term in it. & along with less number of iterations.
  • In case if we talk about Explicit or Implcit, Then Explicit is more effective as it takes the least number of iteration in order to provide the accurate results.
  • SOR or Successive Over Relaxation is more effective iteration as it consumes less time and iterations whether we talk about Explicit ot Implicit.
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.