Study and Analysis of 2D Heat Conduction Equation in Steady and Transient State

1. INTRODUCTION:

The aim of this project is to solve and study the 2D heat conduction equation in the Steady and Transient state. The equation is given by: `frac{delT}{delt} = alpha[frac{del^2T}{delx^2}+ frac{del^2T}{dely^2}]` where `alpha` is the thermal diffusivity.By solving this equation, a heat simulation of the same can be obtained.

In this project, two approaches are followed.

1. Steady state which is independant of time: EQN: `frac{del^2T}{delx^2}+ frac{del^2T}{dely^2} = 0` i.e., `nabla^2(T) = 0` (Laplacian T = 0)
2. Transient state which is dependant on time: EQN: `frac{delT}{delt} = alpha[frac{del^2T}{delx^2}+ frac{del^2T}{dely^2}]`

The steady state is solved using Iterative methods such as Jacobi, Gauss-Siedel and Successive Over-relaxation (SOR) methods and a comparison is done between the results.

The transient state is solved explicitly as well as implicitly and a comparison is done between the results. The implicit method uses 3 iterative solvers as well- Jacobi, Gauss Siedel and SOR.

2. EXPECTED RESULTS:

1. In both analysis, the SOR solver must yield the result in least number of iterations followed by Gauss Siedel followed by Jacobi.
2. Transient state will reach steady state if we time-march long enough.
3. As 'dt' becomes smaller in magnitude, the time for simulation should increase.

3. THEORY:

The partial differential heat equation is initially discretized.

For the Steady State, the equation is an elliptical PDE. Hence we use Central Difference Scheme (CDS) to take information from all sides of the current point. The discretized equation is shown below:

`frac{T(i-1,j) + T(i+1,j) - 2*T(i,j)}{trianglex^2} + frac{T(i,j-1) + T(i,j+1) - 2*T(i,j)}{triangley^2} = 0`.

Simplifying the equation further, we get:

`T(i,j) = (1/K)*([frac{T(i-1,j) + T(i+1,j)}{trianglex^2}] + [frac{T(i,j-1) + T(i,j+1)}{triangley^2}])`

where `K = 2*frac{trianglex^2 + triangley^2}{trianglex^2triangley^2}`

The above equation is used for Jacobi and Gauss Siedel methods. SOR involves taking the result from either Jacobi or Gauss Siedel and make it converge rapidly using the over-relaxation factor (`alpha`) >1. If `alpha`< 1, the convergence will take more number of iterations than the original method used. The SOR formula is given below:

`overset(k+1)(T)underset(n)() = overset(k)(T)underset(n)() + alpha*frac{dT}{di}ni`

where `frac{dT}{di}` is the rate of change of T vector w.r.t. iterations (i), ni is the number of iterations. The `frac{dT}{di}` term is taken from Jacobi or Gauss Siedel (Tgs). The number of iterations is 1 in this project as we calculate the SOR term inside every loop. Hence the equation reduces to `Toverset(k+1)()underset(n)() = overset(k)(T)underset(n)() + alpha*(Tunderset( gs)()-Toverset(k)())`.

For Transient State, a similar discretization is done but with the time derivative term. The discretized equation is:

`frac{overset(n+1)(T)underset(i,j)() -overset(n)(T)underset(i,j)() }{trianglet}= alpha*([frac{T(i-1,j) + T(i+1,j)-2*T(i,j)}{trianglex^2}] + [frac{T(i,j-1) + T(i,j+1)-2*T(i,j)}{triangley^2}])^n`

In Transient Explicit, the method uses already known values of previous time loop and hence the arriving at the new value is easy and less time consuming. In transient Implicit, the method uses guess values initially and converges at the new value taking information from the same time step. The formula for Explicit and Implicit methods are listed below:

EXPLICIT:

`Toverset(n+1)()underset(i,j)() = Toverset(n)()underset(i,j)()*(1-2kunderset(1)()-2kunderset(2)()) + k1*(T(i-1,j)+T(i+1,j))^n+k2*(T(i,j-1)+T(i,j+1))^n`

IMPLICIT:

`Toverset(n+1)()underset(i,j)() = frac{Toverset(n)()underset(i,j)()+ k1*(T(i-1,j)+T(i+1,j))^(n+1)+k2*(T(i,j-1)+T(i,j+1))^(n+1)}{1-2kunderset(1)()-2kunderset(2)()}`

where K1 = `alphafrac{trianglet}{trianglex^2}`

and K2 = `alphafrac{trianglet}{triangley^2}`

4. CFL CONDITION:

In transinet schemes, the Courant Friedrich Lewy (CFL) condition comes into play which describes the relation between the step size and time step in regard to the stability of the solution. According the CFL stability equation for 2D,

`alphafrac{trianglet}{trianglex^2}<=Cmax`

Where `alpha` = Thermal diffusivity and Cmax is the Courant Number which is equal to 0.5 in diffusion processes.

dt is directly proportional to Courant's Number, diffusivity and inversely proportional to square of step size.

Hence the 'dt' calculated must adhere to C< 0.5

5. INITITAL CONDITIONS/ASSUMPTIONS:

Square Domain, nx = ny;

Boundary Consitions:

Left:400K

Top:600K

Right:800K

Bottom:900K

Length of Domain = 1m;

Thermal diffusivity,`alpha` = 1;

Tolerance (For iteration) = 1e-4;

6. MATLAB/OCTAVE CODE:

• Steady state main program

clear all
close all
clc

%steady state eqn: d^2x/dx^2 + d^2y/dy^2 = 0;
%square domain nx = ny;

nx = ny = 10;
tol = 1e-4;

T = zeros(nx,ny);              %Initial conditions;

T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 600;
T(end,:) = 900;

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


ip = input('Enter the method: 1.JACOBI 2.GS 3.SOR 4.COMPARE ALL ');

if(ip==1)
tic;
[T1,jacobi_iter] = jacobi_solver(T,nx,ny)   ;                 % Steady state
sim_time = toc;
elseif(ip==2)
tic;
[T2,GS_iter] = GS_solver(T,nx,ny);
sim_time = toc;
elseif(ip==3)
tic;
[T2,GS_iter] = GS_solver(T,nx,ny);
[T3,SOR_iter] = SOR_GS_solver(T,T2,nx,ny);
sim_time = toc;

elseif(ip==4)
[T1,jacobi_iter,jsim_time] = jacobi_solver(T,nx,ny);    
[T2,GS_iter,gsim_time] = GS_solver(T,nx,ny);
[T3,SOR_iter,ssim_time] = SOR_GS_solver(T,T2,nx,ny);
Y = [jacobi_iter; GS_iter; SOR_iter];  
t = [jsim_time, gsim_time, ssim_time];
figure(1);
subplot(2,1,1);
bar(Y,'g');
xlabel('Methods: 1.Jacobi 2.GS 3.SOR');
ylabel('Number of Iterations');
title('Comparison of no of iterations in different Iterative Methods','fontsize',15);
subplot(2,1,2);
bar(t);
xlabel('Methods: 1.Jacobi 2.GS 3.SOR');
ylabel('Simulation Time ');
title('Comparison of simulation time in different Iterative Methods','fontsize',15);
endif

• Steady state Functions - Jacobi

function [T,jacobi_iter,sim_time] = jacobi_solver(T,nx,ny)
  tic;
  tol = 1e-4;
  L = 1 ;                                   % l  = 1M assuming length of domain;
  error = 9e9;
  dx = dy = L/(nx-1);                            %dx = L/(n-1);

  k = 2*(dx^2+dy^2)/(dx^2*dy^2);
  jacobi_iter = 1;
  Told = T;
  
  while(error>tol)
      for i = 2 : nx-1
        for j = 2 : ny-1
        
        T(i,j) = (  ((Told(i-1,j) + Told(i+1,j))/(dx^2)) + ((Told(i,j-1) + Told(i,j+1))/(dy^2))  ) / k;
        
        end
      end
    error = max(abs(Told(i,j).- T(i,j)));
    jacobi_iter = jacobi_iter+1;
    
    Told = T;
  end
  
[M] = contourf(T);
colorbar;
xlabel('X Axis');
ylabel('Y Axis');
title_text = sprintf('Heat chart at %dth Node Iteration of Jacobi',jacobi_iter);
title(title_text,'fontsize',12);

sim_time = toc;  
endfunction

• Steady state Functions - Gauss Siedel

function [T,GS_iter,sim_time] = GS_solver(T,nx,ny)
  tic;
  L = 1 ;     
  tol = 1e-4;
  error = 9e9;
 dx = dy = L/(nx-1); 
 
  k = 2*(dx^2+dy^2)/(dx^2*dy^2);
  GS_iter = 1;
  Told = T;
  
  while(error>tol)
      for i = 2 : nx-1
        for j = 2 : ny-1
        
        T(i,j) = (  ((T(i-1,j) + Told(i+1,j))/(dx^2)) + ((T(i,j-1) + Told(i,j+1))/(dy^2))  ) / k;
        
        end
      end
    error = max(abs(Told(i,j).- T(i,j)));
    GS_iter = GS_iter+1;
    
    Told = T;
  end
  
  [M] = contourf(T);
  colorbar;
xlabel('X Axis');
ylabel('Y Axis');
title_text = sprintf('Heat chart at %dth Node Iteration of Gauss Siedel',GS_iter);
title(title_text,'fontsize',12);

sim_time = toc;
endfunction

• Steady state Functions - SOR

function [result,SOR_iter,sim_time] = SOR_GS_solver(T,T2,nx,ny)
  tic;
  L = 1 ;   
  tol = 1e-4;
  error = 9e9;
  dx = dy = L/(nx-1); 
  SOR_const = 1.5;         % Scaler1
 
  SOR_iter = 1;
  Told = T;
  
result = zeros(nx,ny);              
result(:,1) = 400;
result(:,end) = 800;
result(1,:) = 600;
result(end,:) = 900;
result(1,1) = (400+600)/2;
result(end,1) = (900+400)/2;
result(1,end) = (800+600)/2;
result(end,end) = (800+900)/2;
  
  while(error>tol)
      for i = 2 : nx-1
        for j = 2 : ny-1
        
        result(i,j) = Told(i,j) +   SOR_const*( T2(i,j) - Told(i,j))   ;
        
        end
      end
    error = max(abs(Told(i,j).- result(i,j)));
    SOR_iter = SOR_iter+1;
    
    Told = result;
  end
  [M] = contourf(result);
  colorbar;
xlabel('X Axis');
ylabel('Y Axis');
title_text = sprintf('Heat chart at %dth Node Iteration of SOR',SOR_iter);
title(title_text,'fontsize',12);

sim_time = toc;
endfunction

• Transient state-Main program

clear all
close all
clc

%steady state eqn: d^2x/dx^2 + d^2y/dy^2 = 0;
%square domain nx = ny;

nx = ny = 10;
tol = 1e-4;

T = zeros(nx,ny);              %Initial conditions;

T(:,1) = 400;
T(:,end) = 800;
T(1,:) = 600;
T(end,:) = 900;

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


ip = input('Enter the transient method: 1.Explicit 2.Implicit Jacobi 3.Implicit GS 4.Implicit SOR 5.Compare Implicit ');

if (ip==1)
tic;
[T4,iter] = tran_solver(T);                                  % Transient Explicit
sim_time = toc;

elseif(ip==2||ip==3||ip==4)
tic;
[T5,node_iter,time_iter] = tran_imp(T,nx,ny,ip) ;            % Transient Implicit 
sim_time = toc;

elseif(ip==5)

tic;
[T5,jnode_iter,jtime_iter] = tran_imp(T,nx,ny,2);
jsim_time = toc;
tic;
[T5,gnode_iter,gtime_iter] = tran_imp(T,nx,ny,3);
gsim_time = toc;
tic;
[T5,snode_iter,stime_iter] = tran_imp(T,nx,ny,4);
ssim_time = toc;

t = [jsim_time; gsim_time; ssim_time];  
N = [jnode_iter, gnode_iter, snode_iter];

figure(1);
subplot(2,1,1);
bar(N,'g');
xlabel('Methods: 1.Jacobi 2.GS 3.SOR');
ylabel('Number of Iterations');
title('Comparison of no of iterations in Implicit Methods','fontsize',15);
subplot(2,1,2);
bar(t);
xlabel('Methods: 1.Jacobi 2.GS 3.SOR');
ylabel('Simulation Time ');
title('Comparison of simulation time in Implicit Methods','fontsize',15);
endif

• Functions - Explicit

function [T,node_iter] = tran_solver(T)
  
   nx = ny = 10;
   L = 1;
   tot_time = 0.3;
   dx = dy = L/(nx-1);
   %dt = CFL_check(dx,dy);
   dt = 3.3e-3;
   
   nt =  tot_time/dt;
  
  alfa = 1;           % Thermal diffusivity
  Told = T;
  node_iter = 1;
  
   for t = 1:nt        %defining time loop
     for i = 2:nx-1
      for j = 2:ny-1          %space loop
       T(i,j) = Told(i,j) + dt*alfa*((Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1)-4*Told(i,j))/(dx^2));
      end
     end
     Told = T;
      node_iter = node_iter+1;
   end
[M] = contourf(T);
colorbar;
xlabel('X Axis');
ylabel('Y Axis');
title_text = sprintf('Heat chart at %dth Node Iteration at Steady-State and dt = %d',node_iter,dt);
title(title_text,'fontsize',12);

endfunction

• Functions - Implicit

  function [T,node_iter,time_iter] = tran_imp(T,nx,ny,ip)
     L = 1;                      %DOmain
     tot_time = 0.3;
     dx = dy = L/(nx-1);
     dt = CFL_check(dx,dy);      %Calculating dt according to CFL
     
     alfa = 1;                  %Heat conduction coefficient
     SOR_const = 1.5;            % Over-relaxation factor
     nt = tot_time/dt;
     
     tol = 1e-4;
     error = 9e9;
     node_iter = 1;
     time_iter = 1;
     Told = T;
     Tprev = T; 
     k1 = (alfa*dt)/dx^2;
     k2 = (alfa*dt)/dy^2;
     Tgs = T;
     
  if(ip ==2)                                 %Jacobi
  for t=1:nt
  error = 9e9;
  while(error>tol) 
   
     for i = 2:nx-1
       for j = 2:ny-1
         T(i,j) = ( Tprev(i,j) + k1*(Told(i+1,j) + Told(i-1,j)) + k2*(Told(i,j+1) + Told(i,j-1)) )/ (1+2*k1+2*k2)   ;
       endfor
     endfor
      error = max(abs(Told(i,j).- T(i,j)));
    Told = T;
    node_iter = node_iter +1;
    end 
  Tprev = T;
  time_iter = time_iter + 1;
  
  endfor
  [M] = contourf(T);
  colorbar;
  xlabel('X Axis');
  ylabel('Y Axis');
  title_text = sprintf('Heat chart at %dth Node Iteration and %d Time Iteration of Jacobi at dt = %d',node_iter,time_iter,dt);
  title(title_text,'fontsize',12);
  
  
  elseif (ip==3)                                                %GS
  [T,node_iter,time_iter] = Imp_GS(T,nx,ny);
  [M] = contourf(T);
  colorbar;
  xlabel('X Axis');
  ylabel('Y Axis');
  title_text = sprintf('Heat chart at %dth Node Itr. and %d Time Itr. of Gauss-Siedel at dt = %d',node_iter,time_iter,dt);
  title(title_text,'fontsize',12);
  
  
  elseif (ip==4)                                                %SOR
  
  
  for t=1:nt
    error = 9e9;
    [Tgs,node_iter1,time_iter] = Imp_GS(T,nx,ny);
    
    while(error>tol) 
   
     for i = 2:nx-1
       for j = 2:ny-1
              T(i,j) =  Tprev(i,j) + SOR_const*(Tgs(i,j) - Tprev(i,j))  ;
        endfor
     endfor
     error = max(abs(Told(i,j).- T(i,j)));
    Told = T;
    node_iter = node_iter +1;
    end 
  
  [M] = contourf(T);
  colorbar;
  xlabel('X Axis');
  ylabel('Y Axis');
  title_text = sprintf('Heat chart at %dth Node Itr. and %d Time Itr. of SOR at dt = %d',node_iter,time_iter,dt);
  title(title_text,'fontsize',12);
  Tprev = T;
  time_iter = time_iter + 1;
  endfor
  
  endif

• Function - GS

function [T,GS_iter,sim_time] = GS_solver(T,nx,ny)
  tic;
  L = 1 ;     
  tol = 1e-4;
  error = 9e9;
 dx = dy = L/(nx-1); 
 
  k = 2*(dx^2+dy^2)/(dx^2*dy^2);
  GS_iter = 1;
  Told = T;
  
  while(error>tol)
      for i = 2 : nx-1
        for j = 2 : ny-1
        
        T(i,j) = (  ((T(i-1,j) + Told(i+1,j))/(dx^2)) + ((T(i,j-1) + Told(i,j+1))/(dy^2))  ) / k;
        
        end
      end
    error = max(abs(Told(i,j).- T(i,j)));
    GS_iter = GS_iter+1;
    
    Told = T;
  end
  
  [M] = contourf(T);
  colorbar;
xlabel('X Axis');
ylabel('Y Axis');
title_text = sprintf('Heat chart at %dth Node Iteration of Gauss Siedel',GS_iter);
title(title_text,'fontsize',12);

sim_time = toc;
endfunction

• Function - CFL condition

function [dt] = CFL_check(dx,dy)
% FOr 2D diffusion, dt <= dx^2/2*alfa;
  
  alfa = 1;               %diffusion coefficient
  dt = (dx^2)/(2*alfa);
  
  endfunction

7. OUTPUT/RESULTS:

7A) STEADY STATE:

7B) INFERENCE FROM STEADY STATE:

1. Contour Plots are plotted for each Iterative method at steady state. All the methods converge at steady state and gives a similar plots.
2. The plots mainly defines the efficieny of the solver used. The SOR solver is the most efficient as it is able to get the output at 26 Iterations, followed by Gauss Siedel solver which gets the result at 96 Iterations. Jacobi solver is the least efficient as results are obtained at 183 Iterations.

A comparison plot of the 3 methods is shown below:

3. The comparison for number of iterations as well as the simulation time has been plotted for all three methods. From the plots, it can be said that SOR takes the least time for simulation as well as converges with the shortest span of Iterations and hence it is the most efficient solver followed by Gauss siedel followed by Jacobi method.

7B) TRANSIENT STATE:

7B1) TRANSIENT EXPLICIT

7B2) INFERENCE FROM TRANSIENT EXPLICIT:

1. Transient state will attain steady state if time marching is done long enough (by increasing time domain). Similar plots of steady state is obtained in transient state after time marching long enough.
2. In explicit methods, the current values at the points are calculated from the values of the previous iteration.
3. Transient Explicit is calculated in Time Domain = 0.3 seconds, at two values of dt (timestep) which is 3.3e-3 and at max CFL condition (Cmax = 0.5) i.e., at dt = 6.17e-3. The solver converges at 91st Iteration when dt = 3.3e-3 and total simulation time obtained is 1.344 seconds.
4. The solution becomes unstable exactly at 3.5e-3 which is before the CFL condition (Cmax = 0.5). Even though the CFL condition is satisfied, there exists a particular point after which oscillations start and solution loses stability.
5. The solution is unstable at dt (6.17e-3) which is calculated from CFL condition. CFL condition is a necessary step to be followed to attain stable solution. The value of dt is dependant on square of dx, directly proportional to Courant's No. and inversely proportional to diffusivity term.
6. As value of dt becomes smaller, the simulation time increases.

Unstable soln dt=1e-4

Unstable sol dt=6.1e-3

7C) TRANSIENT IMPLICIT:

7C1) Jacobi Method: Heat chart plotted at dt = 1e-3, 6.17e-3 and 1e-4

7C2) Gauss Siedel: Heat chart plotted at dt = 1e-3, 6.17e-3 and 1e-4

7C3) SOR METHOD: Heat chart plotted at dt = 1e-3, 6.17e-3 and 1e-4

7C4) INFERENCE FROM TRANSIENT IMPLICIT METHODS:

1. In Implicit methods, the current values at the points are calculated from an initially assumed value and forming equations. Forming equations for an engineering problem may not be feasible because of large number of nodes and hence many equations which cannot be solved simultaneously. Hence Iteration methods are resorted to.
2. Transient state will attain steady state if time marching is done long enough. Plots similar to steady state are obtained.
3. Plots attached are obtained at dt = 1e-3, 1e-4 and at max Courant's No ie 6.17e-3
4. The main observation is that even at CFL boundary (Cmax = 0.5) and past it (dt = 1e-4) the solution does not go unstable and we get tangible results. CFL does not affect Implicit method.
5. To compare the efficiency of the solvers, consider the plots at the same value of dt. At the Cmax = 0.5 i.e., dt = 6.17e-3; At this value of dt, the Jacobi method attains covergence at 913 Iterations and 49 time steps for a Time Domain of 0.3 seconds; Whereas the Gauss Siedel method attains covergence at 498 Iterations and 98 time steps for a Time Domain of 0.6 seconds. The SOR method attains convergence at least number of Node Iterations i.e., 73 and time steps of 99. Hence SOR solver is able to converge with minimum number of iteraions. A comparison plot of the three solvers is shown below:

o

8) CONCLUSION:

1) SOR is the most efficient solver in steady state.

2) The transient explicit method is faster compared to implicit since it calculates current values from previous iterations. But it is bound by CFL condition and cannot take larger values of dt.

3) In transient Implicit, even though SOR converges in the most least number of iterations comparitively, the simulation time taken for SOR is the largest (~58 seconds) followed by Jacobi method (~5.8 seconds) and followed by Gauss Siedel (3.21 seconds); But Implicit method is not bound by CFL condition and is unconditionally stable.