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.