To solve the 2D heat conduction problem for the steady and unsteady state using Jacobi, gauss seidel, and successive overrelaxation iterative techniques.
Theory
A general heat conduction equation is
When considering internal heat generation to be zero,
For an unsteady state the 2-D heat conduction equation is represented as :
For steady state the 2-D heat conduction equation is represented as :
Explicit Method
Explicit methods calculate the state of a system at a later time from the state of the system at the current time.This computation is called Explicit. Thismethod requires stability criteria which when satisfied leads to stable solution. Because of the stabality criteria small time steps hasto be taken to obtain stable solutions.Stability criterias such as CFL number is employed for checking stability.
Implicit Method
Implicit methods find a solution by solving an equation involving both the current state of the system and the later one. When a computation is made for dependent variable are defined by set of equations and a matrix is involved to find the unknownvariable. This computation is called implicit method. This method is unconditionally stable as a result, no stability criteria is required to obtain Stable results.
In order to solve for dependent variables iterative solvers used are
1. Jacobi Iterative solver
2.Gauss seidel iterative solver
3.Successive over relaxation iterative solver.
1. Jacobi Iterative solver
In this method the temperature at (n+1)th term iteration is found using old values of the temperature defined at n'th term iteration.Using the approximation obtained iterations procedure it continued until desired accuracy of the solution is obtained.
2.Gauss seidel iterative solver
In this type iterative solver, temperature at n+1th iterations uses values of n+1 iterations that have already been computed. By this wayold approximated values are overwritten with new approximate values for the iterations performed. Unlike jacobi, Gauss seidel requireone storage vector which can be overwritten with each computation. This helps for the solution to converge faster than jacobimethod.
3.Successive over relaxation iterative solver.
In this type of iterative solver relaxation factor is used in increasing the convergence rate.
Given data
Top boundary temperature = 600 K
Bottom Boundary temperature = 900 K
Left Boundary temperature = 400 K
Right Boundary temperature = 800 K
Error = 1e-4
Steady-State 2D Heat conduction
Initially inputs are defined into variables and the temperature array is created.
Boundary conditions are introduced into the temperature array, and at corners, boundary values are averaged.
Input must be provided by the user, to choose the iterative method.
The if-else loop is used to solve the equation using the selected(input from the user) iterative technique.
while loop is used to run the iterative technique till the required accuracy is achieved.
Two for loop is used to calculate the temperature values at every node at each iterative step.
Finally, the temperature distribution graph is obtained using 'conturf' command.
Program code:
clear all
close all
clc
%Specifying length off domain and no of grid points
l=1;
nx=100;
ny=100;
%Defining space domain
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
%Creating temperatures contour
t=300*ones(nx,ny);
t(:,1)=400;
t(1,:)=600;
t(:,end)=800;
t(end,:)=900;
%Temperature at boundaries
t(1,1)=(600+400)/2;
t(1,end)=(600+800)/2;
t(end,1)=(400+900)/2;
t(end,end)=(800+900)/2;
%Updating old t values
told=t;
%Solving 2D steady state system
[xx,yy]=meshgrid(x,y);
%Assigning tolerance and error values
tolerance=1e-4;
error=1e8;
%Defining K factor
k=2*((1/(dx^2))+(1/(dy^2)));
%Iteration method Input
iteration_method=input('Enter the iteration method number= ')
%Jacobi Method
%Defining iteration variable to store and display no of iterations computated
if iteration_method==1
tic
jacobi_iter=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
term1=(told(i-1,j)+told(i+1,j))/(k*(dx^2));
term2=(told(i,j-1)+told(i,j+1))/(k*(dy^2));
t(i,j)=term1+term2;
end
end
error=max(max(abs(told-t)));
told=t;
jacobi_iter=jacobi_iter+1;
end
time1=toc;
end
%Gauss Siedal Method
if iteration_method==2
tic
gs_iter=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
term1=(t(i-1,j)+told(i+1,j))/(k*(dx^2));
term2=(t(i,j-1)+told(i,j+1))/(k*(dy^2));
t(i,j)=term1+term2;
end
end
error=max(max(abs(told-t)));
told=t;
gs_iter=gs_iter+1;
end
time2=toc;
end
%SOR METHOD
if iteration_method==3
omega=1.8;
tic
sor_iter=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
term1=(t(i-1,j)+told(i+1,j))/(k*(dx^2));
term2=(t(i,j-1)+told(i,j+1))/(k*(dy^2));
t(i,j)=(1-omega)*told(i,j)+omega*(term1+term2);
end
end
error=max(max(abs(told-t)));
told=t;
sor_iter=sor_iter+1;
end
time3=toc;
end
%Plotting the results of temperature
[c,h]=contourf(xx,yy,t);
clabel(c,h);
colorbar;
colormap(jet);
if iteration_method==1
title({'Steady state temperature profile';'Iteration Method= Jacobi Method ';['No of iterations= ',num2str(jacobi_iter)];['Simulation time= ',num2str(time1),'s']})
elseif iteration_method==2
title({'Steady state temperature profile';'Iteration Method = Gauss Seidal Method';['No of iterations= ',num2str(gs_iter)];['Simulation time= ',num2str(time2),'s']})
elseif iteration_method==3
title({'Steady state temperature profile';'Iteration Method= SOR Method';['No of iterations= ',num2str(sor_iter)];['Simulation time= ',num2str(time3),'s']})
end
xlabel('X-axis')
ylabel('Y-axis')
Results:
Jacobi Method
Gauss Seidal Method:
SOR Method:
Transient/Unsteady heat conduction:
Explicit Method
Program code
clear all
close all
clc
l=1;
nx=60;
ny=60;
%Defining space domain
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
%Calculation of time step
%dt=(0.5)/((alpha/dx^2)+(alpha/dy^2));
dt=0.00001;
alpha=1.4;
%Creating temperature contour
t=300*ones(nx,ny);
t(:,1)=400;
t(1,:)=600;
t(:,end)=800;
t(end,:)=900;
%Updating old t values
told=t;
CFL_num=(alpha*dt)*((1/dx^2)+1/(dy^2));
%Solving 2D unsteady state system% [xx,yy]=meshgrid(x,y);
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
iter_method=1;
tic;
for k=1:7000
for i =2:nx-1
for j=2:ny-1
term1=(1-2*k1-2*k2);
term2=k1;
term3=k2;
H=(told(i-1,j)+told(i+1,j));
V=(told(i,j-1)+told(i,j+1));
t(i,j)=(told(i,j)*term1)+(term2*H)+(term3*V);
end
end
told=t;
end
time=toc;
figure(1)
colormap(jet)
contourf(x,y,t,'showtext','on');
colorbar
% clabel(c,h);
title({'Unsteady temperature profile';['Simulation time= ',num2str(time),'s']})
xlabel('X-axis');
ylabel('Y-axis');esult:
Implicit method:
clear all
close all
clc
nx=50;
ny=50;
x=linspace(0,1,nx);
dx=1/(nx-1);
y=linspace(0,1,ny);
dy=1/(ny-1);
alpha=1.4;
%dt=(0.5)/((alpha/dx^2)+(alpha/dy^2));
dt=0.001;
tolerance=1e-4;
error=1e9;
t=300*ones(nx,ny);
t(:,1)=400;
t(1,:)=600;
t(:,end)=800;
t(end,:)=900;
%Temperature at edges
t(1,1)=(600+400)/2;
t(1,end)=(600+800)/2;
t(end,1)=(400+900)/2;
t(end,end)=(800+900)/2;
%Updating old t values
told=t;
CFL_num=(alpha*dt)*((1/dx^2)+1/(dy^2));
%Solving 2D steady state system
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
iter_method=input('Enter the iteration method number= ');
if iter_method==1
%Computing Time loop
jacobi_iter=1;
for nt=1:600
%Jacobi Iterative Solver
tic
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
error=1e9;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
H=(told(i-1,j)+told(i+1,j));
V=(told(i,j-1)+told(i,j+1));
t(i,j)=(told(i,j)*term1)+(term2*H)+(term3*V);
end
end
error=max(max(abs(told-t)));
told=t;
jacobi_iter=jacobi_iter+1;
end
time1=toc;
end
end
%Gauss Siedel Iteration method
if iter_method==2
gs_iter=1;
tic
for nt=1:600
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
error=1e9;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
H=(t(i-1,j)+told(i+1,j));
V=(t(i,j-1)+told(i,j+1));
t(i,j)=(told(i,j)*term1)+(term2*H)+(term3*V);
end
end
error=max(max(abs(told-t)));
told=t;
gs_iter=gs_iter+1;
end
time2=toc;
end
end
%SOR iteration Mtehod
omega=1.2;
if iter_method==3
sor_iter=1;
tic
for nt=1:600
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
error=1e9;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
H=(t(i-1,j)+told(i+1,j));
V=(t(i,j-1)+told(i,j+1));
t(i,j)=(1-omega)*told(i,j)+omega*(term1*told(i,j)+term2*H+term3*V);
end
end
error=max(max(abs(told-t)));
told=t;
sor_iter=sor_iter+1;
end
time3=toc;
end
end
figure(1)
[c,h]=contourf(x,y,t);
clabel(c,h);
colorbar;
colormap(jet);
if iter_method==1
title({'Unsteady state temperature profile';['Iteration Method= Jacobi'];['No of iterations= ',num2str(jacobi_iter)]})
elseif iter_method==2
title({'Unsteady state temperature profile';['Iteration Method= Gauss Siedel'];['No of iterations= ',num2str(gs_iter)]})
elseif iter_method==3
title({'Unsteady state temperature profile';['Iteration Method= SOR'];['No of iterations= ',num2str(sor_iter)]})
end
Results:
Jacobi Method
Gauss Seidel Method
SOR Method:
Conclusions:
From the results obtained it is clear that all the three different techniques gave results almost similar to each other, but what is different is the number of iterations and computational time it consumes.
For steady state heat transfer
For transient or unsteady state heat transfer
From the bar chart it is clear that Jacobi iteration method too many number of iterations to reach the specified accuracy level or the set conditions of convergence.
Successive over relaxation method acheived convergence with lowest number of iterations.
Because it is likewise predicated on future state values, the implicit technique necessitates guess values and iterations for convergence. Which are unknown at the time of calculation It verifies convergence at each time step. The simulation is delayed as a result of this.
Because they are conditionally stable, implicit methods can be utilised to solve issues that require converged answers in a short amount of time. The convergence rate increases as the time step size increases, resulting in a steady state with fewer itearations.
