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 using both implicit and explicit techniques.

OBJECTIVES:-

• To describe the 2D heat conduction equation and the particular case mentioned in the question.
• To write the code to solve the equation and give a detailed description of it.

THEORY:-

As the name suggests the 2D heat conduction equation describes the flow of heat in 2 directions on a plane, particularly for this case we consider the heat flow through a solid slab. Now, given in the question that the left side has a temperature of 400 degrees celsius, 800 on the right and, 600 and 900 at the top and bottom respectively.

This is the 2D heat conduction equation:-

`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`

We are required to solve the equation both implicitly and explicitly.

• IMPLICIT METHOD:- This method uses iterative solvers like Jacobi, Gauss-Siedel, and successive relaxation to find the consecutive iterations, to gradually minimize the error and eventually solve the equation.
• EXPLICIT METHOD:- In this, we don't use any iterative solvers. The equation is discretized and the unknown is found by progressing in the space and time steps, which again gradually but a bit slowly minimizes the error to solve the equation.

Steady-state equation:- This involves only the derivatives with respect to space coordinates.

`alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`

Explanation of the code:-

First, we give an initial condition of a random matrix with all temperatures as 450. Then, we declare the boundary conditions. Since the temperatures at the corner points will be overridden. Therefore, we take the average of the boundary temperatures for the respective corner points.

The tolerance as given in the question is .0001, so we have to make errors less than the tolerance in order to find the solution. We use a counter to find the number of iterations run by each method, to find the best among them.

For each method, we run two loops to cover every element of the temperature matrix.

• Method 1 is Jacobi:- In this, we first discretize the equation, then the entire old matrix is updated to the iterated matrix, and error is calculated. If the error is more than tolerance, the updated matrix again goes for another iteration. This goes on until the error is less than the tolerance.
• Method 2 is Gauss-Siedel:- This method is the same as the Jacobi, the only difference is that every latest element that has been iterated is updated in the old matrix before going to the next iteration. That is why instead of 'Told' in the discretization we use 'T'.
• Method 3 is successive overrelaxation:- Here we use a factor `omega` in order to overrelax the iteration. This reduces the number of iterations and increases the efficiency of the method.

Here is the code for this:-

clear all
close all
clc

nx = 10;
ny = 10;
T = 350*ones(nx,ny);

T(2:nx-1,1) = 400; %left face
T(2:nx-1,10) = 800; %right face
T(1,2:ny-1) = 600; %top face
T(10,2:ny-1) = 900; %bottom face

T(1,1) = 500;
T(1,ny) = 700;
T(nx,ny) = 850;
T(nx,1) = 650;

Told = T;
L = 1;
dx = L/(nx-1);
dy = dx;
error = 1;
tol = .0001;
ct = 1;

%Define method
method = 3

if method == 1    
while error > tol
  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
  error = max(max(abs(T-Told)));
  Told = T;
  ct = ct+1;
end
contourf(flipud(T),'ShowText','ON')
title(sprintf('JACOBI iterations = %d',ct));
colorbar
end


if method == 2    
while error > tol
  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
  error = max(max(abs(T-Told)));
  Told = T;
  ct = ct+1;
end
contourf(flipud(T),'ShowText','ON')
title(sprintf('GS iterations = %d',ct));
colorbar
end


if method == 3
    omega = 1.5;    
while error > tol
  for i = 2:nx-1
      for j = 2:ny-1
          TGS = 0.25*((T(i-1,j)+Told(i+1,j)) + (T(i,j-1)+Told(i,j+1)));
          T(i,j) = (1-omega)*Told(i,j) + omega*TGS;
      end
  end
  error = max(max(abs(T-Told)));
  Told = T;
  ct = ct+1;
end
contourf(flipud(T),'ShowText','ON')
title(sprintf('SOR iterations = %d',ct));
colorbar
end

Transient-state equation:- This includes the time derivative term along with the spatial derivatives.

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

The logic of this method is the same as that of the steady state. The only difference is that here we also run a convergence loop, which checks the value of convergence for every time step. The rest of the methods such as Jacobi, Gauss-Siedel, and SOR remain the same.

Here is the code for the implicit way:-

clear all
close all
clc

nx = 10;
ny = 10;
T = 450*ones(nx,ny);

T(2:nx-1,1) = 400; %left face
T(2:nx-1,10) = 800; %right face
T(1,2:ny-1) = 600; %top face
T(10,2:ny-1) = 900; %bottom face

T(1,1) = 500;
T(1,ny) = 700;
T(nx,ny) = 850;
T(nx,1) = 650;

Told = T;
Tprev = T;
L = 1;
dx = L/(nx-1);
dy = dx;
error = 1;
conv = 1;
tol = .0001;
ct = 1;
t1 = 1;
t = 100;
dt = .001;
alpha = 1;
m1 = alpha*dt/dx^2;
m2 = alpha*dt/dy^2;

method = 3;

if method == 1
    
   
while conv > tol
 while error > tol
  for i = 2:nx-1
      for j = 2:ny-1
          T(i,j) = (1+2*m1+2*m2)^-1*(Tprev(i,j) + m1*(Told(i+1,j)+Told(i-1,j)) + m2*(Told(i,j-1)+Told(i,j+1))); % (1+2m1 + 2m2)
      end
  end
  error = max(max(abs(T-Told)));
  Told = T;
  ct = ct+1;
  figure(1)
  contourf(flipud(T),'ShowText','ON')
  title(sprintf('JACOBI iterations = %d \nTime iterations = %d ',ct,t1));
  colorbar
 end
%  conv = max(max(abs(T-Tprev)));
 Tprev = T;
 t1 = t1+1;
%Reinitialize error for next timestep
error = 9e9;
end

end

if method == 2
    
while conv>tol      
while error > tol
  for i = 2:nx-1
      for j = 2:ny-1
          T(i,j) = (1+2*m1+2*m2)^-1*(Tprev(i,j) + m1*(Told(i+1,j)+T(i-1,j)) + m2*(T(i,j-1)+Told(i,j+1)));
      end
  end
  error = max(max(abs(T-Told)));
  Told = T;
  ct = ct+1;
  figure(1)
  contourf(flipud(T),'ShowText','ON')
  title(sprintf('GS iterations = %d \nTime iterations = %d ',ct,t1));
  colorbar
end
 conv = max(max(abs(T-Tprev)));
 Tprev = T;
 t1 = t1+1;
%Reinitialize error for next timestep
error = 9e9;
end
end


if method == 3
    omega = 1.5;
while conv>tol    
  error = 1;  
while error > tol
  for i = 2:nx-1
      for j = 2:ny-1
          TGS = (1+2*m1+2*m2)^-1*(Tprev(i,j) + m1*(Told(i+1,j)+T(i-1,j)) + m2*(T(i,j-1)+Told(i,j+1)));
          T(i,j) = (1-omega)*Told(i,j) + omega*TGS;
      end
  end
  error = max(max(abs(T-Told)));
  Told = T;
  ct = ct+1;
  figure(1)
  contourf(flipud(T),'ShowText','ON')
  title(sprintf('SOR iterations = %d \nTime iterations = %d ',ct,t1));
  colorbar
end
 conv = max(max(abs(T-Tprev)));
 Tprev = T;
 t1 = t1+1;
%Reinitialize error for next timestep
error = 9e9;
end
end

Here is the code for explicit way:-

clear all
close all
clc

nx = 10;
ny = 10;
T = 450*ones(nx,ny);

T(2:nx-1,1) = 400; %left face
T(2:nx-1,10) = 800; %right face
T(1,2:ny-1) = 600; %top face
T(10,2:ny-1) = 900; %bottom face

T(1,1) = 500;
T(1,ny) = 700;
T(nx,ny) = 850;
T(nx,1) = 650;

Tprev = T;
L = 1;
dx = L/(nx-1);
dy = dx;
conv = 1;
tol = .0001;
ct = 1;
dt = .001;
alpha = 1;
m1 = (alpha*dt)/dx^2;
m2 = (alpha*dt)/dy^2;
cfl = m1+m2;
    
while conv > tol
  for i = 2:nx-1
      for j = 2:ny-1
          T(i,j) = Tprev(i,j)+m1*(Tprev(i-1,j)-2*Tprev(i,j)+Tprev(i+1,j)) + m2*(Tprev(i,j-1)-2*Tprev(i,j)+Tprev(i,j+1)); % negative sign for central difference
      end
  end
  conv = max(max(abs(T-Tprev)));
  Tprev = T;
  ct = ct+1;
end

figure(1)
contourf(flipud(T),'ShowText','ON');   %flipud commands flips the matrix upside down as contourf command inverts it. This is done to make it consistent with BC's
title(sprintf('Time iterations = %d',ct));
colorbar

Here is the output:-