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

KKrishna M Gavali5 Feb 2022
MATLAB

OBJECTIVE:

To solve Steady and Unsteady state 2D Heat Conduction Equation by using the point iterative techniques such as Jacobi, Gauss Seidel & Successive over - relaxation for both implicit and explicit schemes using MATLAB.

Problem Statement:

1.Computation area is to be a square.

2.Assume nx = ny[No. of points on x-direction is equal to no. of points in y-direction.

3.We assune nx = ny = 10

4. Boundary Conditions are

Left = 400k

Right = 600k

Top = 800k

Bottom = 900k

5. Apart from the boundary conditions, all other nodes are given an initial temperature of 300k

6.Various methods to be employed in solving the equation are given below:

1.Various States : Steady-State & Unsteady- state Condition

2. Various Schemes Jacobi, Gauss-Seidel & Succesive over - relaxation

3. Various Methods to be used in Unsteady-state : IMPLICIT & EXPLICIT

By using Point iterative method, solve the 2D Heat Conduction and plot the Temperaature Contour for all the three iterative methods namely Jacobi, Gauss-Siedel &SOR for Steady and Unsteady conditions under IMPLICIT scheme also additionaly solvng the EXPLICIT scheme for Unsteady COndition.

THEORY:

The heat equation is a partial differential equation that descrives how the distribution of some quanttity evolves over time in a solid medium, as it spontaneously flows from places where it is higher towards places where it is lower.

For a function u(x,y,z,t) of three spatial variables (x,y,z) and the time variable t, the heat equation is

`(partial(u))/(partial(t)) = alpha((partial^2(u))/(partial(x^2)) + (partial^2(u))/(partial(y^2))+ (partial^2(u))/(partial(z^2)))`

wheer

alpha is a real coefficient called the diffusivity of the medium.

Types of Heat Conduction schemes

This transfer of Heat energy may occur under two condtions

* Steady - State condition

* Unsteady/Transient state condition

Unsteaady state condtions are a precursor to steady state conditions. No system exists initially under steady state conditions. Some time must pass, after heat transfer is inititated, bedore the system reaches steaady state. During that period of transition the system is under unsteady state conditions.

Governing Equations -

In order to solve the 2D Heat conduction problem, the governing equations for the respective conditions with respect to cartesian coordinates should be dsceretized to enter that in a code which is given under MAIN PROGRAM below.

MAIN PROGRAM :

Below there is a sequence of seven indidvidual programs, starting with IMPLICIT Scheme of Steady & Unsteady - state conditions for three iterative methods Jacobi, Gauss - Siedel & SOR followed by EXXPLICIT Scheme for Unsteaday/Transient Conditions.

IMPLICIT Scheme - Steady - State Condition Program:

Steady-state Equation discretization:

In steady - state conduction, the amount of heat entering any region of an object is equal to the amount of heat coming out

By diiscretizing the equation above it becomes,

Here the value of k is,

If we further simplify this equation for code becomes

0.25(TL + TR + TT + TB)

which can also be written as,

Jacobi Iteration - 

clear all
close all
clc

%soving the steady state 2D conduction 
%Input parameters
%number of grid points
nx = 10;
ny = nx ;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

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

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

%boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx , ny);

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;

%calculating Average temperature at corners
T(1,1) = (T_T + T_L)/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;

%Assigning original values to T
T_initial = T;
T_old = T;

%calculation of 2D steady heat conduction 
iterative_solver = 2;
if iterative_solver == 2;
    jacobi_iteration = 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(i,j-1) + T_old(i,j+1));
            end
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        jacobi_iteration = jacobi_iteration+1;
    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-Label')
title(sprintf('No. of jacobi_iterations(IMPLICIT) = %d', jacobi_iteration));
pause(0.03);

OUTPUT :

Gauss - Siedel Iteration - 

clear all
close all
clc

%soving the steady state 2D conduction 
%Input parameters
%number of grid points
nx = 10;
ny = nx ;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

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

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

%boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx , ny);

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;

%calculating Average temperature at corners
T(1,1) = (T_T + T_L)/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;

%Assigning original values to T
T_initial = T;
T_old = T;

%calculation of 2D steady heat conduction 
iterative_solver = 2;
if iterative_solver == 2;
    Gauss_seidel_iteration = 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(i,j-1) + T_old(i,j+1));
            end
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Gauss_seidel_iteration = Gauss_seidel_iteration+1;
    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-Label')
title(sprintf('No. of Gauss_seidel Iterations(IMPLICIT) = %d', Gauss_seidel_iteration));
pause(0.03);

OUTPUT -

Successive over relaxation 

clear all
close all
clc

%soving the steady state 2D conduction 
%Input parameters
%number of grid points
nx = 10;
ny = nx ;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

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

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

%boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx , ny);

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;

%calculating Average temperature at corners
T(1,1) = (T_T + T_L)/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;

%Assigning original values to T
T_initial = T;
T_old = T;

%calculation of 2D steady heat conduction 
iterative_solver = 3;
if iterative_solver == 3;
    Successive_Over_Relaxation_iteration = 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(i,j-1) + T_old(i,j+1)))+(1-omega)*T_old(i,j);
            end
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Successive_Over_Relaxation_iteration = Successive_Over_Relaxation_iteration+1;
    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-Label')
title(sprintf(' No of Successive Over-Relaxation Iterations(IMPLICIT) = %d', Successive_Over_Relaxation_iteration));
pause(0.03);

OUTPUT -

IMPLICIT Scheme Unsteady/Trabsient - state Condition Program -

Unsteady - state Equation discretization:

Non - steady - state situations appear after an imposed change in temperature at a boundary as well as inside of an object, as a result of a new source or sink of heat suddenly introduced within an object, causing temperatures near the source or sink to change in time.

If we discretize the equation it becomes

If we simplify it further

Alpha is the diffusivity term.

For the better understannding for code.

Jacobi Iteration

clear all
close all
clc

%soving the steady state 2D conduction 
%Input parameters
%number of grid points
nx = 10;
ny = nx ;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

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

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

%boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx , ny);

T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:nx-1) = T_B;

%calculating Average temperature at corners
T(1,1) = (T_T + T_L)/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;

%Assigning original values to T
T_initial = T;
T_old = T;
%calculation of 2D steady heat coduction equation by successive
%over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

%calculation of 2D steady heat conduction 
Jacobi_iteration = 1;
for k = 1:nt
    error = 9e9;
    
    while(error>tolerance)
        for i =  2:nx-1
            for j = 2:ny -1
                Term_1 = 1/(1+(2*k1) + (2*k2));
                Term_2 = k1*Term_1;
                Term_3 = k2*Term_1;
                H = (T_old(i-1,j)) + (T_old(i+1, j));
                V = (T_old(i,j+1)) + (T_old(i,j-1));
                
                T(i,j) = (T_initial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);
            end
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Jacobi_iteration = Jacobi_iteration+1;
    end
    T_initial = T;

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

    xlabel('X - Axis')
    ylabel('Y-Label')
    title(sprintf('No. of Unsteady jacobi_iterations(IMPLICIT) = %d', Jacobi_iteration));
    pause(0.03);
end

OUTPUT :

Gauss - Siedel Iteration

clear all
close all
clc

%soving the steady state 2D conduction 
%Input parameters
%number of grid points
nx = 10;
ny = nx ;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

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

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

%boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx , ny);

T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:nx-1) = T_B;

%calculating Average temperature at corners
T(1,1) = (T_T + T_L)/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;

%Assigning original values to T
T_initial = T;
T_old = T;
%calculation of 2D steady heat coduction equation by successive
%over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

%calculation of 2D steady heat conduction 
Gauss_Seidel_iteration = 1;
for k = 1:nt
    error = 9e9;
    
    while(error>tolerance)
        for i =  2:nx-1
            for j = 2:ny -1
                Term_1 = 1/(1+(2*k1) + (2*k2));
                Term_2 = k1*Term_1;
                Term_3 = k2*Term_1;
                H = (T(i-1,j)) + (T_old(i+1, j));
                V = (T(i,j+1)) + (T_old(i,j-1));
                
                T(i,j) = (T_initial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);
            end
        end
        
        error = max(max(abs(T_old - T)));
        T_old = T;
        Gauss_Seidel_iteration = Gauss_Seidel_iteration+1;
    end
    T_initial = T;

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

    xlabel('X - Axis')
    ylabel('Y-Label')
    title(sprintf('No. of Unsteady Gauss-Seidel Iterations(IMPLICIT) = %d', Gauss_Seidel_iteration));
    pause(0.03);
end

OUTPUT -

Succesive over - relaxation -

clear all
close all
clc

%soving the steady state 2D conduction 
%Input parameters
%number of grid points
nx = 10;
ny = nx ;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

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

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

%boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx , ny);

T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:nx-1) = T_B;

%calculating Average temperature at corners
T(1,1) = (T_T + T_L)/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;

%Assigning original values to T
T_initial = T;
T_old = T;
T_end = T;
%calculation of 2D steady heat coduction equation by successive
%over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

%calculation of 2D steady heat conduction 
Successive_over_relaxation_iteration = 1;
for k = 1:nt
    error = 9e9;
    
    while(error>tolerance)
        for i =  2:nx-1
            for j = 2:ny -1
                Term_1 = 1/(1+(2*k1) + (2*k2));
                Term_2 = k1*Term_1;
                Term_3 = k2*Term_1;
                H = (T(i-1,j)) + (T_old(i+1, j));
                V = (T(i,j+1)) + (T_old(i,j-1));
                
                T(i,j) = (T_initial(i,j)*Term_1) + (H*Term_2) + (V*Term_3);
                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;
        Successive_over_relaxation_iteration = Successive_over_relaxation_iteration+1;
    end
    T_initial = T;

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

    xlabel('X - Axis')
    ylabel('Y-Label')
    title(sprintf('No. of Unsteady SOR iterations(IMPLICIT) = %d', Successive_over_relaxation_iteration));
    pause(0.03)
end

OUTPUT -

EXPLICIT Scheme - Unsteady/Transient state Condition Program

Unsteady - state Equation discretization (EXPLICIT):

If we discretize the equation it becomes

If we simplify it further

For the better understanding for code

Explicit Scheme - 

clear all
close all
clc

%soving the steady state 2D conduction 
%Input parameters
%number of grid points
nx = 10;
ny = nx ;
nt = 1400;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

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

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

%boundary conditions
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;

T = 300*ones(nx , ny);

T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:nx-1) = T_B;

%calculating Average temperature at corners
T(1,1) = (T_T + T_L)/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;

%Assigning original values to T
T_initial = T;
T_old = T;
%calculation of 2D steady heat coduction equation by successive
%over-relaxation method
iterative_solver = 1;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
if iterative_solver == 1
    Steady_Explicit = 1;
    for k = 1:nt
        error = 1;
        while(error>tolerance)
            for i =  2:nx-1
                for j = 2:ny -1
                    Term_1 = T_old(i,j);
                    Term_2 = k1*(T_old(i+1,j) - 2*T_old(i,j) + (T_old(i-1,j)));
                    Term_3 = k2*(T_old(i,j+1) - 2*T_old(i,j) + T_old(i,j-1));
                    T(i,j) = Term_1 + Term_2 + Term_3;

                end
            end

            error = max(max(abs(T_old - T)));
            T_old = T;
            Steady_Explicit = Steady_Explicit+1;
        end
    end
end

%Plottting
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' , Steady_Explicit));
pause(0.03);

OUTPUT -

CONCLUSIONS -

From the above output of various methodds & comparing their No. of iterations, we can conclude that

* Jacobi>Gauss-seidel>Successive over-relaxation according to the No. of iteratoins computed for all methods & Schemes.

* When comparing Steady state and Unsteady state, it is proven that Steady state is the most succesful one as it takes least number of iterations and least time to compute .

*SOR iteration method is the most effective iteration method as it took least No. of iterations for converging solutions in both the schemes compared tto Jacobi & Gauss - siedel.

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.