Comparing Iterative Methods By Solving 2D Heat Conduction Equation
Aim - To compare the fastest steady state simulation vs fastest transient simulation.
Steady state -
Code for steady state is given below.
% steady state solution to 2D heat conduction equation in a unit square slab
clear; close all; clc
nx = 101; % No of Grid points in x-dir
dx = 1/(nx-1);
ny = 101; % No of Grid points in y-dir
dy = 1/(ny-1);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,-y);
u = ones(nx,ny)*298; % Initialising the solution
% Boundary conditions
u(:,1) = 400; % Left wall
u(:,nx) = 800; % Right wall
u(1,:) = 600; % Top wall
u(ny,:) = 900; % Bottom wall
% copying the u
uges = u; % guess values for iteration
% constants in the discritized equation
k1 = dy^2/(2*(dx^2+dy^2));
k2 = dx^2/(2*(dx^2+dy^2));
% Note change the value of variable 'method' to solve using different
% iterative techniques
% 1-jacobi, 2-gaussSeidel, 3-SOR
tol = 1.e-4; % tolerance
err = 1; % initializing error
iter = 0; % initializing iteration number
s = 1.3; % relaxation factor
method = 1;
%% Jacobi
tic;
if method == 1
% convergence loop
while err > tol
% nodal loop
for i = 2:nx-1
for j = 2:ny-1
u(i,j) = k1*(uges(i-1,j)+uges(i+1,j))+k2*(uges(i,j-1)+uges(i,j+1));
end
end
err = max(max(abs(u-uges)));
uges = u; % updating for next iteration
iter = iter + 1;
end
end
Jt = toc;
%% Gauss Seidel
tic;
if method == 2
% convergence loop
while err > tol
% nodal loop
for i = 2:nx-1
for j = 2:ny-1
u(i,j) = k1*(u(i-1,j)+uges(i+1,j))+k2*(u(i,j-1)+uges(i,j+1));
end
end
err = max(max(abs(u-uges)));
uges = u; % updating for next iteration
iter = iter + 1;
end
end
GSt = toc;
%% SOR
tic;
if method == 3
% convergence loop
while err > tol
% nodal loop
for i = 2:nx-1
for j = 2:ny-1
u(i,j) = (1-s) * uges(i,j) + s * ( k1*(u(i-1,j)+uges(i+1,j))+k2*(u(i,j-1)+uges(i,j+1)) );
end
end
err = max(max(abs(u-uges)));
uges = u; % updating for next iteration
iter = iter + 1;
end
end
SORt = toc;
%% plotting
contourf(xx,yy,u,'ShowText','on')
yticklabels({'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'})
for i = 1:3
if method == 1
title(['Steady state, Jacobi, ','Iterations : ',num2str(iter),' Time : ',num2str(Jt)])
elseif method == 2
title(['Steady state, Gauss Seidel, ','Iterations : ',num2str(iter),' Time : ',num2str(GSt)])
elseif method == 3
title(['Steady state, SOR, ','Iterations : ',num2str(iter),' Time : ',num2str(SORt)])
end
endplots for steady state are given below -
Jacobi -
Gauss Seidel -
Successive over relaxation -
From the above figures it can be shown that Successive over relaxation has fastest convegence rate followed by Gauss seidel followed by Jacobi iteration technique. This is because SOR uses a relaxation factor which basically lies between 1-2. The relaxation factor improves the guess value hence faster convergence. If the relaxation factor is 1 then SOR reduces to Gauss Seidel. It can be seen that SOR takes 5097 iterations where as Gauss Seidel takes 8831 iterations, and Jacobi takes higher number of iterations 16234.
Transient - To compare between explicit and implicit method same time step and time interval has been used.
Explicit :
Code for explicit is given below -
% Explicit method to 2D transient heat conduction equation in a unit square slab
clear; close all; clc
nx = 101; % No of Grid points in x-dir
dx = 1/(nx-1);
ny = 101; % No of Grid points in y-dir
dy = 1/(ny-1);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,-y);
tm = 50; % simulating for 50 sec
dt = 0.025; % time step size
nt = tm/dt; % number of time steps
u = ones(nx,ny)*298; % initialising the solution
% Boundary conditions
u(:,1) = 400; % Left wall
u(:,nx) = 800; % Right wall
u(1,:) = 600; % Top wall
u(ny,:) = 900; % Bottom wall
% copying the u
ut = u; % for time increment
% constants in the discritized equation
a = 1.e-3; % Thermal diffusivity
k1 = a*dt/dx^2;
k2 = a*dt/dy^2;
k = 1-2*k1-2*k2;
% Time loop
tic;
for t = 1:nt
% Nodal loop
for i = 2:nx-1
for j = 2:ny-1
u(i,j) = k*ut(i,j) + k1*(ut(i-1,j)+ut(i+1,j)) + k2*(ut(i,j-1)+ut(i,j+1));
end
end
ut = u; % Updating for next time step
end
Et = toc;
%% plotting
contourf(xx,yy,u,'ShowText','on')
yticklabels({'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'})
title(['Explicit, ','TIme step : ',num2str(dt),', Run Time : ',num2str(Et)])plot using explicit method -
Implicit :
Code for implicit method is given below -
% Implicit method to 2D transient heat conduction equation in a unit square slab
clear; close all; clc
nx = 101; % No of Grid points in x-dir
dx = 1/(nx-1);
ny = 101; % No of Grid points in y-dir
dy = 1/(ny-1);
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,-y);
tm = 50; % simulating for 50 sec
dt = 0.025; % time step size
nt = tm/dt; % number of time steps
u = ones(nx,ny)*298; % initialising the solution
% Boundary conditions
u(:,1) = 400; % Left wall
u(:,nx) = 800; % Right wall
u(1,:) = 600; % Top wall
u(ny,:) = 900; % Bottom wall
% copying the u
uges = u; % guess values for iteration
ut = u; % for time increment
% constants in the discritized equation
a = 1.e-3; % Thermal diffusivity
k1 = a*dt/dx^2;
k2 = a*dt/dy^2;
k = 1+2*k1+2*k2;
% Note change the value of variable 'method' to solve using different
% iterative techniques
% 1-jacobi, 2-gaussSeidel, 3-SOR
method = 1;
s = 1.3; % Relaxation factor
tol = 1.e-4; % tolerance
iter = 0; % initializing iteration number
%% Jacobi
tic;
if method == 1
% Time loop
for t = 1:nt
err = 1; % initializing error at each time step
% convergence loop
while err > tol
% Nodal loop
for i = 2:nx-1
for j = 2:ny-1
u(i,j) = (1/k)*( ut(i,j) + k1*(uges(i-1,j)+uges(i+1,j)) + k2*(uges(i,j-1)+uges(i,j+1)) );
end
end
err = max(max(abs(u-uges)));
uges = u; % updating for next iteration
iter = iter + 1;
end
ut = u; % updating for next time step
end
end
Jt = toc;
%% Gauss Seidel
tic;
if method == 2
% Time loop
for t = 1:nt
err = 1; % initializing error at each time step
% convergence loop
while err > tol
% Nodal loop
for i = 2:nx-1
for j = 2:ny-1
u(i,j) = (1/k)*( ut(i,j) + k1*(u(i-1,j)+uges(i+1,j)) + k2*(u(i,j-1)+uges(i,j+1)) );
end
end
err = max(max(abs(u-uges)));
uges = u; % updating for next iteration
iter = iter + 1;
end
ut = u; % updating for next time step
end
end
GSt = toc;
%% SOR
tic;
if method == 3
% Time loop
for t = 1:nt
err = 1; % initializing error at each time step
% convergence loop
while err > tol
% nodal loop
for i = 2:nx-1
for j = 2:ny-1
u(i,j) = (1-s)*uges(i,j) + s *((1/k)*( ut(i,j) + k1*(u(i-1,j)+uges(i+1,j)) + k2*(u(i,j-1)+uges(i,j+1))));
end
end
err = max(max(abs(u-uges)));
uges = u; % updating for next iteration
iter = iter + 1;
end
ut = u; % updating for next time step
end
end
SORt = toc;
%% plotting
contourf(xx,yy,u,'ShowText','on')
yticklabels({'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'})
for i = 1:3
if method == 1
title(['Implicit, Jacobi, ','Iterations : ',num2str(iter),', Run Time : ',num2str(Jt)])
elseif method == 2
title(['Implicit, Gauss Seidel, ','Iterations : ',num2str(iter),', Run Time : ',num2str(GSt)])
elseif method == 3
title(['Implicit, SOR, ','Iterations : ',num2str(iter),', Run Time : ',num2str(SORt)])
end
endPlots using different iterative techniques for implicit method are given below -
Jacobi -
Gauss Seidel -
Succesive over relaxation -
In transient simulation explicit method is more faster than implicit method. This is because explicit method calculates the variable of interest based on the previous time step values. Hence there is no need of convergence loop which reduces the computation time.
In implicit method at each time step a guess solution is assumed and iterated untill the guess solution converges to a predefined tolerance value. If only implicit method is considered than SOR converges faster than Gauss Seidel followed by Jacobi.