Solving the steady and unsteady 2D heat conduction problem(MATLAB)
Objective: Solve the 2D heat conduction equation by using the point iterative techniques
Solution:
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
absolute error criteria = 1e-4
Implementing the following methods for solving implicit equations;
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation
1.Steady-state Equation:
Code for Steady-state Equation:
clear all
close all
clc
%Input Parameter
%Number of Grid Points along x&y axis
nx=10;
ny=10;
x=linspace(0,1,nx);
dx= x(2)-(1);
y=linspace(0,1,ny);
z=fliplr(y);
w=1.1;
iter=1;
t=300*ones([nx,ny]);
%Boundary Condition
t(1,:)=600;
t(end,:)=900;
t(:,1)=400;
t(:,end)=800;
%Creating a copy
t_old=t;
r1=t;
error=19;
tol= 1e-4;
for iterative_solver=1:3
while (error>tol)
for i=2:nx-1
for j=2:ny-1
if(iterative_solver==1)
t(i,j)= jacobi(t,t_old,i,j);
elseif(iterative_solver==2)
t(i,j)= GS(t,t_old,i,j);
elseif(iterative_solver==3)
t(i,j)=SOR(t,t_old,i,j,w);
end
end
end
error= max(max(t-t_old));
t_old=t;
[u1,v1]=meshgrid(x,z);
%Plotting
if(iterative_solver==1)
figure(1);
contourf(u1,v1,t);
[C,h]=contourf(u1,v1,t);
colorbar;
colormap("jet");
clabel(C,h);
xlabel('x-axis')
ylabel('y-axis')
title_text=sprintf('Jacobi iteration number=%d',iter);
title(title_text);
pause(0.003);
elseif(iterative_solver==2)
figure(2);
contourf(u1,v1,t);
[C,h]=contourf(u1,v1,t);
colorbar;
colormap("jet");
clabel(C,h);
xlabel('x-axis')
ylabel('y-axis')
title_text=sprintf('GS iteration number=%d',iter);
title(title_text);
pause(0.003);
else
figure(3)
contourf(u1,v1,t);
[C,h]=contourf(u1,v1,t);
colorbar;
colormap("jet");
clabel(C,h);
xlabel('x-axis')
ylabel('y-axis')
title_text=sprintf('SOR iteration number=%d',iter);
title(title_text);
pause(0.003);
end
iter=iter+1;
end
error=19;
iter=1;
t=r1;
t_old=t;
end
Function Code:
1. For Jacobi Method
function output= jacobi(t,t_old,i,j)
output=0.25*(t_old(i,j-1)+t_old(i,j+1)+t_old(i-1,j)+t_old(i+1,j));
end
2. For Gauss-Seidel Method
function output= GS(t,t_old,i,j)
output= 0.25*(t(i,j-1)+t_old(i,j+1)+t(i-1,j)+t_old(i+1,j));
end
3.For successive over-relaxation (SOR) Method
function output= SOR(t,t_old,i,j,w)
output= (1-w)*t_old(i,j)+ w*(0.25*(t(i,j-1)+t_old(i,j+1)+t(i-1,j)+t_old(i+1,j)));
end
Results:
1.1 Jacobi Method
1.2 Gauss-Seidel Method
1.3 For successive over-relaxation (SOR) Method
From the above results, we can conclude that the SOR method is the fastest method as it solves the equation in 92 iterations.
2. Unsteady or Transient state equation:
The Steady-state equation can be solved with the implicit method but while coming to an unsteady state equation there is a time term, so it can be solved either implicit or explicit.
For Explicit method:
The equation in terms of nodes(i,j);
`T_p^(n+1)=T_p^n+k1*(T_l-2*T_p+T_r )^n+k2*(T_t-2*T_p+T_b )^n`
where `k1=α*(Δt)/(Δx^2 )` `k2=α*(Δt)/(Δy^2 )`
if k1+k2 is greater than 0.5 then it is unstable.
For Implicit method:
The equation in terms of nodes(i,j);
`T_p^(n+1)=Term1*(T_p^n+k1*(T_l+T_r )^(n+1)+k2*(T_t+T_b )^(n+1)`
where `Term1=1/((1+2*k1+2*k2))`
Code for Unsteady state equation:
clear all
close all
clc
%Input Parameter
%Number of Grid Points along x&y axis
nx=10;
ny=10;
x=linspace(0,1,nx);
dx= x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
w=1.2;
t=ones(nx,ny)*298;
%Boundary Condition
t(1,:)=600;
t(end,:)=900;
t(:,1)=400;
t(:,end)=800;
alpha=1.4;
dt=0.001;
%Defining k1 and k2
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);
t_old=t;
t_initial=t;
for iterative_solver=1:3
nt=1:1400;
error=9e9;
tol=1;
tic;
iter=1;
while error>tol
for i=2:nx-1
for j=2:ny-1
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
if iterative_solver==1
t(i,j)=((t_old(i,j)*term1)+(term2*(t(i-1,j)+t(i+1,j))+(term3*(t(i,j-1)+t(i,j+1)))));
elseif iterative_solver==2
t(i,j)=((t_old(i,j)*term1)+(term2*(t(i-1,j)+t_old(i+1,j))+(term3*(t(i,j-1)+t_old(i,j+1)))));
elseif iterative_solver==3
t(i,j)=(1-w)*t_old(i,j)+ w*(((t_old(i,j)*term1)+(term2*(t(i-1,j)+t_old(i+1,j))+(term3*(t(i,j-1)+t_old(i,j+1))))));
end
end
error_n=max(abs(t_old-t));
error=max(error_n);
t_old=t;
iter=iter+1;
end
end
time=toc;
%Plotting
if iterative_solver==1
figure(1)
contourf(x,fliplr(y),t);
[C,h]=contourf(x,fliplr(y),t);
colorbar;
colormap("jet")
clabel(C,h);
text=sprintf('Jacobi Method iteration number=%d; computational time:%0.4f',iter,time);
title(text);
t_old=t_initial;
t=t_initial;
elseif iterative_solver==2
figure(2)
contourf(x,fliplr(y),t);
[C,h]=contourf(x,fliplr(y),t);
colorbar;
colormap("jet")
clabel(C,h);
text=sprintf('GS Method iteration number=%d; computational time:%0.4f',iter,time);
title(text);
t_old=t_initial;
t=t_initial;
else
figure(3)
contourf(x,fliplr(y),t);
[C,h]=contourf(x,fliplr(y),t);
colorbar;
colormap("jet")
clabel(C,h);
text=sprintf('SOR Method iteration number=%d; computational time:%0.4f',iter,time);
title(text);
t_old=t_initial;
t=t_initial;
end
end
Results:
2.1 Jacobi Method
2.2 Gauss-Seidel Method
2.3 For successive over-relaxation (SOR) Method
Conclusion:
The above results show that the SOR method is the fastest while comparing to Gauss-Siedel and Jacobi method as it took 521 iterations to solve.