Solving the steady and unsteady 2D heat conduction problem
This is the coding for solving the steady and unsteady 2D heat conduction problem with the help of three iterative methods
1) Jacobi
2)Gauss-Seidel
3)SOR
The steady-state equation that we are solving is
as this is the steady-state the temperature will show no effect on this system this equation is solved by using three iterative methods in which Gauss-seidel method will take less number of iterations to solve the system.
clear all
close all
clc
%%Defining number of nodes
nx=5;
ny=nx;
l=1;
%%creating grid
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1);
% creating intial array
T=ones(nx,nx);
[xx,yy]=meshgrid(x,y);
% assigning boundary conditions
T(:,1)=400;
T(1,:)=600;
T(:,nx)=800;
T(ny,:)=900;
k=(2*(dx^2+dy^2))/(dx^2*dy^2);
% creating a copy of intial matrix
Told=T;
% introducing error and tolerence
error=5e5;
tol=1e-4;
iterative_solver=input(\'iterative_solver\');
% jacobi method
if iterative_solver==1
jacobi_iter=1;
while (error>tol)
for i=2:nx-1
j=2:ny-1;
T(i,j)= (((Told(i-1,j)+Told(i+1,j))/dx^2)+((Told(i,j+1)+Told(i,j-1))/(dy^2)))/k;
end
error=max(abs(Told(i,j)-T(i,j)));
Told=T;
jacobi_iter= jacobi_iter+1;
end
figure(1)
[c,h]=contourf(x,y,T);
colorbar;
colormap(jet);
clabel(c,h)
Title=sprintf(\'Jacobi iterations for steady state=%d\',jacobi_iter\');
title(Title);
end
%GS method
if iterative_solver==2
gs_iter=1;
while (error>tol)
for i=2:nx-1
j=2:ny-1;
T(i,j)= (((T(i-1,j)+T(i+1,j))/dx^2)+((T(i,j+1)+T(i,j-1))/(dy^2)))/k;
end
error=max(abs(T(i,j)-Told(i,j)));
Told=T;
gs_iter= gs_iter+1;
end
figure(1)
[c,h]=contourf(x,y,T);
colorbar;
colormap(jet);
clabel(c,h)
Title=sprintf(\'gs iterations for steady state=%d\',gs_iter\');
title(Title);
end
% SOR method
if iterative_solver==3
h=1/((nx+1));
w=2/(1+sin(pi*h));
SOR_iter=1;
while (error>tol)
for i=2:nx-1
j=2:ny-1;
T(i,j)= Told(i,j)*(1-w)+(w/k)*((((T(i-1,j)...
+T(i+1,j))/dx^2)+((T(i,j+1)+T(i,j-1))/(dy^2))));
end
error=max(abs(T(i,j)-Told(i,j)));
Told=T;
SOR_iter= SOR_iter+1;
end
figure(1)
[c,h]=contourf(x,y,T);
colorbar;
colormap(jet);
clabel(c,h);
Title=sprintf(\'SOR iterations for steady state=%d\',SOR3_iter\');
title(Title);
end
The unsteady state 2D heat conduction equation is
and here the term alpha represents the diffusivity property for solving this equation in the explicit form it should satisfy the Cournt-friedrids-lewy number this term is used to control the error term and it completely depends on the coding of the present equation and it increases the consistency of the coding. This term is used to tell the boundary limits of the successful execution of the system with minimum error.
Due to the above condition while using an explicit method we should use minimum time step interval otherwise the solution obtained will be in poor accuracy. in this method, the SOR method will work efficiently.
Both implicit and explicit should satisfy the von-Neumann stability condition `modolus of`G should be less than on`
where G= amplification factor
clear all
close all
clc
% inputs of the problem
l=1;
nx=10;
ny=nx;
x=linspace(0,l,nx);
y=linspace(0,l,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
dt=0.01;
alpha=1e-5;
[xx,yy]=meshgrid(x,y);
% B.C of the problem
T=300*ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;
Told=T;
error=9e9;
error1=error;
tol=1e-4;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
iterative_solver=input(\'enter the iterative number\');
jacobi_iter=1;
GS_iter=1;
%coding for jacobi method
if iterative_solver==1
for k=1:20000
while (error>tol)
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(abs(T(i,j)-Told(i,j)));
Told=T;
jacobi_iter=jacobi_iter+1;
end
error=error1;
end
figure(1)
[c,f]=contourf(xx,yy,T);
colormap(jet);
clabel(c,f);
Title=sprintf(\'jacobi iteration for explicit=%d\',jacobi_iter);
title(Title);
end
% coding for GS method
if iterative_solver==2
for k=1:20000
while (error>tol)
for i=2:nx-1
for j=2:ny-1
term1=Told(i,j);
term2=k1*(Told(i+1,j)-2*Told(i,j)+T(i-1,j));
term3=k2*(Told(i,j+1)-2*Told(i,j)+T(i,j-1));
T(i,j)=term1+term2+term3;
end
end
error=max(abs(T(i,j)-Told(i,j)));
Told=T;
GS_iter=GS_iter+1;
end
error=error1;
end
figure(1)
[c,f]=contourf(xx,yy,T);
colormap(jet);
clabel(c,f);
Title=sprintf(\'GS iteration for explicit=%d\',GS_iter);
title(Title);
end
%coding for sor method
if iterative_solver==3
h=1/((nx+1));
w=2/(1+sin(pi*h));
SOR_iter=1;
for k=1:20000
while (error>tol)
for i=2:nx-1
for j=2:ny-1
term1=Told(i,j);
term2=k1*(Told(i+1,j)-2*Told(i,j)+T(i-1,j));
term3=k2*(Told(i,j+1)-2*Told(i,j)+T(i,j-1));
Tf(i,j)=term1+term2+term3;
T(i,j)=w*Tf(i,j)+(1-w)*Told(i,j);
end
end
error=max(abs(T(i,j)-Told(i,j)));
Told=T;
SOR_iter=SOR_iter+1;
end
error=error1;
end
figure(1)
[c,f]=contourf(xx,yy,T);
colormap(jet);
clabel(c,f);
Title=sprintf(\'SOR iteration for explicit=%d\',SOR_iter);
title(Title);
end
clear all
close all
clc
%%Defining number of nodes
nx=80;
ny=nx;
L=1;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
[xx,yy]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%time plot
dt=0.01;
% temperature matrix
T=300*ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;
Told=T;
Tprev=Told;
error=9e9;
error1=error;
tolerence=1e-4;
%constants
alpha=1e-4;
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
%simplifing implicit terms
term1=(1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
jacob_iter=1;
GS_iter=1;
SOR_iter=1;
iterative_solver=input(\'enter the iterative number\');
% coding for jacob solver
if iterative_solver==1
for k=1:20000
while (error>tolerence)
for i=2:nx-1
for j=2:ny-1
Tx=(Told(i-1,j)+Told(i+1,j));
Ty=(Told(i,j-1)+Told(i,j+1));
T(i,j)=(Tprev(i,j)*term1)+(Tx*term2)+(Ty*term3);
end
end
error=max(abs(Told(i,j)-T(i,j)));
Told=T;
jacob_iter=jacob_iter+1;
end
Tprev=T;
error=error1;
end
figure(1);
[c,T]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(c,T);
Title=sprintf(\'jacob iteration for implicit=%d\',jacob_iter);
title(Title);
end
% coding for GS solver
if iterative_solver==2
for k=1:20000
while (error>tolerence)
for i=2:nx-1
for j=2:ny-1
Tx=(T(i-1,j)+Told(i+1,j));
Ty=(T(i,j-1)+Told(i,j+1));
T(i,j)=(Tprev(i,j)*term1)+(Tx*term2)+(Ty*term3);
end
end
error=max(abs(Told(i,j)-T(i,j)));
Told=T;
GS_iter=GS_iter+1;
end
Tprev=T;
error=error1;
end
figure(1);
[c,T]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(c,T);
Title=sprintf(\'GS iteration for implicit=%d\',GS_iter);
title(Title);
end
% coding for SOR solver
if iterative_solver==3
h=1/((nx+1));
w=2/(1+sin(pi*h));
for k=1:20000
while (error>tolerence)
for i=2:nx-1
for j=2:ny-1
Tx=(Told(i-1,j)+Told(i+1,j));
Ty=(Told(i,j-1)+Told(i,j+1));
Tf(i,j)=(Tprev(i,j)*term1)+(Tx*term2)+(Ty*term3);
T(i,j)=(1-w)*Told(i,j)+w*(Tf(i,j));
end
end
error=max(abs(Told(i,j)-T(i,j)));
Told=T;
SOR_iter=SOR_iter+1;
end
Tprev=T;
error=error1;
end
figure(1);
[c,T]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(c,T);
Title=sprintf(\'SOR iteration for implicit=%d\',SOR_iter);
title(Title);
end