Steady and Unsteady (Implicit, Explicit) 2D heat conduction equation using Jacobi, Gauss-Seidel, and SOR(Successive over-relaxation) method
Steady State
Matlab code for steady-state
%Solving 2-D steady and unsteady heat conduction
clear all
close all
clc
Lx = 1; % Length in x-axis
Ly = 1; % Length in y-axis
nx = 50; % number of nodes in x-axis
ny = 50; % number of nodes in y-axis
x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
dx = Lx/(nx-1);
dy = Ly/(ny-1);
%Defining boundary conditions
T = 600*ones(nx,ny);
T(:,1) = 400; % Initial condition for left temperature
T(:,nx) = 800; % Initial condition for right temperature
T(1,:) = 600; % Initial condition for top temperature
T(ny,:) = 900; % Initial condition for botoom temperature
%Define values for corner
T(1,1) = (T(1,2)+T(2,1))/2; %TL
T(1,ny) = (T(1,ny-1)+T(2,ny))/2; %TR
T(nx,1) = (T(nx-1,1)+T(nx,2))/2; %TT
T(nx,ny) = (T(nx-1,ny)+T(nx,ny-1))/2; %TB
%Making a copy of T
meshgrid(x,y);
Told = T;
Tinitial = T;
%Defining K and W
k = 2*(((dx^2) + (dy^2)) /(( dx^2) * (dy^2 )));
w = 1.9; %relaxation factor
%Defining the solver which we want to run
iterative_solver = 1;
error_j = 9e9; %Defining error more than tolerance so that the while loop is executed
error_G = 9e9;
error_SOR = 9e9;
Tol = 1e-4;
%Steady state - Jacobi
if iterative_solver == 1
jacobi_iter = 1;
while (error_j > Tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1/k)*(((Told(i-1,j) + Told(i+1,j))/(dx^2)) + ((Told(i,j-1) + Told(i,j+1))/(dy^2)));
end
end
error_j = max(max(abs(Told - T)));
jacobi_iter = jacobi_iter +1;
Told = T;
figure(1)
contourf(x,y,T,'LineColor','none');
colormap('jet')
set(gca,'YDIR','reverse')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for Jacobi method(steady) = %d',jacobi_iter); %%d is used for integer
title(title_text)
end
end
%Steady state - GS
if iterative_solver == 2
GS_iter = 1;
while (error_G > Tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1/k)*(((T(i-1,j) + Told(i+1,j))/(dx^2)) + ((T(i,j-1) + Told(i,j+1))/(dy^2)));
end
end
error_G = max(max(abs(Told - T)));
GS_iter = GS_iter +1;
Told = T;
figure(2)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for GS method(steady) = %d',GS_iter); %%d is used for integer
title(title_text)
end
end
%Steady state - SOR
if iterative_solver == 3
SOR_iter = 1;
while (error_SOR > Tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (1-w)*Told(i,j) + w*(1/k)*(((T(i-1,j) + Told(i+1,j))/(dx^2)) + ((T(i,j-1) + Told(i,j+1))/(dy^2)));
end
end
error_SOR = max(max(abs(Told - T)));
SOR_iter = SOR_iter +1;
Told = T;
figure(3)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for SOR method(steady) = %d',SOR_iter);
title(title_text)
end
endPlots obtained from the above code are as follows:
Steady-state Jacobi
Steady-state Gauss Seidel
Steady-state SOR
Unsteady-state Implicit
Matlab code for Unsteady-state Implicit
%Solving 2-D steady and unsteady heat conduction
clear all
close all
clc
Lx = 1; % Length in x-axis
Ly = 1; % Length in y-axis
nx = 50; % number of nodes in x-axis
ny = 50; % number of nodes in y-axis
x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
dt = 50;
%Defining boundary conditions
T = 600*ones(nx,ny);
T(:,1) = 400; % Initial condition for left temperature
T(:,nx) = 800; % Initial condition for right temperature
T(1,:) = 600; % Initial condition for top temperature
T(ny,:) = 900; % Initial condition for botoom temperature
%Define values for corner
T(1,1) = (T(1,2)+T(2,1))/2; %TL
T(1,ny) = (T(1,ny-1)+T(2,ny))/2; %TR
T(nx,1) = (T(nx-1,1)+T(nx,2))/2; %TT
T(nx,ny) = (T(nx-1,ny)+T(nx,ny-1))/2; %TB
%Making a copy of T
[X,Y] = meshgrid(x,y);
Told = T;
Tinitial = T;
%Defining k1,k2,alpha and W
alpha = 2e-6;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
w = 4; %relaxation factor
%Defining the solver which we want to run
iterative_solver = 1;
error_j = 999; %Defining error more than tolerance so that the while loop is executed
error_G = 999;
error_SOR = 999;
Tol = 1e-4;
%Unsteady state - Jacobi (implicit)
if iterative_solver == 1
jacobi_iter = 1;
while (error_j > Tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(Told(i-1,j)+Told(i+1,j));
V=(Told(i,j-1)+Told(i,j+1));
T(i,j) = term1*Told(i,j) + term2*H + term3*V;
end
end
error_j = max(max(abs(Told - T)));
jacobi_iter = jacobi_iter +1;
Told = T;
figure(1)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for Jacobi method(unsteady-implicit) = %d',jacobi_iter); %%d is used for integer
title(title_text)
end
end
%Unsteady state - GS (implicit)
if iterative_solver == 2
GS_iter = 1;
while (error_G > Tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j) = term1*Told(i,j) + term2*H + term3*V;
end
end
error_G = max(max(abs(Told - T)));
GS_iter = GS_iter +1;
Told = T;
figure(2)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for GS method(unsteady-implicit) = %d',GS_iter); %%d is used for integer
title(title_text)
end
end
%Unsteady state - SOR (implicit)
if iterative_solver == 3
SOR_iter = 1;
while (error_SOR > Tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j) =(1-w)*Told(i,j)+ w*(term1*Told(i,j) + term2*H + term3*V);
end
end
error_SOR = max(max(abs(Told - T)));
SOR_iter = SOR_iter +1;
Told = T;
figure(3)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for SOR method(unsteady-implicit) = %d',SOR_iter);
title(title_text)
end
endPlots obtained from the above code are as follows:
Unsteady-state Jacobi Implicit
Unsteady-state Gauss-Seidel Implicit
Unsteady-state SOR Implicit
Unsteady-state Explicit
Matlab code for Unsteady-state Implicit
%Solving 2-D steady and unsteady heat conduction
clear all
close all
clc
Lx = 1; % Length in x-axis
Ly = 1; % Length in y-axis
nx = 50; % number of nodes in x-axis
ny = 50; % number of nodes in y-axis
x = linspace(0,Lx,nx);
y = linspace(0,Ly,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
dt = 50;
%Defining boundary conditions
T = 600*ones(nx,ny);
T(:,1) = 400; % Initial condition for left temperature
T(:,nx) = 800; % Initial condition for right temperature
T(1,:) = 600; % Initial condition for top temperature
T(ny,:) = 900; % Initial condition for botoom temperature
%Define values for corner
T(1,1) = (T(1,2)+T(2,1))/2; %TL
T(1,ny) = (T(1,ny-1)+T(2,ny))/2; %TR
T(nx,1) = (T(nx-1,1)+T(nx,2))/2; %TT
T(nx,ny) = (T(nx-1,ny)+T(nx,ny-1))/2; %TB
%Making a copy of T
[X,Y] = meshgrid(x,y);
Told = T;
Tinitial = T;
%Defining k1,k2,alpha and W
alpha = 2e-6;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
w = 1.1; %relaxation factor
nt = 1600;
%Defining the solver which we want to run
iterative_solver = 1;
%Unsteady state - Jacobi (explicit)
if iterative_solver == 1
jacobi_iter = 1;
for z= 1:nt
for i = 2:nx-1
for j = 2:ny-1
term1 = Told(i-1,j) - 2*Told(i,j) + Told(i+1,j);
term2 = Told(i,j-1) - 2*Told(i,j) + Told(i,j+1);
term3 = Told(i,j);
T(i,j) = term3 + k1*term1 + k2*term2;
end
end
jacobi_iter = jacobi_iter +1;
Told = T;
figure(1)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for Jacobi method(unsteady-explicit) = %d',jacobi_iter); %%d is used for integer
title(title_text)
colormap(jet)
end
end
%Unsteady state - GS (explicit)
if iterative_solver == 2
GS_iter = 1;
for z = 1:nt
for i = 2:nx-1
for j = 2:ny-1
term1 = T(i-1,j) - 2*T(i,j) + Told(i+1,j);
term2 = T(i,j-1) - 2*T(i,j) + Told(i,j+1);
term3 = T(i,j);
T(i,j) = term3 + k1*term1 + k2*term2;
end
end
GS_iter = GS_iter +1;
Told = T;
figure(2)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for GS method(unsteady-explicit) = %d',GS_iter); %%d is used for integer
title(title_text)
end
end
%Unsteady state - SOR (explicit)
if iterative_solver == 3
SOR_iter = 1;
for z = 1:nt
for i = 2:nx-1
for j = 2:ny-1
term1 = T(i-1,j) - 2*T(i,j) + Told(i+1,j);
term2 = T(i,j-1) - 2*T(i,j) + Told(i,j+1);
term3 = T(i,j);
T(i,j) = (1-w)*Told(i,j) + w*( term3 + k1*term1 + k2*term2);
end
end
SOR_iter = SOR_iter +1;
Told = T;
figure(3)
contourf(x,y,T,'LineColor','none');
set(gca,'YDIR','reverse')
colormap('jet')
xlabel('X')
ylabel('Y')
title_text = sprintf('Iteration number for SOR method(unsteady-explicit) = %d',SOR_iter);
title(title_text)
end
endPlots obtained from the above code are as follows:
Unsteady-state Jacobi Explicit
Unsteady-state Gauss-Seidel Explicit
Unsteady-state SOR Explicit
Shortforms and their meanings
J = Jacobi method
GS = Gauss-Seidel method
SOR = Successive over-relaxation
w = `omega` = relaxation factor
k1,k2 = CFL = stability constant
alpha = `alpha` = co-efficient of thermal expansion
iter = iteration
nt = number of time steps
nx, ny = number of nodes along x-axis and y-axis
Lx, Ly = length along x and yaxis
Conditions and constants used for the simulation in Matlab code
Lx = Ly = 1m
nx = ny = 50
dt = 50
nt = 1600
`alpha = 2*10^-6`
`omega = 1.9` ( for steady-state)
`omega = 4` (for unsteady-state - implicit)
`omega = 1.1` (for unsteady-state-explicit)
iterative_solver = 1 (for Jacobi)
iterative_solver = 2 (for GS)
iterative_solver = 3 (for SOR)
Explanation of the above-added Matlab codes:
Three different codes were written for steady-state, unsteady-state implicit, and unsteady-state explicit.
Each code gives output for Jacobi, Gauss-Seidel, and SOR method by varying the iterative_solver to 1,2. and 3 respectively.
First, all the initial and given conditions were defined as shown above.
Then the temperatures at boundaries and temperature at inner nodes were defined. The temperature at inner nodes was initially defined to be 600K.
while loop was used for iterative method to get the error below the mentioned/required tolerance value (0.0001 in our case). By varying the value of tolerance we can improve or lower the accuracy of our simulation. Unnecessarily increasing the accuracy without need will lead to higher simulation time and may require more system resources.
Two for loops were initiated for i (x values) and j ( y values) from the 2nd node to the final node.
Told was updated to T after these for loops or the function would run forever in the while loop without updating the Told value.
The figures were plotted using contourf command in Matlab and the details in graphs were managed as per the requirements.
contourf command is used for plotting a contour plot.
The number of iterations was counted by defining a jacobi_iter, GS_iter, and SOR_iter variables.
sprintf command was used to show the number of iterations in the title of the graph.
set(gca,'YDIR','reverse') command was used in the Matlab plot to reverse the y-axis in the contour plot.
Explanation of plots/results:
The following formula was discretized and used in the simulation,
Steady-state -
2D steady-state heat conduction equation
`(del^2T)/(delx^2) + (del^2T)/(dely^2) = 0`
after numerical discretization
`T_(i,j) = 1/k(( T_(i-1,j) + T_(i+1,j))/(trianglex^2)) + 1/k(( T_(i,j-1) + T_(i,j+1))/(trianglex^2))`
where, `k = 2((trianglex^2+triangley^2)/(trianglex^2*triangley^2))`
Unsteady-state (Implicit) -
2D unsteady-state heat conduction equation
`(delT)/(delt)=alpha*[(del^2T)/(delx^2) + (del^2T)/(dely^2) ]`
after numerical discretization
`T_p^(n+1) = 1/(1 + 2*k1 + 2*k2)*[T_p^n + k1(T_L + T_R)^(n+1) + k2(T_T + T_B)^(n+1)]`
where,
`T_p = T_(i,j)`
`T_L = T_(i-1,j)`
`T_R = T_(i+1,j)`
`T_T = T_(i,j-1)`
`T_B = T_(i,j+1)`
`k1 = (alpha*trianglet)/(trianglex^2)`
`k2 = (alpha*trianglet)/(triangley^2)`
Unsteady-state (Explicit) -
2D unsteady-state heat conduction equation
`(delT)/(delt)=alpha*[(del^2T)/(delx^2) + (del^2T)/(dely^2) ]`
after numerical discretization
`T_p^(n+1) = T_p^n + k1(T_L - 2*T_p + T_R)^n + k2(T_T - 2*T_p + T_B)^n`
where,
`T_p = T_(i,j)`
`T_L = T_(i-1,j)`
`T_R = T_(i+1,j)`
`T_T = T_(i,j-1)`
`T_B = T_(i,j+1)`
`k1 = (alpha*trianglet)/(trianglex^2)`
`k2 = (alpha*trianglet)/(triangley^2)`
Results obtained from the above solutions,
The number of iterations for each case is as follows
For steady-state:
| Type of Method | Number of iterations |
| Jacobi method | 3805 |
| Gauss-Seidel method | 2099 |
| Successive over-relaxation method | 171 |
For unsteady-state (implicit)
| Type of Method | Number of iterations |
| Jacobi method | 7061 |
| Gauss-Seidel method | 5572 |
| Successive over-relaxation method | 394 |
For unsteady-state (explicit)
| Type of Method | Number of iterations |
| Jacobi method | 1601 |
| Gauss-Seidel method | 1601 |
| Successive over-relaxation method | 1601 |
It can be observed from the above three cases that,
Jacobi iterations > Gauss-Seidel iterations > Successive over relaxation iterations