Solving the steady and unsteady 2D heat conduction problem
In this project, we will look to solve Steady and Unsteady 2D heat conduction problem using point iterative solvers such as JACOBI, Gauss-Seidel, Successive Over Relaxation (SOR) methods.
For solving the Steady-state equation of 2D heat conduction:
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
Therefore by Central Differencing,
`(del^2T)/(delx^2)=(T(i-1,j)-2T(i,j)+T(i+1,j))/(Deltax^2)`
`(del^2T)/(dely^2)=(T(i,j-1)-2T(i,j)+T(i,j+1))/(Deltay^2)`
For better understanding, we can substitute
T(i-1,j)= T_L
T(i+1,j)=T_R
T(i,j+1)= T_B
T(i,j-1)= T_T
T(i,j)=T_p
Therefore the equation becomes,
`(T(i-1,j)-2T(i,j)+T(i+1,j))/(Deltax^2)+ (T(i,j-1)-2T(i,j)+T(i,j+1))/(Deltay^2)=0`
So,
`Tp=(1/k)((TL+TR)/(Deltax^2)+(T T+TB)/(Deltay^2))`
where k is,
`k=(dx^2+dy^2)/(dx^2.dy^2)`
Following are the code for the Steady state problem solution,
%Steady 2D heat conduction problem using Jacobi Method
%Length of domain
L=1;
%number of nodes
n_x=10;
n_y=n_x;
%Grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%Boundary condition
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
%copy of Temperature
Told=T;
%Valueof K
k=2*((dx^2+dy^2)/((dx^2)*(dy^2)));
tolerance=1e-4;
error=7e7;
%2D conduction equation for steady state analysis
%Using Jacobi Method
tic
jacobi_iteration=1;
while error>tolerance
for i=2:n_x-1
for j=2:n_y-1
term_1= (1/k)*((Told(i-1,j)+Told(i+1,j))/(dx^2));
term_2= (1/k)*((Told(i,j-1)+Told(i,j+1))/(dy^2));
T(i,j)=term_1+term_2;
end
end
%Error
error=max(max(abs(Told-T)));
%Updating the Told
Told=T;
jacobi_iteration=jacobi_iteration+1;
end
Iteration_time=toc;
figure(1)
contourf(x,y,T);
[a,b]= contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
title_name=sprintf('Number of iterations=%d Iteration Time=%d', jacobi_iteration,Iteration_time);
title(title_name);
saveas(gca,'Jacobi Iteration steady state.jpg');
%Steady 2D heat conduction problem using Gauss Seidal
%Length of domain
L=1;
%number of nodes
n_x=10;
n_y=n_x;
%Grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%Boundary condition
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
%copy of Temperature
Told=T;
%Valueof K
k=2*((dx^2+dy^2)/((dx^2)*(dy^2)));
tolerance=1e-4;
error=7e7;
%2D conduction equation for steady state analysis
%Using Guass Seidal Method
tic
GS_iteration=1;
while error>tolerance
for i=2:n_x-1
for j=2:n_y-1
term_1= (1/k)*((T(i-1,j)+T(i+1,j))/(dx^2));
term_2= (1/k)*((T(i,j-1)+T(i,j+1))/(dy^2));
T(i,j)=term_1+term_2;
end
end
%Error
error=max(max(abs(Told-T)));
%Updating the Told
Told=T;
GS_iteration=GS_iteration+1;
end
Iteration_time=toc;
figure(1)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
title_name=sprintf('Iterative solver=%d Iteration time=%d', GS_iteration,Iteration_time);
title(title_name);
saveas(gca,'Gauss seidal steady state.jpg');
%Steady 2D heat conduction problem
%3) SOR METHOD
%Length of domain
L=1;
%number of nodes
n_x=10;
n_y=n_x;
%Grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%Boundary condition
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
%copy of Temperature
Told=T;
%Valueof K
k=2*((dx^2+dy^2)/((dx^2)*(dy^2)));
tolerance=1e-4;
error=7e7;
%Relaxation factor
omega=1.1;
%2D conduction equation for steady state analysis
%Using SOR Method
tic
SOR_iteration=1;
while error>tolerance
for i=2:n_x-1
for j=2:n_y-1
term_1= (1/k)*((T(i-1,j)+T(i+1,j))/(dx^2));
term_2= (1/k)*((T(i,j-1)+T(i,j+1))/(dy^2));
T(i,j)= (1-omega)*Told(i,j)+(omega*(term_1+term_2));
end
end
%Error
error=max(max(abs(Told-T)));
%Updating the Told
Told=T;
SOR_iteration=SOR_iteration+1;
end
Iteration_time=toc;
figure(3)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
title_name=sprintf('SOR Iterative solver=%d Iteration Time=%d', SOR_iteration,Iteration_time);
title(title_name);
saveas(gca,'SOR steady state.jpg');Following are the Output of the JACOBI, Gauss Seidal, SOR methods respectively:
JACOBI METHOD:
The time required for convergence is 0.001731 sec and Iteration carried out are 217 which is the highest number of iterations in Steady State solution.
Gauss-Seidel:
The time required for convergence is 0.00035 sec and Iteration carried out are 117 which is the moderate number of iterations in Steady-State solution.
SOR METHOD:
The time required for convergence is 0.0002363 sec and Iteration carried out are 96 which is the lowest number of iterations in Steady-State solution. Also, the SOR method is the fastest to reach convergence among other iterative solvers. Also in the SOR method, by increasing the relaxation factor omega to 1.2,1.3,... it is seen that the number of iterations reduces from 97 to 62 iterations.
FOR solving the Transient case of 2D heat conduction problem:
`(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))`
by forward differencing for the Temporal terms LHS:
`(delT)/(delt)=(T^(n+1)i,j-T^ni,j)/(Deltat)`
Using central differencing for the Spacial terms.
`(del^2T)/(delx^2)=(T(i-1,j)-2T(i,j)+T(i+1,j))/(Deltax^2``
`(del^2T)/(dely^2)=(T(i,j-1)-2T(i,j)+T(i,j+1))/(Deltay^2)`
Therefore after substituting them letters,
This equation is for Implicit Form.
`Tp^(n+1)=(Tp^n+k1(TL+TR)^(n+1)+k2(T T+TB)^(n+1))/(1+2k1+2k2)`
where k1 & k2 are,
`k1= alpha(Deltat)/(Deltax^2)` `&``k2= alpha(Deltat)/(Deltay^2)`
EXPLICIT form is:
`Tp^(n+1)=(Tp^n+k1(TL+TR)^(n)+k2(T T+TB)^(n))/(1+2k1+2k2)`
The code for the explicit form of Transient case is as follows.
%EXPLICIT METHOD FOR TRANSIENT SOLUTION
clear all
close all
clc
%2D heat conduction
%Length of domain
L=1;
%number of nodes
n_x=10;
n_y=n_x;
%Grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%defining the temperature terms
alpha=1.4;
nt=1000;
dt=0.001;
%Boundary condition
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
%Refining the edges of the T
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
%CFL value
k1= alpha*(dt/(dx^2));
k2= alpha*(dt/(dy^2));
error=5e5;
tolerance=1e-4;
%Copy of Temperature value
Told=T;
method=1;
if method==1
%EXPLICIT CALCULATION
Explicit_iteration=1;
%Time loop
tic
for t=1:nt
%spacial loop
for i=2:n_x-1
for j=2:n_y-1
term_1= Told(i,j);
term_2= k1*(Told(i-1,j)-(2*Told(i,j))+Told(i,j+1));
term_3= k2*(Told(i,j-1)-(2*Told(i,j))+Told(i,j+1));
T(i,j)=term_1+term_2+term_3;
end
end
%Updating the old value of T
Told=T;
Explicit_iteration=Explicit_iteration+1;
end
end
iteration_time=toc;
figure(1)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
colorbar
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
title_text=sprintf('Number of iteration=%d Iteration Time=%d',Explicit_iteration,iteration_time);
title(title_text);
saveas(gca,'Explicit-Transient.jpg');The result of the explicit scheme is :
The iterations required to converge the simulation is 1001 and time is 0.0071 sec. The explicit scheme also when decreased or increased the value of dt from 0.001 sec, the results become unstable and we don't obtain any solution for it. Also, the explicit scheme takes less time than any other implicit scheme methods.
Code for the IMPLICIT SCHEMES ARE AS FOLLOWS:
%Implicit method
clear all
close all
clc
%2D heat conduction using JACOBI METHOD
%Length of domain
L=1;
%number of nodes
n_x=10;
n_y=n_x;
%Grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%defining the temperature terms
alpha=1.4;
nt=1000;
dt=0.001;
%Boundary condition
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
%Refining the edges of the T
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
%CFL value
k1= alpha*(dt/(dx^2));
k2= alpha*(dt/(dy^2));
cfl = k1+k2;
error=7e7;
tolerance=1e-4;
%Copy of Temperature value
Told=T;
T_prev_dt=Told;
%Simplication of JACOBI equation
term_1= (1+2*k1+2*k2)^-1;
term_2= k1*term_1;
term_3= k2*term_1;
method=1;
%then JACOBI method
jacobi_itr=1;
%Time Loop
tic
for t=1:nt
error = 7e7;
%CONVERGENCE LOOP
while(error>tolerance)
%Nodal loop
for i=2:n_x-1
for j=2:n_y-1
V=(Told(i-1,j))+(Told(i+1,j));
H=(Told(i,j-1))+(Told(i,j+1));
T(i,j)=(term_1*(T_prev_dt(i,j)))+(term_2*V)+(term_3*H);
end
end
error=max(max(abs(Told-T)));
Told=T;
jacobi_itr=jacobi_itr+1;
end
T_prev_dt=T;
end
Iteration_time=toc;
figure(1)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
clabel(a,b);
colormap(jet);
xlabel('X axis');
ylabel('Y axis');
colorbar
set(gca,'YDIR','reverse');
title_txt=sprintf('Number of iterations=%d Iteration Time=%d',jacobi_itr,Iteration_time);
title(title_txt);
saveas(gca,'Implicit Jacobi Method.jpg');
%Implicit method
clear all
close all
clc
%2D heat conduction using GS METHOD
%Length of domain
L=1;
%number of nodes
n_x=10;
n_y=n_x;
%Grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%defining the temperature terms
alpha=1.4;
nt=1000;
dt=0.001;
%Boundary condition
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
%Refining the edges of the T
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
%CFL value
k1= alpha*(dt/(dx^2));
k2= alpha*(dt/(dy^2));
CFL = k1+k2;
error=7e7;
tolerance=1e-4;
%Copy of Temperature value
Told=T;
T_prev_dt=Told;
%Simplication of GS equation
term_1= (1+2*k1+2*k2)^-1;
term_2= k1*term_1;
term_3= k2*term_1;
%then GUASS SAIDEL method
GS_itr=1;
%Time Loop
tic
for t=1:nt
error = 7e7;
%CONVERGENCE LOOP
while(error>tolerance)
%Nodal loop
for i=2:n_x-1
for j=2:n_y-1
V=(T(i-1,j))+(Told(i+1,j));
H=(T(i,j-1))+(Told(i,j+1));
T(i,j)=(term_1*(T_prev_dt(i,j)))+(term_2*V)+(term_3*H);
end
end
error=max(max(abs(Told-T)));
Told=T;
GS_itr=GS_itr+1;
end
T_prev_dt=T;
end
Iteration_time=toc;
figure(1)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
clabel(a,b);
colormap(jet);
xlabel('X axis');
ylabel('Y axis');
colorbar
set(gca,'YDIR','reverse');
title_txt=sprintf('Number of Iterations=%d Iteration Time=%d', GS_itr,Iteration_time);
title(title_txt);
saveas(gca,'Implicit GS METHOD.jpg');
%Implicit method
clear all
close all
clc
%2D heat conduction using GS METHOD
%Length of domain
L=1;
%number of nodes
n_x=10;
n_y=n_x;
%Grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%defining the temperature terms
alpha=1.4;
nt=1000;
dt=0.001;
%Boundary condition
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
%Refining the edges of the T
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
%CFL value
k1= alpha*dt/(dx^2);
k2= alpha*dt/(dy^2);
CFL = k1+k2;
%Relaxation Factor
omega=1.1;
error=7e7;
tolerance=1e-4;
%Copy of Temperature value
Told=T;
T_prev_dt=Told;
%Simplication of GS equation
term_1= (1+2*k1+2*k2)^-1;
term_2= k1*term_1;
term_3= k2*term_1;
%then GUASS SAIDEL method
SOR_itr=1;
%Time Loop
tic
for t=1:nt
error = 7e7;
%CONVERGENCE LOOP
while(error>tolerance)
%Nodal loop
for i=2:n_x-1
for j=2:n_y-1
H=(T(i-1,j))+(Told(i+1,j));
V=(T(i,j-1))+(Told(i,j+1));
GS=(term_1*(T_prev_dt(i,j)))+(term_2*H)+(term_3*V);
T(i,j)=((1-omega)*Told(i,j))+(omega*(GS));
end
end
error=max(max(abs(Told-T)));
Told=T;
SOR_itr=SOR_itr+1;
end
T_prev_dt=T;
end
Iteration_time=toc;
figure(1)
contourf(x,y,T);
[a,b]=contourf(x,y,T);
clabel(a,b);
colormap(jet);
xlabel('X axis');
ylabel('Y axis');
colorbar
set(gca,'YDIR','reverse');
title_txt=sprintf('Number of Iterations=%d Iteration Time=%d', SOR_itr,Iteration_time);
title(title_txt);
saveas(gca,'Implicit SOR method.jpg');The result of the implicit scheme is as follows:
JACOBI METHOD:
The time required for convergence is 0.0198 sec and Iteration carried out are 3483 which is the highest number of iterations in Transient-State solution.
GAUSS SEIDEL METHOD:
The time required for convergence is 0.0180 sec and Iteration carried out are 2874 which is the moderate number of iterations in Transient-State solution.
SOR METHOD:
The time required for convergence is 0.0172 sec and Iteration carried out are 2467 which is the lowest number of iterations in Transient-State solution. Also, it is the fastest method in Implicit schemes.
CONCLUSIONS:
1) The steady-state solution is the fastest schemes to compute the solution. This is also because they are absent we the temperature terms in the equation which leads to faster computation.
2) Explicit scheme forms to be the fastest among the transient method. It requires way less number of iterations and computation time to converge compared to Implicit schemes.
3) Relaxation factor (omega) must be maintained within a range of 1.1-1.3. More or less than the range it may become over-relaxed or under-relaxed.