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

MMohd Salim17 Jul 2022
MATLAB

Steady-state 2D heat conduction equation:

`(delta^2T)/(delta^2x)+(delta^2T)/(delta^2y)=0`

This equation can be solved using three iterative methods such as Jacobi, Gauss-Seidel and SOR.

MATLAB PROGRAM

FUNCTIONS FOR JACOBI

function [T jacobi_iter]  = Steady_2D_Heat_Conduction_jacobi(dx,dy,T,error_tol,k,error,nx,ny)
                     
Told = T;

jacobi_iter = 1;
while (error>error_tol)
    % Space Loop
    for i=2:nx-1
         for j=2:ny-1
             term1 = (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
             term2 = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
          T(i,j) = term1+term2;
         end
    end
    
%error = max(abs(T(i,j)-Told(i,j)));
error = max(max(abs(T-Told)));
Told = T;
jacobi_iter = jacobi_iter+1;
end
end

FUNCTIONS FOR GAUSS-SEIDAL

function [T GS_iter] = Steady_2D_Heat_Conduction_GS(dx,dy,T,error_tol,k,error,nx,ny)
                     
Told = T;

GS_iter = 1;
while (error>error_tol)
    % Space Loop
    for i=2:nx-1
         for j=2:ny-1
             term1 = (T(i-1,j)+T(i+1,j))/(k*(dx^2));
             term2 = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
          T(i,j) = term1+term2;
         end
    end
    
%error = max(abs(T(i,j)-Told(i,j)));
error = max(max(abs(T-Told)));
Told = T;
GS_iter = GS_iter+1;
end
end

FUNCTIONS FOR Successive Over Relaxation

function [T SOR_iter] = Steady_2D_Heat_Conduction_SOR(dx,dy,T,error_tol,k,error,nx,ny)
                     
Told = T;
r = 1.5;
SOR_iter = 1;
while (error>error_tol)
    % Space Loop
    for i=2:nx-1
         for j=2:ny-1
             term1 = (T(i-1,j)+T(i+1,j))/(k*(dx^2));
             term2 = (T(i,j-1)+T(i,j+1))/(k*(dy^2));
          T(i,j) = (1-r)*Told(i,j) + r*(term1+term2);
         end
    end
    
%error = max(abs(T(i,j)-Told(i,j)));
error = max(max(abs(T-Told)));
Told = T;
SOR_iter = SOR_iter+1;
end
end

Main Script

clc
clear all
close all

%Inputs
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error_tol =  1e-4;
error = 1e5;

%IC & BC
T = 300*ones(nx,ny);
T(1,:) = 600;
T(nx,:) = 900;
T(:,1) = 400;
T(:,ny) = 800;
T(1,1) = mean([400,600]);
T(1,ny) = mean([600,800]);
T(nx,1) = mean([400,900]);
T(nx,ny) = mean([900,800]);
Told = T;
k=2*((dx^2+dy^2)/(dx^2*dy^2));

%Jacobi Solver

iterative_solver = input('Enter the desired number with respect to solver:');
if iterative_solver==1
     tic;
     [T, jacobi_iter] = Steady_2D_Heat_Conduction_jacobi(dx,dy,T,error_tol,k,error,nx,ny); %Function 
     time_jacobi_in_sec=toc;
     figure(1)
     contourf(x,y,T)
     clabel(contourf(x,y,T))
     colorbar
     colormap(jet)
     set(gca, 'ydir', 'reverse')     
     title(sprintf('No. of Jacobi Iterations = %d', jacobi_iter)) 
     
  % Gauss-Seidal Solver   
elseif iterative_solver==2
     tic;
     [T GS_iter] = Steady_2D_Heat_Conduction_GS(dx,dy,T,error_tol,k,error,nx,ny); %Function call
     time_GS_in_sec = toc;
     figure(1)
     contourf(x,y,T)
     clabel(contourf(x,y,T))
     colorbar
     colormap(jet)
     set(gca, 'ydir', 'reverse')     
     title(sprintf('No. of Gauss-Siedel = %d', GS_iter))  
     
     % SOR Solver
     
 elseif iterative_solver==3
     tic;
     [T, SOR_iter] = Steady_2D_Heat_Conduction_SOR(dx,dy,T,error_tol,k,error,nx,ny); %Function Call
     time_GS_in_sec = toc;
     figure(1)
     contourf(x,y,T)
     clabel(contourf(x,y,T))
     colorbar
     colormap(jet)
     set(gca, 'ydir', 'reverse')     
     title(sprintf('No. of Successive Over-Relaxation Iterations = %d',SOR_iter))
end
xlabel('X-Axis')
ylabel('Y-Axis')

OUTPUT

JACOBI SOLVER OUTPUT

GAUSS-SEIDAL SOLVER OUTPUT

SUCCESSIVE OVER RELAXATION SOLVER OUTPUT

IMPLICIT Scheme Unsteady/Transient-state Condition Program:

The unsteady/transied 2d heat conduction equation is given by

MATLAB CODE

%% 2D HEAT CONDUCTION ON A FLAT PLATE - TRANSIENT STATE
clear all
close all
clc

%INPUTS
nx = 10;
ny = 10;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);

% IC & BC
T = 300*ones(nx,ny);

TL = 400;
TR = 800;
TT = 600;
TB = 900;

T(:,1) = TL;
T(:,end) = TR;
T(1,:) = TT;
T(end,:) = TB;

T(1,1) = 0.5*(TL+TT);
T(1,nx) = 0.5*(TR+TT);
T(ny,1) = 0.5*(TL+TB);
T(nx,ny) = 0.5*(TR+TB);

tol = 1e-4;
Terror = 1;

% CFL 
dt = 1;
alpha = 1e-5;
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
omega=2/(1+sin(pi*dx));
CFL = k1+k2;

% Copy Previous time matrix
T_prev = T;
T_old = T;
% Initial matrix
T_initial = T;

% IMPLICIT SOLVER
Titer = 0;
Siter = 0;
tspan = 0;
Serror = 1;


iterative_solver = input('input the appropriate number to select required the method:');

%% JACOBI METHOD

if iterative_solver==1    
% TIME LOOP
% CONVERGENCE criterion for time
while(Terror>tol) 
    Titer = Titer+1;
    Serror = 1;
% SPACE LOOP - JACOBIAN   
while(Serror>tol) 
    
for j = 2:ny-1
    for i = 2:nx-1
    % Current time level but old iteration neighbours
    TL = T_old(i,j-1); TR = T_old(i,j+1); TT = T_old(i-1,j); TB = T_old(i+1,j); 
    H = TL + TR; V = TT + TB; k1k2 = (1+2*k1+2*k2);
    T(i,j) = (T_prev(i,j) + H*k1 + V*k2)/k1k2;
    end
end
Siter = Siter+1;
Serror(Siter) = max(max(abs(T-T_old)));
T_old = T;
end
tspan = tspan + (dt*Titer);
Terror(Titer) = max(max(abs(T-T_prev)));
T_prev = T;
end

%% GAUSS-SEIDAL METHOD
elseif iterative_solver==2
    
    while(Terror>tol) 
    Titer = Titer+1;
    Serror = 1;
% SPACE LOOP - GAUSS-SEIDAL
while(Serror>tol) 
    
for j = 2:ny-1
    for i = 2:nx-1
    % Current time level but old iteration neighbours
          term1 = (1+(2*k1)+(2*k2))^-1;
          term2= k1*term1;
          term3= k2*term1;
          H=(T(i-1,j)+T_old(i+1,j));
          v=(T(i,j-1)+T_old(i,j+1));
          T(i,j)= T_prev(i,j)*term1+term2*H+term3*v;
    end
end

Siter = Siter+1;
Serror(Siter) = max(max(abs(T-T_old)));
T_old = T;


end

tspan = tspan + (dt*Titer);
Terror(Titer) = max(max(abs(T-T_prev)));
T_prev = T;
    end
 
%% SUCCESSIVE 0VER RELAXATION METHOD

elseif iterative_solver==3
    
    while(Terror>tol) 
    Titer = Titer+1;
    Serror = 1;
% SPACE LOOP - SUCCESSIVE OVER RELAXATION
while(Serror>tol) 
    
for j = 2:ny-1
    for i = 2:nx-1
    % Current time level but old iteration neighbours
          term1 = (1+(2*k1)+(2*k2))^-1;
          term2= k1*term1;
          term3= k2*term1;
          H=(T(i-1,j)+T_old(i+1,j));
          v=(T(i,j-1)+T_old(i,j+1));
          T(i,j)= omega*(T_prev(i,j)*term1+term2*H+term3*v)+(1-omega)*T_prev(i,j);
    end
end

Siter = Siter+1;
Serror(Siter) = max(max(abs(T-T_old)));
T_old = T;


end

tspan = tspan + (dt*Titer);
Terror(Titer) = max(max(abs(T-T_prev)));
T_prev = T;
end
end
%Plotting
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 Unsteady Iterations(EXPLICIT) = %d', Titer));

OUTPUT OF CODE

JACOBI OUTPUT

GAUSS-SEIDAL OUTPUT

SUCCESIVE OVER RELAXATION OUTPUT

EXPLICIT Scheme Unsteady/Transient-state Condition Program:

MATLAB CODE:

%% 2D HEAT CONDUCTION ON A FLAT PLATE - TRANSIENT STATE
clear all
close all
clc

%% 1) EXPLICIT METHOD
nx = 10;
ny = 10;

x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);


T = 300*ones(nx,ny);

TL = 400;
TR = 800;
TT = 600;
TB = 900;

T(:,1) = TL;
T(:,end) = TR;
T(1,:) = TT;
T(end,:) = TB;

T(1,1) = 0.5*(TL+TT);
T(1,nx) = 0.5*(TR+TT);
T(ny,1) = 0.5*(TL+TB);
T(nx,ny) = 0.5*(TR+TB);

tol = 1e-4;
error = 1;

% CFL 

dt = 1e-2;
CFL = 0.1;
alpha = (CFL*dx^2)/dt;

k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;

% Copy 
% Previous time matrix
T_prev = T;
% Current time matrix
T_initial = T;

% EXPLICIT SOLVER
iter = 0;
tspan = 0;
% CONVERGENCE criterion for time
while(error>tol) % FOR LOOP (Repeats until nt)

for j = 2:ny-1
    for i = 2:nx-1
        % Previous time level neighbours
    TL = T_prev(i,j-1); TR = T_prev(i,j+1); TT = T_prev(i-1,j); TB = T_prev(i+1,j);
    H = TL + TR; V = TT + TB; 
    k1k2 = 1-2*k1-2*k2;
    T(i,j) = T_prev(i,j)*k1k2 + H*k1 + V*k2;
    end
end

iter = iter+1;
tspan = tspan + (dt*iter);
error(iter) = max(max(abs(T-T_prev)));
T_prev = T;

end

%Plotting
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('Unsteady 2D heat conduction - Explicit Method, Iterations = %d',iter));

OUTPUT

RESULTS/DISCUSSION:

  • In Steady state equation, SOR methods yields the reults very quickly when compare to other two methods.whereas Jacobi takes more number of iterations and more time to solve the unknows.
    • Jacobi > Gauss-seidel > Successive over-relaxation according to the No. of iterations computed
    • So, the SOR method is most efficient when compared to the other two methods to solve the steady-state equation.
  • When comparing Steady state and Unsteady state, it is proven that Steady state is the most successful one as it takes least number of iterations and least time to compute.
  • In transient state equation, explicit methods yield the results quickly when compared to the implicit method, but is recommended under the situation when CFL is less than 0.25.
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.