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

Aim:To Solve the Steady and transient 2D heat conduction equation by using the point iterative techniques .

Objective:

To Solve the Steady 2D heat conduction equation by using the point iterative techniques .

To Solve the transient 2D heat conduction equation by using the point iterative techniques .

To write code for implicit used for solving transient 2D heat conduction equation by implementing the following methods for solving implicit equations.

1-Jacobi

2-Gauss-seidel

3-Successive over-relaxation

To write code for explicit schemes used for solving transient 2D heat conduction equation.

To compare the convergence rate [Number of iterations] of the steady state and transient simulations.

Theory:

The 2-D heat conduction equation is a partial differential equation which governs the heat transfer through a medium by thermal conduction. The equation is defined as:

`(del T)/(delt)-alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))`

For the steady state, the temperature does not vary with time, hence the time differential is equal to 0">00. The equation then gets modified to: `(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

Here, `alpha` is the thermal diffusivity which is characteristic for a medium.The domain over which the above equation

is solved is a square domain with unit length.The domain over which the above equation is solved is a square domain

with unit length. The domain over which the above equation is solved is a

square domain with unit length. It is divided into small elements to implement

the discretized schemes. For the purpose of this project, it is divided into an equal number of elements in

the nx and ny directions (nx=ny) The domain also has specified boundary conditions in order

to obtain a unique solution to the heat conduction equation. The boundary conditions are defined as:

Left boundary: T(x=0,y)=400K

Bottom boundary: T(x,y=0)=900K

Right boundary: T(x=L,y)=800K

Top boundary: T(x,y=L)=600K

The values at the corner nodes are taken as the average of the values from the surrounding edges. The equation is

solved for steady and transient state separately. The finite difference schemes solve for the temperatures at each node

throughout the domain using the data from the surrounding nodes and the boundary conditions applied on the domain.

Steady State:

The steady state equation is discretized using the central difference scheme in both x and y directions. Suppose i and jare the node numbers in the x and y directions respectively, then

This is the final equation that needs to be solved. We use three iterative solvers for the above equation, namely - Jacobi, Gauss-Seidel and Successive Over-Relaxation (SOR).

Each of the solvers solves the equation iteratively until the iterative error falls below the specified tolerance value, which is 10^-4 in this case.

The difference between the methods is the manner in which the values for the successive iterations are updated through the values from the previous iterations.

For the Jacobi solver, all the values in the RHS are taken from the previous iteration. The superscript k in the following equation denotes the iteration number.

For the Gauss-Seidel solver, the values that have already been updated in the current iteration are used and the rest of the values are taken from the previous iteration

The SOR solver is a modification of the Gauss-Seidel solver in which the values are updated according to the relaxation paramater. The relaxation parameter can take a range of values 0<ω<2 which

determine whether the solver is over-relaxed or under-relaxed.

The above solvers are implemented in MATLAB for the defined domain, and the steady state results are compared.

Matlab main code to call functions written for the solution of 2-D steady conduction equation:

% code to solve 2-D steady state conduction equation using various iterative solver
clc;
close all;
clear all;

% Inputs given for the domain 
H=1;    % Height of the domain in m
W=1;    % Width of the domain in m
nx=31;  % No of grid points along x 
ny=nx;  % No of grid points along y 

x=linspace(0,W,nx);      % Cordinates of all grids along width in Matrix
h=W/(nx-1);              % Grid spacing along width of the domain.
y=linspace(0,H,ny);      % Cordinates of all grids along hight in Matrix
dy=H/(ny-1);             % Grid spacing along the hight of the domain.

% Solution Initialization 

T=zeros(nx,ny);
T(1,2:end-1) = 900;    %Bottom
T(end,2:end-1) = 600;  %Top
T(2:end-1,1) = 400;    %Left
T(2:end-1,end) = 800;  %Right
T(1,1) = 650; T(1,end) = 850; T(end,end) = 700; T(end,1) = 500; %Corners
T(2:end-1,2:end-1)=300
w=1.1;

 % Calling function to compute time marching solution using Jacobi methods

tic   % simulation time computation
[T_j, err_j,count_j] = jacobi_steady_state(T,nx,ny);
time_j=toc;

%Calling function to compute time marching solution using Gauss siedal methods

tic
[T_g, err_g, count_g] = gauss_siedel_steady_state(T,nx,ny);
time_g=toc;

%  Calling function to compute time marching solution using SOR iterative methods
tic
[T_s, err_s, count_s] = SOS_steady_state(T,nx,ny,w);
time_s=toc;


% Plotting temperature on X and Y directions computed using various
% iterative methods

% Jacobi
figure(1)
contourf(x,y,T_j,'ShowText','on');
c = colorbar; colormap(jet); c.FontSize = 11;
set(get(c,'label'),'string','Contour: Temperature (K)','FontSize',14);
xlabel('X (m)','FontSize',14);
ylabel('Y (m)','FontSize',14);
title('Steady state temperature contour using jacobi ','FontSize',11)

% Gauss siedel
figure(2)
contourf(x,y,T_g,'ShowText','on');
c = colorbar; colormap(jet); c.FontSize = 11;
set(get(c,'label'),'string','Contour: Temperature (K)','FontSize',14);
xlabel('X (m)','FontSize',14);
ylabel('Y (m)','FontSize',14);
title('Steady state temperature contour using  gauss_siedel','FontSize',11)

% Successive over Relaxation
figure(3)
contourf(x,y,T_s,'ShowText','on');
c = colorbar; colormap(jet); c.FontSize = 11;
set(get(c,'label'),'string','Contour: Temperature (K)','FontSize',14);
xlabel('X (m)','FontSize',14);
ylabel('Y (m)','FontSize',14);
title('Steady state temperature contour using SOR ','FontSize',11)

% Comparing No of Iteration required to converge soution 
figure(4)
Iterative_method=categorical({'jacobi','gauss siedal','SOR'});
No_of_Iteration=[count_j count_g count_s];
bar(Iterative_method,No_of_Iteration)
xlabel('Iterative method')
ylabel('No of Iteration')
title('No of iteration for convergence in various Iterative Method ')
grid on

% Comparing Time required to converge soution in each iterative techniques
figure(5)
Iterative_method=categorical({'jacobi','gauss siedal','SOR'});
Time_of_convergence=[time_j time_g time_s];
bar(Iterative_method,No_of_Iteration)
xlabel('Iterative method')
ylabel('Time of convergence')
title('Time taken for convergence in various Iterative Method ')
grid on

Function code:

% Function to obtain output temperation error and No of Iteration using
% jacobi method 

function [T err count]=jacobi_steady_state(T,nx,ny)
count=0;        % initial no of iteration
err=1;          % initializing Error
tol=0.0001;      % Tolerance

% while loop to check error
while (err >tol)
% Store the temperature values from previous iteration
Told=T;
% Space loop
    for i=2:nx-1
        for j=2:ny-1
            T(i,j) = 0.25*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));
        end
    end
     % updating the error and updating no of iteration after each while loop
    count=count+1
    err = max(max(abs(T-Told)));
end

% Function to obtain output temperation error and No of Iteration using
% Gauss siedel method

function [T,err,count]=gauss_siedel_steady_state(T,nx,ny)
count=0;     % initial no of iteration
err=1;       % initializing Error
tol=0.0001;  % Tolerance

  % While Loop
while (err>tol)
Told=T;      % Store the temperature values from previous iteration
% Space loop
    for i=2:nx-1
        for j=2:ny-1
            T(i,j) = 0.25*(T(i-1,j) + Told(i+1,j) + T(i,j-1) + Told(i,j+1));
        end
    end
   % updating the error and updating no of iteration after each while loop
    count=count+1
    err = max(max(abs(T-Told)));
end

% Function to obtain output temperation error and No of Iteration using
% SOR method


function [T err count]=SOS_steady_state(T,nx,ny,w)
count=0;     % initial no of iteration
err=1;       % initializing Error
tol=0.0001;  % Tolerance

w=1.1;
while (err>tol)
 % Store the temperature values from previous iteration
 Told=T; 

    % Space loop
    for i=2:nx-1
        for j=2:ny-1
            T(i,j) = (1-w)*Told(i,j) + w*0.25*(T(i-1,j) + Told(i+1,j)+ T(i,j-1) + Told(i,j+1));
        end
    end
   % updating the error and updating no of iteration after each while loop
    count=count+1
    err = max(max(abs(T-Told)));
   Told=T; 
end

Temperature contours:

Bar graph for time and No of iteration comparision:

Transient simulation :

1-Explicit scheme:

In the transient equation, the time derivative does not disappear and we use the entire heat conduction equation. We discretize the time derivative using forward difference scheme. The transient scheme

can be implicit or explicit depending on the time-step at which the spatial derivatives are chosen. First, we will focus on the explicit scheme where all the spatial derivatives are taken from the previous

time step. The heat conduction equation, after discretization, modifies to:

In this case, since Δx=Δy=h

Here subscript n denotes the time step at which the value is taken .let `K1=(alpha DeltaT)/h^2` then

For the purpose of stability, the coefficient of the term `T_(i,j)^(n+1` in the RHS cannot be negative. This comes from the Courant-Friedrich-Lewy (CFL) condition for stability.

The term `K1=(alpha DeltaT)/h^2``is called the CFL number. From this condition, we can calculate the maximum time-step that can be chosen for solving the transient heat conduction equation explicitly.

`1 − 4 K1 ≥ 0` `K 1 ≤ 1 /4`

`(α Δ t) /h^2 < 1 /4`

`Δ t ≤ h^2/(4 α)`

The above condition must be satisfied in order to obtain convergence of the solution. As a borderline case, `Deltat=h^2/(4*alpha)`is chosen for the explicit solver.

Matlab main code to call functions written for the solution of 2-D Transient conduction equation in Explicit Scheme:

% code to solve 2-D transient conduction equation in Explicit scheme using
% various iterative solver
clc;
close all;
clear all;
% Inputs given for the domain 
H=1;    % Height of the domain in m
W=1;    % Width of the domain in m
nx=31;  % No of grid points along x 
ny=nx;  % No of grid points along y 

x=linspace(0,W,nx);      % Cordinates of all grids along width in Matrix
h=W/(nx-1);              % Grid spacing along width of the domain.
y=linspace(0,H,ny);      % Cordinates of all grids along hight in Matrix
dy=H/(ny-1);             % Grid spacing along the hight of the domain.

t_start=0;               %start time in s
t_end=0.35;             % End time of simulation in s
alpha=1.4;
dt=h^2/(4*alpha);                   % Time step size
t = t_start:dt:t_end;    % Time discritization
nt=length(t);            % Total no of time step

tol=0.0001;              % Tolrance value for error
K1 = alpha*dt/(h^2);     % CFL number since dx=dy=h

% Solution Initialization 

T=zeros(nx,ny);
T(1,2:end-1) = 900; %Bottom
T(end,2:end-1) = 600; %Top
T(2:end-1,1) = 400; %Left
T(2:end-1,end) = 800; %Right
T(1,1) = 650; T(1,end) = 850; T(end,end) = 700; T(end,1) = 500; %Corners


error = 1;
count=0;
Told = T;
%time marching solution

tic;

for k=1:nt    %Time loop

        % Storing Temperature from previous time step
       while(error > tol)
        
     % splace loop
    for i=2:nx-1
        for j=2:ny-1
             T(i,j)=(1-4*K1)*Told(i,j)+K1*(Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1));
        end
    end
    
    error = max(max(abs(Told - T)));
    Told = T;
    count=count+1;
     
       end
 
end

 % updating  no of iteration after each while loop
  
time_c=toc;
% Plotting temperature on X and Y directions computed using Explicit scheme

figure(1)
contourf(x,y,T,'ShowText','on'); ax = gca; ax.FontSize = 11;
c = colorbar; colormap(jet); c.FontSize = 11;
set(get(c,'label'),'string','Contour: Temperature (K)','FontSize',14);
xlabel('X (m)','FontSize',14); ylabel('Y (m)','FontSize',14);
  str = {['Temperature contour at t = ',num2str(t_end),' using Explicit scheme']};
title(str,'FontSize',11);

Temperature contour:

TRANSIENT STATE (IMPLICIT):

For the implicit solver, the terms T(i,j) on the RHS are taken at the time-step n+1 as opposed to n in the explicit solver. This changes the discretization of the transient equation to

substituting K1=`(alpha*Deltat)/h^2`

As we can see above, the coefficient of the term `T^n(i,j)` on the RHS can never be negative. Hence, the implicit scheme is unconditionally stable for any value of K1 Since the values on the RHS are also

unknowns, we use iterative solvers to implicitly solve the heat conduction equation. The iterative solvers - Jacobi, Gauss-Seidel, and SOR, are used similarly as described in the steady-state section.

For Jacobi method,

For Gauss-siedal methad,

For SOR method,

The superscript n denotes the time-step and k denotes the iteration within that time-step. The terms missing the superscript k are the values taken from the previous time-step which are constant throughout the current time-step and do not iterate ("T_prev_dt"). Since the implicit scheme is unconditionally stable, a time-step of 0.01 seconds is chosen which is much higher than the time-step for the explicit scheme.

Matlab main code to call functions written for the solution of 2-D Transient conduction equation in Implicit Scheme:

% code to solve 2-D transient conduction equation in Implicit scheme using
% various iterative solver
clc;
close all;
clear all;
% Inputs given for the domain 
H=1;    % Height of the domain in m
W=1;    % Width of the domain in m
nx=31;  % No of grid points along x 
ny=nx;  % No of grid points along y 
x=linspace(0,W,nx);      % Cordinates of all grids along width in Matrix
h=W/(nx-1);              % Grid spacing along width of the domain.
y=linspace(0,H,ny);      % Cordinates of all grids along hight in Matrix
dy=H/(ny-1);             % Grid spacing along the hight of the domain.
t_start=0;               %start time in s
t_end=12000;             % End time of simulation in s
dt=10;                   % Time step size
t = t_start:dt:t_end;    % Time discritization
nt=length(t);            % Total no of time step
tol=0.0001;              % Tolrance value for error
alpha=1.4;
K1 = alpha*dt/(h^2);     % CFL number since dx=dy=h
 
% Solution Initialization 
T=zeros(nx,ny);
T(1,2:end-1) = 900;    %Bottom
T(end,2:end-1) = 600;  %Top
T(2:end-1,1) = 400;    %Left
T(2:end-1,end) = 800;  %Right
T(1,1) = 650; T(1,end) = 850; T(end,end) = 700; T(end,1) = 500; %Corners
T_prev_dt=T;
Told=T;
 % Calling function to compute time marching solution using Jacobi methods
tic;   % simulation time computation
      
[T_j, err_j count_j] = jacobi_implicit(T,T_prev_dt,nt,nx,ny,K1);
time_Jacobi=toc;
%Calling function to compute time marching solution using Gauss siedal methods
tic;
[T_g, err_g count_g] = gauss_siedel_implicit(T,T_prev_dt,nt,nx,ny,K1);
time_gauss=toc;
%  Calling function to compute time marching solution using SOR iterative methods
w=1.1;
tic;
[T_SOR, err_SOR count_SOR] = SOR_implicit(T,T_prev_dt,nt,nx,ny,K1,w)
time_SOR=toc;
% Plotting temperature on X and Y directions computed using various
% iterative methods
% Jacobi
figure(1)
contourf(x,y,T_j,'ShowText','on')
c = colorbar; colormap(jet); c.FontSize = 11;
set(get(c,'label'),'string','Contour: Temperature (K)','FontSize',14);
xlabel('X (m)','FontSize',14);
ylabel('Y (m)','FontSize',14);
title('Steady state temperature contour using jacobi ','FontSize',11)
% Gauss siedel
figure(2)
contourf(x,y,T_g,'ShowText','on')
c = colorbar; colormap(jet); c.FontSize = 11;
set(get(c,'label'),'string','Contour: Temperature (K)','FontSize',14);
xlabel('X (m)','FontSize',14);
ylabel('Y (m)','FontSize',14);
title('Steady state temperature contour using gauss siedel ','FontSize',11)
% Successive over Relaxation
figure(3)
contourf(x,y,T_SOR,'ShowText','on')
c = colorbar; colormap(jet); c.FontSize = 11;
set(get(c,'label'),'string','Contour: Temperature (K)','FontSize',14);
xlabel('X (m)','FontSize',14);
ylabel('Y (m)','FontSize',14);
title('Steady state temperature contour using SORl ','FontSize',11)
% Comparing No of Iteration required to converge soution 
figure(4)
Iterative_method=categorical({'jacobi','gauss siedal','SOR'});
No_of_Iteration=[count_j count_g count_SOR];
bar(Iterative_method,No_of_Iteration)
xlabel('Iterative method')
ylabel('No of Iteration')
title('No of iteration for convergence in various Iterative Method used in Implicit scheme')
grid on
% Comparing Time required to converge soution in each iterative techniques
figure(5)
Iterative_method=categorical({'jacobi','gauss siedal','SOR'});
Time_of_convergence=[time_Jacobi time_gauss time_SOR];
bar(Iterative_method,No_of_Iteration)
xlabel('Iterative method')
ylabel('Time of convergence')
title('Time taken for convergence in various Iterative Method used in Implicit scheme')
grid on
 

Function code:

% Function to obtain output temperation error and No of Iteration using
% jacobi method 
function [T, err count] = jacobi_implicit(T,T_prev_dt,nt,nx,ny,K1)

count=0;        % initial no of iteration
err=1;          % initializing Error
tol=0.0001;      % Tolerance

% Time Loop
for k=1:nt
     % Storing Temperature from previous time step
       T_prev_dt=T;
    while (err>=tol)
    % Store the temperature values from previous iteration
   Told = T;
    % Space loop
    for i=2:nx-1
        for j=2:ny-1
            
           T(i,j) = (1/(1+4*K1))*T_prev_dt(i,j) + (K1/(1+4*K1))*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));
        end
    end
    % updating the error and updating no of iteration after each while loop
   
   err = max(max(abs(Told-T)));
   count=count+1;
 
   
    end
   
end

% Function to obtain output temperation error and No of Iteration using
% Gauss siedel method

function [T, err count] = gauss_siedel_implicit(T,T_prev_dt,nt,nx,ny,K1)

count=0;     % initial no of iteration
err=1;       % initializing Error
tol=0.0001;  % Tolerance

% Time Loop
for k=1:nt
    % Storing Temperature from previous time step
       T_prev_dt=T;
       % While Loop
    while (err>=tol)
    % Store the temperature values from previous iteration
    Told = T;
    % Space loop
    for i=2:nx-1
        for j=2:ny-1
            
           term1 = 1/(1+(4*K1));
               term2 = K1*term1;
            
               h = (T((i-1),j)+Told((i+1),j));
               v = (T(i,(j-1))+Told(i,(j+1)));
               T(i,j) =(T_prev_dt(i,j)*term1)+(h+v)*term2;
        end
    end
    % updating the error and updating no of iteration after each while loop
    count=count+1
   err = max(max(abs(Told-T)));
   
    end
   
    end
   

% Function to obtain output temperation error and No of Iteration using
% SOR method


function [T err count]=SOR_implicit(T,T_prev_dt,nt,nx,ny,K1,w);
count=0;     % initial no of iteration
err=1;       % initializing Error
tol=0.0001;  % Tolerance

% Time Loop
for k=1:nt
    % Storing Temperature from previous time step
       T_prev_dt=T;
   
    while (err>=tol)
    % Store the temperature values from previous iteration
    Told = T;
    % Space loop
    for i=2:nx-1
        for j=2:ny-1
            term1 = 1/(1+(4*K1));
               term2 = K1*term1;
             
               h = T(i-1,j)+Told(i+1,j);
               v = (T(i,j-1)+Told(i,j+1));
               T(i,j) =(T_prev_dt(i,j)*term1)+(h+v)*term2;
               T(i,j) = Told(i,j)*(1-w)+w*T(i,j);
            
          
        end
    end
     % updating the error and updating no of iteration after each while loop
    err = max(max(abs(Told-T)));
     count=count+1;
    end
   
    end

Temperature contours:

Bar graph for time and No of iteration comparision:

Conclusion:

In steady state simultion convergence rate of SOR iterative method is highest while in unsteady case Explicit scheme conveges solution faster.

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.

Time Taken in explicit scheme is least due to uncoupled Equation.But because of the stability criteria it limits the time step size incresing cost of simulation.

Steady state simulation takes less time and no of iteration to converge the solution.