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

Week 5 Assignment 1. Steady state analysis & Transient State Analysis for 2D heat conduction equation .

Steady state analysis & Transient State Analysis

What is steady and unsteady State Heat Transfer?

Heat transfer is the transfer of thermal energy from a body, at a high temperature, to another at a lower temperature. This transfer of thermal energy may occur under steady or unsteady state conditions. Under Steady state conditions the temperature within the system does not change with time. Conversely, under unsteady state conditions the temperature within the system does vary with time.

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.

Clearly, no system can remain under unsteady state conditions perpetually. The temperature of the system will eventually reach the temperature of the heat source, and once this happens, the system will be at steady state. Even if the amount of heat being transferred into the system is increased, at some point the system reaches its critical temperature and the energy transferred into it the starts causing phase changes within the system rather than temperatures increases.

2D heat conduction equation is given by:

Transient state form:-

where α is the thermal diffusivity of material.

Steady State form:-

Choose the type of solver to be computed:-

1. Jacobi steady state,
2. GS steady state,
3. SOR steady state,
4. Jacobi transient explicit,
5. GS transient explicit,
6. SOR transient explicit,
7. Jacobi transient implicit,
8. GS transient implicit,
9. SOR transient implicit.

Corresponding number to be entered in the command window.

Matlab code for steady and transient state analysis:-

clear all
close all
clc

% Given Conditions 
nx = 10;
ny = 10;

% Domain is unit square along x and y
 x = linspace(0,1,nx);
 y = linspace(0,1,ny);
 dx = 1/(nx-1);
 dy = 1/(ny-1);
  % Time step dt 
 dt = 1e-2;
 nt = 2/dt;
 % let the initial temp be 100k
 T = 100*ones(nx,ny);
 
 % BC
 T(:,1)  =   400; % Left
 T(:,nx) =   800; % Right
 T(1,:)  =   600; % Top
 T(ny,:) =   900; % Bottom
 
 % Tolerance and error
 error = 9e9;
 tol = 1e-4;
 
 % STEADY STATE ANALYSIS
 % EXPLICIT METHOD
 T_old = T;
 T_initial = T;
 type_of_solver= input('Choose the type of solver to be computed:- 1-Jacobi steady state, 2- GS steady state, 3-SOR steady state, 4-Jacobi transient explicit, 5-GS transient explicit, 6-SOR transient explicit, 7-Jacobi transient implicit, 8- GS transient implicit, 9-SOR transient implicit. number =');

 alpha = 1.25;
 K = 2*(dx^2+dy^2)/((dx^2)*(dy^2));
% Transient state analysis
%Explcit
beta = 0.25;
m1 = beta*dt/(dx^2);
m2 = beta*dt/(dy^2);
m = 1-(2*m1)-(2*m2);
%Implicit
L1 = (beta*dt/(dx^2));
L2 = (beta*dt/(dy^2));
L = 1+(2*L1)+(2*L2);
   
% Steady State Method
 % JACOBI
 
 tic;
 if type_of_solver == 1
     jacobi_iter=1;
     while (error>tol)
         for i = 2:nx-1
             for j = 2:ny-1
                 T(i,j) = (K^-1)*((T_old(i-1,j)+T_old(i+1,j))/(dx^2)+(T_old(i,j+1)+T_old(i,j-1))/(dy^2));
             end
         end
         error = max(max(abs(T-T_old)));
         T_old = T;
         jacobi_iter = jacobi_iter+1;
         toc
         simulation_steady_jacobi = toc;
     end
     figure(1)
     [X,Y] = contourf(x,y,T);
     clabel(X,Y)
     colormap(jet)
     colorbar
     xlabel('X-Axis')
     ylabel('Y-Axis')
     title_text =sprintf('Jacobi steady state Iteration =%d, simulation time =%0.3fs', jacobi_iter, simulation_steady_jacobi);
     title(title_text)
     pause(0.0005)
 end
 
 % Gauss Siedel Method

 tic;
 if type_of_solver ==2 
     gauss_iter = 1;
     while (error>tol)
         for i = 2:nx-1
             for j = 2:ny-1
                 T(i,j) = (K^-1)*((T(i-1,j)+T_old(i+1,j))/(dx^2)+(T_old(i,j+1)+T(i,j-1))/(dy^2));
             end
         end
         error = max(max(abs(T-T_old)));
         T_old = T;
         gauss_iter = gauss_iter+1;
         toc
         simulation_steady_gauss = toc;
     end
     figure(2)
     [X,Y] = contourf(x,y,T);
     clabel(X,Y)
     colormap(jet)
     colorbar
     xlabel('X-Axis')
     ylabel('Y-Axis')
     title_text =sprintf('GS steady state Iteration =%d, Simulation time =%0.3fs', gauss_iter, simulation_steady_gauss);
     title(title_text)
     pause(0.0005)
 end
 
 % Successive over relaxation
 
 tic;
 if type_of_solver == 3;
     sor_iter=1;
     while (error>tol)
         for i = 2:nx-1
             for j = 2:ny-1
                 k1 = (T(i-1,j)+T_old(i+1,j))/(dx^2);
                 k2 = (T_old(i,j+1)+T(i,j-1))/(dy^2);
                 T_gauss(i,j) = (k1+k2)/K;
                 T(i,j) = T_old(i,j)*(1-alpha)+alpha*(T_gauss(i,j));
             end
         end
         error = max(max(abs(T-T_old)));
         T_old = T;
         sor_iter = sor_iter+1;
         toc
         simulation_steady_sor=toc;
     end
     
     figure(3)
     [X,Y] = contourf(x,y,T);
     clabel(X,Y)
     colormap(jet)
     colorbar
     xlabel('X-Axis')
     ylabel('Y-Axis')
     title_text =sprintf('SOR steady stateIteration =%d, Simulation time =%0.3fs', sor_iter, simulation_steady_sor);
     title(title_text)
     pause(0.0005)
 end
 % TRANSIENT EXPLICIT METHOD

% JACOBI

tic
jacobi_iter = 1;
if type_of_solver == 4
    for s = 1:nt
        while (error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    T(i,j) = T_old(i,j)*m+m1*(T_old(i+1,j)+T_old(i-1,j))+m2*(T_old(i,j+1)+T_old(i,j-1));
                end
            end
            error = max(max(abs(T-T_old)));
            jacobi_iter = jacobi_iter+1
            T_old = T;
        end
    end
    toc
    transient_explicit_jacobi = toc;
    
    % Plot
    figure(4)
    [X,Y] = contourf(x,y,T);
    clabel(X,Y)
    colormap(jet)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title_text =sprintf('Jacobi transient explicit Iteration =%d, simulation time =%0.3fs', jacobi_iter, transient_explicit_jacobi);
     title(title_text)
     pause(0.0005)
end

% Gauss Siedel 

tic
if type_of_solver == 5
    gauss_iter = 1;
    for s = 1:nt
        while (error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    T(i,j) = T_old(i,j)*m+m1*(T_old(i+1,j)+T(i-1,j))+m2*(T_old(i,j+1)+T(i,j-1));
                end
            end
            error = max(max(abs(T-T_old)));
            gauss_iter = gauss_iter+1
            T_old = T;
        end
    end
    toc
    transient_explicit_gauss = toc;
    
    % Plot
    figure(5)
    [X,Y] = contourf(x,y,T);
    clabel(X,Y)
    colormap(jet)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title_text =sprintf('GS transient explicit Iteration =%d, simulation time =%0.3s', gauss_iter, transient_explicit_gauss);
     title(title_text)
     pause(0.0005)
end

% SOR

tic
if type_of_solver == 6
    sor_iter = 1;
    for s = 1:nt
        while (error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    T(i,j) = T_old(i,j)*(1-alpha)+alpha*(T_old(i,j)*m+m1*(T_old(i+1,j)+T(i-1,j))+m2*(T_old(i,j+1)+T(i,j-1)));
                end
            end
            error = max(max(abs(T-T_old)));
            sor_iter = sor_iter+1
            T_old = T;
        end
    end
    toc
    transient_explicit_sor = toc;
    
    % Plot
    figure(6)
    [X,Y] = contourf(x,y,T);
    clabel(X,Y)
    colormap(jet)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title_text =sprintf('SOR steady state Iteration =%d, simulation time =%0.3fs', sor_iter, transient_explicit_sor);
     title(title_text)
     pause(0.0005)
end
%Transient Implicit Method
% JACOBI

tic
jacobi_iter_implicit = 1;
if type_of_solver == 7
    for s = 1:nt
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                   T_left = T_old(i-1,j);
                   T_right = T_old(i+1,j);
                   T_top = T_old(i,j+1);
                   T_bottom = T_old(i,j-1);
                   T(i,j) = ((T_initial(i,j)+L1*(T_left+T_right)+L2*(T_top+T_bottom))/L);                   
                end
            end
            error = max(max(abs(T-T_old)));
            T_old = T;
            jacobi_iter_implicit = jacobi_iter_implicit+1
            
        end
    end
    toc
    transient_implicit_jacobi = toc;
    
    % Plot
    figure(7)
    [X,Y] = contourf(x,y,T);
    clabel(X,Y)
    colormap(jet)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title_text =sprintf('Jacobi transient implicit Iteration =%d, simulation time =%0.3fs', jacobi_iter_implicit, transient_implicit_jacobi);
    title(title_text)
    pause(0.0005)
end

% GAUSS SIEDEL

tic
gauss_iter_implicit = 1;
if type_of_solver == 8
    for s = 1:nt
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                   T_left = T(i-1,j);
                   T_right = T(i+1,j);
                   T_top = T(i,j+1);
                   T_bottom = T(i,j-1);
                   T(i,j) = ((T_initial(i,j)+L1*(T_left+T_right)+L2*(T_top+T_bottom))/L );                   
                end
            end
            error = max(max(abs(T-T_old)));
            T_old = T;
            gauss_iter_implicit = gauss_iter_implicit+1
            
        end
    end
    toc
    transient_implicit_gauss = toc;
    
    % Plot
    figure(8)
    [X,Y] = contourf(x,y,T);
    clabel(X,Y)
    colormap(jet)
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title_text =sprintf('GS transient implicit iteration =%d, simulation time =%0.3fs', gauss_iter_implicit, transient_implicit_gauss);
    title(title_text)
    pause(0.0005)
end
    
% SOR

tic
 sor_iter_implicit = 1;
if type_of_solver == 9
  sor_iter_implicit = 1;
    for s = 1:nt
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                   T_left = T(i-1,j);
                   T_right = T(i+1,j);
                   T_top = T(i,j+1);
                   T_bottom = T(i,j-1);
                   T(i,j) = (1-alpha)*(T_old(i,j))+alpha*((T_initial(i,j)+L1*(T_left+T_right)+L2*(T_top+T_bottom))/L);                   
                end
            end
            error = max(max(abs(T-T_old)));
            T_old = T;
            sor_iter_implicit = sor_iter_implicit+1
            
        end
    end
    toc
    transient_implicit_sor = toc;
    
    % Plot
    figure(9)
    [X,Y] = contourf(x,y,T);
    clabel(X,Y)
    colormap(jet) 
    colorbar
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title_text =sprintf('SOR transient implicit Iteration =%d, simulation time =%0.3fs', sor_iter_implicit, transient_implicit_sor);
    title(title_text)
    pause(0.0005)
end

Results-

The following graphs are the results for the above iterative techniques.

Plots for steady state

1) Jacobi steady state

2) GS steady state

3) SOR steady state.

Plots for Transient state analysis.

Transient Explicit

4) Jacobi transient explicit

5) GS transient explicit

6) SOR transient explicit.

Transient Implicit

7) Jacobi transient Implicit.

8) GS transient Implicit

9) SOR transient implicit

Observations and Conclusions:-

• It can be observed that Transient (explicit) takes more time to compute compared to steady state methods.
• Transient Explicit method has more number of iterations when compared t steady state methods. this is because of the addition of time derivative which increases the number of time steps in the given solution.
• Implicit solutions are unconditionally stable. This means any value of dt, no matter how large will result in a stable solution.
• In all cases it can be observed that Simulation time for Jacobi> Gauss siedle> SOR.
• Jacobi takes more time as it dosenot update the database after each iteration.
• Gauss Siedle takes the latest value and computes which makes it faster than Jacobi method.
• Transient Implicit methods take the least amount of time.
• SOR method takes the minimum time of all as it is an enhanced version of Gauss Siedle and is devised by applying extrapolation to the Gauss-Seidel method.This extrapolation takes the form of a weighted average between the previous iterate and the computed iterate successively for each component.