Challenge

Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

JJaya Prakash2 Sept 2022

AIM:

To solve the 2D heat conduction equation for steady and unsteady state by using the point iterative technique to using method of t Jacobi,

gauss seidel, and successive over-relaxation for both implicit and explict method.

Explanation:

Heat Conduction:

Heat conduction is the flow of thermal energy through a substance from a region of higher temperature

to a region of lower temperature,which occurs by atomic or molecular interactions.

Given Inputs and Boundary conditions:

Domain is assumed to be a shape of a square of length 1m.

nodes in x and y direction be 20. (n = 20)

Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
Initial temperature inside the domain = 298 K
Top - Left corner = 500 K
Top - Right corner = 700 K
Bottom - Left corner = 650 K
Bottom - Right corner = 850 K
Absolute error criteria = 1e-4

Temperature distribution as per given data:

Equations:

1. Steady state equation

2. Unsteady state equation

1. Steady state equation:

Discretizing the above equation using the Central difference scheme (CDS), we get

Simplifying the above eqation,

Steady state - Jacobi Method:

Steady state - Gauss-Seidel Method:

Steady state - Successive over Relaxation Method:

2. Unsteady state equation - Explicit Method:

Unsteady State heat conduction equation is given as,

Discretizing the above time term using the forward differencing scheme (FDS) and for space term using central difference scheme (CDS).

we get,

Unsteady state Implicit Method:

Unsteady state Implicit - Jacobi- Method:

Unsteady state Implicit - Gauss-Seidel Method:

Unsteady state Implicit - SOR Method:

Steady state MATLAB CODE:

clear all

close all

clc

%input array

n = 10;

x = linspace(0,1,n);

y = linspace(0,1,n);

dx = x(2) - x(1);

dy = dx;

T = 300*ones(n,n);

% Given boundary condtion apply.

T(1,:) = 600; % Top Boundary Temp.... in Kelvin

T(end,:) = 900; % Bottom Boundary Temp..... in Kelvin

T(:,1) = 400; % Left Boundary Temp..... in Kelvin

T(:,end) = 800; % Right Boundary Temp.... in Kelvin

T(1,1) = 500; % Top - Left corner Temp.... in Kelvin

T(1,end) = 700; % Top - Right corner Temp.... in Kelvin

T(end,1) = 650; % Bottom - Left corner Temp.... in Kelvin

T(end,end) = 850; % Bottom - Right corner Temp.... in Kelvin

%copying

Told_1 = T;

Told_2 = T;

Told_3 = T;

T_J = T;

T_G = T;

T_S = T;

%apply rand - error and iter values

jacobi_error = 10;

gauss_sidel_error = 10;

SOR_error = 10;

error_req = 1e-4;

jacobi_iter = 0;

gauss_sidel_iter = 0;

SOR_iter = 0;

w = 2/(1+sin(pi*dx)); % SOR Method

% appy condtion for jacobi

for method = 1:3

tic

if method == 1

while jacobi_error > error_req

for i = 2:(n-1)

for j = 2:(n-1)

T_J(i,j) = 0.25*(Told_1(i+1,j)+Told_1(i-1,j)+Told_1(i,j+1)+Told_1(i,j-1));

end

jacobi_iter = jacobi_iter + 1;

jacobi_error = max(max(abs(Told_1-T_J)));

Told_1 = T_J;

time_jc = toc;

end

figure(1)

contourf(x,y,T_J,'Showtext','on')

colormap(jet)

colorbar

xlabel('X Axis')

ylabel('Y Axis')

title(['Result=' num2str(jacobi_iter) ' Jacobi iterations'])

pause(0.3)

set(gca,'ydir','reverse')

% appy condtion for gauss_sidel

tic

if method == 2

while gauss_sidel_error > error_req

for i = 2:(n-1)

for j = 2:(n-1)

T_G(i,j) = 0.25*(Told_2(i+1,j)+T_G(i-1,j)+Told_2(i,j+1)+T_G(i,j-1));

end

gauss_sidel_iter = gauss_sidel_iter + 1;

gauss_sidel_error = max(max(abs(Told_2-T_G)));

Told_2 = T_G;

time_gs = toc;

end

figure(2)

contourf(x,y,T_G,'Showtext','on')

colormap(jet)

colorbar

xlabel('X Axis')

ylabel('Y Axis')

title(['Result =' num2str(gauss_sidel_iter) ' Gauss Sidel iterations'])

pause(0.3)

set(gca,'ydir','reverse')

% appy condtion for SOR_error

tic

if method == 3

while SOR_error > error_req

for i = 2:(n-1)

for j = 2:(n-1)

T_S(i,j) = (1-w)*Told_3(i,j) + w*(0.25*(Told_3(i+1,j)+T_S(i-1,j)+Told_3(i,j+1)+T_S(i,j-1)));

end

SOR_iter = SOR_iter + 1;

SOR_error = max(max(abs(Told_3-T_S)));

Told_3 = T_S;

time_sor = toc;

end

figure(3)

contourf(x,y,T_S,'Showtext','on')

colormap(jet)

colorbar

xlabel('X Axis')

ylabel('Y Axis')

title(['Result=' num2str(SOR_iter) ' SOR iterations'])

pause(0.3)

set(gca,'ydir','reverse')

end

Un Steady state MATLAB CODE:

Explicit:

clc

clear all

close all

n = 10;

x = linspace(0,1,n);

y = linspace(0,1,n);

dx = x(2) - x(1);

dy = dx;

dt = 2e-4;

a = 1.2;

k1 = a*dt/(dx^2);

k2 = a*dt/(dy^2);

CFL = k1+k2;

T = 298*ones(n,n);

T(1,:) = 600; % Top Boundary Temp. in Kelvin

T(end,:) = 900; % Bottom Boundary Temp. in Kelvin

T(:,1) = 400; % Left Boundary Temp. in Kelvin

T(:,end) = 800; % Right Boundary Temp. in Kelvin

T(1,1) = 500; % Top - Left corner Temp. in Kelvin

T(1,end) = 700; % Top - Right corner Temp. in Kelvin

T(end,1) = 650; % Bottom - Left corner Temp. in Kelvin

T(end,end) = 850; % Bottom - Right corner Temp. in Kelvin

T_old = T;

T_US = T;

error_req = 1e-4;

error_mag = 1;

Ex_US_iter = 0;

tic

for nt = 1:1000

for i = 2:(n-1)

for j = 2:(n-1)

term1 = T_old(i,j);

term2 = k1*(T_old(i-1,j)-2*T_old(i,j)+T_old(i+1,j));

term3 = k2*(T_old(i,j-1)-2*T-old(i,j)+T_old(i,j+1));

T_US (i,j) = term1 + term2 + term3;

end

Ex_US_iter = Ex_US_iter+1;

error_mag = max(max(abs(T_old - T_US)));

T_old = T_US;

time_req = toc;

end

figure(1)

contourf(x,y,T_US,'Showtext','on')

colormap(jet)

colorbar

xlabel('X Axis')

ylabel('Y Axis')

title(['Result of Explicit Method for unsteady state=' num2str( Ex_US_iter) ' Iteration'])

Unsteady state Implicit Jacobi:

clear all
close all
clc

n = 10;

x = linspace(0,1,n);

y = linspace(0,1,n);

dx = x(2) - x(1);
dy = dx;
nt = 1000;
dt = 0.02;
a = 1200;

k1 = (a*dt)/(dx^2);
k2 = (a*dt)/(dy^2);

CFL = k1+k2;

T = 300*ones(n,n);
T(1,:) = 600; % Top Boundary Temp. in Kelvin
T(end,:) = 900; % Bottom Boundary Temp. in Kelvin
T(:,1) = 400; % Left Boundary Temp. in Kelvin
T(:,end) = 800; % Right Boundary Temp. in Kelvin
T(1,1) = 500; % Top - Left corner Temp. in Kelvin
T(1,end) = 700; % Top - Right corner Temp. in Kelvin
T(end,1) = 650; % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850; % Bottom - Right corner Temp. in Kelvin
jacobi_iter = 0;
jacobi_error = 9e9;

tol = 1e-4;

Told = T;

T_previous = Told;

term1 = (1+2*k1+2*k2)^(-1);
term2 = (k1*term1);
term3 = (k2*term1);

tic
for k = 1:nt
while jacobi_error > tol
for i = 2:(n-1)
for j = 2:(n-1)
H = Told(i-1,j)+Told(i+1,j);
V = Told(i,j-1)+Told(i,j+1);
T(i,j) = (T_previous(i,j)*term1)+(H*term2)+(V*term3);
end
end

jacobi_iter = jacobi_iter + 1;

jacobi_error = max(max(abs(Told-T)));

Told = T;

time_jc = toc;

end

T_previous = T;

end

figure(1)
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
xlabel('X Axis')
ylabel('Y Axis')
title(['Unsteady state Implicit Jacobi = ' num2str( jacobi_iter) ' Iteration'])

UNSTEADY STATE IMPLICIT GAUSS SEIDEL:

clear all
close all
clc

n = 10;

x = linspace(0,1,n);
y = linspace(0,1,n);

dx = x(2) - x(1);

dy = dx;

nt = 1000;

dt = 0.02;

a = 1200;

k1 = (a*dt)/(dx^2);
k2 = (a*dt)/(dy^2);

CFL = k1+k2;

T = 300*ones(n,n);

T(1,:) = 600; % Top Boundary Temp. in Kelvin
T(end,:) = 900; % Bottom Boundary Temp. in Kelvin
T(:,1) = 400; % Left Boundary Temp. in Kelvin
T(:,end) = 800; % Right Boundary Temp. in Kelvin
T(1,1) = 500; % Top - Left corner Temp. in Kelvin
T(1,end) = 700; % Top - Right corner Temp. in Kelvin
T(end,1) = 650; % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850; % Bottom - Right corner Temp. in Kelvin

gauss_iter = 0;

gauss_error = 9e9;

tol = 1e-4;

Told = T;

T_previous = Told;

term1 = (1+2*k1+2*k2)^(-1);

term2 = (k1*term1);

term3 = (k2*term1);

tic
for k = 1:nt
while gauss_error > tol
for i = 2:(n-1)
for j = 2:(n-1)
H = T(i-1,j)+T(i+1,j);
V = T(i,j-1)+T(i,j+1);
T(i,j) = (T_previous(i,j)*term1)+(H*term2)+(V*term3);
end
end

gauss_iter = gauss_iter + 1;

gauss_error = max(max(abs(Told-T)));

Told = T;
time_gs = toc;
end
T_previous = Told;
end
figure(1)
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
xlabel('X Axis')
ylabel('Y Axis')
title(['UNSTEADY STATE IMPLICIT GAUSS SEIDEL = ' num2str( gauss_iter) ' Iteration'])

UNSTEADY STATE IMPLICIT SOR:

clear all
close all
clc
n = 10;

x = linspace(0,1,n);
y = linspace(0,1,n);
dx = x(2) - x(1);
dy = dx;
nt = 1000;
dt = 0.02;
a = 1200;

k1 = (a*dt)/(dx^2);
k2 = (a*dt)/(dy^2);
CFL = k1+k2;
T = 300*ones(n,n);
T(1,:) = 600; % Top Boundary Temp. in Kelvin
T(end,:) = 900; % Bottom Boundary Temp. in Kelvin
T(:,1) = 400; % Left Boundary Temp. in Kelvin
T(:,end) = 800; % Right Boundary Temp. in Kelvin
T(1,1) = 500; % Top - Left corner Temp. in Kelvin
T(1,end) = 700; % Top - Right corner Temp. in Kelvin
T(end,1) = 650; % Bottom - Left corner Temp. in Kelvin
T(end,end) = 850; % Bottom - Right corner Temp. in Kelvin
SOR_iter = 0;
SOR_error = 9e9;
tol = 1e-4;
Told = T;
T_previous = Told;
term1 = (1+2*k1+2*k2)^(-1);
term2 = (k1*term1);
term3 = (k2*term1);
w = 2/(1+sin(pi*dx)); % Relaxation parameter for SOR Method
%w = 1.2;
%%
tic
for k = 1:nt
while SOR_error > tol
for i = 2:(n-1)
for j = 2:(n-1)
H = T(i-1,j)+T(i+1,j);
V = T(i,j-1)+T(i,j+1);
T(i,j) = w*((T_previous(i,j)*term1)+(H*term2)+(V*term3)) - (w-1)*T(i,j);
end
end
SOR_iter = SOR_iter + 1;
SOR_error = max(max(abs(Told-T)));
Told = T;
time_sor = toc;
end
T_previous = T;
end
figure(1)
contourf(x,y,T,'Showtext','on')
colormap(jet)
colorbar
xlabel('X Axis')
ylabel('Y Axis')
title(['UNSTEADY STATE FOR IMPLICIT SOR = ' num2str( SOR_iter) ' Iteration'])

Result:

Steady State:

Result for Jacobi:

Result for GS:

Result for SOR:

UNSTEADY STATE EXPLICIT:

UNSTEADY STATE IMPLICIT:

Unsteady state implicit Jacobi:

Unsteady state implicit Gauss – seidel:

Unsteady state implicit SOR:

EXPLANATION:

Here, We find the Iteration and errors from the each results

Sr. No.	Method	Iterations	Errors
1.	Steady State (Jacobi)	207	9.48e-05
2.	Steady State (Gauss-Seidel)	111	9.586e-05
3.	Steady State (SOR)	29	8.57e-05
4.	Unsteady state Explicit	-	0
5.	Unsteady state Implicit (Jacobi)	206	9.83e-05
6.	Unsteady state Implicit (Gauss-Seidel)	111	9.32e-05
7.	Unsteady state Implicit (SOR)	29	7.71e-05

From the table, we can conclude that the Jacobi method took lots of iteration to converge and the Successive over-relaxation method took the least, whereas the Gauss – seidel iteration lies between the Jacobi and SOR.

CONCLUSION:

Thus, in MATLAB, the 2D linear heat conduction equation is been solved by using the point iterative technique using the methods Jacobi, Gauss – seidel, Successive over Relaxation (SOR) for both implicit and explicit scheme.

The number of iteration and error values for each method is calculated and displayed in the above table.

The temperature distribution plots are plotted for each method.

Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.