Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

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 assume 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 & Successive over-relaxation
      
   3. Various Methods to be used in Unsteady-state: IMPLICIT & EXPLICIT

By using Point iterative method, solve the 2D Heat Conduction Equation and plot the Temperature Contour for all the three iterative methods namely Jacobi, Gauss-Siedel & SOR for Steady and Unsteady conditions under IMPLICIT scheme also additionally solving the EXPLICIT scheme for Unsteady Condition.

THEORY:

The heat equation is a partial differential equation that describes how the distribution of some quantity (such as heat) 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 of three spatial variables (see Cartesian coordinate system) and the time variable , the heat equation is

dou T / dou t = alpha*(dou T^2 / dou x ^2 + dou T ^2/dou y ^2 + dou T^2/dou z^2)

where,

alpha = real coefficient called the diffusivity of the medium,

Types of Heat Conduction schemes

This transfer of Heat energy may occur under two conditions

• Steady-state condition
• Unsteady/Transient state condition

Unsteady state conditions are a precursor to steady state conditions. No system exists initially under steady state conditions. Some time must pass, after heat transfer is initiated, before the system reaches steady state. During that period of transition the system is under unsteady state conditions.

MAIN PROGRAM:

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

MPLICIT 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 (if this were not so, the temperature would be rising or falling, as thermal energy was tapped or trapped in a region).

Jacobi Iterations

clear all
close all
clc

%solving the steady state 2d heat conduction (Jacobi method)

%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 > tolerance

error = 9e9;
tolerance = 1e-4;

%defining 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 corner
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 equation by jacobi method
iterative_solver = 1;

if iterative_solver == 1
    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_old(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

    %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_text = (sprintf('Number of Jacobi Iterations (implicit) = %d',jacobi_iteration));
    title(title_text)
pause(0.03)
end

OUTPUT

Gauss Siedel Iterations

clear all
close all
clc

%solving the steady state 2d heat conduction (gauss seidel method)

%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 > tolerance

error = 9e9;
tolerance = 1e-4;

%defining 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 corner
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 equation by jacobi method
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

    %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_text = (sprintf('Number of gauss seidel Iterations (implicit) = 
    %d',gauss_seidel_iteration));
    title(title_text)
    pause(0.03)
end

OUTPUT

Succesive over relaxation iterations

clear all
close all
clc

%solving the steady state 2d heat conduction (succesive over relaxation method)

%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 > tolerance

error = 9e9;
tolerance = 1e-4;
omega = 1.2;

%defining 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 corner
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 equation by succesive over relaxation method
iterative_solver = 3;
if iterative_solver == 3
    succesive_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;
        succesive_over_relaxation_iteration = succesive_over_relaxation_iteration + 1;
    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_text = (sprintf('Number of SOR Iterations (implicit) = %d',succesive_over_relaxation_iteration));
    title(title_text)
    pause(0.03)
end
               

IMPLICIT Scheme Unsteady/Transient-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.

Jacobi Iterations :-

clear all
close all
clc

%solving the unsteady state 2d heat conduction by jacobi method(exp scheme)

%Input parameters
%No 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;

%defining 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_old = T;
T_initial = T;

%calculation of 2d steady heat conduction by jacobi method
k1 = 1.1*(dt/dx^2);
k2 = 1.1*(dt/dy^2);

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-axis')
    title(sprintf('No of transient state Jacobi Iterations (implicit) = %d',jacobi_iteration))
  
    pause(0.03)
end     

Gauss Siedel Iterations :-

clear all
close all
clc

%solving the unsteady state 2d heat conduction by gauss seidel method(exp scheme)

%Input parameters
%No 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;

%defining 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_old = T;
T_initial = T;

%calculation of 2d steady heat conduction by jacobi method
k1 = 1.1*(dt/dx^2);
k2 = 1.1*(dt/dy^2);

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-axis')
    title(sprintf('No of transient state gauss seidel Iterations (implicit) = %d',gauss_seidel_iteration))
  
    pause(0.03)
end     

succesive over relaxation iterations :-

clear all
close all
clc

%solving the unsteady state 2d heat conduction by sor method(implicit scheme)

%Input parameters
%No 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;

%defining 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_old = T;
T_initial = T;
T_end = T;

%calculation of 2d steady heat conduction by sor method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

sor_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_end;
        sor_iteration = sor_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-axis')
    title(sprintf('No of transient state SOR Iterations (implicit) = %d',sor_iteration))
  
    pause(0.03)
end     

EXPLICIT Scheme - Unsteady/Transient state Condition Program

Unsteady-state Equation discretization(EXPLICIT):

Explicit schemes :-

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction (EXPLICIT SCHEME)

%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;

%Defining 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 orginal values to T
T_intial = T;
T_old = T;


%Calculation of 2D steady heat conduction 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

%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 steady Iterations(EXPLICIT) = %d', Steady_Explicit));
pause(0.03);

CONCLUSIONS:

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

• Jacobi > Gauss-seidel > Successive over-relaxation according to the No. of iterations computed for all methods & Schemes.
  
  
• 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.
  
  
• SOR iteration method is the most effective iteration method as it took least No. of iterations for coverging solutions in both the schemes (EXPLICIT & IMPLICIT) compared to Jacobi & Gauss-siedel.