2D Heat Conduction flow analysis using Matlab

Theory:

To solve the Steady-state analysis & Transient State Analysis for the 2D heat conduction equation by using the point iterative techniques by implementing the following methods:

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation (SOR)

Your absolute error criteria is 1e-4.

Assumptions:

1. The domain is a unit square.
2. nx = ny [Number of points along the x-direction is equal to the number of points along the y direction]
3. Boundary conditions of temperatures in Kelvin for steady and transient case: Left:400K; Top:600K; Right:800K; Bottom:900K.
4. alpha=0.25.
5. dt=0.005.
6. error=9e9
7. tolerance(tol)=1e-5 (for the program and analysis, this value is used)

Note: There are a total of 8 programs in all as per all the given conditions of our problem. Though the best way was to write the programs as sub-functions and call it from the main program as per the need of the user, the programs here are written individual as a separate one for each one of them for clear viewing of the programs their solution and their out animation. Some references have been taken from different sources available for a better understanding of the problem.

*************************************************START OF PROGRAMS*******************************************************************

_____________________________________________________________________________________________________________________________________________________

Program 1: Transient State analysis for 2D Heat conduction flow using Jacobi method Implicitly.

clear all
close all
clc

nx=10; %no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005; %assumption

Told=T; %copying our original array

alpha=0.25; % Assumption

k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);

%jacobi method

error=9e9;
tol=1e-5;

itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=Told(i-1,j);
Tt=Told(i,j+1);
Tb=Told(i,j-1);
term1= (1+2*k1+2*k2);
T(i,j)=(Told(i,j)+k1*(Tr+Tl)+k2*(Tt+Tb))/term1;
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet) %used to colormap our output
contourf(x,y,T) % used for getting a filled colorplot
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf); % to get an animation movie of our output
itno=itno+1;
colorbar % shows the reference color scale with our plot
pause(0.001)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'Transient_state_jacobi_method_implicit.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 823

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

_____________________________________________________________________________________________________________________________________________________

Program 2: Transient State analysis for 2D Heat conduction flow using Jacobi method Explicitly.

clear all
close all
clc

nx=10;%no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005;%assumption

Told=T; %copying our original array

alpha=0.25; % Assumption

k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);

%jacobi method

error=9e9;
tol=1e-5;

itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=Told(i-1,j);
Tt=Told(i,j+1);
Tb=Told(i,j-1);
term1= (1-2*k1-2*k2);
T(i,j)=Told(i,j)*term1+k1*(Tr+Tl)+k2*(Tt+Tb);
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet)
contourf(x,y,T)
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf);
itno=itno+1;
colorbar
pause(0.03)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'Transient_state_jacobi_method_explicit.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 598

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

_____________________________________________________________________________________________________________________________________________________

Program 3: Steady-State analysis for 2D Heat conduction flow using Jacobi method.

clear all
close all
clc

nx=10; %no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005; %assumption

Told=T; %copying our original array

%jacobi method

error=9e9;
tol=1e-5;
k=2*(dx^2+dy^2)/(dx^2*dy^2);
itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=Told(i-1,j);
Tt=Told(i,j+1);
Tb=Told(i,j-1);
T(i,j)=((Tr+Tl)/(dx^2)+(Tt+Tb)/(dy^2))/k;
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet)
contourf(x,y,T)
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf);
itno=itno+1;
colorbar
pause(0.03)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'steady_state_jacobi_method.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 253

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

_____________________________________________________________________________________________________________________________________________________

Program 4: Transient State analysis for 2D Heat conduction flow using Gauss Siedel method Implicitly.

clear all
close all
clc

nx=10; %no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005; %assumption

Told=T; %copying our original array

alpha=0.25; % Assumption

k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);

%Gauss-Siedel method

error=9e9;
tol=1e-5;

itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=T(i-1,j);
Tt=Told(i,j+1);
Tb=T(i,j-1);
term1= (1+2*k1+2*k2);
T(i,j)=(Told(i,j)+k1*(Tr+Tl)+k2*(Tt+Tb))/term1;
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet)
contourf(x,y,T)
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf);
itno=itno+1;
colorbar
pause(0.001)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'transient_state_gauss_method_implicit.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 718

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

_____________________________________________________________________________________________________________________________________________________

Program 5: Transient State analysis for 2D Heat conduction flow using Gauss Siedel method Explicitly.

clear all
close all
clc

nx=10; %no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005; %assumption

Told=T; %copying our original array

alpha=0.25; % Assumption

k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);

%Gauss-Siedel method

error=9e9;
tol=1e-5;

itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=T(i-1,j);
Tt=Told(i,j+1);
Tb=T(i,j-1);
term1= (1-2*k1-2*k2);
T(i,j)=(Told(i,j)*term1+k1*(Tr+Tl)+k2*(Tt+Tb));
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet)
contourf(x,y,T)
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf);
itno=itno+1;
colorbar
pause(0.01)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'transient_state_gauss_method_explicit.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 490

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

_____________________________________________________________________________________________________________________________________________________

Program 6: Steady-State analysis for 2D Heat conduction flow using Gauss Siedel method.

clear all
close all
clc

nx=10; %no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005; %assumption

Told=T; %copying our original array

%Gauss-Siedel method method

error=9e9;
tol=1e-5;
k=2*(dx^2+dy^2)/(dx^2*dy^2);
itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=T(i-1,j);
Tt=Told(i,j+1);
Tb=T(i,j-1);
T(i,j)=((Tr+Tl)/(dx^2)+(Tt+Tb)/(dy^2))/k;
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet)
contourf(x,y,T)
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf);
itno=itno+1;
colorbar
pause(0.03)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'steady_state_gauss_method.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 134

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

_____________________________________________________________________________________________________________________________________________________

Program 7: Transient State analysis for 2D Heat conduction flow using the Successive Over Relaxation (SOR) method. (Note: Under relaxation is taken i.e alpha<1, hence the time and iterations to converge will be more. Had Over Relaxation been taken, our time to converge and iterations would have been the least as compared to other methods)

clear all
close all
clc

nx=10; %no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005; %assumption

Told=T; %copying our original array

alpha=0.25; % Assumption

k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);

%SOR Method

error=9e9;
tol=1e-5;

itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=Told(i-1,j);
Tt=Told(i,j+1);
Tb=Told(i,j-1);
term1= (1-2*k1-2*k2);
T(i,j)=(1-alpha)*(Told(i,j))+alpha*((Told(i,j)*term1)+(k1*(Tr+Tl))+k2*(Tt+Tb));
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet)
contourf(x,y,T)
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf);
itno=itno+1;
colorbar
pause(0.01)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'transient_state_SOR_method.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 2183

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

_____________________________________________________________________________________________________________________________________________________

Program 8: Steady-State analysis for 2D Heat conduction flow using the Successive Over Relaxation (SOR) method. (Note: Under relaxation is taken i.e alpha<1, hence the time and iterations to converge will be more. Had Over Relaxation been taken, our time to converge and iterations would have been the least as compared to other methods)

clear all
close all
clc

nx=10; %no. of grids in x and y direction
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);

T=ones(nx,ny); %initialising the temperature array to defauult value

%implementing the initial boundary conditions

T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;

dt=0.005; %assumption

Told=T; %copying our original array
alpha=0.25; % Assumption

%SOR Method

error=9e9;
tol=1e-5;
k=2*(dx^2+dy^2)/(dx^2*dy^2);
itno=1;

while(error>tol) % Start of time loop

for i=2:nx-1 %start of convergence loop for a multi row/column matrix
for j=2:ny-1
Tr=Told(i+1,j);
Tl=Told(i-1,j);
Tt=Told(i,j+1);
Tb=Told(i,j-1);
T(i,j)=(1-alpha)*Told(i,j)+alpha*((((Tr+Tl)/dx^2)+((Tt+Tb)/dy^2))/k);
end
end
error=max(max(abs(Told-T)));
Told=T;

colormap(jet)
contourf(x,y,T)
title(sprintf(\'Iteration number = %d\',itno))
M(itno)=getframe(gcf);
itno=itno+1;
colorbar
pause(0.0001)
end

movie(M);%making a movie of array M
videofile=VideoWriter(\'Steady_state_SOR_method.avi\',\'uncompressed avi\');% creating the videofile as pendulum.avi and in uncompressed format
open(videofile)%opening the video file
writeVideo(videofile,M)%writing the video file by the frames stored in array M
close(videofile)

Output:

Final no. of iterations: 940

Temperature Range: 400-900K

X (length) = 0 to 1 metre

_____________________________________________________________________________________________________________________________________________________

***************************************************END OF PROGRAMS*******************************************************************

Comparison and conclusions:

• Transient state Jacobi method (implicit v/s explicit):

  	Implicit	Explicit
  No. of Iterations	823	598
  Video File time (seconds)	28	20

• Transient State Gauss Siedel Method(implicit v/s Explicit):

  	Implicit	Explicit
  No. of Iterations	718	490
  Video File time (seconds)	23	17

• Steady-state Jacobi method v/s Steady-state Gauss-Siedel method:

  	Jacobi	Gauss-Siedel
  No. of Iterations	253	134
  Video File time (seconds)	09	05

• Transient State SOR method v.s Steady-state SOR method:

  	Transient SOR	Steady State SOR
  No. of Iterations	2183	940
  Video File time (seconds)	71	32

• From the above data, the following conclusions are made. (Keeping in mind that under relaxation SOR method is used and not over-relaxation)

1. The transient state implicit Gauss-Siedel method is faster than the Transient state implicit Jacobi method. (As expected)
2. The transient state explicit Gauss-Siedel method is faster than the Transient state explicit Jacobi method. (As expected)
3. The steady-state Gauss-Siedel method is faster than the Steady-state Jacobi method. (As expected)
4. The steady-state under relaxed SOR method is way faster than the Transient state under relaxed SOR method (As expected, because the transient state is an unsteady one and takes time to first converge and then stabilize)
5. In terms of time, i.e faster method for the transient state: Transient State OVER relaxed SOR > Transient state Gauss-Siedel Explicit > Transient state Jacobi Explicit > Transient State Gauss-Siedel Implicit > Transient state Jacobi Implicit.
6. In terms of time, i.e faster method for the steady-state: Steady-state OVER relaxed SOR> Steady-state Gauss-Siedel > Steady-state Jacobi > Steady-state under relaxed SOR