Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
The heat Equation is a PDE that describes how the distribution of quantities like heat changes with time in a solid medium as it flows from higher to lower places.
Heat equation:

α = diffusivity of medium
Types of heat schemes:
- Steady-state condition
- Transient state condition
Equations used for the codes:
- Steady-state equation

- unsteady-state equation

- unsteady-state Explicit equation:

JACOBI ITERATION:
close all
clear all
clc
%jacobi method
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;
error=9e9;
tolerance=1e-4;
%boundry condn
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 300*ones(nx,ny);
T(2:ny-1, 1) = T_L;
T(2:ny-1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;
%avg temp
T(1,1)= (T_T+T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1)= (T_L+T_B)/2;
T_initial=T;
T_old=T;
%calculating jacobi
iterative_solver = 1;
if iterative_solver ==1
jacobi_iteration=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
T(i,j) = 0.25*(T_old(i-1,j) + T_old(i+1,j)+T_old(i, j-1)+T_old(i,j+1));
end
end
error= max(max(abs(T_old-T)));
T_old=T;
jacobi_iteration=jacobi_iteration+1;
end
end
%plotting
figure(1);
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('Xaxis')
ylabel('yaxis')
title(sprintf('no:of jacobi=%d',jacobi_iteration));
pause(0.03);

GAUSS SIDEL IMPLICIT METHOD:
close all
clear all
clc
%gauss siedel method
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;
error=9e9;
tolerance=1e-4;
%boundry condn
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 300*ones(nx,ny);
T(2:ny-1, 1) = T_L;
T(2:ny-1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;
%avg temp
T(1,1)= (T_T+T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1)= (T_L+T_B)/2;
T_initial=T;
T_old=T;
%calculating gauss sidel method
iterative_solver = 2;
if iterative_solver ==2
Gaussseidel_iteration=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
T(i,j) = 0.25*(T(i-1,j) + T_old(i+1,j)+T(i, j-1)+T_old(i,j+1));
end
end
error= max(max(abs(T_old-T)));
T_old=T;
Gaussseidel_iteration= Gaussseidel_iteration+1;
end
end
%plotting
figure(1);
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('Xaxis')
ylabel('yaxis')
title(sprintf('no:of Gaussseidel=%d', Gaussseidel_iteration));
pause(0.03);

succesive_relaxation implicit:
close all
clear all
clc
%succesive_relaxation method
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;
error=9e9;
tolerance=1e-4;
omega=1.2;
%boundry condn
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 300*ones(nx,ny);
T(2:ny-1, 1) = T_L;
T(2:ny-1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;
%avg temp
T(1,1)= (T_T+T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1)= (T_L+T_B)/2;
T_initial=T;
T_old=T;
%calculating succesive_relaxation method
iterative_solver = 3;
if iterative_solver ==3
succesive_relaxation_iteration=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
T(i,j) = omega*(0.25*(T(i-1,j) + T_old(i+1,j)+T(i, j-1)+T_old(i,j+1)))+(1-omega)*T_old(i,j);
end
end
error= max(max(abs(T_old-T)));
T_old=T;
succesive_relaxation_iteration= succesive_relaxation_iteration+1;
end
end
%plotting
figure(1);
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('Xaxis')
ylabel('yaxis')
title(sprintf('no:of succesive_relaxation=%d', succesive_relaxation_iteration));
pause(0.03);

UNSTEADY:
Jacobi method:
close all
clear all
clc
%jacobi method
nx=10;
ny=nx;
nt=1000;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;
error=9e9;
tolerance=1e-4;
dt=1e-3;
%boundry condn
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 300*ones(nx,ny);
T(2:ny-1, 1) = T_L;
T(2:ny-1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;
%avg temp
T(1,1)= (T_T+T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1)= (T_L+T_B)/2;
T_intial=T;
T_old=T;
%calculating jacobi
k1= 1.1*(dt/(dx^2));
k2= 1.1*(dt/(dy^2));
jacobi_iteration=1;
for k=1:nt
error=9e9;
while(error > tolerance)
for i=2:nx-1
for j=2:ny-1
Term_1= 1/(1+(2*k1)+(2*k2));
Term_2= k1*Term_1;
Term_3= k2*Term_1;
H= (T_old(i-1,j))+(T_old(i+1,j));
V= (T_old(i,j+1))+(T_old(i,j-1));
T(i,j)= (T_intial(i,j)*Term_1)+(H*Term_2)+(V*Term_3);
end
end
error = max(max(abs(T_old -T)));
T_old=T;
jacobi_iteration=jacobi_iteration+1;
end
T_intial=T;
%plotting
figure(1);
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('Xaxis')
ylabel('yaxis')
title(sprintf('no:of unsteady jacobi=%d',jacobi_iteration));
pause(0.03);
end

gauss seidel method:
close all
clear all
clc
%gausssidel method
nx=10;
ny=nx;
nt=1400;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;
error=9e9;
tolerance=1e-4;
dt=1e-3;
%boundry condn
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 300*ones(nx,ny);
T(2:ny-1, 1) = T_L;
T(2:ny-1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;
%avg temp
T(1,1)= (T_T+T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1)= (T_L+T_B)/2;
T_intial=T;
T_old=T;
%calculating jacobi
k1= 1.1*(dt/(dx^2));
k2= 1.1*(dt/(dy^2));
gausssidel_iteration=1;
for k=1:nt
error=9e9;
while(error > tolerance)
for i=2:nx-1
for j=2:ny-1
Term_1= 1/(1+(2*k1)+(2*k2));
Term_2= k1*Term_1;
Term_3= k2*Term_1;
H= (T(i-1,j))+(T_old(i+1,j));
V= (T(i,j+1))+(T_old(i,j-1));
T(i,j)= (T_intial(i,j)*Term_1)+(H*Term_2)+(V*Term_3);
end
end
error = max(max(abs(T_old -T)));
T_old=T;
gausssidel_iteration=gausssidel_iteration+1;
end
T_intial=T;
%plotting
figure(1);
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('Xaxis')
ylabel('yaxis')
title(sprintf('no:of unsteady gausssidel=%d',gausssidel_iteration));
pause(0.03);
end

succesive_relaxation(SOR):
close all
clear all
clc
%succesive_relaxation method
nx=10;
ny=nx;
nt=1400;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;
error=9e9;
tolerance=1e-4;
dt=1e-3;
omega=1.1;
%boundry condn
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 300*ones(nx,ny);
T(2:ny-1, 1) = T_L;
T(2:ny-1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;
%avg temp
T(1,1)= (T_T+T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1)= (T_L+T_B)/2;
T_intial=T;
T_old=T;
T_end=T;
%calculating jacobi
k1= 1.1*(dt/(dx^2));
k2= 1.1*(dt/(dy^2));
succesive_relaxation_iteration=1;
for k=1:nt
error=9e9;
while(error > tolerance)
for i=2:nx-1
for j=2:ny-1
Term_1= 1/(1+(2*k1)+(2*k2));
Term_2= k1*Term_1;
Term_3= k2*Term_1;
H= (T(i-1,j))+(T_old(i+1,j));
V= (T(i,j+1))+(T_old(i,j-1));
T(i,j)= (T_intial(i,j)*Term_1)+(H*Term_2)+(V*Term_3);
T_end(i,j)= ((1-omega)*T_old(i,j))+(T(i,j)*omega);
end
end
error = max(max(abs(T_old -T)));
T_old=T_end;
succesive_relaxation_iteration=succesive_relaxation_iteration+1;
end
T_intial=T_end;
%plotting
figure(1);
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('Xaxis')
ylabel('yaxis')
title(sprintf('no:of unsteady succesive_relaxation=%d',succesive_relaxation_iteration));
pause(0.03);
end

EXPLICIT method - unsteady/Transient condition:
close all
clear all
clc
%succesive_relaxation method(EXPLICIT)
nx=10;
ny=nx;
nt=1400;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=dx;
error=9e9;
tolerance=1e-4;
dt=1e-3;
omega=1.1;
%boundry condn
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = 300*ones(nx,ny);
T(2:ny-1, 1) = T_L;
T(2:ny-1, nx) = T_R;
T(1, 2:nx - 1) = T_T;
T(ny, 2:nx - 1) = T_B;
%avg temp
T(1,1)= (T_T+T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1)= (T_L+T_B)/2;
T_intial=T;
T_old=T;
%calculating jacobi
iterative_solver=1;
k1= 1.1*(dt/(dx^2));
k2= 1.1*(dt/(dy^2));
if iterative_solver==1
steady_explicit=1;
for k=1:nt
error=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
Term_1=T_old(i,j);
Term_2=k1*(T_old(i+1,j)-2*T_old(i,j)+(T_old(i-1,j)));
Term_3=k2*(T_old(i,j+1)-2*T_old(i,j)+T_old(i,j-1));
T(i,j)=Term_1+Term_2+Term_3;
end
end
error=max(max(abs(T_old-T)));
T_old=T;
steady_explicit= steady_explicit+1;
end
end
end
%plotting
figure(1);
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('Xaxis')
ylabel('yaxis')
title(sprintf('no:of unsteady iteration(EXPLICIT)=%d', steady_explicit));
pause(0.03);

RESULTS:
| TYPE OF SCHEME | METHODS | ITERATIONS |
| Implicit steady-state analysis | Jacobi | 208 |
| Implicit steady-state analysis |
Gauss seidel |
112 |
| Implicit steady-state analysis | SOR | 46 |
| Implicit Transient-state analysis | Jacobi | 3880 |
| Implicit Transient-state analysis | SOR | 3290 |
| Implicit Transient-state analysis | Gauss seidel | 3612 |
- Jacobi > Gauss sidel>SOR according to number of iterations
- Steady-state takes less time to compute and fewer iterations