Solving 2D Heat Conduction Equation Using Iterative Methods
Aim - To solve 2D heat conduction equation using iterative techniques for steady and transient analysis.
Governing equation is given by -
`((partial u)/(partial t)) = alpha*(((partial^2 u)/(partial x^2)) + ((partial^2 u)/(partial y^2)))`
In the above equation, temperature is denoted by `u` . For the above equation we need 1 initial condition and 4 boundary conditions i.e. 2 for x and 2 for y.
From the given data we can write -
`u(0,y,t) = 400 k` Left face
`u(1,y,t) = 800 k` Right face
`u(x,0,t) = 900 k` Bottom face
`u(x,1,t) = 600 k` Top face
Since the equation is solved numerically, solution is found at grid points only i.e. `u(x,y,t) = u_(i,j)^n` .
Initial temperature at all the interior grid points is assumed to be `298 k` i.e `u(x,y,t_o) = 298 k` .
Assumptions -
`alpha = 0.001`
Relaxation factor `omega = 1.5`
1 ) Steady state analysis : Steady state solution can be approached by two methods , one method is to remove the time derivative term in the equation and the second method is to run the simulation for larger time interval. 1st method is choosen for the analysis.
The governing equation reduces to -
`(partial^2 u)/(partial x^2) + (partial^2 u)/(partial y^2) = 0`
Using finite difference it is written as -
`(u_(i-1,j) - 2*u_(i,j) + u_(i+1,j))/(Delta x)^2 + (u_(i,j-1) - 2*u_(i,j) + u_(i,j+1))/(Delta y)^2 = 0`
On rearranging
`u_(i,j) = k1*(u_(i-1,j) + u_(i+1,j)) + k2*(u_(i,j-1) + u_(i,j+1))`
Where
`k1 = (Delta y)^2/(2*((Delta x)^2 + (Delta y)^2))` and `k2 = (Delta x)^2/(2*((Delta x)^2 + (Delta y)^2))`
Note - superscript k represents iteration index
Jacobi -
`u_(i,j)^(k+1) = k1*(u_(i-1,j)^k+u_(i+1,j)^k) + k2*(u_(i,j-1)^k+u_(i,j+1)^k)`
Gauss Seidel -
`u_(i,j)^(k+1) = k1*(u_(i-1,j)^(k+1)+u_(i+1,j)^k) + k2*(u_(i,j-1)^(k+1)+u_(i,j+1)^k)`
Successive over-relaxation -
`u_(i,j)^(k+1) = (1-omega)*u_(i,j)^k + omega*u_(GS)`
Code for steady state is pasted below. In order to run using different iterative techniques change the value of method.
% 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.5; % 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 using different iterative techniques are presented below.
Jacobi -
Gauss Seidel -
Successive over-relaxation -
2 ) Transient analysis : For discretizing, forward time central spacing i.e. FTCS has been used. There are 2 methods to conduct transient analysis. Simulation is carried out for 200 seconds for both the cases.
a ) Explicit method - In this method `u_(i,j)^(n+1)` is directly evaluated based on the previous time step values.
The discretized equation is given by -
`((u_(i,j)^(n+1) - u_(i,j)^n)/(Delta t)) = alpha * ( (u_(i-1,j)^n - 2*u_(i,j)^n + u_(i+1,j)^n)/(Delta x)^2 + (u_(i,j-1)^n - 2*u_(i,j)^n + u_(i,j+1)^n)/(Delta y)^2 )`
On rearranging
`u_(i,j)^(n+1) = k*u_(i,j)^n + k1*(u_(i-1,j)^n + u_(i+1,j)^n) + k2*(u_(i,j-1)^n+u_(i,j+1)^n)`
where
`k1 = alpha*(Delta t)/(Delta x)^2` , `k2 = alpha*(Delta t)/(Delta y)^2` and `k = 1-2*k1-2*k2`
Code for the explicit method 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 = 200; % simulating for 200 sec
dt = 1.e-2; % 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 : ',num2str(Et)])Plot using explicit method.
b ) Implicit method -
The discretized equation is written as -
`((u_(i,j)^(n+1) - u_(i,j)^n)/(Delta t)) = alpha * ( (u_(i-1,j)^(n+1) - 2*u_(i,j)^(n+1) + u_(i+1,j)^(n+1))/(Delta x)^2 + (u_(i,j-1)^(n+1) - 2*u_(i,j)^(n+1) + u_(i,j+1)^(n+1))/(Delta y)^2 )`
On rearranging
`u_(i,j)^(n+1) = (1/k)*(u_(i,j)^n + k1*(u_(i-1,j)^(n+1) + u_(i+1,j)^(n+1))+k2*( u_(i,j-1)^(n+1) + u_(i,j+1)^(n+1)))`
where
`k1 = alpha * (Delta t)/(Delta x)^2` , `k2 = alpha * (Delta t)/(Delta y)^2` and `k = 1+2*k1+2*k2`
Note - superscript k represents iteration index
Jacobi -
`u_(i,j)^(n+1,k+1) = (1/k)*(u_(i,j)^n + k1*(u_(i-1,j)^(n+1,k) + u_(i+1,j)^(n+1,k))+k2*( u_(i,j-1)^(n+1,k) + u_(i,j+1)^(n+1,k)))`
Gauss Seidel -
`u_(i,j)^(n+1,k+1) = (1/k)*(u_(i,j)^n + k1*(u_(i-1,j)^(n+1,k+1) + u_(i+1,j)^(n+1,k))+k2*( u_(i,j-1)^(n+1,k+1) + u_(i,j+1)^(n+1,k)))`
Successive over-relaxation -
`u_(i,j)^(n+1,k+1) = (1-omega)*u_(i,j)^(n+1,k) + omega*u_(GS)`
Code for the Implicit method is given below. In order to run using different iterative techniques change the value of method.
% 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 = 200; % simulating for 200 sec
dt = 1; % 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.5; % 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),', Time : ',num2str(Jt)])
elseif method == 2
title(['Implicit, Gauss Seidel, ','Iterations : ',num2str(iter),', Time : ',num2str(GSt)])
elseif method == 3
title(['Implicit, SOR, ','Iterations : ',num2str(iter),', Time : ',num2str(SORt)])
end
endPlots using different iterative techniques are presented below.
Jacobi -
Gauss Seidel -
Successive over-relaxation -
Conclusion - It has been found that solution for either steady or transient is not going to change with different iterative techniques. Only the rate of convergence is changing. It can be seen that rate of convergence for SOR > Gauss Seidel > Jacobi.