Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

Aim: Solving the 2D heat conduction problem for steady and transient state using various iterative techniques (Jacobi, Gauss seidel and Successive over relaxation methods) using Matlab

Given that:

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

 Your absolute error criteria is 1e-4

Steady State flow:

The generalised 2D temperature equation is,

Discretize the equation by central differencing schemes

Since we assumed nx =ny There fore delta x^2 = delta y^2

Then the equation becomes

T(i,j) = 1/4 (T(i-1,j) + T(i+1,j) + T(i, j-1) + T(i, j+1))

Matlab code for Steady State:

clear all
close all
clc

% Solve 2nd order 1-dimensional PDE
% d^2u/dx^2 = 0

% Inputs
nx = 20; % no of grids in x direction
ny = nx; % no of grids in y direction

tolerance = 1e-5;
error = 9e9; % initial guess error
T =5; % time limit
dt = 1e-2 % time step
x = linspace(0,1,nx);
y = linspace(0,1,ny);

% grid spacing
dx = x(2)-x(1);
dy = y(2)-y(1);

[xx, yy] = meshgrid(x,y);

% matrix defination
t = 300*ones(nx,ny);
t(end,:) = 900;
t(1,:) = 600;
t(:,1) = 400;
t(:,end) = 800;

k = (2*(dx^2 + dy^2)/(dx^2*dy^2));

t_old = t;
t_p = t;

% solving method
method = 'which method'
solver = input(method); % 1 for jacobi
                                    % 2 for gauss siedel
                                    % 3 for SOR

% Jacobi solver

if solver ==1
iteration_1 =1;
tic;
%converge loop
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
TL = t_old(i-1,j);
TR = t_old(i+1,j);
TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (TL + TR)/dx^2;
V = (TB + TT)/dy^2;

t(i,j) = (H + V)/k;
end
end
time =toc;
error = max (max(abs(t_old-t)));
t_old =t;
iteration_1 = iteration_1+1;

[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('Steady state_jacobi = %d',iteration_1);
title_text2 = sprintf('Simulation time +%f s',time);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end
end

% Gauss siedel Method
if solver ==2
iteration_2 = 1;
tic;

%converge loop
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
% TL = t_old(i-1,j);
TR = t_old(i+1,j);
% TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (t(i-1,j)+ TR)/dx^2;
V = (t(i,j-1) + TT)/dy^2;

t(i,j) = (H + V)/k;
end
end
time2 =toc;
error = max (max(abs(t_old-t)));
t_old =t;
iteration_2 = iteration_2+1;
time =toc;

[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('Steady state_Gauss siedel = %d',iteration_2);
title_text2 = sprintf('Simulation time +%f s',time);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end
end

% Sucessive over relaxation method (SOR)
alpha = 0.3; %Thermal diffusivity
omega = 1.2; % SOR value, should be grater than 1 and less than 2
if solver ==3
iteration_3 =1;
tic;
% convergence loop
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
% TL = t_old(i-1,j);
TR = t_old(i+1,j);
% TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (t(i-1,j)+ TR)/dx^2;
V = (t(i,j-1) + TT)/dy^2;
tgs(i,j) = (H + V)/k;

t(i,j) = (1-omega)*(t_old(i,j))+omega*tgs(i,j);

end
end
time3 =toc;
error = max (max(abs(t_old-t)));
t_old =t;
iteration_3 = iteration_3+1;
time =toc;

[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('Steady state_SOR = %d',iteration_2);
title_text2 = sprintf('Simulation time +%f s',time);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end
end

Plots:

Jacobi method:

Gauss siedel method:

SOR method:

 As from the figure, the no of iteration required to converge the solution is more in jacobi iteration scheme followed by Gauss Siedle iteration and SOR scheme.

Computational time required for a computer program to run the iteration is more in jacobi followed by Gauss siedel and SOR scheme.

Transient State (implicit):

Implicit method uses the parameters which are inter dependent with each other in the same time step. 

This method is good because of its stable characteristics.

Matlab Code:

clear all
close all
clc

% Inputs
nx = 20; % no of grids in x direction
ny = nx; % no of grids in y direction
alpha = 1.4;
omega = 1.2;

tolerance = 1e-5;
error = 9e9; % initial guess error
T =5; % time limit
dt = 1e-2 % time step
x = linspace(0,1,nx);
y = linspace(0,1,ny);

% grid spacing
dx = x(2)-x(1);
dy = y(2)-y(1);

[xx, yy] = meshgrid(x,y);

% matrix defination
t = 300*ones(nx,ny); % initial temperature guess
t(end,:) = 900;
t(1,:) = 600;
t(:,1) = 400;
t(:,end) = 800;
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
k = 1+(2*k1)+(2*k2);

t_old = t;
t_p = t;

% solving method
method = 'which method'
solver = input(method); % 1 for jacobi
% 2 for gauss siedel
% 3 for SOR

% Jacobi solver
if solver ==1
iteration_1 =1;
tic;
%converge loop
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
TL = t_old(i-1,j);
TR = t_old(i+1,j);
TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (TL + TR);
V = (TB + TT);

t(i,j) = ((t_p(i,j)+H*k1 + V*k2)/k);
end
end
time =toc;
error = max (max(abs(t_old-t)));
t_old =t;
iteration_1 = iteration_1+1;

[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('transient state implicit jacobi = %d',iteration_1);
title_text2 = sprintf('Simulation time +%f s',time);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end
end

% Gauss siedel Method
if solver ==2
iteration_2 = 1;
tic;

%converge loop
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
% TL = t_old(i-1,j);
TR = t_old(i+1,j);
% TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (t(i-1,j)+ TR);
V = (t(i,j-1) + TT);

t(i,j) = ((t_p(i,j)+H*k1 + V*k2)/k);
end
end
time2 =toc;
error = max (max(abs(t_old-t)));
t_old =t;
iteration_2 = iteration_2+1;
time =toc;

[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('Transient state implicit Gauss siedel = %d',iteration_2);
title_text2 = sprintf('Simulation time +%f s',time);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end
end

% Sucessive over relaxation method (SOR)
if solver ==3
iteration_3 =1;
tic;
% convergence loop
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
TL = t_old(i-1,j);
TR = t_old(i+1,j);
TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (t(i-1,j)+ TR);
V = (t(i,j-1) + TT);
tgs(i,j) = ((t_p(i,j)+H*k1 + V*k2)/k);

t(i,j) = (1-omega)*(t_old(i,j))+omega*tgs(i,j);

end
end
time3 =toc;
error = max (max(abs(t_old-t)));
t_old =t;
iteration_3 = iteration_3+1;
time =toc;

[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('Transient state implict SOR = %d',iteration_3);
title_text2 = sprintf('Simulation time +%f s',time);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end
end

Plots:

Jacobi Method:

Gauss siedel :

SOR method:

Transient State 2D Explicit:

In the explicit method the next values is estimated based on the already existed parameters, which means n+1 nodes is calculated by nth step parameters.

The time derivative is discretised by FDS first order while space derivative is discretised by 2nd order CDS:

Matlab code:

clear all
close all
clc

% Inputs
nx = 10; % no of grids in x direction
ny = nx; % no of grids in y direction
alpha = 0.3;
omega = 1.2;

tolerance = 1e-4;
error = 9e9; % initial guess error
T =5; % time limit
dt = 0.01 % time step
nt = T/dt;
x = linspace(0,1,nx);
y = linspace(0,1,ny);

% grid spacing
dx = x(2)-x(1);
dy = y(2)-y(1);

[xx, yy] = meshgrid(x,y);

% matrix defination
t = 300*ones(nx,ny); % initial temperature guess
t(end,:) = 900;
t(1,:) = 600;
t(:,1) = 400;
t(:,end) = 800;
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
k = 1-(2*k1)-(2*k2);

t_old = t;
t_p = t;

% solving method
method = 'which method'
solver = input(method); % 1 for jacobi
% 2 for gauss siedel
% 3 for SOR

% Jacobi solver
if solver ==1
iteration_1 =1;
nt = T/dt;
tic;
for m = 1:nt
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
TL = t_old(i-1,j);
TR = t_old(i+1,j);
TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (TL + TR);
V = (TB + TT);

t(i,j) = t_p(i,j)*k + k1*H + k2*V;
end
end

error = max (max(abs(t-t_old)));
t_old =t;
iteration_1 = iteration_1+1;
end
t_p = t;
end
time1 =toc;
[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('transient state explicit jacobi = %d',iteration_1);
title_text2 = sprintf('Simulation time +%f s',time1);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end

% Gauss siedel Method
if solver ==2
iteration_2 = 1;
nt = T/dt;
tic;
for m = 1:nt
while (error>tolerance)
for i = 2:nx-1
for j =2:ny-1
%TL =t_old(i-1,j);
TR = t_old(i+1,j);
%TB = t_old(i,j-1);
TT = t_old(i,j+1);
H = t(i-1,j) + TR;
V = t(i, j-1) + TT;
t(i,j) = t_p(i,j)*k + H*k1 + V*k2;
end
end
error = max(max(abs(t-t_old)));
t_old =t;
iteration_2 =iteration_2 +1;
end
t_p =t;
end
time2 = toc;
[x,y] = contourf(xx, yy,t)
colormap (jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('Transient state explicit Gauss siedel = %d',iteration_2);
title_text2 = sprintf('Simulation time +%f s',time2);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end

% Sucessive over relaxation method (SOR)
if solver ==3
iteration_3 =1;
tic;
for m = 1:nt
while(error>tolerance)
for i = 2:(nx-1) % loop of equations
for j = 2:(ny-1)
% TL = t_old(i-1,j);
TR = t_old(i+1,j);
% TB = t_old(i,j-1);
TT = t_old(i,j+1);

H = (t(i-1,j)+ TR);
V = (t(i,j-1) + TT);
tgs(i,j) = ((t_p(i,j)+H*k1 + V*k2)/k);

t(i,j) = (1-omega)*(t_old(i,j))+omega*tgs(i,j);

end
end
error = max (max(abs(t-t_old)));
t_old =t;
iteration_3 = iteration_3+1;
end
t_p =t;
end
time3 =toc;
[x y] = contourf(xx,yy,t)
colormap(jet)
colorbar
clabel(x,y)
pause(0.03)
title_text1 = sprintf('Transient state Explict SOR = %d',iteration_3);
title_text2 = sprintf('Simulation time +%f s',time3);
title_text3 = sprintf('final error = %f',error)
title({title_text1; title_text2; title_text3})
end

Plots:

Jacobi Method:

Gauss siedel:

SOR method:

For Explicit Scheme, we cannot use the same time steps and length of mesh domain, as it yeilds to higher CFL number, which leads to unstable solution.

Conclusions:

Comparing with different iterations approaches, we can conclude that

The computational time and iterations are in the order of  SOR<GS<Jacobi for both steady and unsteady states.

Steady state are less complicated to analyse because time marching properties change has been avoided.

SOR method proves t be the cost and time effective method when compared to Jacobi and Gauss seidel iterative methods.