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

VVaibhav Carpenter17 Feb 2023
MATLAB

Aim:- To write the code which solves the 2D heat conduction equation by using the point iterative techniques by implementing the following methods.

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

Introduction :- We know the general 2D heat conduction equations-

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

where T is temperature,t stands for time.

For steady-state Time derivative will be zero hence the equation will be;

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

since `(d^2t)/dx^2 + (d^2T)/dy^2 = 0`, it has elliptic nature properties which we an discretize them using central differencing scheme

`(T_(i+1,j)-2T_(i,j)+T_(i-1,j))^n/(Deltax^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))^n/(Deltay^2)=0`

`T_(i,j)*2(Deltax^2+Deltay^2)/(Deltax^2Deltay^2)=(T_(i-1,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2)`

Let -`k=2(Deltax^2+Deltay^2)/(Deltax^2Deltay^2)`

`T_(i,j)=(1/k)(T_(i-1,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2)`

now what happens when represent them as jacobi

`T_(i,j)^(n+1)=(1/k)((T_(i-1,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2))^n`

for Gauss seidel

`T_(i,j)^(n+1)=(1/k)((T_(i-1,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2))^(n+1)`

and SOR is enhanced gauss-seidel method where nth time step and n+1 time step is all same.

`T_(i,j)^(n+1)=T_(i,j)^n(1-omega)+omega*(1/k)((T_(i-1,j)+T_(i+1,j))/(Deltax^2)+(T_(i,j+1)+T_(i,j-1))/(Deltay^2))^(n+1)`

where omega is the relaxation factor as 1.4.

now here we begin the coding part which is done in MATLAB where these all terms are considered by applying for loop and while loop.

Steady state heat conduction :-

clear all
close all
clc

% Initial conditions

nx = 10; % number of nodes 
ny = 10; % number of nodes
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha=0.3; % thermal diffusivity

% BC
t=300*ones(nx,ny);
t(:,1) = 400; % Left Temperature
t(1,:) = 900; % Top Temperature
t(end,:) = 600; % Bottom Temperature
t(:,end) = 800; % Right Temperature

% Editing Corner Values 
t(1,1) = (400+900)/2;
t(end,1) = (400 + 600)/2;
t(1,end) = (900+800)/2;
t(end,end)= (800+600)/2;
told = t;

k = 2*(dx^2+dy^2)/(dx^2*dy^2);
error =100;
iterative_solver=3;
tol=1e-4;

% jacobi method ( steady state form)

if iterative_solver==1
   jacobi_iter = 1;
 while(error>tol)
    for i =2:(nx-1)
        for j =2:(ny-1)
            t(i,j)=(((told(i-1,j)+told(i+1,j))/(dx^2))+((told(i,j+1)+told(i,j-1))/(dy^2)))/k;      
        end
    end  
    error=max(max(abs(told-t)));
    told=t;
    jacobi_iter=jacobi_iter+1;    
 end
 contourf(x,y,t,'ShowText','on')
 colormap(jet)
 colorbar
 %set(gca,'YDir','reverse')
 xlabel('X axis')
 ylabel('t')
 text1 = sprintf('number of iteration = %d', jacobi_iter);
 title({'Steady state heat conduction profile by jacobi'; text1})
end

 % guass seidel method
 
if iterative_solver==2
   gs_iter = 1;
 while(error>tol)
    for i =2:(nx-1)
        for j =2:(ny-1)
            t(i,j)=(((t(i-1,j)+told(i+1,j))/(dx^2))+((told(i,j+1)+t(i,j-1))/(dy^2)))/k;      
        end
    end  
    error=max(max(abs(told-t)));
    told=t;
    gs_iter=gs_iter+1;    
 end
 contourf(x,y,t,'ShowText','on')
 colormap(jet)
 colorbar
 xlabel('X axis')
 ylabel('t')
 text1 = sprintf('number of iteration = %d', gs_iter);
 title({'Steady state heat conduction profile by Guess Seidel'; text1})
end

% successive over relaxation iterative method

if iterative_solver==3                                                                                                                                                                      
   sor_iter = 1;
   omega=1.4;
     while(error>tol)
        for i =2:(nx-1)
             for j =2:(ny-1)
                 term1 = (t(i-1,j) + told(i+1,j))/(dx^2);
                 term2 = (told(i,j+1) + t(i,j-1))/(dy^2);
                 t(i,j) = ((told(i,j)*(1-omega))) + (omega*(term1+term2))/k;
             end
        end  
        error=max(max(abs(told-t)));
        told=t;
        sor_iter= sor_iter+1;
    end
    contourf(x,y,t,'ShowText','on')
    colormap(jet)
    colorbar
    xlabel('X axis')
    ylabel('t')
    text1 = sprintf('number of iteration = %d', sor_iter);
    title({'Steady state heat conduction profile by SOR'; text1})
end

Explanation:-

  • First provide intial conditions and boundry conditions.
  • nx and ny is number of nodes in x and y direction,dx and dy is distance b/w nodes, x and y is the distance of nodes inside.
  • t is temperature provide as a boundry condition as in the question.
  • Iterative_solver provide for run methods respectively by providing input same have methods.
  • First 'if' loop belongs to run methods.
  • then 'while' loop is for Iterative_solver used for represent method then jacobi_iter is use for start caculation of iteration.
  • 'for' loop provide i and j values as its 2D heat equation sto be solved.
  • In this loop equation to be solved as required as jacobi,guass seidal,successive over relaxation methods.
  • equations is provide before the matlab codes.
  • the close loops after closing loops error to calculated by using difference of told and t.
  • then updating told=t
  • then contourf command used for getting contours.
  • same procedure in later cases and the equation is provide before the matlab codes.

Output:-

By jacobi method

By guass Seidal method

By successive over relaxation method

From above results we can see that the SOR method gives a better in terms of iterations carried out by them.

Transient (Unsteady) state simulation:

Moving on to the Transient State simulation where we are going to work on Implicit and explicit schemes.

1. Implicit schemes

2. Explicit scheme

Equations for transient cases use

for jacobi

`T_(i,j)^(n+1)=(T_(i,j)^n+k1(T_(i+1,j)+T_(i+1,j))^n+k2(T_(i,j+1)+T_(i,j-1))^n)/(1+2k1+2k2)`

for guass seidal

`T_(i,j)^(n+1)=(T_(i,j)^n+k1(T_(i+1,j)+T_(i+1,j))^(n+1)+k2(T_(i,j+1)+T_(i,j-1))^(n+1))/(1+2k1+2k2)`

for successive over relaxation

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

Transient state heat conduction by implicit method:-

clear all
close all
clc

% Initial conditions

nx = 10; % number of nodes 
ny = 10; % number of nodes
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha=0.3; % thermal diffusivity

% BC
t=300*ones(nx,ny);
t(:,1) = 400; % Left Temperature
t(1,:) = 900; % Top Temperature
t(end,:) = 600; % Bottom Temperature
t(:,end) = 800; % Right Temperature

% Editing Corner Values 
t(1,1) = (400+900)/2;
t(end,1) = (400 + 600)/2;
t(1,end) = (900+800)/2;
t(end,end)= (800+600)/2;
told = t;
tprev=t;
dt = 1e-2;

k_1 = alpha*((dt)/(dx)^2);
k_2 = alpha*((dt)/(dy)^2);

error = 100;
iterative_solver=1;
tol=1e-4; 

% jacobi method ( steady state form)
nt=100; % time step
if iterative_solver==1
   jacobi_iter = 1;    
    for z=1:nt
        error=100;
        while(error>tol)
            for i =2:(nx-1)
                  for j =2:(ny-1)
                     term1=(1+(2*k_1)+(2*k_2))^-1;
                     term2=k_1*term1;
                     term3=k_2*term1;
                     H=(told(i-1,j)+told(i+1,j));
                     V=(told(i,j-1)+told(i,j+1));
                     t(i,j)=((tprev(i,j)*term1)+(H*term2)+(V*term3));     
                   end
             end
             error=max(max(abs(told-t)));
             told=t;
             jacobi_iter=jacobi_iter+1;
         end
         tprev=t;
    end
    contourf(x,y,t,'ShowText','on')
    colormap(jet)
    colorbar
    %set(gca,'YDir','reverse')
    xlabel('X axis')
    ylabel('t')
    text1 = sprintf('number of iteration = %d', jacobi_iter);
    title({'Transient state heat conduction profile by jacobi'; text1})
end



 % guass seidel method
 nt=100;
if iterative_solver==2
    gs_iter = 1; 
     for z=1:nt
         error=100;
        while(error>tol)
            for i =2:(nx-1)
                for j =2:(ny-1)
                     term1=(1+(2*k_1)+(2*k_2))^-1;
                     term2=k_1*term1;
                     term3=k_2*term1;
                     H=(t(i-1,j)+told(i+1,j));
                     V=(told(i,j-1)+t(i,j+1));
                     t(i,j)=((tprev(i,j)*term1)+(H*term2)+(V*term3));
                end
           end  
           error=max(max(abs(told-t)));
           told=t;
           gs_iter=gs_iter+1;
        end
        tprev=t;
     end
     contourf(x,y,t,'ShowText','on')
     colormap(jet)
     colorbar
     xlabel('X axis')
     ylabel('t')
     text1 = sprintf('number of iteration = %d', gs_iter);
     title({'Transient state heat conduction profile by Guess Seidel'; text1})
end

% successive over relaxation iterative method
nt=100;
if iterative_solver==3                                                                                                                                                                      
   sor_iter = 1;
   omega=1.4;
    for z=1:nt
        error=100;
        while(error>tol)
            for i =2:(nx-1)
                for j =2:(ny-1)
                    term1=(1+(2*k_1)+(2*k_2))^-1;
                    term2=k_1*term1;
                    term3=k_2*term1;
                    H=(t(i-1,j)+told(i+1,j));
                    V=(told(i,j-1)+t(i,j+1));
                    t(i,j) = (told(i,j)*(1-omega))+(omega*((tprev(i,j)*term1)+(H*term2)+(V*term3)));             
                end
            end  
            error=max(max(abs(told-t)));
            told=t;
            sor_iter= sor_iter+1;
        end
        tprev=t;
    end
    contourf(x,y,t,'ShowText','on')
    colormap(jet)
    colorbar
    xlabel('X axis')
    ylabel('t') 
    text1 = sprintf('number of iteration = %d', sor_iter);
    title({'Transient state heat conduction profile by SOR'; text1})
end

Output:-

By Jacobi method

By Guass Seidal Method:-

By Successive Over relaxation method:-

Transient state heat conduction by explicit method:-

equation for explicit method

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

clear all
close all
clc

% Initial conditions

nx = 10; % number of nodes 
ny = 10; % number of nodes
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha=0.3; % thermal diffusivity
dt=1e-2;


% BC
t=300*ones(nx,ny);
t(:,1) = 400; % Left Temperature
t(1,:) = 900; % Top Temperature
t(end,:) = 600; % Bottom Temperature
t(:,end) = 800; % Right Temperature

% Editing Corner Values 
t(1,1) = (400+900)/2;
t(end,1) = (400 + 600)/2;
t(1,end) = (900+800)/2;
t(end,end)= (800+600)/2;
told = t;

%k = 2*(dx^2+dy^2)/(dx^2*dy^2);
error = 100;
iterative_solver=1;
tol=1e-4;
k_1 = alpha*((dt)/(dx)^2);
k_2 = alpha*((dt)/(dy)^2);


if iterative_solver==1
   explicit = 1;
 while(error>tol)
    for i =2:(nx-1)
        for j =2:(ny-1)
            term1=(told(i-1,j)-2*told(i,j)+told(i+1,j));
            term2=(told(i,j-1)-2*told(i,j)+told(i,j+1));
            t(i,j)=(t(i,j))+(k_1*term1)+(k_2*term2);
            %t(i,j)=(((told(i-1,j)+told(i+1,j))/(dx^2))+((told(i,j+1)+told(i,j-1))/(dy^2)))/k;      
        end
    end  
    error=max(max(abs(told-t)));
    told=t;
    explicit=explicit+1;    
 end
 contourf(x,y,t,'ShowText','on')
 colormap(jet)
 colorbar
 %set(gca,'YDir','reverse')
 xlabel('X axis')
 ylabel('t')
 text1 = sprintf('Time = %d', explicit);
 title({'Transient state heat conduction profile Explicit '; text1})
end

Output:-

Conclusion:-

For explicit methods produce greater accuracy with less computational period require to run it.

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.