Steady and Unsteady 2D Heat Conduction

STEADY STATE ANALYSIS FOR 2D HEAT CONDUCTION

1. JACOBI METHOD

In the jacobi method the value of the (x_i)^k  obtained in the k^th iteration remain unchanged until the entire (k+1)^th iteration has been calculated. The value during iteration are used are from the previous iteration so the number of iteration is higher when compared to other methods.

2. GAUSS SEIDAL METHOD

With the Jacobi method, the values of (x_i)^k obtained in the k th iteration remain unchanged until the entire (k+1)th iteration has been calculated. With the Gauss-Seidel method, we use the new values (x_i)^k+1 as soon as they are known. For example, once we have computed (x₁)^k+1 from the first equation, its value is then used in the second equation to obtain the new (x_i)^k+1 and so on.

3. SUCCESSIVE OVER RELAXATION METHOD 

         a) Using GAUSS SEIDAL

This method is an advanced version of the Gauss-Seidel method. It uses the result of the Gauss-Seidel method and scales the convergence rate by a factor \'omega\' called the relaxation factor. The term \'over-relaxation\' means that the relaxation factor is greater than 1 in this method.

The SOR iteration formula can be given by

Tk+1(i,j)=Tk(i,j)*(1−ω)+ω*T_gs(i,j)

where T_gs is the result of the Gauss-Seidel method. The above formula means that when the value of omega increases from  1, we are moving from no update to update using the Gauss-Seidel method. When omega > 1, we\'re scaling-up the convergence rate.

         b) Using JACOBI

This method is an advanced version of the jacobi method. It uses the result of the jacobi method and scales the convergence rate by a factor \'omega\' called the relaxation factor. The term \'over-relaxation\' means that the relaxation factor is greater than 1 in this method.

The SOR iteration formula can be given by

Tk+1i,j=Tki,j(1−ω)+ω T_j(i,j)

where T_j is the result of the Jacobi method. The above formula means that when the value of omega increases from  1, we are moving from no update to update using the Jacobi method. When omega > 1, we\'re scaling-up the convergence rate.

MATLAB PROGRAM FOR STEADY STATE :

% Code to solve a 2D HEAT CONDUCTION

clear all
close all
clc

% inputs
nx = 20;             % number of points along x 
ny = nx;             % number of along x & y are same
T = 298*ones(nx,ny); % assume the initial temp as 298

% Assigning BC
T(1,1:end) = 900;   % bottom condition
T(end,1:end) = 600; % top condition
T(1:end,1) = 400;   % left condition
T(1:end,end) = 800; % right condition

% creating a copy of T
Told = T;

x = linspace(0,1,nx);  
y = linspace (0,1,ny); 
dx = x(2)-x(1);        
dy = y(2)-y(1);
[xx, yy] = meshgrid(x,y); % create 2D vector
tol = 1e-4;               % tolerence limit
error = 9e9;              % error vaule set

% calculation simplification
k = (2*(dx^2+dy^2)/(dx^2*dy^2));

% To choose method of iteration to be performed
prompt= \'Which method to be used(1. Jacobi 2.GS 3.SOR+GS 4.SOR+Jacobi)?\';   
iterative_solver=input(prompt);

% Jacobi method

if iterative_solver == 1
  tic
  jacobi_iter=1;
  while(error > tol)
     for i=2:nx-1
         for j=2:ny-1
          T(i,j) = (1/k)*(((Told(i-1,j)+Told(i+1,j))/dx^2) + ((Told(i,j-1)+Told(i,j+1))/dy^2));
         end
     end
     error=max(max(abs(Told-T)));
     Told=T;
     jacobi_iter=jacobi_iter+1;
     figure(1)  
     [c,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     title_text=sprintf(\'Jacobi iteration number=%d @steady state\',jacobi_iter);
     title(title_text)
     xlabel(\'X-axis\')
     ylabel(\'y-axis\')
     clabel(c,d)
     pause(0.0001)
  end
Jacobi_time = toc
end

% Gauss Seidal method
if iterative_solver == 2
    tic
    GS_iter=1;
     while(error > tol)
     for i=2:nx-1
         for j=2:ny-1
          T(i,j) = (1/k)*(((T(i-1,j)+Told(i+1,j))/dx^2) + ((Told(i,j+1)+T(i,j-1))/dy^2));
         end
     end
     error=max(max(abs(Told-T)));
     Told=T;
     GS_iter=GS_iter+1;
     figure(1)  
     [c,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     title_text=sprintf(\'Gauss iteration number=%d @steady state\',GS_iter);
     title(title_text)
     xlabel(\'X-axis\')
     ylabel(\'y-axis\')
     clabel(c,d)
     pause(0.03)
     end
Gauss_time = toc;
end   
    
% Succesive over relaxation with Gauss seidal
if iterative_solver == 3
  tic
  SOR_GS_iter=1;
  T_GS = zeros();
  l = 1.724;
     while(error > tol)
        for i=2:nx-1
           for j=2:ny-1
             T_GS(i,j) = (1/k)*(((T(i-1,j)+Told(i+1,j))/dx^2) + ((Told(i,j+1)+T(i,j-1))/dy^2));
             T(i,j) = Told(i,j)*(1-l) + l*T_GS(i,j); 
           end
        end
     error=max(max(abs(Told-T)));
     Told=T;
     SOR_GS_iter = SOR_GS_iter+1;
     figure(1)  
     [c,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     title_text=sprintf(\'SOR by Gauss iteration number=%d @steady state\',SOR_GS_iter);
     title(title_text)
     xlabel(\'X-axis\')
     ylabel(\'y-axis\')
     clabel(c,d)
     pause(0.000001)
     end
SOR_GS_time = toc
end


% Successive over relaxation with jacobi
if iterative_solver == 4
   tic
   SOR_J_iter=1;
   T_J = zeros();
   l1 = 1.01;
     while(error > tol)
     for i=2:nx-1
         for j=2:ny-1
             T_J(i,j) = (1/k)*(((Told(i-1,j)+Told(i+1,j))/dx^2) + ((Told(i,j+1)+Told(i,j-1))/dy^2));
         T(i,j) = Told(i,j)*(1-l1) + l1*T_J(i,j);
         end
     end
     error=max(max(abs(Told-T)));
     Told=T;
     SOR_J_iter=SOR_J_iter+1;
     figure(1)  
     [c,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     title_text=sprintf(\'SOR by Jacobi iteration number=%d @steady state\',SOR_J_iter);
     title(title_text)
     xlabel(\'X-axis\')
     ylabel(\'y-axis\')
     clabel(c,d)
     pause(0.00000001)
     end
SOR_J_time = toc
end

PLOT OF VARIOUS ITERATIONS FOR STEADY STATE :

1. JACOBI METHOD

2 . GAUSS SEIDAL METHOD

3. SUCCESSIVE OVER RELAXATION WITH GAUSS SEIDAL METHOD

4. SUCCESSIVE OVER RELAXATION WITH JACOBI METHOD

RESULTS:

JACOBI TIME PERIOD

GAUSS SEIDAL TIME PERIOD

SOR using GAUSS TIME PERIOD 

SOR using JACOBI TIME PERIOD

CONCLUSION:

For Steady State solver, we compare the results of 4 iterative methods i.e. Jacobian, Gauss-Seidel and Successive Over-Relaxation methods which are being presented in the table above. Total 20 no. of grids and the assumed temperatures for all 4 sides of thin fluid plate were applied, we got the converged temperature distribution over the plate. The parameters like no. of iterations and the converging time was during the simulation with the help of tic-toc command in the programme. Following the stats or results, we clearly see that Successive Over-Relaxation iterative method is the one which takes least no. of iterations and Converging time and hence, the Computational time which plays an important role in CFD. SOR with gauss seidal takes a total of 57 iterations and 12.182 seconds & SOR with jacobi takes a total of 819 iterations and 137.013 seconds whereas the Jacobian iterative method takes the most with 827 iterations and 137.013 seconds to converge. Gauss-Seidel iterative method stays in the middle with 441 no. of iterations and 83.071 seconds as the converging time.

TRANSIENT STATE ANALYSIS FOR 2D HEAT CONDUCTION USING IMPLICIT METHOD

1. JACOBI METHOD

2. GAUSS SEIDAL METHOD

3. SUCCESSIVE OVER RELAXATION Using GAUSS SEIDAL

MATLAB PROGRAM FOR UNSTEADY STATE :

% implicit using Jacobi, Gauss-Seidel, and SOR iterative solvers for transient state

clear all
close all
clc
                                    
nx =20;                     % number of points along x 
ny = nx;                    % number of along x & y are same
dt = 1e-2;                  % time step size
x = linspace(0,1,nx);       % x nodes
y = linspace(0,1,ny);       % y nodes
dx = x(2)-x(1);             % gird size along x
dy = y(2)-y(1);             % grid size along y
alpha = 1.4;                % thermal diffusivity 
[xx yy] = meshgrid(x,y)

%Initialization
T = 298*ones(nx, ny);           

% Boundary conditions
T(1,1:end) = 900;   % bottom boundary condition
T(end,1:end) = 600; % top boundary condition
T(1:end,1) = 400;   % left boundary condition
T(1:end,end) = 800; % right boundary condition

% creating a copy of T
Told = T;          
T_prev_dt = T;
tol = 1e-4;
error = 9e9;

% for ease of calculation
k1 = alpha*dt/dx^2;        
k2 = alpha*dt/dy^2;
term1 = (1+(2*k1)+(2*k2));
  
prompt= \'Which method to be used(1.Implicit with jacobi  2.Implicit with Gauss sidel 3. Implicit with SOR)?\';
iterative_solver=input(prompt);

% Jacobi iteration using implicit method for unsteady equation
jacobi_iter = 1;
if(iterative_solver == 1)
    tic
   while(error>tol)
       
      for i = 2:(nx-1)
          for j = 2:(ny-1)
              h = Told(i-1,j)+Told(i+1,j);
              v = Told(i,j-1)+Told(i,j+1);
              T(i,j)= ((T_prev_dt(i,j) +h*k1 +v*k2)/term1);
           end
      end
            
     error=max(max(abs(Told-T)));
     jacobi_iter=jacobi_iter+1;
            
     % updating old values
     Told = T;
     % creating a contour plot
     [C,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     clabel(C,d);
     xlabel(\'X axis\');
     ylabel(\'Y axis\');
     title_text=sprintf(\'Jacobi iteration number=%d @ unsteady state\',jacobi_iter);
     title(title_text);       % display the title for plot
     pause(0.0001)
   end
     T_prev_dt = T; 
     jacobi_iteration_time=toc;
 end

% Gauss seidal iteration using implicit method for unsteady equation
GS_iter = 1;
if(iterative_solver == 2)
    tic
   while(error>tol)
      for i = 2:(nx-1)
         for j = 2:(ny-1)
             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) +h*k1 +v*k2)/term1);
          end
      end
            
     error=max(max(abs(Told-T)));
     GS_iter=GS_iter+1;
            
     % updating old values
     Told = T;
     % creating a contour plot
     [C,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     clabel(C,d);
     xlabel(\'X axis\');
     ylabel(\'Y axis\');
     title_text=sprintf(\'Gauss Seidel iteration number=%d @ unsteady state\',GS_iter);
     title(title_text);      % display the title for plot
     pause(0.0001)
   end
     T_prev_dt = T;   
     GS_iteration_time = toc;
end

% Succesive over relaxation with Gauss seidal iteration using implicit method for unsteady equation
SOR_GS_iter = 1;
omega = 1.473;
if(iterative_solver == 3)
    tic
   while(error>tol)
      
       for i = 2:(nx-1)
           for j = 2:(ny-1)
               h = T(i-1,j)+Told(i+1,j);
               v = T(i,j-1)+Told(i,j+1);
               T(i,j)= Told(i,j)*(1-omega)+omega*((T_prev_dt(i,j) +h*k1 +v*k2)/term1);
           end
        end
            
     error=max(max(abs(Told-T)));
     SOR_GS_iter=SOR_GS_iter+1;
            
     % updating old values
     Told = T;
     % creating a contour plot
     [C,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     clabel(C,d);
     xlabel(\'X axis\');
     ylabel(\'Y axis\');
     title_text=sprintf(\'SOR + GS iteration number=%d @ unsteady state\',SOR_GS_iter);
     title(title_text);       % display the title for plot
     pause(0.0001)
   end
     T_prev_dt = T;
     SOR_GS_iteration_time=toc;
end

PLOT OF VARIOUS ITERATION FOR UNSTEADY STATE:

1. JACOBI METHOD

2. GAUSS SEIDAL METHOD

3. SUCCESSIVE OVER RELAXATION USING GAUSS SEIDAL METHOD

RESULTS

JACOBI TIME PERIOD

GAUSS SEIDAL TIME PERIOD

SOR using GAUSS TIME PERIOD

CONCLUSIONS

For UnSteady State Solver, a similar trend gets followed as that of Steady-State solver in terms of parameters like no. of iterations and the total converging time.  About the implicit method, SOR iterative method takes the least iterations and convergence time i.e.39 and 9.4569 seconds respectively. The Jacobian method takes the most time 184 iterations and 33.9841 seconds which way too much computational time when compared with the SOR method. The Gauss-Seidel Method stays in between with 101 iterations and 21.191 seconds. In the above implicit solution plots, we notice that the Jacobi method took the most number of iterations to converge followed by SOR (omega = 1.473) and then Gauss-Seidel. In the SOR method, increasing the relaxation factor only caused a slower convergence rate. However, the number of iterations decreased for omega<1 (under-relaxation).

TRANSIENT STATE ANALYSIS FOR 2D HEAT CONDUCTION USING EXPLICIT METHOD

% Solving 2D transient state heat conduction equation using explicit method

clear all
close all
clc
                                    
nx = 20;                    % number of grid points along x and y
ny = nx;
dt = 1e-4;                  % time step size
x = linspace(0,1,nx);       % x nodes
y = linspace(0,1,ny);       % y nodes
dx = x(2)-x(1);             % gird size along x
dy = y(2)-y(1);             % grid size along y
alpha = 1.4;                % thermal diffusivity 
[xx,yy] = meshgrid(x,y);

%Initialization
T = 298*ones(nx, ny);           

T(1,1:end) = 900;   % bottom boundary condition
T(end,1:end) = 600; % top boundary condition
T(1:end,1) = 400;   % left boundary condition
T(1:end,end) = 800; % right boundary condition

% creating a copy of T
Told = T;          
% for ease of calculation
k1 = alpha*dt/dx^2;        
k2 = alpha*dt/dy^2;
  
% iteration loop
tic
for k = 1:1000  
    % nodal loop
    for i = 2:(nx-1)
        for j = 2:(ny-1)
         % explicit scheme for 2D heat equation
         T(i,j)=(1-2*k1-2*k2)*Told(i,j)+k1*(Told(i+1,j)+Told(i-1,j))+k2*(Told(i,j+1)+Told(i,j-1));    
         end
    end
     % updating old values
     Told = T;
     % creating a contour plot
     [C,d]=contourf(xx,yy,T); % to fill contour plot
     colormap(jet)            % sets the colormap for the figure specified by FIG
     colorbar                 % display colorbar to the current axes in the default
     clabel(C,d);
     xlabel(\'X axis\');
     ylabel(\'Y axis\');
     title_text=sprintf(\'Explicit iteration number=%d @ unsteady state\',k);
     title(title_text);       % display the title for plot
     pause(0.00001)
end
Explicit_time=toc    

PLOT OF THE EXPLICIT METHOD

RESULTS

CONCLUSION

The implicit schemes are slower and more complicated than the explicit scheme, since we\'re using guess values and iterating for convergence in each time step. However, they are unconditionally stable for any choice of dt with respect to dx. They offer better stability and larger time step sizes. The explicit method the number of iteration is fixed by keeping the time loop and error> tol is not used since in explicit method depends on the time loop not the conditions which are set.