Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
2D Heat Conduction Equation
`(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))`
For solving this equation numerically a two dimensional fluid domain is created with both length and width of 1m. The domain is divided into 10 number of nodes in each x and y axis. The coefficient of thermal diffusivity is taken to be 1.4 by hit and trial method. And error is kept as 1e-4 for the convergence of the solution. The relaxation factor for succesive over relaxtion solver is 1.18. The maximum CFL number for the stability of the solution is kept 0.5. The boundary conditions are given in the challenge as-
Temperature at left boundary = 400K
Temperature at right boundary = 800K
Temperature at top boundary = 600K
Temperature at bottom boundary = 900K
Initial a temperature of 300K is used as the guess value for the iterative solvers.
Unsteady 2D Heat conduction equation
Implicit
The numerical equation for the unsteady state can derived as follows-
`(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))`
The temperature derivative with respect to time is approximated by using implicit first order forder difference scheme as
`(delT)/(delt) = (T^(n+1)(i,j) - T^n(i,j))/(Deltat)`
And the second order partial derivatives of the temperature with respect to x and y are approximated using second order central difference scheme so that the information from all the neighbouting point can be taken into account for the parabolic solution of convective part as -
`(del^2T)/(delx^2) = (T^(n+1)(i+1,j) - 2T^(n+1)(i,j) + T^(n+1)(i-1,j))/(Deltax^2)`
`(del^2T)/(dely^2) = (T^(n+1)(i,j+1) - 2T^(n+1)(i,j) + T^(n+1)(i,j-1))/(Deltay^2)`
where,
T(i,j) = current point to be calculated
T(i+1,j) = Point right to the the current point
T(i-1,j) = Point left to the the current point
T(i,j+1) = Point above to the the current point
T(i,j-1) = Point below to the the current point
And n is for the current iteration and n+1 is the next iteration.
Therefore the equation is -
`T^(n+1)(i,j) = (T^(n)(i,j) + k1(T^(n+1)(i+1,j) + T^(n+1)(i-1,j)) + K2(T^(n+1)(i,j+1) + T^(n+1)(i,j-1)))/(1 + 2K1 +2K2)`
where `K1 = (alpha(Deltat))/(Deltax^2)` and `K2 = (alpha(Deltat))/(Deltay^2)`are CFL number along x and y axis respectivel`
K1 + K2 = CFL =0.5
two loops are used one for iterating time and one for space marching over the complete domain
All of the three iterative solvers are used for finding the solution of unsteady implicit equation
Explicit
For explicit solution the temperature at the next iteration is formulated in terms of previously known temperature (previous iterations) for convective part.
Therefore-
`(delT)/(delt) = (T^(n+1)(i,j) - T^n(i,j))/(Deltat)`
`(del^2T)/(delx^2) = (T^(n)(i+1,j) - 2T^(n)(i,j) + T^(n)(i-1,j))/(Deltax^2)`
`(del^2T)/(dely^2) = (T^(n)(i,j+1) - 2T^(n)(i,j) + T^(n)(i,j-1))/(Deltay^2)`
The final equation is given as
`T^(n+1)(i,j) = (1 - 2K1 - 2K2)T^n(i,j) + K1(T^n(i+1,j) + T^n(i-1,j)) + K2(T^n(i,j+1) + T^n(i,j-1))`
All of the three given solvers are used for computing the solution.
Steady 2D Heat conduction equation
For steady equations the temperature is independent of time therefore the time derivative of temperature will vanish from the equation. The 2D heat conduction equation for steady state is given as:-
`(del^2T)/(delx^2) + (del^2T)/(delx^2) = 0`
For the above equation the spatial derivatives of temperature will be discritized as discussed above.
`(del^2T)/(delx^2) = (T^(n+1)(i+1,j) - 2T^(n+1)(i,j) + T^(n+1)(i-1,j))/(Deltax^2)`
`(del^2T)/(dely^2) = (T^(n+1)(i,j+1) - 2T^(n+1)(i,j) + T^(n+1)(i,j-1))/(Deltay^2)`
By substituting these, the equation for steady 2D heat conduction equation is given as -
`T(i,j) = (K1(T(i+1,j) + T(i-1,j)) + K2(T(i,j+1) + T(i,j-1)))/(2(K1 + K2))`
where K1 = `1/(Deltax^2)`and K2 = `1/(Deltay^2)`
Since there is no time derivative all the temperatures in the above equations are at the same instant of time therefore there is no need for time marching, only space marching is required to find the solution of steady 2D heat conduction equation. Since all the temperatures are unknown at the same time, iterative solvers are also used to find the solution of the equation.
All of the three solvers mentioned in the problem are used to find the solution.
The code for the problem is given as below:-
clear all
close all
clc
Lx = 1; % length of rectangle
Ly = 1; % width of rectangle
nx = 10; % number of nodes in x
ny = 10; % number of nodes in y
nt = 1400; % number of iterations for time or time loop
x = linspace(0,Lx,nx); % Meshing in x
dx = x(2) - x(1); % length of element along x
y = linspace(0,Ly,ny); % Meshing in y
dy = y(2) - y(1); % length of element alongy
absolute_error = 1e-4; % absolute error criteria for convergence
alpha = 1.4; % coefficient of thermal diffusivity
w = 1.18; % relaxation factor for successive over relaxation iterative solver
CFL_max = 0.5; % criteria for stability
dt = (CFL_max*((dx^2)*(dy^2)))/(alpha*((dx^2)+(dy^2))); % time step
k1 = (alpha*dt)/(dx^2); % CFL along x
k2 = (alpha*dt)/(dy^2); % CFL along y
k1_steady = (dt)/(dx^2); % CFL along x
k2_steady = (dt)/(dy^2); % CFL along y
T_initial = 300*ones(ny,nx); % Initial guess for iterative method
T_initial(1,:) = 600; % Top boundary condition
T_initial(ny,:) = 900; % Bottom boundary condition
T_initial(:,1) = 400; % Left boundary condition
T_initial(:,nx) = 800; % Right boundary condition
T_initial(1,1) = 500; % Average boundary condition for left top corner
T_initial(1,nx) = 700; % Average boundary condition for right top corner
T_initial(ny,1) = 650; % Average boundary condition for left bottom corner
T_initial(ny,nx) = 850; % Average boundary condition for right bottom corner
T_previous_jacobi_unsteady_implicit = T_initial; % temperature at previous time step
T_old_jacobi_unsteady_implicit = T_initial; % temperature before iteration at the same time step
T_jacobi_unsteady_implicit = T_initial; % temperature after the iteration at the same time step
%Jacobi iterative solver for unsteady implicit 2D heat conduction equation
jacobi_iteration_unsteady_implicit = 1;
for t = 1:nt % loop for time marching
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_jacobi_unsteady_implicit(i,j) = (T_previous_jacobi_unsteady_implicit(i,j) + k1*(T_old_jacobi_unsteady_implicit(i+1,j) + T_old_jacobi_unsteady_implicit(i-1,j)) + k2*(T_old_jacobi_unsteady_implicit(i,j+1) + T_old_jacobi_unsteady_implicit(i,j-1)))/(1+(2*k1)+(2*k2)); % Discritized 2D heat condunction equation in implicit form using first order forward difference scheme
error(i,j) = abs(T_old_jacobi_unsteady_implicit(i,j) - T_jacobi_unsteady_implicit(i,j)); % error between the temperature before and after iteration
end
end
if any(error(:) > absolute_error)
T_old_jacobi_unsteady_implicit = T_jacobi_unsteady_implicit;
jacobi_iteration_unsteady_implicit = jacobi_iteration_unsteady_implicit + 1;
else
error_num = -1;
end
end
T_previous_jacobi_unsteady_implicit = T_jacobi_unsteady_implicit;
figure(1)
contourf(x,y,T_jacobi_unsteady_implicit)
clabel(contourf(x,y,T_jacobi_unsteady_implicit))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Unsteady Implicit Jacobi Iterations = %d', jacobi_iteration_unsteady_implicit));
pause(0.03)
end
T_previous_GS_unsteady_implicit = T_initial; % temperature at previous time step
T_old_GS_unsteady_implicit = T_initial; % temperature before iteration at the same time step
T_GS_unsteady_implicit = T_initial; % temperature after the iteration at the same time step
%Gauss-Seidel iterative solver for unsteady implicit 2D heat conduction equation
Gauss_seidel_iteration_unsteady_implicit = 1;
for t = 1:nt % loop for time marching
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_GS_unsteady_implicit(i,j) = (T_previous_GS_unsteady_implicit(i,j) + k1*(T_GS_unsteady_implicit(i+1,j) + T_GS_unsteady_implicit(i-1,j)) + k2*(T_GS_unsteady_implicit(i,j+1) + T_GS_unsteady_implicit(i,j-1)))/(1+(2*k1)+(2*k2)); % Discritized 2D heat condunction equation in implicit form using first order forward difference scheme
error(i,j) = abs(T_old_GS_unsteady_implicit(i,j) - T_GS_unsteady_implicit(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_old_GS_unsteady_implicit = T_GS_unsteady_implicit;
Gauss_seidel_iteration_unsteady_implicit = Gauss_seidel_iteration_unsteady_implicit + 1;
else
error_num = -1;
end
end
T_previous_GS_unsteady_implicit = T_GS_unsteady_implicit;
figure(2)
contourf(x,y,T_GS_unsteady_implicit)
clabel(contourf(x,y,T_GS_unsteady_implicit))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Unsteady Implicit Gauss-siedel Iterations = %d', Gauss_seidel_iteration_unsteady_implicit));
pause(0.03)
end
T_previous_SOR_unsteady_implicit = T_initial; % temperature at previous time step
T_old_SOR_unsteady_implicit = T_initial; % temperature before iteration at the same time step
T_SOR_unsteady_implicit = T_initial; % temperature after the iteration at the same time step
%Successive Over Relaxation iterative solver for unsteady implicit 2D heat conduction equation
SOR_iteration_unsteady_implicit = 1;
for t = 1:nt % loop for time marching
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_SOR_unsteady_implicit(i,j) = w*(T_previous_SOR_unsteady_implicit(i,j) + k1*(T_SOR_unsteady_implicit(i+1,j) + T_SOR_unsteady_implicit(i-1,j)) + k2*(T_SOR_unsteady_implicit(i,j+1) + T_SOR_unsteady_implicit(i,j-1)))/(1+(2*k1)+(2*k2)) + (1 - w)*T_SOR_unsteady_implicit(i,j); % Discritized 2D heat condunction equation in implicit form using first order forward difference scheme
error(i,j) = abs(T_old_SOR_unsteady_implicit(i,j) - T_SOR_unsteady_implicit(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_old_SOR_unsteady_implicit = T_SOR_unsteady_implicit;
SOR_iteration_unsteady_implicit = SOR_iteration_unsteady_implicit + 1;
else
error_num = -1;
end
end
T_previous_SOR_unsteady_implicit = T_SOR_unsteady_implicit;
figure(3)
contourf(x,y,T_SOR_unsteady_implicit)
clabel(contourf(x,y,T_SOR_unsteady_implicit))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Unsteady Implicit Successive Over Relaxation Iterations = %d', SOR_iteration_unsteady_implicit));
pause(0.03)
end
T_previous_jacobi_unsteady_explicit = T_initial; % temperature at previous time step or before iteration
T_jacobi_unsteady_explicit = T_initial; % temperature after the iteration at the same time step
%Jacobi iterative solver for unsteady explicit 2D heat conduction equation
jacobi_iteration_unsteady_explicit = 1;
for t = 1:nt
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_jacobi_unsteady_explicit(i,j) = ((1 - (2*k1) - (2*k2))*T_previous_jacobi_unsteady_explicit(i,j)) + k1*(T_previous_jacobi_unsteady_explicit(i+1,j) + T_previous_jacobi_unsteady_explicit(i-1,j)) + k2*(T_previous_jacobi_unsteady_explicit(i,j+1) + T_previous_jacobi_unsteady_explicit(i,j-1))
error(i,j) = abs(T_previous_jacobi_unsteady_explicit(i,j) - T_jacobi_unsteady_explicit(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_previous_jacobi_unsteady_explicit = T_jacobi_unsteady_explicit;
jacobi_iteration_unsteady_explicit = jacobi_iteration_unsteady_explicit + 1;
else
T_previous_jacobi_unsteady_explicit = T_jacobi_unsteady_explicit;
error_num = -1;
end
end
figure(4)
contourf(x,y,T_jacobi_unsteady_explicit)
clabel(contourf(x,y,T_jacobi_unsteady_explicit))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Unsteady explicit Jacobi Iterations = %d', jacobi_iteration_unsteady_explicit));
pause(0.3)
end
T_previous_GS_unsteady_explicit = T_initial; % temperature at previous time step or before iteration
T_GS_unsteady_explicit = T_initial; % temperature after the iteration at the same time step
%Gauss Seidal iterative solver for unsteady explicit 2D heat conduction equation
GS_iteration_unsteady_explicit = 1;
for t = 1:nt
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_GS_unsteady_explicit(i,j) = (1 - (2*k1) - (2*k2))*T_GS_unsteady_explicit(i,j) + k1*(T_GS_unsteady_explicit(i+1,j) + T_GS_unsteady_explicit(i-1,j)) + k2*(T_GS_unsteady_explicit(i,j+1) + T_GS_unsteady_explicit(i,j-1))
error(i,j) = abs(T_previous_GS_unsteady_explicit(i,j) - T_GS_unsteady_explicit(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_previous_GS_unsteady_explicit = T_GS_unsteady_explicit;
GS_iteration_unsteady_explicit = GS_iteration_unsteady_explicit + 1;
else
error_num = -1;
end
end
figure(5)
contourf(x,y,T_GS_unsteady_explicit)
clabel(contourf(x,y,T_GS_unsteady_explicit))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Unsteady explicit Gauss-seidel Iterations = %d', GS_iteration_unsteady_explicit));
pause(0.03)
end
T_previous_SOR_unsteady_explicit = T_initial; % temperature at previous time step or before iteration
T_SOR_unsteady_explicit = T_initial; % temperature after the iteration at the same time step
%Successive Over Relaxation iterative solver for unsteady explicit 2D heat conduction equation
SOR_iteration_unsteady_explicit = 1;
for t = 1:nt
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_SOR_unsteady_explicit(i,j) = w*((1 - (2*k1) - (2*k2))*T_SOR_unsteady_explicit(i,j) + k1*(T_SOR_unsteady_explicit(i+1,j) + T_SOR_unsteady_explicit(i-1,j)) + k2*(T_SOR_unsteady_explicit(i,j+1) + T_SOR_unsteady_explicit(i,j-1))) + (1-w)*T_SOR_unsteady_explicit(i,j)
error(i,j) = abs(T_previous_SOR_unsteady_explicit(i,j) - T_SOR_unsteady_explicit(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_previous_SOR_unsteady_explicit = T_SOR_unsteady_explicit;
SOR_iteration_unsteady_explicit = SOR_iteration_unsteady_explicit + 1;
else
error_num = -1;
end
end
figure(6)
contourf(x,y,T_SOR_unsteady_explicit)
clabel(contourf(x,y,T_SOR_unsteady_explicit))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Unsteady explicit Successive Over Relaxation Iterations = %d', SOR_iteration_unsteady_explicit));
pause(0.03)
end
T_old_Jacobi_steady = T_initial; % temperature at previous time step or before iteration
T_Jacobi_steady = T_initial; % temperature after the iteration at the same time step
%Jacobi iterative solver for steady explicit 2D heat conduction equation
Jacobi_iteration_steady = 1;
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_Jacobi_steady(i,j) = (k1_steady*(T_old_Jacobi_steady(i+1,j) + T_old_Jacobi_steady(i-1,j)) + k2_steady*(T_old_Jacobi_steady(i,j+1) + T_old_Jacobi_steady(i,j-1)))/(2*(k1_steady + k2_steady))
error(i,j) = abs(T_old_Jacobi_steady(i,j) - T_Jacobi_steady(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_old_Jacobi_steady = T_Jacobi_steady;
Jacobi_iteration_steady = Jacobi_iteration_steady + 1;
else
error_num = -1;
end
figure(7)
contourf(x,y,T_Jacobi_steady)
clabel(contourf(x,y,T_Jacobi_steady))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Steady Jacobi Iterations = %d', Jacobi_iteration_steady));
pause(0.03)
end
T_old_GS_steady = T_initial; % temperature at previous time step or before iteration
T_GS_steady = T_initial; % temperature after the iteration at the same time step
%Gauss Seidal iterative solver for steady explicit 2D heat conduction equation
GS_iteration_steady = 1;
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_GS_steady(i,j) = (k1_steady*(T_GS_steady(i+1,j) + T_GS_steady(i-1,j)) + k2_steady*(T_GS_steady(i,j+1) + T_GS_steady(i,j-1)))/(2*(k1_steady + k2_steady))
error(i,j) = abs(T_old_GS_steady(i,j) - T_GS_steady(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_old_GS_steady = T_GS_steady;
GS_iteration_steady = GS_iteration_steady + 1;
else
error_num = -1;
end
figure(8)
contourf(x,y,T_GS_steady)
clabel(contourf(x,y,T_GS_steady))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Steady Guass-seidel Iterations = %d', GS_iteration_steady));
pause(0.03)
end
T_old_SOR_steady = T_initial; % temperature at previous time step or before iteration
T_SOR_steady = T_initial; % temperature after the iteration at the same time step
%Succesive Over Relaxation iterative solver for steady explicit 2D heat conduction equation
SOR_iteration_steady = 1;
error_num = 1;
while(error_num > absolute_error) % loop for convergence
for i = 2:nx-1 % loop for space marching
for j = 2:ny-1
T_SOR_steady(i,j) = w*((k1_steady*(T_SOR_steady(i+1,j) + T_SOR_steady(i-1,j)) + k2_steady*(T_SOR_steady(i,j+1) + T_SOR_steady(i,j-1)))/(2*(k1_steady + k2_steady))) + (1-w)*T_SOR_steady(i,j)
error(i,j) = abs(T_old_SOR_steady(i,j) - T_SOR_steady(i,j));
end
end
if any(error(:) > absolute_error)
error_num = 1;
T_old_SOR_steady = T_SOR_steady;
SOR_iteration_steady = SOR_iteration_steady + 1;
else
error_num = -1;
end
figure(9)
contourf(x,y,T_SOR_steady)
clabel(contourf(x,y,T_SOR_steady))
colorbar
colormap(jet)
set(gca, 'ydir', 'reverse')
xlabel('X-Axis')
ylabel('Y-Axis')
title(sprintf('Steady Successive Over Relaxation Iterations = %d', SOR_iteration_steady));
pause(0.03)
end
figure(10)
y = [jacobi_iteration_unsteady_implicit Gauss_seidel_iteration_unsteady_implicit SOR_iteration_unsteady_implicit jacobi_iteration_unsteady_explicit GS_iteration_unsteady_explicit SOR_iteration_unsteady_explicit Jacobi_iteration_steady GS_iteration_steady SOR_iteration_steady]
bar(y)
The output of the program is given below:-
Observations:-
From the above bar graph it is observed that the succesive over relaxation method is faster in convergence than the jacobi and gauss-seidel method. The number of iterations reduces from jacobi to gauss siedel to succesive over relaxation method. The order of the rate of convergene is-
Succesive Over Relaxation > Gauss siedel Method > Jacobi Method
Moreover if we compare the steady case with the unsteady case then the number of iteration in steady case is less than the number of iterations in unsteady case, since for unsteady case both space marching and time marching is needed to be solved.
Therefore, from the results it is shown that the succesive over relaxation method with proper relaxation factor (w) converges faster than the gauss siedel and jacobi method.