Steady and Unsteady Heat Conduction Analysis
clear all
close all
clc
%given paramaters
Lx=1; %length of the domain in the x-direction
Ly=1; %length of the domain in the y-direction
%Constructing the computational domain
nx=25; %no. of discretized elements in the x-direction
ny=nx; %no. of discretized elements in the y-direction
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1); %Element size or grid size along x-direction
dy=y(2)-y(1);
[X Y]=meshgrid(x,y); %constructing the mesh
%applying Boundary conditions
T=ones(nx,ny);
T(:,1)=400; %left side Boundary condition
T(:,nx)=800; %Right side Boundary condition
T(1,:)=900; %bottom side Boundary condition
T(ny,:)=600; %top side boundary condition
T_old=T;
a=1.5; %relaxation parameter
convergence_criteria=1e-04;
%initializing the loop
error=8e9;
iteration=1;
iterative_solver=3;
tic
while error>convergence_criteria
for i=2:nx-1
for j=2:ny-1
%the equation for the steady state heat conduction is
%d^2(T)/d(x)^2+d^2(T)/d(y)^2=0;
%discretizing this equation using the central difference scheme for the approximation of the 2nd order derivative
%we get the equation as,
if iterative_solver==1 %jacobi Method
T(i,j)=0.25*(T_old(i+1,j)+T_old(i-1,j)+T_old(i,j+1)+T_old(i,j-1));
elseif iterative_solver==2 %Gauss-Seidal Method
T(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1));
elseif iterative_solver==3 %SOR Method
T(i,j)=(1-a)*T_old(i,j)+a*0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1));
end
end
end
error=max(max(abs(T_old-T)));
%update the old values
T_old=T;
iteration=iteration+1;
toc
t=toc
C=contourf(X,Y,T,\'edgecolor\',\'none\')
set(gca,\'YDIR\')
colormap(jet)
clabel(C,\'Fontsize\',10)
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
title({\'2-D Steady State Heat Conduction analysis using SOR Method\';sprintf(\'No. of
iterations:%d\', iteration);sprintf(\'Total Simulation time:%d\', toc)})
end
The above code was written to solve the steady state Heat conduction equation which is given by,
d^2(T)/d(x)^2+d^2(T)/d(y)^2=0
This governing equation was discretized using the 2nd order approximation formula for the 2nd order derivative of the tempurature. The computational domain was considered as a unit square with boundary conditions given as,
1) Top side of the square: 600K
2) Bottom side of the square: 900K
3) Right side of the square:800K
4) Left side of the square:400K
The computational domain was discretized into further smaller squares of length dx=dy=0.0417 units. Then the boundary conditions were applied to the domain. To solve the governing equation in the constructed domain, three iterative methods were used namely 1) Jacobi Method; 2) Gauss-Seidal Method and 3) Successive Over-Relaxation Method. The results and plots of the following methods used to solve the 2-D steady Heat Conducion equation are as follows:
Fig.1-Analysis of 2-D steady State Heat Conduction using the Jacobi iterative solver.
Fig.2-Analysis of 2-D steady State Heat Conduction using the Gauss-Seidal iterative solver.
Fig.3-Analysis of 2-D steady State Heat Conduction using the SOR iterative solver.
From the above three figures there are two points worth noting,
1) The SOR method with a relaxation factor of a=1.5 converged with the least number of iterations i.e, 253 and the Jacobi method takes the most number of iterations for convergence which is 1335.
2) Due to less number of iterations, the SOR method with an over-relaxation parameter of a=1.5 is the fastest method (18.79 sec)to solve the 2-D steady state Heat Conduction equation.
.................................................................................................................................
clear all
close all
clc
%given paramaters
Lx=1; %length of the domain in the x-direction
Ly=1; %length of the domain in the y-direction
%Constructing the computational domain
nx=20; %no. of discretized elements in the x-direction
ny=nx; %no. of discretized elements in the y-direction
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1); %Element size or grid size along x-direction
dy=y(2)-y(1); %Element size or grid size along y-direction
[X Y]=meshgrid(x,y); %constructing the mesh
%applying Boundary conditions
T=ones(nx,ny);
T(:,1)=400; %left side Boundary condition
T(:,nx)=800; %Right side Boundary condition
T(1,:)=900; %bottom side Boundary condition
T(ny,:)=600; %top side boundary condition
T_old=T;
%Unsteady state Heat conduction analysis
gamma=1; %Thermal diffusivity
dt=0.0001; %time step
simulation_time=0.1;
nt=simulation_time/dt; %no. of time steps
k=(gamma*dt)/dx^2; %Fourier number calculation
iteration=1;
tic
for m=1:nt
for i=2:nx-1
for j=2:ny-1
%the equation for the Unsteady state heat condution is
%dT/dt=gamma(d^2(T)/d(x)^2+d^2(T)/d(y)^2);
%discretizing this equation using the central difference scheme for the
approximation of the 2nd order derivative and first order approximation
for the 1st order time derivative
%we get the equation as,
T(i,j)=(1-4*k)*T_old(i,j)+k*(T_old(i+1,j)+T_old(i-
1,j)+T_old(i,j+1)+T_old(i,j-1));
end
end
%update the old values
T_old=T;
iteration=iteration+1;
toc
t=toc;
C=contourf(X,Y,T,\'edgecolor\',\'none\');
pause(0.1)
set(gca,\'YDIR\');
colormap(jet);
clabel(C,\'Fontsize\',10);
title({\'2-D Unsteady State Heat Conduction analysis using explicit
Method\';sprintf(\'No. of iterations:%d\', iteration);sprintf(\'Total Simulation
time:%d\', toc)})
end
The above code was written to solve the Unsteady state Heat conduction equation which is given by,
d^2(T)/d(x)^2+d^2(T)/d(y)^2=dT/dt
To solve this equation using explicit approach, the governing equation was first discretized using the 2nd order approximation formula for the 2nd order spatial derivative of the tempurature and first order forward approximation formula for the temporal derivative of tempurature. The computational domain was considered as a unit square with boundary conditions given as,
1) Top side of the square: 600K
2) Bottom side of the square: 900K
3) Right side of the square:800K
4) Left side of the square:400K
The computational domain was discretized into further smaller squares of length dx=dy=0.0417 units. Then the boundary conditions were applied to the domain. The countour plot for the explicit analyis of the heat conduction equation is shown below.
Fig.4-Analysis of 2-D Unsteady State Heat Conduction using the explicit method.
The impotant points worth noting during this analysis are,
1) The stability of an explicit method is highly dependent on the Courant-Freidrich-Lewy condition which provides a dimensionless parameter called the Courant Number (or Fourier Number). This number is expressed as,
k=(gamma*dt)/dx^2;
For any scheme of disretization used, the courant number always has a limiting value beyond which the solution diverges and becomes unstable. In this case an FTCS scheme has been used which has a theoretical stabiity criteria for,
k<=0.5;
This analysis was run with a courant number of k=0.0576 and the time step dt=0.0001 sec. Keeping all the parameters the same if the time-step was increase 10 times to dt=0.001 sec, the solution diverged because the courant number bcame more that its limiting value i.e, k=0.576>0.5.
2) The number of iterations for which this simulation ran was 1001 and the total simulation time was 0.1 sec. When the total simulation time (nearly to 0.35 sec) was increased, the solution reached steady state values which mimicked the contour plots of Fig. 1, 2 and 3 in the steady state case. An example has bee shown below with an increased simulation time of 0.5 sec.
Fig.5-Analysis of 2-D Unsteady State Heat Conduction using the explicit method with increased simulation time.
.................................................................................................................................
clear all
close all
clc
%given paramaters
Lx=1; %length of the domain in the x-direction
Ly=1; %length of the domain in the y-direction
%Constructing the computational domain
nx=10; %no. of discretized elements in the x-direction
ny=nx; %no. of discretized elements in the y-direction
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
dx=x(2)-x(1); %Element size or grid size along x-direction
dy=y(2)-y(1); %Element size or grid size along y-direction
[X Y]=meshgrid(x,y); %constructing the mesh
%applying Boundary conditions
T=ones(nx,ny);
T(:,1)=400; %left side Boundary condition
T(:,nx)=800; %Right side Boundary condition
T(1,:)=900; %bottom side Boundary condition
T(ny,:)=600; %top side boundary condition
T_old=T;
T_previous=T;
%Given convergence criteria
%for convergence it has been given that the difference between the new and the previous value of tempurature at a grid point should be less than 1e-04
%error=T-T_old;
error=9e9;
%Unsteady state Heat conduction analysis
%Jacobi Method
%Gauss-seidal method
%SOR-Method
gamma=1; %Thermal diffusivity
dt=0.01; %time step
total_simulation_time=10;
nt=total_simulation_time/dt; %no. of time steps
k1=(gamma*dt)/dx^2; %Fourier number calculation
k2=(gamma*dt)/dy^2;
term1=(1+(2*k1)+(2*k2))^(-1);
term2=k1*term1;
term3=k2*term1;
iterative_solver=3;
iteration=1;
tic
for m=1:nt
while error>1e-04
for i=2:nx-1
for j=2:ny-1
%the equation for the Unsteady state heat condution is
%dT/dt=gamma(d^2(T)/d(x)^2+d^2(T)/d(y)^2);
%discretizing this equation using the central difference scheme for the
approximation of the 2nd order derivative and first order approximation
for the 1st order time derivative
%we get the equation as,
if iterative_solver==1 %Jacobi Method
H1=(T_old(i+1,j)+T_old(i-1,j));
V1=(T_old(i,j+1)+T_old(i,j+1));
T(i,j)=(T_previous(i,j)*term1)+(term2*H1)+(term3*V1);
elseif iterative_solver==2 %Gauss-Seidal Method
H2=T(i+1,j)+T(i-1,j);
V2=T(i,j-1)+T(i,j+1);
T(i,j)=(T_previous(i,j)*term1)+(term2*H2)+(term3*V2);
elseif iterative_solver==3 %SOR Method
a=1.2; %relaxation parameter
H2=T(i+1,j)+T(i-1,j);
V2=T(i,j-1)+T(i,j+1);
T(i,j)=(1-a)*T_old(i,j)+a*(T_previous(i,j)*term1+term2*H2+term3*V2);
end
end
end
error=max(max(abs(T_old-T)));
%update the old values
T_old=T;
iteration=iteration+1;
end
T_previous=T;
toc
t=toc
C=contourf(X,Y,T,\'edgecolor\',\'none\');
set(gca,\'YDIR\');
colormap(jet);
clabel(C,\'Fontsize\',10);
title({\'2-D Unsteady State Heat Conduction analysis using SOR Method\';sprintf(\'No.
of iterations:%d\', iteration)})
end
The above code was written to solve the Unsteady state Heat conduction equation which is given by,
d^2(T)/d(x)^2+d^2(T)/d(y)^2=dT/dt
To solve this equation using implicit approach, the governing equation was first discretized using the 2nd order approximation formula for the 2nd order spatial derivative of the tempurature and first order forward approximation formula for the temporal derivative of tempurature. The computational domain was again considered as a unit square with boundary conditions given as,
1) Top side of the square: 600K
2) Bottom side of the square: 900K
3) Right side of the square:800K
4) Left side of the square:400K
The computational domain was discretized into further smaller squares of length dx=dy=0.1111 units. Then the boundary conditions were applied to the domain. The countour plot for the implicit analyis of the heat conduction equation is shown below.
Fig.6-Implicit Analysis of 2-D Unsteady State Heat Conduction using the Jacobi method.
Fig.7-Implicit analysis of 2-D Unsteady State Heat Conduction using the Gauss-Seidal method.
Fig.8-Implicit Analysis of 2-D Unsteady State Heat Conduction using the SOR method.
Some of the important points worth noting are,
1) The stability of the implicit method is independent of the courant number. The unsteady Heat conduction equation is discretized using the FTCS scheme and according to the von-neumann stability analysis, its critical courant number is <=0.5. However, this simulation runs at a courant number of 0.81 and still provides a stable convergent solution.
2)Due to the corresponding independence on the courant number, higher time-steps can be taken for running the simulation.
3)Similar to the results of the steady state heat conduction, the SOR method with a relaxation parameter of a=1.2 simulates the problem in the domain the fastest and with the least number of iterations i.e, 19.
.................................................................................................................................