Solving the steady and unsteady 2D heat conduction problem
Aim:To Solve the 2D heat conduction equation by using the point iterative techniques.
Inputs we have assumed here
L=1
nx=20
ny=nx
error criteria is 1e-4
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
Techniques or methods we are using :
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation
Steady-state solution:
`(∂^2T)/(∂x^2) + (∂^2T)/(∂y^2) = 0`
Discretizing the above equation:
`(∂^2T)/(∂x^2) = (T(x-h) -2T(x)+T(x+h))/h^2`
`(∂^2T)/(∂y^2) = (T(y-h) -2T(y)+T(y+h))/h^2`
`(∂T)/(∂t) = (T^(n)-T)/(△t)`
MATLAB Coding:
% solving the steady 2D heat conduction problem
close all
clear all
clc
%inputs
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=30;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=ones(nx,ny);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
% jacobi method
for iterative_solver=1;
tic
if iterative_solver==1;
jacobi_iterations=1;
while (error >tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=0.25*(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;
jacobi_iterations = jacobi_iterations+1;
end
end
F=toc;
%plotting
figure(1);
[P,Q]=contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('Jacobi method,Iteration number=%d,computation time F=%fs',jacobi_iterations,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf('jacobi=%d/n',jacobi_iterations);
pause(0.03);
end
% for GS method
for iterative_solver=2;
% update values
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=30;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=298*ones(nx,ny);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
tic;
if iterative_solver==2;
gs_iterations=1;
while (error >tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=0.25*(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;
gs_iterations = gs_iterations+1;
end
end
F=toc;
%plotting
figure(2);
[P,Q]=contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('Gauss seidal method,Iteration number=%d,computation time F=%fs',gs_iterations,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf('gs=%d/n',gs_iterations);
pause(0.03);
end
%SOR
for iterative_solver=3;
%update values
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=30;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=ones(nx,ny);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
tic;
if iterative_solver==3;
SOR_iterations=1;
while (error >tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=((1-omega)*Told(i,j)+(omega*0.25*(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;
SOR_iterations = SOR_iterations+1;
end
end
F=toc;
%plotting
figure(3);
[P,Q]=contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('SOR method,Iteration number=%d,computation time F=%fs',SOR_iterations,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf('SOR=%d/n',SOR_iterations);
pause(0.03);
endOutput:
Results:
By comparing above plots that is jacobi,gs & sor method, sor method is solving faster with less iterations and compution time when compared to two other methds.
SOR is taking 93 iterations and computation time =0.009769s
Gauss seidel is taking 954 iterations and computation time=0.04521s
Jacobi method is taking 1884 iterations and computation time=0.11530s
Transient method:
`(T(x-h) -2T(x)+T(x+h))/h^2 + (T(y-h) -2T(y)+T(y+h))/h^2 = (1/alpha)(T^n-T)/(△t)`
Discretizing:
`((T(x-h) -2T(x)+T(x+h))/h^2 + (T(y-h) -2T(y)+T(y+h))/h^2)^n= (1/alpha)(T^n-T)/(△t)`
For explicit method
(Tn+1)i,j=λ(Tn)i−1,j+(Tn)i,j−1
lambda is also called the CFL number, CFL number is used to check the stability of the solution.
alpha is the convective coefficient.
MATLAB coding:
% solving the steady 2d transient unsteady state heat conduction problem
% explicit
clear all
close all
clc
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=51;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1e-3; %convective coefficient
t=200;
%CFL number
M = alpha*dt/(dx^2);
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=298*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
nt=t/dt;
%2d transient unsteady state heat conduction equation in explicit method
tic;
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
Told=T;
end
time=toc;
%plottings
figure(1);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('Explicit method,computation time =%fs,M=%f',time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.03);Output:
For implicit
Tn+1p=Tnp+k1(Tn+1L+Tn+1R)+K2(T
% solving the steady 2d transient unsteady state heat conduction problem
% implicit
clear all
close all
clc
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=20;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1; %convective coefficient
t=200;
%CFL number
M = alpha*dt/(dx^2);
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=298*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
nt=t/dt;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
T_prev_dt=Told;
%2d transient unsteady state heat conduction equation in explicit method
nt=t/dt;
iterative_solver=1;
tic;
jacobi_iterations=1;
for v=1:nt
if iterative_solver==1
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=100;
while(error >tol)
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)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
jacobi_iterations=jacobi_iterations+1;
end
T_prev_dt=T;
end
time=toc;
%plottings
figure(1);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('Jacobi method,Iteration number=%d,computation time =%fs,M=%f',jacobi_iterations,time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.03);
end
% %for gs method
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=20;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1; %convective coefficient
t=200;
% %CFL number
M = alpha*dt/(dx^2);
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=298*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
% %calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
nt=t/dt;
% %creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
T_prev_dt=Told;
% %2d transient unsteady state heat conduction equation in explicit method
iterative_solver=2;
tic;
gs_iterations=1;
for w=1:nt
if iterative_solver==2
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=100;
while(error >tol)
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)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
gs_iterations=gs_iterations+1;
end
T_prev_dt=T;
end
F=toc;
% %plottings
figure(2);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('Gauss siedal,Iteration number=%d,computation time =%fs,M=%f',gs_iterations,time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.03);
end
% %for SOR method
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=20;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1; %convective coefficient
t=200;
% %CFL number
M = alpha*dt/(dx^2);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
%define boundary conditions
T_L=400;% T_R=800;
T_T=600;
T_B=900;
T=298*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
nt=t/dt;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
T_prev_dt=Told;
%2d transient unsteady state heat conduction equation in explicit method
iterative_solver=3;
tic;
SOR_iterations=1;
for Z=1:nt
if iterative_solver==3
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=100;
while(error >tol)
for i=2:nx-1
for j=2:ny-1
H=(omega*T(i-1,j)+(1-omega)*Told(i-1,j)+Told(i+1,j));
V=(omega*T(i,j-1)+(1-omega)*Told(i,j-1)+Told(i,j+1));
T(i,j)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
SOR_iterations=SOR_iterations+1;
end
T_prev_dt=T;
end
F=toc;
%plottings
figure(3);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('SOR,Iteration number=%d,computation time =%fs,M=%f',SOR_iterations,time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.03);
endOutput:
Jacobi method
Gauss-Seidel method
SOR
Result:
As we can see from the results of explicit and implicit methods,explicit method is much faster than the implicit method.
The computational time for the explicit method is 0.171388s.
The fastest method in implicit is the sor with 2588 iterations.
But from the explicit and implicit method, implicit is much accurate as it uses various iterative solvers to converge the solution and above it uses all values from the entire domain.