2D steady and transient heat conduction solvers with detailed interpretation of results

AIM: To solve 2D heat conduction equation using steady and unsteady solvers with different iterative techniques

ASSUMPTIONS:

1. Assume that the domain is a unit square.

2. Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]

3. Boundary conditions for steady and transient case

1. Left:400K

2. Top:600K

3. Right:800K

4. Bottom:900K

4. Your absolute error criteria is 1e-4

SOLUTION:

1. Program to solve the steady state 2D heat conduction equation by various means, is given below

clear all
close all
clc

L = 1;
nx = ny = 10;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;

% Need to be changed depending upon the required method of solving
% iterative_solver = 1 for Jacobi_iter
% iterative_solver = 2 for Gauss siedel
% iterative_solver = 3 for SOR
iterative_solver = 3 ;


% Defining boundary and initial conditions
% While plotting, the 1st element of the Temperature matrix is assigned at the point x(1),y(1) & similarly other points.
% As a result, the obtained plot is VERTICALLY INVERSE.
% So that, Top and bottom temperatures are initially assigned vice versa, for the purpose of obtaining expected result as the output .
T_left = 400;
T_right = 800;
T_bottom = 600;
T_top = 900;
T_matrix = ones(nx,ny);
T_matrix(nx,2:end-1) = T_bottom;
T_matrix(1,2:end-1) = T_top;
T_matrix(2:end-1,1) = T_left;
T_matrix(2:end-1,ny) = T_right;

% Assigning each corner with the average temperature of the two intersecting edges. 
T_matrix(1,1) = (T_top + T_left)/2;
T_matrix(1,ny) = (T_top + T_right)/2;
T_matrix(nx,1) = (T_left + T_bottom)/2;
T_matrix(nx,ny) = (T_right + T_bottom)/2; 

% Copy of the Temperature matrix, for the purpose of iterative solving
Told = T_matrix;

k_1 = (dx^2)/(2*(dx^2+dy^2));
k_2 = (dy^2)/(2*(dx^2+dy^2));

graphics_toolkit("gnuplot")

   % STEADY STATE, JACOBI
   % 217 ITERATIONS
   
if iterative_solver == 1
    jacobi_iter = 1;
 tic   
 while (error > tol)
  for i = 2:nx-1
   for j = 2:ny-1
   T_matrix(i,j) = k_1*(Told(i,j-1)+Told(i,j+1))+k_2*(Told(i-1,j)+Told(i+1,j));
   end
  end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     jacobi_iter = jacobi_iter + 1; 
 end
 time = toc   
     
     contourf(x, y, T_matrix);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', jacobi_iter);
     title(title_text)
     file_text = sprintf("steady_%d.png",iterative_solver);
     saveas(gca,file_text);
end

    
    
    % STEADY STATE, GS
    % 116 ITERATIONS
    
if iterative_solver == 2
  GS_iter = 1;
  
 tic 
 while (error > tol)
  for i = 2:nx-1
   for j = 2:ny-1
   T_matrix(i,j) = k_1*(T_matrix(i,j-1)+T_matrix(i,j+1))+k_2*(T_matrix(i-1,j)+T_matrix(i+1,j));
   end
  end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     GS_iter = GS_iter + 1;   
 end
  time = toc  
    
     contourf(x, y, T_matrix);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', GS_iter);
     title(title_text) 
     file_text = sprintf("steady_%d.png",iterative_solver);
     saveas(gca,file_text); 
end
    
    
    
     % STEADY STATE SOR
     % 31 ITERATIONS
  
if iterative_solver == 3
   opt_omega = 2/(1+sin(pi*dx));
   SOR_iter = 1;
 tic  
 while (error > tol)
  for i = 2:nx-1
   for j = 2:ny-1
   T_matrix(i,j) =opt_omega*(k_1*(T_matrix(i,j-1)+T_matrix(i,j+1))+k_2*(T_matrix(i-1,j)+T_matrix(i+1,j)))+(1-opt_omega)*Told(i,j);
   end
  end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     SOR_iter = SOR_iter + 1; 
 end
  time = toc   
      
     contourf(x, y, T_matrix);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', SOR_iter);
     title(title_text)
     file_text = sprintf("steady_%d.png",iterative_solver);
     saveas(gca,file_text);
end

2. Program to solve the 2D unsteady state heat conduction equation by various means, is given below

clear all
close all
clc

% Assigning required variables
L = 1;
nx = ny = 10;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;


% Defining boundary and initial conditions.
% While plotting, the 1st element of the Temperature matrix is assigned at the point x(1),y(1) & similarly other points.
% As a result, the obtained plot is VERTICALLY INVERSE.
% So that, Top and bottom temperatures are initially assigned vice versa, for the purpose of obtaining expected result as the output .
T_left = 400;
T_right = 800;
T_bottom = 600;
T_top = 900;
T_matrix = ones(nx,ny);
T_matrix(nx,2:end-1) = T_bottom;
T_matrix(1,2:end-1) = T_top;
T_matrix(2:end-1,1) = T_left;
T_matrix(2:end-1,ny) = T_right;

% Assigning each corner with the average temperature of the two intersecting edges. 
T_matrix(1,1) = (T_top + T_left)/2;
T_matrix(1,ny) = (T_top + T_right)/2;
T_matrix(nx,1) = (T_left + T_bottom)/2;
T_matrix(nx,ny) = (T_right + T_bottom)/2; 

% Copy of the Temperature matrix, for the purpose of iterative solving
Told = T_matrix;

% dt refers to the time step
 dt = 32;
 
% To shift between various methods of solution

% iterative_solver = 1 [Explicit]
% iterative_solver = 2 [Implicit (jacobi inside convergence loop)]
% iterative_solver = 3 [Implicit (Gauss siedel inside convergence loop)]
% iterative_solver = 4 [Implicit (Gauss siedel with successive over relaxation inside convergence loop)]
% iterative_solver = 5 [Implicit Crank Nicolson ]
% iterative_solver = 6 [Implicit Crank Nicolson (Gauss siedel inside convergence loop)]
% iterative_solver = 7 [Implicit Crank Nicolson (Gauss siedel with successive over relaxation inside convergence loop)]
 
 iterative_solver = 7;

% alpha is the thermal diffusivity of the material 
 alpha = 1e-4;
 k_3 = alpha*dt/(dx^2);
 k_4 = alpha*dt/(dy^2);
 
 
 graphics_toolkit("gnuplot")
 

   % UNSTEADY EXPLICIT 
   % For values of dt, greater than 33, the solution tends to become UNSTABLE
if iterative_solver == 1
tic
% Time marching solution 
% With dt = 32, the steady state is obtained at less no.of iterations @ k = 200, with the negligible error of 1.7071e-004.
% Therefore the solution resembles the solution of the steady state solution
 for k = 1:200
  for i = 2:nx-1
   for j = 2:ny-1
  T_matrix(i,j) = Told(i,j)+k_3*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j))+k_4*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
   end
  end
     Told = T_matrix;  
 end
 time = toc
     contourf(x, y, T_matrix);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     file_text = sprintf("unsteady_ex_%d.png",iterative_solver);
     saveas(gca,file_text);
end
    
    
  
  
                 % STABILITY OF THE SOLUTION IS INCREASED IN THE IMPLICIT METHOD
                 % LARGE VALUES OF TIME STEP CAN BE USED    
    
     % UNSTEADY IMPLICIT JACOBI
     % For dt = 32, 2107 ITERATIONS
     
     term1 = (1+2*k_3+2*k_4)^-1;
     term2 = k_3*term1;
     term3 = k_4*term1;
    
if iterative_solver == 2
tic

% Copy of the Temperature matrix for the purpose of equation solving
  T_prev_dt = Told;
  jacobi_iter = 1;

% Time loop  
for k = 1:200  
 % Convergence loop
  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_matrix(i,j) = (T_prev_dt(i,j)*term1)+(term2*H)+(term3*V);
    end
   end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     jacobi_iter = jacobi_iter + 1; 
 end
    
    T_prev_dt = T_matrix;
    error = 9e9;
end
time = toc
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', jacobi_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end
    
    
    
   
    
     % UNSTEADY IMPLICIT GAUSS SIEDEL
     % For dt = 32, 1518 ITERATIONS
     
     term1 = (1+2*k_3+2*k_4)^-1;
     term2 = k_3*term1;
     term3 = k_4*term1;
if iterative_solver == 3
tic
% Copy of the Temperature matrix for the purpose of equation solving
  T_prev_dt = Told;
  GS_iter = 1;
  
% Time loop
 for k = 1:200
 % Convergence loop
  while (error > tol)
   for i = 2:nx-1
    for j = 2:ny-1
   H = (T_matrix(i-1,j)+Told(i+1,j));
   V = (T_matrix(i,j-1)+Told(i,j+1));
  T_matrix(i,j) = (T_prev_dt(i,j)*term1)+(term2*H)+(term3*V);
    end
   end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     GS_iter = GS_iter + 1; 
  end
   
    T_prev_dt = T_matrix;
    error = 9e9;
    
 end
time = toc 
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', GS_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end
    
    
    
    
    
     % UNSTEADY IMPLICIT GAUSS SIEDEL WITH OVER RELAXATION
     % For dt = 32, 1141 ITERATIONS
     
     term1 = (1+2*k_3+2*k_4)^-1;
     term2 = k_3*term1;
     term3 = k_4*term1;
     
if iterative_solver == 4
tic
  opt_omega = 2/(1+sin(pi*dx));
  
% Copy of the Temperature matrix for the purpose of equation solving
  T_prev_dt = Told;
  GS_iter = 1;
  
% Time loop
 for k = 1:200
 % Convergence loop
  while (error > tol)
   for i = 2:nx-1
    for j = 2:ny-1
   H = (T_matrix(i-1,j)+Told(i+1,j));
   V = (T_matrix(i,j-1)+Told(i,j+1));
  T_matrix(i,j) = opt_omega*((T_prev_dt(i,j)*term1)+(term2*H)+(term3*V))+(1-opt_omega)*T_prev_dt(i,j);
    end
   end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     GS_iter = GS_iter + 1; 
  end
   
    T_prev_dt = T_matrix;
    error = 9e9;
    
 end 
 time = toc
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', GS_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end
   
     % UNSTEADY IMPLICIT CRANK NICOLSON
     % dx = dy, Therefore k_3 = K_4
     % ORDER OF ACCURACY is increased in Crank Nicolson method
     % For dt = 32, 1489 ITERATIONS
     term1 = (1-2*k_3)/(1+2*k_3);
     term2 = k_3/(2*(1+2*k_3));
    
  if iterative_solver == 5
tic  
  % Copy of the Temperature matrix for the purpose of equation solving
  T_prev_dt = Told;
  CN_iter = 1;

% Time loop
for k = 1:200

 % Convergence loop
 while (error > tol)
  for i = 2:nx-1
   for j = 2:ny-1
   H = (T_prev_dt(i+1,j)+T_prev_dt(i-1,j)+T_prev_dt(i,j+1)+T_prev_dt(i,j-1));
   V = (Told(i+1,j)+Told(i-1,j)+Told(i,j+1)+Told(i,j-1));
  T_matrix(i,j) = (T_prev_dt(i,j)*term1)+(term2*H)+(term2*V);
   end
  end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     CN_iter = CN_iter + 1; 
 end
   
     T_prev_dt = T_matrix;
     error = 9e9;
end 
 time = toc  
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', CN_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end
    
    
     
    
    
    % UNSTEADY IMPLICIT CRANK NICOLSON  [GS(Immediate update of the terms inside the convergence loop)]
    % By implementing the same idea of GAUSS SIEDEL technique inside the convergence loop, the required no.of iterations get reduced, thus rate of convergence is increased
    % For dt = 32, 1153 ITERATIONS
  
  % k_3 = k_4
     term1 = (1-2*k_3)/(1+2*k_3);
     term2 = k_3/(2*(1+2*k_3));
     
if iterative_solver == 6
tic
  % Copy of the Temperature matrix for the purpose of equation solving
  T_prev_dt = Told;
  CN1_iter = 1;
  
% Time loop  
 for k = 1:200

% Convergence Loop
  while (error > tol)
   for i = 2:nx-1
    for j = 2:ny-1
   H = (T_prev_dt(i+1,j)+T_prev_dt(i-1,j)+T_prev_dt(i,j+1)+T_prev_dt(i,j-1));
   V = (Told(i+1,j)+T_matrix(i-1,j)+Told(i,j+1)+T_matrix(i,j-1));
   T_matrix(i,j) = (T_prev_dt(i,j)*term1)+(term2*H)+(term2*V);
    end
   end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     CN1_iter = CN1_iter + 1; 
  end

    T_prev_dt = T_matrix;
    error = 9e9; 
 end 
 time = toc
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', CN1_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end
   
  
 



    % UNSTEADY IMPLICIT CRANK NICOLSON WITH OVER RELAXATION [GS(Immediate update of the terms inside the convergence loop)]
    % By implementing the same idea of GAUSS SIEDEL technique inside the convergence loop, the required no.of iterations get reduced, thus increases the rate of convergence
    % Over relaxation factor, further increases the rate of convergence
    % For dt = 32, 918 ITERATIONS
    
    
    % k_3 = k_4 
     term1 = (1-2*k_3)/(1+2*k_3);
     term2 = k_3/(2*(1+2*k_3));
     
if iterative_solver == 7
tic
  % Copy of the Temperature matrix for the purpose of equation solving
  T_prev_dt = Told;
  CN1_iter = 1;
  
% Time loop  
 for k = 1:200
opt_omega = 2/(1+sin(pi*dx));
% Convergence Loop
  while (error > tol)
   for i = 2:nx-1
    for j = 2:ny-1
   H = (T_prev_dt(i+1,j)+T_prev_dt(i-1,j)+T_prev_dt(i,j+1)+T_prev_dt(i,j-1));
   V = (Told(i+1,j)+T_matrix(i-1,j)+Told(i,j+1)+T_matrix(i,j-1));
   T_matrix(i,j) = opt_omega*((T_prev_dt(i,j)*term1)+(term2*H)+(term2*V))+(1-opt_omega)*T_prev_dt(i,j);
    end
   end
     error = max(max(abs(Told-T_matrix)));
     Told = T_matrix; 
     CN1_iter = CN1_iter + 1; 
  end

    T_prev_dt = T_matrix;
    error = 9e9; 
 end 
 time = toc
     contourf(x, y, Told);
     [c,h] = contourf(x, y, T_matrix);
     clabel(c,h);
     colorbar
     title_text = sprintf('iteration number = %d', CN1_iter);
     title(title_text)
     file_text = sprintf("unsteady_im_%d.png",iterative_solver);
     saveas(gca,file_text);
end 

STEADY STATE EQUATION - It is solved using iterative techniques(jacobi, GS, SOR) with the condition that the error criteria is within the tolerance limit. SOR method reaches the solution with less no.of iterations and also less simulation time.

With 10 grid points on both axes

A. jacobi - 217 iterations with 0.513s simulation time

B. GS - 116 iterations with 0.285s simulation time

C. SOR - 31 iterations with 0.101s simulation time

UNSTEADY STATE EQUATION:

1. EXPLICIT - It is solved using the concept of time marching, and doesn\'t require formation of any convergence loop.

1. Large values for time step can\'t be used and it results in blowing up of the solution.
2. Stability is based on certain criteria and it depends upon the type of equation involved. For 2D heat conduction equation with equal grid spacing on both the axes, the stability criteria is r = alpha*(dt/dx^2) & r <= 1/4
3. With dt = 32, the steady state is obtained at less no.of iterations @ k = 200, with the negligible error of 1.7071e-004. If dt is chosen to be greater than 32 then the solution starts to blow up.
4. We have chosen dt as 32 and the no.of time loops as 200, therefore the solution resembles the solution of the steady state solution, due to the reason that the steady state has been attained.
5. In other words, it could be stated that it takes approximately 6400s (dt*k = 32*200) to reach steady state, for the assumed set of boundary and initial conditions & the parameters chosen (thermal diffusivity).
6. Solving an explicit equation requires less time, than solving the implicit equation. The simulation time taken to solve our given equation is 0.7465s.

2. IMPLICIT - Influence of both time loop and convergence loop is experienced. The rate of convergence is modified by using the above mentioned three iterative methods(Jacobi, GS, SOR).

1. The solution remains stable even for large values of dt
2. r = alpha*(dt/dx^2) , here the value of 'r' doesn't affect the stability of the solution.
3. Time accuracy decreases as the time step increases, in the case of solving the implicit form of the unsteady equation
4. Solving implicit equations iteratively involves higher no.of iterations & higher simulation time than the explicit one.
5. Comparing no.of iterations taken by different methods of solving, by increasing the rate of convergence. The value of dt (time step) is taken to be 32, no.of loops of time as 200 and no.of grid points on both the axes as 10, for all the methods involved so that the effect on the no.of iterations can be observed specifically.

A. [Implicit (Jacobi inside convergence loop)] 2107 ITERATIONS with 6.2424s simulation time

B. [Implicit (Gauss Siedel inside convergence loop)] 1518 ITERATIONS with 4.666s simulation time

C. [Implicit (Gauss Siedel with successive overrelaxation inside convergence loop)] 1141 ITERATIONS with 4.188s simulation time

D. [Implicit Crank Nicolson ] 1489 ITERATIONS with 7.1441s simulation time

E. [Implicit Crank Nicolson (Gauss Siedel inside convergence loop)] 1153 ITERATIONS with 5.4719s simulation time

F. [Implicit Crank Nicolson (Gauss Siedel with successive overrelaxation inside convergence loop)] 918 ITERATIONS with 4.97s simulation time

f. For the methods involving successive over relaxation, that is C and F, the solution becomes unstable for large values of dt.

STEADY VS UNSTEADY STATE ANALYSIS

Steady state

1. Fastest steady-state solution is obtained by the SOR iterative method with less no.of iterations and less simulation time. It took 31 iterations with 0.101s simulation time to solve 2D heat conduction equation with the given initial & boundary conditions.
2. The optimum value of the relaxation factor increases the rate of convergence, thus leading to less no.of iterations.

Unsteady state

1. Fastest transient simulation is obtained by solving the explicit form of the unsteady equation. It took 0.7465 seconds.
2. The explicit form of unsteady equation can be solved by the approach of time marching and there is no need of solving a set of equations, which makes the calculation part easier and less time-consuming.
3. Choosing the value of time step is very important because large values of dt results in solution blow up. The stability of the solution of the 2D heat conduction equation is dependent on the value of the variable 'r' which is equal to, r = alpha*(dt/dx^2).
4. The value of 'r' should be less than 0.25. The following figure is the blownup explicit solution with r value equal to 0.40500
5. So whenever we are solving the explicit form, it is safe to have a minimum value of timestep to avoid blowup of the solution and to maintain a good range of accuracy
6. Whereas the implicit form of unsteady equations require high simulation time but are stable even for large values of dt, and the value of 'r' won't affect the stability of the solution.
7. This is because implicit schemes use the entire domain to calculate each timestep. The drawback of implicit form is, even though it is stable for large values of dt, the time accuracy decreases as the value of time step increases.

The reason why both the steady and transient plots look similar - While solving the unsteady equation, we had chosen dt = 32 and no.of time loops = 200. So that the final time at which solution gets calculated will be at 6400th second at which the steady state is already obtained with the negligible error of 1.7071e-004. And that's why the solution of the unsteady equation resembles the solution of the steady state equation.