Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Solving the steady and unsteady 2D heat conduction problem
Objective:
Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation
Your absolute error criteria is 1e-4
Explaination:
Steady state analysis:
`(delx^2)/(delt^2)+(dely^2)/(delt^2)`
Unsteady state/ Transient State Analysis:
Initial value for inbetween nodes
Initial value for inbetween nodes= 675K
Inputs:
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
Source code:
Steady state:
`(del^2T)/(del x^2) +(del^2T)/(del y^2) = 0`
% Code to solve a second order 2 D PDE
% d^u/dx^2 = 0
% boundary conditions:
% x=0 u=5
%y=0 u=8
clear all
close all
clc
% inputs
nx = 10 % number of grid points along x axis
ny = 10 % number of grid points along y axis
t=ones(nx)*300
t(:,1) = 400; % Boundary condition - Left
t(:,ny) = 800; % Boundary condition - Right
t(1,:) = 600; % Boundary condition - Top
t(nx,:) = 900; % Boundary condition - Bottom
t(1,1)= 500
t(1,ny)= 650
t(nx,1)= 700
t(nx,ny)= 850
error= 9e+9 % random value higher than the tolerance
tol=1e-4; % Tolerance limit
omega = 1.2; %SOR omega value
iterative_solver = 1;
% Creating a copy
told=t;
t_initial= t;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx= x(2)-x(1);
dy= y(2)-y(1);
k=0.25
%-----------------------------------------------------------------------
% Jacobi method
if iterative_solver==1
jacobi_counter =1;
while(error> tol)
for i= 2:nx-1 % Space loop along x axis
for j= 2:ny-1 % Space loop along y axis
t(i,j)=k*(told(i-1,j)+told(i+1,j)+told(i,j-1)+told(i,j+1));
end
end
error = max(max(abs(told-t)));
told=t; % Updating values
jacobi_counter =jacobi_counter+1;
end
contourf(t)
xlabel('X axis')
ylabel('U')
colorbar
title_text = sprintf('Iteration number = %d',jacobi_counter);
title(title_text)
pause(0.03)
end
%-----------------------------------------------------------------------
% GS method
if iterative_solver==2
gs_counter =1;
while(error> tol)
for i= 2:nx-1 % Space loop along x axis
for j= 2:ny-1 % Space loop along y axis
t(i,j)=k*(t(i-1,j)+told(i+1,j)+t(i,j-1)+told(i,j+1));
end
end
error = max(max(abs(told-t)));
told=t; % Updating values
gs_counter =gs_counter+1;
end
contourf(t)
xlabel('X axis')
ylabel('U')
title_text = sprintf('Iteration number = %d',gs_counter);
title(title_text)
pause(0.03)
end
%-----------------------------------------------------------------------
% SOR-GS method
if iterative_solver==3
sor_counter =1;
while(error> tol)
for i= 2:nx-1 % Space loop along x axis
for j= 2:ny-1 % Space loop along y axis
t_gs=k*(t(i-1,j)+told(i+1,j)+t(i,j-1)+told(i,j+1));
t(i,j)=omega*t_gs+(1-omega)*told(i,j);
end
end
error = max(max(abs(told-t)));
told=t; % Updating values
sor_counter =sor_counter+1;
end
contourf(t)
xlabel('X axis')
ylabel('U')
title_text = sprintf('Iteration number = %d',sor_counter);
title(title_text)
pause(0.03)
end
%-----------------------------------------------------------------------Changing the iterative_solver value will change the method.
1 - Jacobi
2 - Gauss Seidal
3 - SOR - Gauss Seidal
Output:
Jacobi
Gauss Seidal:
SOR - Gauss Seidal:
Transient State:
`(del^2T)/(del t^2)= alpha((del^2T)/(del x^2) +(del^2T)/(del y^2) )`
`(T^(n+1) -T^(n))/(del t) = alpha((T - 2`(T^(n+1) ) = T^(n)+(del t)*alpha((T - 2T+T)/(deltax^2)+(T-2T+T)/(deltay^2))`T+T)/(deltax^2)+(T-2T+T)/(deltay^2))`
`(T^(n+1) ) = T^(n)+K1((T - 2T+T))+K2((T-2T+T))`
`K1= (alpha*delta t)/(deltax^2)`
`K2= (alpha*delta t)/(deltay^2)`
Source code:
Explicit method:
% Code to solve a second order 2 D PDE - Explicit method
% d^u/dx^2 = 0
% boundary conditions:
% x=0 u=5
%y=0 u=8
clear all
close all
clc
% inputs
nx = 10 % number of grid points along x axis
ny = 10 % number of grid points along y axis
T=ones(nx)*300
T(:,1) = 400; % Boundary condition - Left
T(:,ny) = 800; % Boundary condition - Right
T(1,:) = 600; % Boundary condition - Top
T(nx,:) = 900; % Boundary condition - Bottom
T(1,1)= 500
T(1,ny)= 650
T(nx,1)= 700
T(nx,ny)= 850
error= 9e+9 % random value higher than the tolerance
dt=1e-3
nt=1400
alpha = 1.4;
iterative_solver = 1;
% If iterative_solver = 1 ,then Explicit method
% If iterative_solver = 2,then Implicit Gauss Seidal method
% If iterative_solver = 3,then Implicit Jacobi method
% If iterative_solver = 4,then Implicit SOR
% Creating a copy
Told=T;
T_initial= T;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx= x(2)-x(1);
dy= y(2)-y(1);
K1= (alpha*dt)/(dx^2)
K2= (alpha*dt)/(dy^2)
cfl = K1+K2;
%--------------------------------------------------------------------------
%Explicit method
if iterative_solver==1
counter =1;
for v=1:nt
for i=2:(nx-1)
for j=2:ny-1
term1=Told(i,j);
term2=K1*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j));
term3=K2*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1));
T(i,j)=term1+term2+term3;
end
end
error = max(max(abs(Told-T)));
Told=T; % Updating values
counter =counter+1;
end
contourf(T)
xlabel('X axis')
ylabel('U')
colorbar
title_text = sprintf('Explicit method - Iteration number = %d',counter);
title(title_text)
pause(0.03)
end
%-----------------------------------------------------------------------
Implicit method:
% Code to solve a second order 2 D PDE
% d^u/dx^2 = 0
% boundary conditions:
% x=0 u=5
%y=0 u=8
clear all
close all
clc
% inputs
nx = 10 % number of grid points along x axis
ny = 10 % number of grid points along y axis
T=ones(nx)*300
T(:,1) = 400; % Boundary condition - Left
T(:,ny) = 800; % Boundary condition - Right
T(1,:) = 600; % Boundary condition - Top
T(nx,:) = 900; % Boundary condition - Bottom
T(1,1)= 500
T(1,ny)= 650
T(nx,1)= 700
T(nx,ny)= 850
error= 9e+9 % random value higher than the tolerance
dt=1/50
nt=800
tol=1e-4; % Tolerance limit
alpha = 1.4;
omega = 1.4 % SOR
iterative_solver = 4;
% If iterative_solver = 2,then Implicit Gauss Seidal method
% If iterative_solver = 3,then Implicit Jacobi method
% If iterative_solver = 4,then Implicit SOR
% Creating a copy
Told=T;
T_initial= T;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx= x(2)-x(1);
dy= y(2)-y(1);
K1= (alpha*dt)/(dx^2)
K2= (alpha*dt)/(dy^2)
cfl = K1+K2;
%-----------------------------------------------------------------------
%Implicit GS method
if iterative_solver==2
term1= (1+2*K1+2*K2)^(-1)
term2= K1*term1;
term3= K2*term1;
counter =1;
for v=1:nt
error = 10;
while(error> tol)
for i= 2:nx-1 % Space loop along x axis
for j= 2:ny-1 % Space loop along y axis
Hor=(T(i-1,j)+Told(i+1,j));
Ver=(T(i,j-1)+Told(i,j+1));
T(i,j)= (T_initial(i,j)*term1)+(Hor*term2)+(Ver*term3);
end
end
error = max(max(abs(Told-T)));
Told=T; % Updating values
counter =counter+1;
end
T_initial=T;
end
contourf(T)
xlabel('X axis')
ylabel('U')
title_text = sprintf('GS method - Iteration number = %d',counter);
title(title_text)
pause(0.03)
end
%-----------------------------------------------------------------------
%Implicit Jacobi method
if iterative_solver==3
term1= (1+2*K1+2*K2)^(-1)
term2= K1*term1;
term3= K2*term1;
counter =1;
for v=1:nt
error = 10;
while(error> tol)
for i= 2:nx-1 % Space loop along x axis
for j= 2:ny-1 % Space loop along y axis
Hor=(Told(i-1,j)+Told(i+1,j));
Ver=(Told(i,j-1)+Told(i,j+1));
T(i,j)= (T_initial(i,j)*term1)+(Hor*term2)+(Ver*term3);
end
end
error = max(max(abs(Told-T)));
Told=T; % Updating values
counter =counter+1;
end
T_initial=T;
end
contourf(T)
xlabel('X axis')
ylabel('U')
title_text = sprintf('Jacobi method - Iteration number = %d',counter);
title(title_text)
pause(0.03)
end
%-----------------------------------------------------------------------
%Implicit SOR
if iterative_solver==4
counter =1;
for v=1:nt
error=10;
while(error> tol)
for i= 2:nx-1 % Space loop along x axis
for j= 2:ny-1 % Space loop along y axis
term1= (1+2*K1+2*K2)^(-1)
term2= K1*term1;
term3= K2*term1;
Hor=T(i-1,j)+Told(i+1,j);
Ver=T(i,j-1)+Told(i,j+1);
T(i,j)=(1-omega)*Told(i,j) + omega*(T_initial(i,j)*term1 +(Hor*term2)+(Ver*term3));
end
end
error = max(max(abs(Told-T)));
Told=T; % Updating values
counter =counter+1;
end
T_initial=T
end
contourf(T);
xlabel('X axis')
ylabel('U')
title_text = sprintf('Iteration number = %d',counter);
title(title_text)
pause(0.03)
end
%-----------------------------------------------------------------------
Gauss-Seidel Method:
Jacobi method:
SOR:
Observations:
As we can see that the explicit method is much faster than the implicit method.
But from the results we can see that the implicit method is much accurate compared to the explicit method.
Reason:
In the implicit method, for each time step we make multiple iterations to get a the error lesser than the tolerance limit because the the RHS term values are not fully known.
In the explicit method, for each time step, it makes only one iteration because the RHS term values are already known.