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

Aim:

    Solving 2D heat conduction equations using Implicit & explicit methods.

Objective:

    For steady state condition, only implicit method needs to implement while for transient condition, both methods need to be implemented to solve the 2D heat conduction equation. For implicit method, point iterative methods are requried to use. The point iterative methods are jacobi, gauss seidel & successive over relaxation with either of the first two methods. 

    The temperature distribution plots for each program needs to be presented with the report with detailed explanation.  

Design Inputs:

    Boundary conditions:

       Top Boundary = 600 K

       Bottom Boundary = 900 K

       Left Boundary = 400 K

       Right Boundary = 800 K

      - absolute error criteria is 1e-4

      - Domain length = 1 meter

      - Assumption: Domain is in square size. 

      `

  - 2D heat equation for steady state condition

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

  - 2D heat conduction equation for transient condition

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

1. Steady state analysis & Transient State Analysis

    The heat conduction eqautions has been solved using explicit & implicit method. In steady state condition, the equations has been solved using only implicit method while in transient condition, equations will be solved using both methods.

    In explicit method, the current node values has been calculated based on the values calculated in previous time step. The Courant Numbers or the CFL value is very important in explicit method to maintain the stability while implementing it. The values should be below 0.5 to maintain the stability of the process. The CFL values has been determined which will be presented in the respective section. 

    The implicit method try to seek the value of the next node and the values of each variable on that node has not been determined while implementing implicit method on the current node. To determine the values on the next node in mesh, iterative method has been used. Three methods has been implemented to compare their peroformance in terms of numner of iterations & time taken by each method to evaluate their effectiveness. 

Following three methods, with their short introduction about working, has been used:

a) Jacobi Iterative method

    This method is used to solve linear equation by point-wise iteration method. Jacobi method uses the values determined in the previus time step to solve the equations in current  iteration.

b) Gauss Seidel (GS) Iterative method

    Gauss seidel method is also similar to jacobi iterative method in terms of performing iteration by point-wise method. But in contrast to jacobi, this method uses the latest updated values to solve the linear equation.

c) Successive Over Relaxation (SOR) iterative method

    SOR is not an independent method like previous two, but this method has been variant of GS method to improve the performance of the iterative method. This method can be implemented to solve the linear equation with either of the two previously defined method wi th addition to relaxation factor which helps in improving the performance of the method.

1.1 Steady State Condition

    From the given boundary conditions, temperature will be calculated at the centre of the domain using following equation:

    `T_p = 1/k (((T_L + T_R) / (del^2x)) + ((T_T + T_B)/ (del ^2 y)))`

    The code will be executed using implicit method & the supportive function for each iterative method will be attached after the main code.

Main Code:

% Solving 2D heat conduction equation
% Implicit Scheme
% Steady state condition

clear all 
close all 
clc

% inputs
% Grid points
nx = 10;
ny = nx;
l = 1; % length of square domain

% Creating mesh
x = linspace(0,l,nx);
y = linspace(0,l,ny);
dx = x(2) - x(1); % step size in x
dy = dx; % step size in y
[xx,yy] = meshgrid(x,y); % mesh grid

% Initial value of the error
error = 9e9; 
tol = 1e-4; % tolerance value

% Calculating constant "k"
k = 2*((1/dx^2)+ (1/dy^2));

% Boundary Conditions
T = ones(nx,ny);
% inserting temperature values at different side of the domain
T(:,1) = 400; % Left boundary 
T(1,:) = 600; % Top boundary
T(:,end) = 800; % Right boundary
T(end,:) = 900; % Bottom boundary

% Calculating 4 corner nodes temperature
T(1,1) = (600 + 400)/2; 
T(1,end) = (800 + 600)/2;
T(end,1) = (400 + 900)/2;
T(end,end) = (800 + 900)/2;

% Told = Temprature during previous iteration
Told = T; 

% Defining loop to execute solvers 
for iter_solver = 1:3
    if iter_solver == 1 % for Jacobi Iterative method

        % Calling jacobi iterative method function
        jacobi_sol = jacobi_iter(error, tol,k, nx, ny,T,Told, dx, dy, xx,yy); 

    elseif iter_solver == 2 % for Gauss Seidel Iterative method
        
        % Calling Gauss Seidel iterative method function
        gs_sol = gs_iter(error, tol,k, nx, ny,T,Told, dx, dy, xx,yy ); 
        
    elseif iter_solver == 3 % for SOR Iterative method

        % Calling Successive Over Relaxation iterative method function
        sor_sol = sor_iter(error, tol,k, nx, ny,T,Told, dx, dy, xx,yy ); % executing SOR Iterative method
    end
end


% Jacobi Iterative Method Function:

function [out] = jacobi_iter(error, tol,k, nx, ny,T,Told, dx, dy, xx,yy)

% Jacobi iterative method function for the 2D heat conduction equation
% solving code for steady state condition

tic; % start timer

jacobi_iter = 1;
while (error > tol) % program will stop when this criteria will meet
    for i = 2: nx-1 % 2:nx-1 because first & last node has been computed before
        for j = 2:ny-1
            
            % dividing eqaution in two terms
            t1 = ((Told(i-1,j) + Told(i+1,j))/(dx^2))/k;
            t2 = ((Told(i,j-1) + Told(i,j+1))/(dy^2))/k;
            T(i,j) = t1 + t2;
        end
    end
    error = max(max(abs(Told - T))); 
    Told = T;
    
    % given incrementation to continue iterative process
    jacobi_iter = jacobi_iter + 1;
    % stop timer & saving time into variable
    sol_time_jacobi = toc; 
    
    % plotting temperature distribution in the domain
    figure(1)
    % plotting contour graph
    [a,b] = contourf(xx,yy,T); 
    % labelling each curve in contour graph from the data points of a & b
    clabel(a,b)
    % indicating chart color pattern
    colormap(autumn)
    % turining on color bar
    colorbar
    
    % setting y-direction axis in reverse order
    set(gca,'ydir','reverse')
    % iserting title to mention number of iteration & time taken by solver
    title(sprintf('Jacobi Iteration = %d\n Iteration Time = %f sec', jacobi_iter,sol_time_jacobi ))
    % X & Y axis labels
    xlabel('X Axis')
    ylabel('Y Axis')
    
end

out = sol_time_jacobi;

end

% Guass-Seidel Iterative Method Function:

function [out] = gs_iter(error, tol,k, nx, ny,T,Told, dx, dy, xx,yy)

% Gauss Seidel iterative method function for the 2D heat conduction equation
% solving code for steady state condition

tic;

gs_iter = 1;
while (error > tol)
    for i = 2: nx-1
        for j = 2:ny-1
            t1 = (T(i-1,j) + T(i+1,j))/(k*dx^2);
            t2 = (T(i,j-1) + T(i,j+1))/(k*dy^2);
            T(i,j) = t1 + t2;
        end
    end
    error = max(max(abs(Told - T)));
    Told = T;
    gs_iter = gs_iter + 1;
    sol_time_gs = toc; 
  
end

out = sol_time_gs;

% Plotting temperature distribution in domain
figure(2)
[a,b] = contourf(xx,yy,T);
clabel(a,b)
colormap(summer)
colorbar

set(gca,'ydir','reverse')
title(sprintf('Gauss Seidel Iteration = %d\n Iteration Time = %f sec', gs_iter, sol_time_gs))
xlabel('X Axis')
ylabel('Y Axis')

end

% Successive Over Relaxation Method function using Gauss-Seidel Iterative Method:

function [out] = sor_iter(error, tol,k, nx, ny,T,Told, dx, dy, xx,yy)

% Successive Over Relaxation (SOR) iterative method function for the 2D heat conduction equation
% solving code for steady state condition

% Gauss-Seidel method has been used.


omega = 1.2;

tic; % start timer
sor_iter = 1;
while (error > tol)
    for i = 2: nx-1
        for j = 2:ny-1
            t1 = (T(i-1,j) + T(i+1,j))/(k*dx^2);
            t2 = (T(i,j-1) + T(i,j+1))/(k*dy^2);
            T_GS = t1 + t2;
            T(i,j) = (1- omega)*(Told(i,j)) + omega*(T_GS);
        end
    end
    error = max(max(abs(Told - T)));
    Told = T;
    sor_iter = sor_iter + 1;
    sol_time_sor = toc; % stop timer & saving time into variable
   
end

out = sol_time_sor;

% Plotting temperature distribution in domain
figure(3)
[a,b] = contourf(xx,yy,T);
clabel(a,b)
colormap(winter)
colorbar

set(gca,'ydir','reverse')
title(sprintf('SOR Iteration = %d\n Iteration Time = %f sec', sor_iter, sol_time_sor ))
xlabel('X Axis')
ylabel('Y Axis')

end

Results:

     As it can be seen from the results images that number of iterations taken by Jacobi method is highest comapre to other two which is 217 while the gauss seidel method has solved the same problem in 100 less iterations and also in half the time. Jacobi method has taken 0.004 seconds while Gauss seidel method has taken 0.002 seconds. But the SOr method has solved the 2D heat conduction problem in just 78 iterations which is almost 1/3 of jacobi iterative method. Thats why SOR is very efficient method compare to other two.

1.2 Transient Condition - Implicit Method:

Main Code

% Solving 2D heat conduction equation
% Implicit Scheme
% Transient/Unsteady state condition

clear all
close all
clc

% inputs
% Grid points
nx = 10;
ny = nx;
l = 1; % length of square domain
nt = 1400; % number of time steps


% Creating mesh
x = linspace(0,l,nx); % Number of nodes in X - direction
y = linspace(0,l,ny); % Number of nodes in Y - direction
dx = x(2) - x(1); % step size in x
dy = dx; % step size in y
[xx,yy] = meshgrid(x,y); % mesh grid

% Initial value of the error
error = 9e9; 
tol = 1e-4; % tolerance value
dt = 0.01; % time step
alpha = 1.4; % diffusivity


% Calculating Courant numbers (CFL values)
k1 = alpha*(dt/(dx^2));
k2 = alpha*(dt/(dy^2));

% Boundary Conditions
T = ones(nx,ny);
% inserting temperature values at different side of the domain
T(:,1) = 400; % Left boundary 
T(1,:) = 600; % Top boundary
T(:,end) = 800; % Right boundary
T(end,:) = 900; % Bottom boundary

% Calcualating 4 corner nodes' temperature values
T(1,1) = (600 + 400)/2;
T(1,end) = (800 + 600)/2;
T(end,1) = (400 + 900)/2;
T(end,end) = (800 + 900)/2;

Tprev = T; % value of previous time step
Told = T; % value of previous iteration


% Defining loop to execute solvers 
for iter_solver = 1:3
    if iter_solver == 1 % for Jacobi Iterative method

        % Calling jacobi iterative method function
        jacobi_sol = jacobi_iter_transient(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt); 

    elseif iter_solver == 2 % for Gauss Seidel Iterative method
        
        % Calling Gauss Seidel iterative method function
        gs_sol = gs_iter_transient(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt ); % executing Gauss Seidel Iterative method

        % Calling Successive Over Relaxation iterative method function
    elseif iter_solver == 3 % for SOR Iterative method

        % Calling Successive Over Relaxation iterative method function
        sor_sol = sor_iter_transient(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt ); % executing SOR Iterative method
    end
end



% Jacobi Iterative Method Function:

 
function [out] = jacobi_iter_transient(error, tol,k1, k2,nx, ny,T,Told,Tprev, xx,yy, nt)

% Jacobi iterative method function for the 2D heat conduction equation
% solving code for transient condition.

% The function is part of the transient_heat_conduction_2d_implicit code.

tic; % start timer

jacobi_iter = 1;
for k = 1:nt
    while (error > tol)
        for i = 2: nx-1
            for j = 2:ny-1
                t1 = 1/(1 + 2*k1 + 2*k2);
                t2 = k1*t1;
                t3 = k2*t1;
                H = Told(i-1,j) + Told(i+1,j);
                V = (Told(i,j-1) + Told(i,j+1));
                T(i,j) = Tprev(i,j)*t1 + H*t2 + V*t3;
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        jacobi_iter = jacobi_iter + 1;
        sol_time_jacobi = toc; % stop timer & saving time into variable
    end
Tprev = T;
out = sol_time_jacobi;
end

figure(1)
[a,b] = contourf(xx,yy,T);
clabel(a,b)
colormap(autumn)
colorbar

set(gca,'ydir','reverse')
title(sprintf('Jacobi Iteration = %d\n Iteration Time = %f sec', jacobi_iter,sol_time_jacobi ))
xlabel('X Axis')
ylabel('Y Axis')

end


% Guass-Seidel Iterative Method Function:

function [out] = gs_iter_transient(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt)

% Gauss Seidel iterative method function for the 2D heat conduction equation
% solving code for transient condition.

% The function is part of the transient_heat_conduction_2d_implicit code.

tic;
gs_iter_transient = 1;
for k = 1:nt
    while (error > tol)
        for i = 2: nx-1
            for j = 2:ny-1
                t1 = 1/(1 + 2*k1 + 2*k2);
                t2 = k1*t1;
                t3 = k2*t1;
                H = Told(i-1,j) + Told(i+1,j);
                V = (Told(i,j-1) + Told(i,j+1));
                T(i,j) = Tprev(i,j)*t1 + H*t2 + V*t3;
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        gs_iter_transient = gs_iter_transient + 1;
        sol_time_gs = toc; 
    end
    
Tprev = T;
out = sol_time_gs;

end

figure(2)
[a,b] = contourf(xx,yy,T);
clabel(a,b)
colormap(summer)
colorbar

set(gca,'ydir','reverse')
title(sprintf('Gauss Seidel Iteration = %d\n Iteration Time = %f sec', gs_iter_transient, sol_time_gs))
xlabel('X Axis')
ylabel('Y Axis')

end


% Successive Over Relaxation Method function using Gauss-Seidel Iterative Method:

function [out] = sor_iter_transient(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt)

% Successive Over Relaxation (SOR) iterative method function for the 2D
% heat conduction equation solving code for transient condition.

% The function is part of the transient_heat_conduction_2d_implicit code.

% Gauss Seidel method has been used in support to SOR iterative method. 

omega = 1.7; % relaxatioan factor

tic; % start timer
sor_iter_transient = 1;
for k = 1:nt
    while (error > tol)
        for i = 2: nx-1
            for j = 2:ny-1
                t1 = 1/(1 + 2*k1 + 2*k2);
                t2 = k1*t1;
                t3 = k2*t1;
                H = T(i-1,j) + Told(i+1,j);
                V = (T(i,j-1) + Told(i,j+1));
                T(i,j) = (1-omega)*(Told(i,j))+ omega*(Tprev(i,j)*t1 + H*t2 + V*t3);
            end
        end
        error = max(max(abs(Told - T)));
        Told = T;
        sor_iter_transient = sor_iter_transient + 1;
        sol_time_sor = toc; % stop timer & saving time into variable
    end
    
Tprev = T;
out = sol_time_sor;

end

figure(3)
[a,b] = contourf(xx,yy,T);
clabel(a,b)
colormap(winter)
colorbar

set(gca,'ydir','reverse')
title(sprintf('SOR Iteration = %d\n Iteration Time = %f sec', sor_iter_transient, sol_time_sor ))
xlabel('X Axis')
ylabel('Y Axis')

end

Results:

    In transient condition, the number of iterations taken by Jacobi & Gauss Seidel method is similar but the Gauss seidel mehod has solved the equations in less than half the time taken by Jacobi method. The SOR method has again solved the equations in less iterations compare to other two. 

1.3 Transient - Explicit

Main Code:

% 2D heat conduction equations
% Explicit Scheme
% Transient/Unsteady state condition

clear all
close all
clc

% inputs
% Grid points
nx = 10;
ny = nx;
l = 1; % length of square domain
nt = 1400;


% Creating mesh
x = linspace(0,l,nx);
y = linspace(0,l,ny);
dx = x(2) - x(1); % step size in x
dy = dx; % step size in y
[xx,yy] = meshgrid(x,y); % mesh grid

% Initial value of the error
error = 9e9; 
tol = 1e-4; % tolerance value
dt = 0.001;
alpha = 1.4; % diffusivity


% Calculating Courant Numbers (CFL numbers)
k1 = alpha*(dt/(dx^2))
k2 = alpha*(dt/(dy^2))

% Boundary Conditions
T = ones(nx,ny);
% inserting temperature values at different side of the domain
T(:,1) = 400; % Left boundary 
T(1,:) = 600; % Top boundary
T(:,end) = 800; % Right boundary
T(end,:) = 900; % Bottom boundary


T(1,1) = (600 + 400)/2;
T(1,end) = (800 + 600)/2;
T(end,1) = (400 + 900)/2;
T(end,end) = (800 + 900)/2;

% Tprev = T; % value of previous time step
Told = T; % value of previous iteration


% % Defining loop to execute solvers 
% for iter_solver = 1:3
%     if iter_solver == 1 % for Jacobi Iterative method
% %         tic; % start timer
%         % executing jacobi iterative method
%         jacobi_sol = jacobi_iter_transient_explicit(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt,dy); 
% %         sol_time_jacobi = toc; % stop timer & saving time into variable
%     elseif iter_solver == 2 % for Gauss Seidel Iterative method
% %         tic;
%         gs_sol = gs_iter_transient_explicit(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt ); % executing Gauss Seidel Iterative method
% %         sol_time_gs = toc;
%     else
% %         iter_solver = 3; % for SOR Iterative method
% %         tic;
%         sor_sol = sor_iter_transient_explicit(error, tol,k1, k2, nx, ny,T,Told,Tprev, xx,yy, nt ); % executing SOR Iterative method
% %         sol_time_sor = toc;    
%     end
% end


tic; % start timer
explicit_method = 1;
while (error > tol)
    for k = 1:nt
        for i = 2: nx-1
            for j = 2:ny-1
%                 t1 = 1/(1 + 2*k1 + 2*k2);
%                 t2 = k1*t1;
%                 t3 = k2*t1;
%                 H = Told(i-1,j) + Told(i+1,j);
%                 V = (Told(i,j-1) + Told(i,j+1));
%                 T(i,j) = Tprev(i,j)*t1 + H*t2 + V*t3;
                  T(i,j) = Told(i,j) + k1*(Told(i-1,j) - 2*Told(i,j) + Told(i+1,j) + k2*(Told(i,j-1) - 2*Told(i,j) + Told(i,j+1)));
            end
        end
%         error = max(max(abs(Told - T)));
        error = max(abs(Told - T));
        Told = T;
        explicit_method = explicit_method + 1;
%         sol_time_explicit = toc; % stop timer & saving time into variable
    end
% Tprev = T;
% out = sol_time_explicit;
end

sol_time_explicit = toc; % stop timer & saving time into variable 

figure(1)
[a,b] = contourf(xx,yy,T);
clabel(a,b)
colormap(autumn)
colorbar

set(gca,'ydir','reverse')
title(sprintf('Explicit Method Iteration = %d\n Iteration Time = %f sec', explicit_method,sol_time_explicit ))
xlabel('X Axis')
ylabel('Y Axis')

Results:

CFL Numbers

Temperature Distribution

    The CFL values has been presented for the current problem which is below 0.5 and stood at 0.1134. 

    Compare to implicit method, the explicit method has taken almost 7 times the iteration of Jacobi method & 7 times the Gauss seidel method. The time taken to solved the equations is also quite high compare to implicit method. 

    Thats why, it can be said that the implicit method is very efficient compare to explicit method. But in implicit methods, the SOR is the most efficient.