Solving the steady and unsteady 2D heat conduction problem
2D Heat Conduction equation is a partial differential equation that explains how heat gets distributed over a while. Based on the time, the transfer of heat occurs may occur in two ways
- Steady-state condition
- Transient State condition
For any system to achieve a steady state, it has to first undergo an unsteady process. Only after a period of unsteadiness, steady-state is achieved.
The steady-state 2D heat conduction equation is
The Transient state 2D heat conduction equation is
The above described Partial Differential Equations needs to be discredited, a central difference scheme in space is used to discretize this equation.
Discretization of the above PDE leads to a set of equations, which can be solved either Implicitly or Explicitly.
In laymen's terms, Explicit means clear; this is a direct way or an easy way t solve these equations. It is less complex when compared to its counterpart. In this case, we use terms from (n)th iteration or timestep to find terms in (n+1)th iteration or timestep. This seems normal.
In an Implicit method, we use terms from both (n)th time step and (n+1)th time step to form coupled equations and solve them to get the terms in (n+1)th time step.
This explanation of Implicit and Explicit methods is not accurate, but it is easy to understand. For an accurate definition, one may refer to standard Textbooks or videos from youtube.
Now, these Implicit and Explicit methods can further be solved using various iterative methods. In this project, we are limiting to only 3types of methods
- Jacobi Method
- Gauss Sidel Method (GS)
- Successive Over-Relaxation Method (SOR)
The user of the code is given the freedom to select whichever method he/she may want to use.
Temperature Function :
%% Temperature matrix
function T = temperature_matrix(ttop,tbot,tleft,tright,nodex,nodey)
T=298*ones(nodex);
T([1:nodex],nodey)=ttop;
T([1:nodex],1)=tbot;
T(1,[1:nodey])=tleft;
T(nodex,[1:nodey])=tright;
T(1,1)=0.5*[T(2,1)+T(1,2)];
T(nodex,1)=0.5*[T(nodex-1,1)+T(nodex,2)];
T(1,nodey)=0.5*[T(1,nodey-1)+T(2,nodey)];
T(nodex,nodey)=0.5*[T(nodex,nodey-1)+T(nodex-1,nodey)];
end
Program - This program is written in the form of a tree, using loops and conditions, where the selection of Type of Heat flow is at the top level followed by the Explicit/Implicit Method ( for Transient flow ) further followed by the type of iterative method.
clear all
close all
clc
nodex=10;
nodey=10;
x=linspace(0,1,nodex);
delx=x(2)-x(1);
y=linspace(0,1,nodey);
dely=y(2)-y(1);
denominator=2*[(delx)^2+(dely)^2];
T =temperature_matrix(600,900,400,800,nodex,nodey);
Told=T;
Tprev=T;
T_SOR=T;
Threshold=1e-4;
error=1;
display('Type 1 for Steady and 2 for Transient ')
Flow=input('');
% Steady Flow
if Flow==1
display('Type 100 for Jacobi; 101 for Gauss and 102 for SOR')
Method=input('Choose the method-');
%Steady - Jacboi
iteration=1;
if Method==100
while error > Threshold
for i=2:nodex-1
for j=2:nodey-1
T(i,j)=[(Told(i+1,j)+Told(i-1,j))*(dely)^2+(Told(i,j+1)+Told(i,j-1))*(delx)^2]/denominator;
error=max(max(abs(T-Told)));
end
end
Told=T;
iteration=iteration+1;
[xx,yy]=meshgrid(x,y);
figure(1)
surfc(xx,yy,T)
colorbar
shading interp
pause(0.01)
title_text=sprintf('Steady state-Jacobi iteration number= %d',iteration);
title(title_text)
figure(2)
contourf(x,y,T)
[a b]=contourf(x,y,T);
clabel(a,b);
colorbar
pause(0.01)
title_text=sprintf('Steady state-Jacobi iteration number= %d',iteration);
title(title_text)
end
Jacobi_iterations=iteration
end
% Steady - Gauss Sidel
if Method==101
while error > Threshold
for i=2:nodex-1
for j=2:nodey-1
T(i,j)=[(Told(i+1,j)+T(i-1,j))*(dely)^2+(Told(i,j+1)+T(i,j-1))*(delx)^2]/denominator
error=max(max(abs(Told-T)));
end
end
Told=T;
iteration=iteration+1;
[xx,yy]=meshgrid(x,y);
figure(1)
surfc(xx,yy,T)
colorbar
shading interp
pause(0.01)
title_text=sprintf('Steady state-GS iteration number= %d',iteration);
title(title_text)
figure(2)
contourf(x,y,T)
[a b]=contourf(x,y,T);
clabel(a,b);
colorbar
pause(0.01)
title_text=sprintf('Steady state-GS iteration number= %d',iteration);
title(title_text)
end
GS_iterations=iteration
end
% Steady - SOR - GS
if Method==102
SOR_Factor=input('Relaxation Factor = ')
while error > Threshold
for i=2:nodex-1
for j=2:nodey-1
T(i,j)=Told(i,j)+SOR_Factor*[[[(Told(i+1,j)+T(i-1,j))*(dely)^2+(Told(i,j+1)+T(i,j-1))*(delx)^2]/denominator]-Told(i,j)];
error=max(max(abs(Told-T)));
end
end
Told=T;
iteration=iteration+1;
[xx,yy]=meshgrid(x,y);
figure(1)
surfc(xx,yy,T)
colorbar
shading interp
pause(0.01)
title_text=sprintf('Steady state-SOR iteration number= %d',iteration);
title(title_text)
figure(2)
contourf(x,y,T)
[a b]=contourf(x,y,T);
clabel(a,b);
colorbar
pause(0.01)
title_text=sprintf('Steady state-SOR iteration number= %d',iteration);
title(title_text)
end
SOR_iterations=iteration;
end
end
% Transient Flow
if Flow==2
display('Type 1 for Explicit; 2 for Implicit')
Method=input('Choose the method-');
delt=1e-3;
nt=1400;
alpha=1.1;
CFLx=(alpha*delt)/(delx^2);
CFLy=(alpha*delt)/(dely^2);
%For explicit schemes CFLx+CFLy < 0.5
% Explicit
if Method==1
iteration=1
for t=1:nt
error=1
while error > Threshold
for i=2:nodex-1
for j=2:nodey-1
Term1=1/(1+2*CFLx+2*CFLy);
Term2=CFLx*(Told(i-1,j)+Told(i+1,j));
Term3=CFLy*(Told(i,j-1)+Told(i,j+1));
T(i,j)=Term1*(Told(i,j)+Term2+Term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
iteration=iteration+1;
end
[xx,yy]=meshgrid(x,y);
figure(1)
surfc(xx,yy,T)
shading interp
colorbar
pause(0.01)
title_text=sprintf('Unsteady-Explicit iteration number= %d',iteration);
title(title_text)
figure(2)
contourf(x,y,T)
[a b]=contourf(x,y,T);
clabel(a,b);
colorbar
pause(0.001)
title_text=sprintf('Unsteady-Explicit iteration number= %d',iteration);
title(title_text)
end
end
% Implicit
if Method==2
display('Type 100 for Jacobi; 101 for Gauss and 102 for SOR')
Method=input('Choose the method-');
% Transient - Implicit - JAcobi
if Method==100
iteration=1
for t=1:nt
error=1;
while error > Threshold
for i=2:nodex-1
for j=2:nodey-1
Term1=1/(1+2*CFLx+2*CFLy);
Term2=CFLx*(Told(i-1,j)+Told(i+1,j));
Term3=CFLy*(Told(i,j-1)+Told(i,j+1));
T(i,j)=Term1*(Tprev(i,j)+Term2+Term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
iteration=iteration+1;
end
Tprev=T;
[xx,yy]=meshgrid(x,y);
figure(1)
surfc(xx,yy,T)
shading interp
colorbar
pause(0.01)
title_text=sprintf('Unsteady-Jacobi iteration number= %d',iteration);
title(title_text)
figure(2)
contourf(x,y,T)
[a b]=contourf(x,y,T);
clabel(a,b);
colorbar
pause(0.001)
title_text=sprintf('Unsteady-Jacobi iteration number= %d',iteration);
title(title_text)
end
end
% Transient - Implicit - GS
if Method==101
iteration=1
for t=1:nt
error=1;
while error > Threshold
for i=2:nodex-1
for j=2:nodey-1
Term1=1/(1+2*CFLx+2*CFLy);
Term2=CFLx*(T(i-1,j)+Told(i+1,j));
Term3=CFLy*(T(i,j-1)+Told(i,j+1));
T(i,j)=Term1*(Tprev(i,j)+Term2+Term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
iteration=iteration+1;
end
Tprev=T;
[xx,yy]=meshgrid(x,y);
figure(1)
surfc(xx,yy,T)
shading interp
colorbar
pause(0.01)
title_text=sprintf('Unsteady-GS iteration number= %d',iteration);
title(title_text)
figure(2)
contourf(x,y,T)
[a b]=contourf(x,y,T);
clabel(a,b);
colorbar
pause(0.001)
title_text=sprintf('Unsteady-GS iteration number= %d',iteration);
title(title_text)
end
end
% Transient - Implicit - SOR
if Method==102
Tprev=Told;
SOR_Factor=input('Relaxation Factor = ')
iteration=1
for t=1:nt
error=1;
while error > Threshold
for i=2:nodex-1
for j=2:nodey-1
Term1=1/(1+2*CFLx+2*CFLy);
Term2=CFLx*(T(i-1,j)+Told(i+1,j));
Term3=CFLy*(T(i,j-1)+Told(i,j+1));
T(i,j)=Term1*(Tprev(i,j)+Term2+Term3);
T_SOR(i,j)=((1-SOR_Factor)*Tprev(i,j))+(T(i,j)*SOR_Factor);
end
end
error=max(max(abs(Told-T)));
Told=T;
iteration=iteration+1;
end
Tprev=T;
[xx,yy]=meshgrid(x,y);
figure(1)
surfc(xx,yy,T_SOR)
shading interp
colorbar
pause(0.01)
title_text=sprintf('Unsteady-SOR iteration number= %d',iteration);
title(title_text)
figure(2)
contourf(x,y,T_SOR)
[a b]=contourf(x,y,T_SOR);
clabel(a,b);
colorbar
pause(0.001)
title_text=sprintf('Unsteady-SOR iteration number= %d',iteration);
title(title_text)
end
end
end
end
Output for Steady Flow:
1. Jacobi
2. Gauss-Sidel
3. Successive Over Relaxation (SOR Factor =1.4)
From the above graphs, it is clear that the number of iterations required for the solution to converge decreases as we move from Jacobi to Gauss Sidel to the SOR method.
Output for Transient Flow:
1. Explicit Method
2. Implicit-Jacobi Method
3. Implicit- Gauss Sidel Method
4. Implicit- Successive Over Relaxation (Relaxation Factor = 1.4)
Here, we can observe that the number of iterations for Explicit is far less when compared to Implicit Methods. Further observation into Implicit Methods provides similar results as seen in Steady flow. But the number of iterations for Transient flow is very high when compared to Steady-state flow. One thing that seems wrong is that the number of iterations required for implicit methods of Gauss Sidel and SOR methods is the same, which does not make any sense.