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

NNirmal Kumar V23 Feb 2021
MATLAB

AIM:

      To write a code to Solve the 2D heat conduction equation by using the point iterative techniques such as

  • Jacobi
  • Gauss-seidel
  • Successive over-relaxation (SOR)

ASSUMPTIONS:

  • The domain is a unit square.
  • Number of points in x direction is equal to number of points along y direction.
  • Boundary condition for steady and transient case.
    • Top Boundary = 600 K
    • Bottom Boundary = 900 K
    • Left Boundary = 400 K
    • Right Boundary = 800 K
  • Absolute error criteria is 1e-4
  • Initial temperature is taken as 298k.

THEORY:

The general two dimensional heat conduction equation is

             `(∂T)/(∂t)+ alpha *((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))=0`  

              Where, T is temperature,

                          t is time.

                          Alpha is thermal diffusivity.

There are two types of heat conduction equation,

  • Steady state
  • Unsteady or transient

STEADY STATE:

The rate of change of temperature remains 0. Therefore, the time derivative term is 0.                         

Therefore the equation becomes,

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

Discretising the above equation by central differencing scheme,

             

MATLAB CODE:

MAIN CODE:

%solving steady stste 2d heat conduction equation
%d^2t/dx^2 + d^2t/dy^2 = 0

clear all
close all
clc

%inputs
nx = 10;
ny = 10;

%grids in x and y direction
x = linspace(0,1,nx);
y = linspace(0,1,ny);

%grid spacing
dx = x(2)-x(1);
dy = y(2)-y(1);

%2d array creation
[xx,yy] = meshgrid(x,y);

%initial temperature
t = 298*ones(nx,ny);

%applying boundary conditions
%left
t(:,1) = 400;
%right
t(:,end) = 800;
%top
t(1,:) = 600;
%bottom
t(end,:) = 900;

%four corners
t(1,1) = (t(2,1)+t(1,2))/2;
t(end,end) = (t(end,9)+t(9,end))/2;
t(1,end) = (t(1,9)+t(2,end))/2;
t(end,1) = (t(9,1)+t(end,2))/2;

%creating copy of t
told = t;

%error and tolerance
error = 9e9;
tol = 1e-4;

%calculation of k
k = 2*((dx^2+dy^2)/((dx^2)*(dy^2)));

%defining the type of solver
iterative_solver = 3;

%jacobi
if iterative_solver == 1
    t_jacobi = jacobi_steady(x,y,nx,ny,t,told,dx,dy,tol,k);
    
%gauss siedel
elseif iterative_solver == 2
    t_gs = gauss_siedel_steady(x,y,nx,ny,t,told,dx,dy,tol,k);
    
%sor
else iterative_solver == 3
    t_sor = sor_steady(x,y,nx,ny,t,told,dx,dy,tol,k);
    
end

FUNCTION CODE:

 JACOBI:

function t = jacobi_steady(x,y,nx,ny,t,told,dx,dy,tol,k)  
%time loop
jacobi_iter = 1;
tic
%convergence loop
error = 9e9;
    while(error>tol)
        %spatial 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
        error = max(max(abs(told - t)));
        %updating the value of t
        told=t;
        jacobi_iter = jacobi_iter+1;
    end
    t_time = toc;
    %plotting
    contourf(x,y,t)
    [a,b] = contourf(x,y,t);
    colormap(jet)
    clabel(a,b)
    colorbar
    set(gca,'ydir','reverse')
    %title and labelling
    xlabel('x axis')
    ylabel('y axis')
    title_text = sprintf('number of jacobi iterations : %d \nconvergence time:%d ',jacobi_iter,t_time);
    title(title_text)
end

 GAUSS SIEDEL:

function t = gauss_siedel_steady(x,y,nx,ny,t,told,dx,dy,tol,k)
%time loop
gs_iter = 1;
tic
%convergence loop
error = 9e9;
    while(error>tol)
        %spatial loop
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = (t(i-1,j)+told(i+1,j))/(k*(dx^2));
                term2 = (t(i,j-1)+told(i,j+1))/(k*(dy^2));
                t(i,j) = term1+term2;
            end
        end
        %error
        error = max(max(abs(told - t)));
        %updating the value of t
        told=t;
        gs_iter = gs_iter+1;
    end   
    t_time = toc;
    %plotting
    contourf(x,y,t);
    [a,b] = contourf(x,y,t);
    colormap(jet);
    clabel(a,b);
    colorbar
    set(gca,'ydir','reverse')
    %title and labelling
    xlabel('x axis')
    ylabel('y axis')
    title_text = sprintf('number of gauss siedel iterations : %d \nconvergence time :%d',gs_iter,t_time);
    title(title_text)
end

 SUCCESSIVE OVER RELAXATION:

function t = sor_steady(x,y,nx,ny,t,told,dx,dy,tol,k)
omega = 1.2;
%time loop
   sor_iter = 1;
   tic
   error = 9e9;
   %convergence loop
    while(error>tol)
        %spatial loop
        for i = 2:nx-1
            for j = 2:ny-1
                term1 = (t(i-1,j)+told(i+1,j))/(k*(dx^2));
                term2 = (t(i,j-1)+told(i,j+1))/(k*(dy^2));
                t(i,j) = ((1-omega)*told(i,j))+(omega*(term1+term2));
            end
        end
        %error
        error = max(max(abs(told - t)));
        %updating the value of t
        told=t;
        sor_iter = sor_iter+1;
    end
    t_time = toc;
    %plotting
    contourf(x,y,t)
    [a,b] = contourf(x,y,t);
    colormap(jet)
    clabel(a,b)
    colorbar
    set(gca,'ydir','reverse')
    %title and labelling
    xlabel('x axis')
    ylabel('y axis')
    title_text = sprintf('number of sor iterations : %d \nconvergence time :%d',sor_iter,t_time);
    title(title_text)
end

OUTPUT:

JACOBI:

 

GAUSS SIEDEL:

SOR:

 From the three iteration methods, the jacobi method takes more number of iterations and time to make convergences while, successive over relaxation takes less time and less number of iterations to converge and if we increase the omega value in sor method the number of iteration decreases.

TRANSIENT STATE:

Both terms are considered in this condition,

           `(∂T)/(∂t)+ alpha *((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))=0`              

The above transient equation can be solved in two ways by using three iterative methods,

  • Implicit method
  • Explicit method.

EXPLICIT METHOD:

           

for explicit scheme we can use the direct methods so there is no need to use iterative methods.

MATLAB CODE:

MAIN CODE:

%solving transient 2d heat conduction equation
%dt/dt = alpha*(d^2t/dx^2 + d^2t/dy^2 )

clear all
close all
clc

%inputs
nx = 10;
ny = 10;

%grids in x and y direction
x = linspace(0,1,nx);
y = linspace(0,1,ny);

%grid spacing
dx = x(2)-x(1);
dy = y(2)-y(1);

%2d array creation
[xx,yy] = meshgrid(x,y); 

%initial temperature
t = 298*ones(nx,ny);

%applying boundary conditions
%left
t(:,1) = 400;
%right
t(:,end) = 800;
%top
t(1,:) = 600;
%bottom
t(end,:) = 900;

%four corners
t(1,1) = (t(2,1)+t(1,2))/2;
t(end,end) = (t(end,9)+t(9,end))/2;
t(1,end) = (t(1,9)+t(2,end))/2;
t(end,1) = (t(9,1)+t(end,2))/2;

%creating copy of t
told = t;

%thermal diffusivity
alpha = 1.4;
dt = 0.001;
nt = 1400;

%error and tolerance
error = 9e9;
tol = 1e-4;

%calculation of k
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);

%calling the function
t_explicit = transient_explicit(x,y,nt,nx,ny,t,told,k1,k2);

FUNCTION CODE:

function t = transient_explicit(x,y,nt,nx,ny,t,told,k1,k2)
    %explicit calculation
    counter = 1;
    %time loop
    tic
    for kt = 1:nt
            %spatial loop
            for i = 2:nx-1
                for j = 2:ny-1
                    term1 = told(i,j);
                    term2 = k1*(told(i-1,j)-2*told(i,j)+told(i+1,j));
                    term3 = k2*(told(i,j-1)-2*told(i,j)+told(i,j+1));
                    t(i,j) = term1+term2+term3;
                end
            end
            %updating the value of t
            told=t;
            counter = counter+1;
    end
    t_time = toc;
    %plotting
    contourf(x,y,t)
    [a,b]= contourf(x,y,t);
    colormap(jet)
    clabel(a,b)
    colorbar
    set(gca,'ydir','reverse')
    %labelling and title
    xlabel('x axis')
    ylabel('y axis')
    title_text = sprintf('number of iterations explicit method : %d \n convergence time : %d',counter,t_time);
    title(title_text)
end

OUTPUT:

 

The explicit scheme takes 1401 iterations and 0.00607 sec to converge. when we increase or decrease the value of dt the cfl number goes more than 0.5 and the solution becomes unstable.

IMPLICIT METHOD:

 

SOR MODIFIED EQUATION:

 

MATLAB CODE:

MAIN CODE:

%solving transient 2d heat conduction equation
%dt/dt = alpha*(d^2t/dx^2 + d^2t/dy^2 )

clear all
close all
clc

%inputs
nx = 10;
ny = 10;

%grids in x and y direction
x = linspace(0,1,nx);
y = linspace(0,1,ny);

%grid spacing
dx = x(2)-x(1);
dy = y(2)-y(1);

%2d array creation
[xx,yy] = meshgrid(x,y);

%initial temperature
t = 298*ones(nx,ny);

%applying boundary conditions
%left
t(:,1) = 400;
%right
t(:,end) = 800;
%top
t(1,:) = 600;
%bottom
t(end,:) = 900;

%four corners
t(1,1) = (t(2,1)+t(1,2))/2;
t(end,end) = (t(end,9)+t(9,end))/2;
t(1,end) = (t(1,9)+t(2,end))/2;
t(end,1) = (t(9,1)+t(end,2))/2;

%creating copy of t
told = t;
tprev = t;

%thermal diffusivity
alpha = 1.4;
dt = 0.001;
nt = 1400;

%error and tolerance
error = 9e9;
tol = 1e-4;

%calculation of k
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);

%defining the type of solver
iteration_solver = 1;

if iteration_solver == 1
    %jacobi
    t_jacobi = transient_jacobi_implicit(x,y,nt,nx,ny,t,told,tol,k1,k2,tprev);
    
elseif iteration_solver == 2
    %gauss siedel
    t_gs= transient_gs_implicit(x,y,nt,nx,ny,t,told,tol,k1,k2,tprev);

else iteration_solver == 3
    %sor
    t_sor = transient_sor_implicit(x,y,nt,nx,ny,t,told,tol,k1,k2,tprev);
end

FUCTION CODE:

 JACOBI:

function t = transient_jacobi_implicit(x,y,nt,nx,ny,t,told,tol,k1,k2,tprev)
    jacobi_iter = 1;
    %time loop
    tic
    for kt = 1:nt
        error = 9e9;
        %convergence loop
        while(error>tol)
            %spatial loop
            for i = 2:nx-1
                for j = 2:ny-1
                    term1 = (1+(2*k1)+(2*k2))^(-1);
                    term2 = k1*term1;
                    term3 = k2*term1;
                    h = (told(i-1,j)+told(i+1,j));
                    v = (told(i,j-1)+told(i,j+1));
                    t(i,j) = (term1*tprev(i,j))+(term2*h)+(term3*v);
                end
            end
            %error
            error = max(max(abs(told - t)));
            %updating the value of t
            told=t;
            jacobi_iter = jacobi_iter+1; 
        end
        %updating the previous value of t
        tprev = t;
    end
    t_time = toc;
    %plotting
    contourf(x,y,t)
    [a,b] = contourf(x,y,t);
    colormap(jet)
    clabel(a,b)
    colorbar
    set(gca,'ydir','reverse')
    %title and labelling
    xlabel('x axis')
    ylabel('y axis')
    title_text = sprintf('number of jacobi iterations implicit method : %d \nconvergence time : %d',jacobi_iter,t_time);
    title(title_text)
end

 GAUSS SIEDEL:

function t = transient_gs_implicit(x,y,nt,nx,ny,t,told,tol,k1,k2,tprev)
    gs_iter = 1;
    tic
    %time loop
    for kt = 1:nt
        error = 9e9;
        %convergence loop
        while(error>tol)
            %spatial loop
            for i = 2:nx-1
                for j = 2:ny-1
                    term1 = (1+(k1*2)+(2*k2))^(-1);
                    term2 = k1*term1;
                    term3 = k2*term1;
                    h = (t(i-1,j)+told(i+1,j));
                    v = (t(i,j-1)+told(i,j+1));
                    t(i,j) = (term1*tprev(i,j))+(term2*h)+(term3*v);
                end
            end
            %error
            error = max(max(abs(told - t)));
            %updating the old value of t
            told=t;
            gs_iter = gs_iter+1;
        end
        %updating the previous value of t
        tprev = t;
    end
    t_time = toc;
    %plotting
    contourf(x,y,t)
    [a,b] = contourf(x,y,t);
    colormap(jet)
    clabel(a,b)
    colorbar
    set(gca,'ydir','reverse')
    %title and labelling
    xlabel('x axis')
    ylabel('y axis')
    title_text = sprintf('number of gs iterations implicit method : %d \nconvergence time : %d',gs_iter,t_time);
    title(title_text)
end

 SOR:

function t = transient_sor_implicit(x,y,nt,nx,ny,t,told,tol,k1,k2,tprev)
omega = 1.2;
   sor_iter = 1;
   tic
   %time loop
       for kt = 1:nt
           error = 9e9;
           %convergence loop
            while(error>tol)
                %spatial loop
                for i = 2:nx-1
                    for j = 2:ny-1
                        term1 = (1+(k1*2)+(2*k2))^(-1);
                        term2 = k1*term1;
                        term3 = k2*term1;
                        h = (t(i-1,j)+told(i+1,j));
                        v = (t(i,j-1)+told(i,j+1));
                        t_gs = ((tprev(i,j)*term1)+(h*term2)+(v*term3));
                        t(i,j) = ((1-omega)*told(i,j))+(omega*t_gs);
                    end
                end
                %error
                error = max(max(abs(told - t)));
                %updating the value of t
                told=t;
                sor_iter = sor_iter+1;      
            end
            %updating the previous value of t
            tprev = t;
       end   
       t_time = toc;
       %plotting
       contourf(x,y,t)
       [a,b] = contourf(x,y,t);
       colormap(jet)
       clabel(a,b)
       colorbar
       set(gca,'ydir','reverse')
       %title and labelling
       xlabel('x axis')
       ylabel('y axis')
       title_text = sprintf('number of sor iterations implicit method : %d \nconvergence time : %d',sor_iter,t_time);
       title(title_text)
end

OUTPUT:

 JACOBI:

GAUSS SIEDEL:

SOR:

 

The jacobi method takes more time and more number of iterations to converge. while the sor method takes less number of iterations and less time to converge and as we increase the relaxation factor it tooks further less time to converge.

CONCLUSION:

 FOR STEADY STATE:

              Among the three iterative methods ,the successive over relaxation iterative method takes the least no of iterations and converging time. sor takes total of 76 iterations and 0.01084 seconds whereas the jacobi iteration takes the most with 208 iterations and 0.00308 seconds to converge and gauss siedel lies in between them by taking 112 iterations and 0.00270 seconds to converge.

FOR UNSTEADY STATE:

EXPLICIT METHOD:

           In this method if the CFL number exceeds 0.5 the solution becomes unstable.

           This method takes 1401 iterations and 0.00607 seconds to converge which is the most among the transient methods.

 IMPLICIT METHOD:

           Among the three iterative methods ,the successive over relaxation iterative method takes the least no of iterations and converging time. sor takes total of 2993 iterations and 0.01554 seconds whereas the jacobi iteration takes the most with 3669 iterations and 0.01785 seconds to converge and gauss siedel lies in between them by taking 3120 iterations and 0.01605 seconds to converge.

           Implicit linear convection is unconditionally stable, and also lacks a CFL condition. This happens because implicit schemes use the entire domain to calculate each timestep. 

overall, the steady state takes less iterations to converge than transient state. so, the sor method is the best method which takes less time to converge in both the above scenarios.

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.