Numerical solution of 2 dimensional 'Heat Conduction Equation' (in steady state and transient forms) using 3 different iterative solvers
Project Aim: To solve the 2 dimensional form of heat conduction equation numerically by using point iterative techniques. The forms of the equation that would be solved are:
1. Steady state form : `(del^2T)/(delx^2) + (del^2T)/(dely^2)=0` or `gradT = 0`
2. Transient form : `(delT)/(delt) + alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))=0`
where : x and y are the co-ordinate axes,
t = time,
T = Temperature,
`alpha =` thermal diffusivity
Assumptions:
1. Domain assumed is a unit square (i.e. x = y = 1m).
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. Absolute error criteria is 1e-4.
5. No internal heat generation is considered in the equation.
Expected outputs:
1. To write correct code for both implicit and explicit schemes (This is for the Transient Analysis only).
2. Point iterative techniques used are
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation (with assumed value of 'w' the relaxation factor).
Solution procedure for Part 1 :
Part1: Solving the steady state form : `(del^2T)/(delx^2) + (del^2T)/(dely^2)=0` or `gradT = 0`
Discretization equation used for numerical solution is :
`T(i,j) = 1/k((T(i-1,j)+T(i+1,j))/(Deltax^2)+(T(i,j-1)+T(i,j+1))/(Deltay^2))`
where `k = 2xx (Deltax^2+Deltay^2)/(Deltax^2xxDeltay^2`
Step 1 : Declaration of fixed variables.
Step 2 : Creation of 2D grid.
Step 3 : Creation of Initial 2D Temperature matrix.
Step 4 : Assigning BC to the 2D Temperature matrix.
Step 5 : Creating a copy of T.
Step 6 : Formulation of iterative solvers.
The Matlab code based on the above mentioned steps with detailed description in the form of comments is given below:
% Code to solve 2D Steady State Heat Conduction Equation without internal
% heat generation using 3 iterative techniques
% d^2T*dx^2 + d^2T/dy^2 = 0
clear all
close all
clc
% Step 1 : Declaration of fixed variables
lx = 1;
ly = 1;
% Step 2 : Creation of 2D grid
nx = 11;
ny = nx;
dx = lx/(nx-1);
dy = ly/(ny-1);
k = (2*(dx^2+dy^2))/(dx^2*dy^2);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% Step 3 : Creation of Initial 2D Temperature matrix
T = 300 .* ones(nx,ny);
% Step 4 : Assigning BC to the 2D Temperature Matrix
T(:,1) = 400; % Left
T(:,end) = 800; % Right
T(1,:) = 600; % Top
T(end,:) = 900; % Bottom
T(1,1) = 500; % Top left
T(1,end) = 700; % Top right
T(end,end) = 850; %Bottom right
T(end,1) = 650; % Bottom left
% Step 5 : Creating a copy of T
Told = T;
Tsor = T;
error = 9e9;
% Step 6 : Formulation of iterative solvers
iterative_solver = 3; % 1 = jacobi , 2 = Gauss seidel, 3 = SOR
% Jacobi method
threshold_error = 1e-4;
if iterative_solver == 1
jacobi_iter = 1;
while (error > threshold_error)
% For loop for iterations (space loop)
for i = 2:nx-1
for j = 2:ny-1
term1 = (1/k)*((Told(i-1,j) + Told(i+1,j))/dx^2);
term2 = (1/k)*((Told(i,j-1) + Told(i,j+1))/dy^2);
T(i,j) = term1 + term2 ;
end
end
error = max(abs(Told(i,j)-T(i,j))); % Computing the error
Told = T;% Updating the velocities after computing the error. This value will be used for the next iteration
jacobi_iter = jacobi_iter + 1;
end
sprintf('Number of iterations for Jacobi Method = %d',jacobi_iter)
contourf(x,y,T,'ShowText','on')
grid on;
xlabel('length');
ylabel('height');
title('Temperature variation in 2D grid');
set(gca,'YDIR','reverse');
end
% Gauss Seidel method
if iterative_solver == 2
Gauss_seidel_iter = 1;
while (error > threshold_error)
% For loop for iterations (space loop)
for i = 2:nx-1
for j = 2:ny-1
term1 = (1/k)*((T(i-1,j) + Told(i+1,j))/dx^2);
term2 = (1/k)*((T(i,j-1) + Told(i,j+1))/dy^2);
T(i,j) = term1 + term2 ;
end
end
error = max(abs(Told(i,j)-T(i,j))); % Computing the error
Told = T;% Updating the velocities after computing the error. This value will be used for the next iteration
Gauss_seidel_iter = Gauss_seidel_iter + 1;
end
sprintf('Number of iterations for Gauss Seidel Method = %d',Gauss_seidel_iter)
contourf(x,y,T,'ShowText','on')
grid on;
xlabel('length');
ylabel('height');
title('Temperature variation in 2D grid');
set(gca,'YDIR','reverse');
end
% Successive Over Relaxation (SOR) method
if iterative_solver == 3
w = 1.1;
SOR_iter = 1;
while (error > threshold_error)
% For loop for iterations (space loop)
for i = 2:nx-1
for j = 2:ny-1
term1 = (1/k)*((T(i-1,j) + Told(i+1,j))/dx^2);
term2 = (1/k)*((T(i,j-1) + Told(i,j+1))/dy^2);
T(i,j) = term1 + term2;
Tsor(i,j) = (1-w)*Told(i,j) + w*T(i,j);
end
end
error = max(abs(Told(i,j)- Tsor(i,j))); % Computing the error
Told = Tsor;% Updating the Temperatures after computing the error. This value will be used for the next iteration
SOR_iter = SOR_iter + 1;
end
sprintf('Number of iterations for SOR Method = %d',SOR_iter)
contourf(x,y,Tsor,'ShowText','on')
grid on;
xlabel('length');
ylabel('height');
title('Temperature variation in 2D grid');
set(gca,'YDIR','reverse');
end
Results and Explanation for Part 1 :
Results :
for iterative solver = 1 i.e. for Jacobi method, Number of iterations are 207
for iterative solver = 2 i.e. for Gauss Seidel method, Number of iterations are 109
for iterative solver = 1 i.e. for successive over relaxation method, Number of iterations are 99
Explanation :
Since the successive over relaxation method extrapolates the results obtained from Gauss Seidel method (in the above written code), it coverges very quickly as compared to Gauss Seidel method. Also, as Gauss Seidel method uses the latest values at the previous nodes, it takes lesser number of iterations to converge as compared to the Jacobi method which takes old values at nodes during calculations at each node.
Solution procedure for Part 2 :
Part 2 : Solving the transient form : `(delT)/(delt) + alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))=0`
I. With explicit iterative approach :
Discretization equation used for numerical solution is :
`T(i,j)^(n+1) =T(i,j)^(n) + alpha Deltat ( (T(i+1,j)^n - 2 T(i,j)^n +T(i-1,j)^n)/(Deltax^2)+ (T(i,j+1)^n - 2 T(i,j)^n +T(i,j-1)^n)/(Deltay^2) )`
Step 1 : Declaration of fixed variables.
Step 2 : Creation of 2D grid.
Step 3 : Creation of Initial 2D Temperature matrix.
Step 4 : Assigning BC to the 2D Temperature matrix.
Step 5 : Declaration of other constants required for the iterations.
Step 6 : Formulations of time and space loops.
The Matlab code based on the above mentioned steps with detailed description in the form of comments is given below:
% Code to solve 2D Transient State Heat Conduction Equation with explicit
% iterations
% d^2T/dt^2 + alpha * (d^2T*dx^2 + d^2T/dy^2) = 0
% alpha = thermal diffusivity
close all
clc
% Step 1 : Declaration of fixed variables
lx = 1; % Length of X axis
ly = 1; % Length of Y axis
% Step 2 : Creation of 2D grid
nx = 11;
ny = nx;
dx = lx/(nx-1);
dy = ly/(ny-1);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% Step 3 : Creation of Initial 2D Temperature matrix
T_n = 300 .* ones(nx,ny);
% Step 4 : Assigning BC to the 2D Temperature Matrix
T_n(:,1) = 400; % Left edge
T_n(:,end) = 800; % Right edge
T_n(1,:) = 600; % Top face edge
T_n(end,:) = 900; % Bottom edge
T_n(1,1) = 500; % Top left corner
T_n(1,end) = 700; % Top right corner
T_n(end,end) = 850; %Bottom right corner
T_n(end,1) = 650; % Bottom left corner
%Step 5 : Other constants required for the iterations
alpha = 1.4;
dt = 1.7e-3;
no_of_time_steps = 1400 ;
total_time = dt * no_of_time_steps;
k1 = alpha * (dt/dx^2);
k2 = alpha * (dt/dy^2);
T_n_plus_one = T_n;
% Step 6 : Formulations of time and space loops
iterations = 0;
for nt = 1:no_of_time_steps
for i = 2:nx-1
for j = 2:ny-1
term1 = k1 * (T_n(i-1,j) - 2*(T_n(i,j)) + T_n(i+1,j));
term2 = k2 * (T_n(i,j-1) - 2*(T_n(i,j)) + T_n(i,j+1));
T_n_plus_one(i,j) = T_n(i,j) + term1 + term2;
T_n(i,j) = T_n_plus_one(i,j);
iterations = iterations + 1;
end
end
end
sprintf('Number of iterations = %d',iterations)
contourf(x,y,T_n_plus_one,'ShowText','on')
grid on;
xlabel('length');
ylabel('height');
title('Temperature variation in 2D grid');
set(gca,'YDIR','reverse'); Results and Explanation for Part 2 (I) :
Result :
After running the code the ouput shows:
'Number of iterations = 113400'
Explanation :
Number of iterations = No. of itertaions for j x No. of itertaions for i x No. of itertaions for time steps
= 9 x 9 x 1400
= 113400
II. With implicit iterative approach :
Discretization equation used for numerical solution is :
`T(i,j)^(n+1) =(1/(1+2k1+k2)) ( T(i,j)^(n) + k1 (T(i-1,j)^(n+1) +T(i+1,j)^(n+1)) +k2(T(i,j+1)^(n+1) +T(i,j-1)^(n+1) )`
where : `k1 = (alphaDeltat)/(Deltax^2`
`k2 = (alphaDeltat)/(Deltay^2`
Step 1 : Declaration of fixed variables.
Step 2 : Creation of 2D grid.
Step 3 : Creation of Initial 2D Temperature matrix.
Step 4 : Assigning BC to the 2D Temperature Matrix.
Step 5 : Creating a copy of T_n.
Step 6 : Other constants required for the iterations.
The Matlab code based on the above mentioned steps with detailed description in the form of comments is given below:
% Code to solve 2D Transient State Heat Conduction Equation with implicit solver
% d^2T/dt^2 + alpha * (d^2T*dx^2 + d^2T/dy^2) = 0
% alpha = thermal diffusivity
clear all
close all
clc
% Step 1 : Declaration of fixed variables
lx = 1; % Length of X axis
ly = 1; % Length of Y axis
% Step 2 : Creation of 2D grid
nx = 15;
ny = nx;
dx = lx/(nx-1);
dy = ly/(ny-1);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
% Step 3 : Creation of Initial 2D Temperature matrix
T_n = 300 .* ones(nx,ny);
% Step 4 : Assigning BC to the 2D Temperature Matrix
T_n(:,1) = 400; % Left edge
T_n(:,end) = 800; % Right edge
T_n(1,:) = 600; % Top face edge
T_n(end,:) = 900; % Bottom edge
T_n(1,1) = 500; % Top left corner
T_n(1,end) = 700; % Top right corner
T_n(end,end) = 850; %Bottom right corner
T_n(end,1) = 650; % Bottom left corner
% Step 5 : Creating a copy of T_n
Told = T_n;
T_n_plus_one = T_n;
Tsor = T_n;
% Step 6 : Other constants required for the iterations
alpha = 1.2;
dt = 1.7e-3;
no_of_time_steps = 1400 ;
total_time = dt * no_of_time_steps;
k1 = alpha * (dt/dx^2);
k2 = alpha * (dt/dy^2);
jacobi_iter = 0;
Gauss_seidel_iter = 0;
SOR_iter = 0;
% Step 7 : Formulation of the loops
time_iterations = 0;
iterative_solver = 3;
for nt = 1:no_of_time_steps
threshold_error = 1e-4;
%Convergence loops
% Jacobi Method
if (iterative_solver == 1)
error_jacobi = 1;
while (error_jacobi > threshold_error)
for i = 2:nx-1
for j = 2:ny-1
term1 =(1+(2*k1)+(2*k2))^-1;
term2 = k1 * term1;
term3 = k2 * term1;
H = (Told(i-1,j) + Told(i+1,j));
V = (Told(i,j-1) + Told(i,j+1));
T_n_plus_one(i,j) = (term1*T_n(i,j))+(term2*H)+(term3*V);
end
end
jacobi_iter = jacobi_iter + 1;
error_jacobi = max(max(abs(Told - T_n_plus_one)));
Told = T_n_plus_one;
end
T_n = T_n_plus_one;
contourf(x,y,T_n_plus_one,'ShowText','on')
grid on;
xlabel('length');
ylabel('height');
title('Temperature variation in 2D grid');
set(gca,'YDIR','reverse');
end
% Gauss Seidel method
if (iterative_solver == 2)
error_Gauss_seidel = 1;
while (error_Gauss_seidel > threshold_error)
% For loop for iterations (space loop)
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
term2 = k1 * term1;
term3 = k2 * term1;
H = (T_n_plus_one(i-1,j) + Told(i+1,j));
V = (T_n_plus_one(i,j-1) + Told(i,j+1));
T_n_plus_one(i,j) = (term1*T_n(i,j))+(term2*H)+(term3*V);
end
end
Gauss_seidel_iter = Gauss_seidel_iter + 1;
error_Gauss_seidel = max(max(abs(Told - T_n_plus_one)));
Told = T_n_plus_one;
end
T_n = T_n_plus_one;
contourf(x,y,T_n_plus_one,'ShowText','on')
grid on;
xlabel('length');
ylabel('height');
title('Temperature variation in 2D grid');
set(gca,'YDIR','reverse');
end
% Successive Over Relaxation (SOR) method
if (iterative_solver == 3)
w = 1.3;
error_SOR = 1;
while (error_SOR > threshold_error)
% For loop for iterations (space loop)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1/(1+(2*k1)+(2*k2));
term2 = k1 * term1;
term3 = k2 * term1;
H = T_n_plus_one(i-1,j) + Told(i+1,j);
V = T_n_plus_one(i,j-1) + Told(i,j+1);
T_n_plus_one(i,j) = (term1*T_n(i,j))+(term2*H)+(term3*V);
Tsor(i,j) = (1-w)*Told(i,j) + w*(T_n_plus_one(i,j));
end
end
error_SOR = max(max(abs(Told - Tsor))); % Computing the error
Told = Tsor;% Updating the Temperatures after computing the error. This value will be used for the next iteration
SOR_iter = SOR_iter + 1;
end
T_n = Tsor;
contourf(x,y,Tsor,'ShowText','on')
grid on;
xlabel('length');
ylabel('height');
title('Temperature variation in 2D grid');
set(gca,'YDIR','reverse');
end
endResults and Explanation for Part 2 (II) :
Result :
After running the code the ouput shows:
Number of iterations for the Jacobi solver was found to be maximum and that for the SOR solver was found to be minimum.
Explanation :
Since the successive over relaxation method extrapolates the results obtained from Gauss Seidel method (in the above written code), it coverges very quickly as compared to Gauss Seidel method. Also, as Gauss Seidel method uses the latest values at the previous nodes, it takes lesser number of iterations to converge as compared to the Jacobi method which takes old values at nodes during calculations at each node. Error obtained in all the three cases is below the threshold error. So, it means all the three solvers converged.
Conclusions :
Temperature profile initially is as shown in figure:
Temperature variation after execution of the codes for all the possible cases is as shown:
1) This project demonstrates the understanding of how a simple second order Partial Differential Equation can be solved using Steady state and transient state consideration by considering large number of time steps (1400) to reach the steady state solution.
2) This projects demonstrates the proper application of various basic iterative solvers i.e. Jacobi, Gauss Seidel and Successive Over Relaxation(SOR) and provides the basic codes structure to develop one's own solver.