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

AIM: The main goal of this challenge is to solve Steady and Unsteady state 2D Heat Conduction Equation by using the iterative techniques such as Jacobi method, Gauss-Seidel method & Successive over-relaxation method for both implicit and explicit schemes using MATLAB.

The assumption for solving the problem:

1. Assume that the domain is to be a square.
2. Assume nx = ny [Number of points along the x-direction is equal to number of points along y-direction]
3. We are assuming nx = ny = 10
4. Boundary Conditions for steady and unsteady/transient case,
   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. We will implementing the fillowing methods,

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

7. Absolute error criteria [Tolerance] is 1e-4.

THEORY:

Considering General Three dimensional Heat Conduction Equation which has both time derivative and the spatial derivative terms.

The above equation has two terms where, term1 represents time derivative and term2 is spatial derivative, where alpha is thermal diffusivity.

Now by considering the assumption that the above equation is solved for a thin plate, dimensions in the Z-direction are neglected, thereby the above equation is converted into the following form which is as follows:

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

The above equation can be calculated by considering two approaches:

1.Steady State equation:

Under these conditions, the rate of change of temperature remains 0. Therefore in this condition time, the derivative term is equated to zero.

Steady-state 2D heat conduction equation:

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

2.Transient/Unsteady-State equation:

Under these conditions, both terms are considered therefore the equation is as follows:

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

Above two equations can be solved in 2 ways by using 3 iterative methods such as Jacobi, Gauss-Seidel and SOR.

1.Explicit Method

2.Implicit Method

Here in this challenge Transient state equation is solved by using both ways.

CODE STRUCTURE:

BOUNDRY CONDITION MAP:

Here,

`T_L` = Left boundary

`T_T` = Top boundary

`T_P` = Middle/Point boundary

`T_B` = Bottom boundary

`T_R` = Right boundary

From the above Boundary conditions, we can calculate the central temperature of the domain `(T_P)=1/K*((T_L+T_R/dx^2)+(T_T+T_B/dy^2))`

For 2D heat conduction equation we are using the following method:

• JACOBI – IN this method, the 2D heat conduction equation solves using calculating the value of Temperature at inner nodes using several iterations, in which each value of temperature will be previously assigned temperature values at nodes.

• GAUSS-SEIDEL – In this method, the 2D heat conduction equation solves, we are calculating the value of temperatures at inner nodes using several iterations, in which every value of temperature will be updated successively with respect to the previous iteration.

• SUCCESSIVE OVER-RELAXATION- In this method, the 2D heat conduction equation solves. We are calculating the value of temperatures at inner nodes using several iterations, in which every value of temperature will be updated successively with respect to the previous iteration. And also “Relaxation factor” multiplies to the algebraic equation to alter the path of the iteration.

  When RF > 1 over relaxation

  When RF < 1 under relaxation

  The value of Relaxation factor can be taken as 1.1 or 1.2

  The Relaxation factor plays a vital role in determining the rate of convergence for CFD simulation.

Solving the Steady-state equation:

Considering Two dimensional Steady State Equation

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

after discretising above equation,

JACOBI ITERATION [STEADY-STATE]:

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Jacobi Method)

%Input conditions
%Number of nodes
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;

%Input 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 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
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 Jacobi Iterations(IMPLICIT) = %d', jacobi_iteration));
pause(0.03);

OUTPUT:

GAUSS SEIDAL ITERATION [STEADY-STATE]:

close all
clear 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 criteria > 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 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 Gauss-Seidel 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
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 Gauss-Seidel Iterations(IMPLICIT) = %d', Gauss_seidel_iteration));
pause(0.03);

OUTPUT:

SUCCESSIVE OVER RELAXATION ITERATION [STEADY-STATE]:

close all
clear all
clc

%Solving the Steady state 2D heat conduction (Successive 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 criteria > 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 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 = 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-Axis')
title(sprintf('No. of Successive Over-Relaxation Iterations(IMPLICIT) = %d', Successive_Over_Relaxation_iteration));
pause(0.03);

OUTPUT:

Solving Transient state Equation:

Transient state 2D Steady state heat conduction equation is solved in both ways i.e explicitly and implicitly. Following is the 2D heat conduction equation

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

Solving the above equation Explicitly:

Explicit Approach is used in numerical analysis for obtaining numerical approximations to the solution of time-dependent ordinary or partial differential equations as in required in computer simulations of a physical process. Explicit methods calculate the state of a system at a later time from the state of the system at the current time.

In this method, each differential equation approach contains only one unknown and can be solved explicitly for this one unknown using this method.

Here in this equation time derivative terms are calculated using the first order differential equations and whereas the term 2 has second-order derivative terms and elliptical in nature so information needs to be taken from all the points equally, So, therefore, a central differencing scheme is applied for the discretising the term2.

the discretised equation for 2D heat conduction equation for the explicit form is as follows:

JACOBI ITERATION [UNSTEADY-STATE]:

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Jacobi Method(IMPLICIT 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;

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

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation 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_intial(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_intial = 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 Unsteady Jacobi Iterations(EXPLICIT) = %d', Jacobi_iteration));
    pause(0.03)
end

Output

GAUSS SEIDAL ITERATION [UNSTEADY-STATE]:

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Gauss Seidel Method(IMPLICIT 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;

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

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation 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_intial(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_intial = 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 Unsteady Gauss-Seidel Iterations(IMPLICIT) = %d', Gauss_Seidel_iteration));
    pause(0.03)
end

OUTPUT:

SUCCESSIVE OVER RELAXATION ITERATION [UNSTEADY-STATE]:

close all
clear all
clc

%Solving the Unsteady state 2D heat conduction by Successive over-relaxation Method(IMPLICIT 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_old = T;
T_intial = T;
T_end = T;

%Calculation of 2D steady heat conduction EQUATION by Successive over-relaxation method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));

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_intial(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;
        Successive_over_relaxation_iteration = Successive_over_relaxation_iteration + 1;  
        
    end
    
    T_intial = 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('No. of Unsteady SOR Iterations(IMPLICIT) = %d', Successive_over_relaxation_iteration));
    pause(0.0000003)
end

OUTPUT:

EXPLICIT SCHEME - UNSTEADY/TRANSEINT STATE CONDITION

EXPLICIT ITERATION:

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

OUTPUT:

CONCLUSION:

Conclusion:

1.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.

So the SOR method is most efficient when compared to the other two methods to solve the steady-state equation.

2. For transient equations, the time required to solve is more when compared to steady-state equations.So, therefore it is recommended to solve the 2D heat conduction equation with steady state conditions.

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.

When transient state equation is solved implicitly Jacobi takes less time when compare to the other two methods whereas SOR is more stable compared to gauss seidal.