Solving the steady and unsteady 2D heat conduction problem
The main objective of this challenge is to write a code for simulating transient and steady 2D heat conduction equation. As well as Jacobian, Gauss-Seidel and Successive overrelaxation methods are involved for solving the conduction equation. Then the Implicit and Explicit method too taken into consideration. At last the most efficient method for solving 2D heat conduction both steady and transient case is concluded by taking both Numbers of Iterations and time required to solve the simulation.
We know that general 2D heat conduction equation is
`((delT)/(delt))+ alpha*((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
where T is temperature
and t is time.
For steady-state Time derivative will be zero hence the equation will be;
`alpha*((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
Boundary Conditions:
Domain to be solved is a unit square
The number of divisions(mesh) is equal in both x- and y-direction.
Boundary conditions for steady and transient case
- Left:400K
- Top:600K
- Right:800K
- Bottom:900K
Initial Conditions (Assumption):
The domain is to be at a temperature of 298K
Alpha is assumed to be 1e-3 and
Tolerance is assumed to be 1e-4
BASE CODE:
In this case, the average temperature of neighboring elements is taken for four corner elements temperature. Base code is used and functions are called separately for three solvers to reduce the mess of codes.
clc
clear all
close all
%Defining mesh criteria
nx=20;
ny=20;
length_xaxis=1;
length_yaxis=1;
x=linspace(0,length_xaxis,nx);
y=linspace(0,length_yaxis,ny);
yflip=fliplr(y);
dx=abs(x(2)-x(1));
dy=abs(y(2)-y(1));
%Defining Intial and boundary conditions
%Initial condition: Plate is assumed to be at 298K
t=ones(nx,ny)*298;
%Boundary condtions:
%Left=400K
%Top=600K
%Right=800K
%Bottom=900K
t(1,1:nx)=600;
t(nx,1:nx)=900;
t(1:nx,1)=400;
t(1:nx,ny)=800;
t(1,1)=(t(2,1)+t(1,2))/2;
t(end,end)=(t(end,end-1)+t(end-1,end))/2;
t(1,end)=(t(1,end-1)+t(2,end))/2;
t(end,1)=(t(end,2)+t(end-1,1))/2;
told=t;
%Error and tolerance
error=1;
tolerance=1e-4;
%Inputs for calculation of k
%k=alpha*delta(t)/(delta(x)^2)
alpha=1e-3;
deltax=dx;
deltay=dy;
k=2*(dx^2+dy^2)/(dx^2*dy^2);
alpha=2/(1+sin(pi*dx));
%Function calling for different method of solving the problem
iterative_solver=3;
%Jacobi for explicit method
if iterative_solver==1
t=jacobian_steady(x,error,tolerance,k,told,t,nx,ny,yflip,dx,dy);
elseif iterative_solver==2
t=gauss_seidal_steady(x,error,tolerance,k,told,t,nx,ny,yflip,dx,dy);
else
t=sor_steady(x,error,tolerance,k,told,t,nx,ny,yflip,dx,dy,alpha);
end
FUNCTION CODE FOR JACOBIAN SOLVER:
function t=jacobian_steady(x,error,tolerance,k,told,t,nx,ny,yflip,dx,dy)
[X,Y]=meshgrid(x,yflip);
jacobian_iteration_number=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
t1=1/k;
t2=told(i-1,j)+told(i+1,j);
t3=t2/dx^2;
t4=told(i,j-1)+told(i,j+1);
t5=t4/dy^2;
t(i,j)=t1*(t3+t5);
end
end
error=max(abs(t(i,j)-told(i,j)));
told=t;
jacobian_iteration_number=jacobian_iteration_number+1;
end
jacobian_iteration_number
plot=contourf(X,Y,t);
colormap(jet);
title('Plot for Steady State Simulation using Jacobian method')
xlabel('X-axis')
ylabel('Y-axis')
endFUNCTION CODE FOR GAUSS SEIDAL SOLVER:
function t=gauss_seidal_steady(x,error,tolerance,k,told,t,nx,ny,yflip,dx,dy)
[X,Y]=meshgrid(x,yflip);
gauss_seidal_iteration_number=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
t1=1/k;
t2=t(i-1,j)+told(i+1,j);
t3=t2/dx^2;
t4=t(i,j-1)+told(i,j+1);
t5=t4/dy^2;
t(i,j)=t1*(t3+t5);
end
end
error=max(abs(t(i,j)-told(i,j)));
told=t;
gauss_seidal_iteration_number=gauss_seidal_iteration_number+1;
end
gauss_seidal_iteration_number
plot=contourf(X,Y,t);
colormap(jet);
title('Plot for Steady State Simulation using Gauss Seidal method')
xlabel('X-axis')
ylabel('Y-axis')
endFUNCTION CODE FOR SUCCESSIVE OVER RELAXATION SOLVER:
function t=sor_steady(x,error,tolerance,k,told,t,nx,ny,yflip,dx,dy,alpha)
[X,Y]=meshgrid(x,yflip);
sor_iteration_number=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
t1=1/k;
t2=t(i-1,j)+told(i+1,j);
t3=t2/dx^2;
t4=t(i,j-1)+told(i,j+1);
t5=t4/dy^2;
t6=(1-alpha)*told(i,j);
t7=alpha*t1*(t3+t5);
t(i,j)=t6+t7;
end
end
error=max(abs(t(i,j)-told(i,j)));
told=t;
sor_iteration_number=sor_iteration_number+1;
end
sor_iteration_number
plot=contourf(X,Y,t);
colormap(jet);
title(sprintf('Plot for Steady State Simulation using Successive Over Relaxation method. Number of Iterations=%d',sor_iteration_number))
xlabel('X-axis')
ylabel('Y-axis')
endTRANSIENT 2D HEAT CONDUCTION EQUATION
BASE CODE
clc
clear all
close all
tic
%Defining mesh criteria
nx=10;
ny=10;
length_xaxis=1;
length_yaxis=1;
x=linspace(0,length_xaxis,nx);
y=linspace(0,length_yaxis,ny);
yflip=fliplr(y);
dx=abs(x(2)-x(1));
dy=abs(y(2)-y(1));
%defining total time for transient simulation
start_time=0;
end_time=5000;
time_interval=0.1;
time_range=[start_time:time_interval:end_time];
%Defining Intial and boundary conditions
%Initial condition: Plate is assumed to be at 298K
t=ones(nx,ny)*298;
%Boundary condtions:
%Left=400K
%Top=600K
%Right=800K
%Bottom=900K
t(1,1:nx)=600;
t(nx,1:nx)=900;
t(1:nx,1)=400;
t(1:nx,ny)=800;
t(1,1)=(t(2,1)+t(1,2))/2;
t(end,end)=(t(end,end-1)+t(end-1,end))/2;
t(1,end)=(t(1,end-1)+t(2,end))/2;
t(end,1)=(t(end,2)+t(end-1,1))/2;
told=t;
%Error and tolerance
error=1;
tolerance=1e-4;
%Inputs for calculation of k1 and k2
%k=alpha*delta(t)/(delta(x)^2)
alpha=1e-3;
deltat=time_interval;
deltax=dx;
deltay=dy;
k1=alpha*deltat/(deltax^2);
k2=alpha*deltat/(deltay^2);
%Function calling for different method of solving the problem
iterative_solver=4;
%Jacobi for explicit method
if iterative_solver==1
t=jacobian_explicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
elseif iterative_solver==2
t=jacobian_implicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
elseif iterative_solver==3
t=gauss_seidal_explicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
elseif iterative_solver==4
t=gauss_seideal_implicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
elseif iterative_solver==6
t=sor_gauss_seidal_explicit(x,y,end_time,error,tolerance,k1,k2,told,t,nx,ny,yflip,dx)
elseif iterative_solver==7
t=sor_jacobian_implicit(x,y,end_time,error,tolerance,k1,k2,told,t,nx,ny,yflip,dx)
else
t=sor_jacobian_explicit(x,y,end_time,error,tolerance,k1,k2,told,t,nx,ny,yflip,dx)
end
time=toc
FUNCTION CODE FOR JACOBIAN EXPLICIT
function t=jacobian_explicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
[X,Y]=meshgrid(x,yflip);
for z=1:end_time
jacobian_iteration_number=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
t1=told(i,j);
t2=told(i-1,j)+told(i+1,j)-(2*told(i,j));
t3=k1*t2;
t4=told(i,j-1)+told(i,j+1)-(2*told(i,j));
t5=k1*t4;
t(i,j)=t1+t3+t5;
error=max(abs(t(i,j)-told(i,j)));
told=t;
end
end
jacobian_iteration_number=jacobian_iteration_number+1;
end
end
time=toc;
jacobian_iteration_number
plot=contourf(X,Y,t);
colormap(jet);
title(sprintf('Plot for transient simulation using Jacobian Explicit method \n Time=%d',time))
xlabel('X-axis')
ylabel('Y-axis')
end
FUNCTION CODE FOR JACOBIAN IMPLICIT
function t=jacobian_implicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
[X,Y]=meshgrid(x,yflip);
constant_matrix=told;
t3=1+(2*k1)+(2*k2);
t4=1/t3;
for k=1:end_time
Jacobian_Implicit_iteration_number=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
t1=told(i-1,j)+told(i+1,j);
t2=told(i,j-1)+told(i,j+1);
t5=k1*t1*t4;
t6=k2*t2*t4;
t7=constant_matrix(i,j)*t4;
t(i,j)=t6+t7+t5;
error=max(max(abs(t-told)));
end
end
told=t;
constant_matrix=t;
Jacobian_Implicit_iteration_number=Jacobian_Implicit_iteration_number+1
end
end
time=toc;
plot=contourf(X,Y,t);
colormap(jet);
title(sprintf('Plot for transient simulation using Jacobian Implicit method \n Time=%d',time))
xlabel('X-axis')
ylabel('Y-axis')
end
FUNCTION CODE FOR GAUSS SEIDAL EXPLICIT
function t=gauss_seidal_explicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
[X,Y]=meshgrid(x,yflip);
for z=1:end_time
Gauss_Seidal_Explicit_iteration_number=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
t1=told(i,j);
t2=t(i-1,j)+told(i+1,j)-(2*told(i,j));
t3=k1*t2;
t4=t(i,j-1)+told(i,j+1)-(2*told(i,j));
t5=k2*t4;
t(i,j)=t1+t3+t5;
error=max(abs(t(i,j)-told(i,j)));
told=t;
end
end
Gauss_Seidal_Explicit_iteration_number=Gauss_Seidal_Explicit_iteration_number+1;
end
end
time=toc
plot=contourf(X,Y,t);
colormap(jet);
title(sprintf('Plot for transient simulation using Gauss Seidel Explicit method \n Time=%d',time))
xlabel('X-axis')
ylabel('Y-axis')
end
FUNCTION CODE FOR GAUSS SEIDAL IMPLICIT
function t=gauss_seideal_implicit(x,y,time_range,error,tolerance,k1,k2,told,t,nx,ny,yflip,end_time)
[X,Y]=meshgrid(x,yflip);
constant_matrix=told;
t3=1+(2*k1)+(2*k2);
t4=1/t3;
for k=1:end_time
Gauss_Seidal_Implicit_iteration_number=1;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
t1=t(i-1,j)+told(i+1,j);
t2=t(i,j-1)+told(i,j+1);
t5=k1*t1*t4;
t6=k2*t2*t4;
t7=constant_matrix(i,j)*t4;
t(i,j)=t6+t7+t5;
error=max(abs(t(i,j)-told(i,j)));
end
end
Gauss_Seidal_Implicit_iteration_number= Gauss_Seidal_Implicit_iteration_number+1;
told=t;
constant_matrix=t;
end
end
time=toc
plot=contourf(X,Y,t);
colormap(jet);
title(sprintf('Plot for transient simulation using Gauss Seidal Implicit method \n Time=%d',time))
xlabel('X-axis')
ylabel('Y-axis')
endFUNCTION FOR SUCCESSIVE OVER RELAXATION JACOBIAN EXPLICIT CODE:
function t=sor_jacobian_explicit(x,y,end_time,error,tolerance,k1,k2,told,t,nx,ny,yflip,dx)
[X,Y]=meshgrid(x,yflip);
opt_omega = 2/(1+sin(pi*dx));
previous_dt = told;
t=told;
GS_iteration_SOR = 1;
% t3=1+(2*k1)+(2*k2);
% t4=1/t3;
for k = 1:end_time
while (error > tolerance)
for i = 2:nx-1
for j = 2:ny-1
H = (told(i-1,j)+told(i+1,j)-(2*told(i,j)));
V = (told(i,j-1)+told(i,j+1)-(2*told(i,j)));
GS = told(i,j)+(k1*(H+V));
t(i,j) = opt_omega*(GS)+(1-opt_omega)*told(i,j);
end
end
told = t;
GS_iteration_SOR = GS_iteration_SOR + 1;
error = max(max(abs(told-t)));
end
previous_dt = t;
error = 1;
end
time=toc;
plot=contourf(X,Y,t);
colormap(jet);
title(sprintf('Plot for SOR Jacobian Implicit method \n Time=%d',time))
xlabel('X-axis')
ylabel('Y-axis')
Jacobian_Explicit_Iteration_Number=GS_iteration_SOR
end
FUNCTION FOR SUCCESSIVE OVER RELAXATION JACOBIAN IMPLICIT CODE:
function t=sor_jacobian_implicit(x,y,end_time,error,tolerance,k1,k2,told,t,nx,ny,yflip,dx)
[X,Y]=meshgrid(x,yflip);
opt_omega = 2/(1+sin(pi*dx));
previous_dt = told;
t=told;
GS_iteration_SOR = 1;
t3=1+(2*k1)+(2*k2);
t4=1/t3;
for k = 1:end_time
while (error > tolerance)
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) = opt_omega*((previous_dt(i,j)/t3)+(t4*k1*H)+(t4*k2*V))+(1-opt_omega)*previous_dt(i,j);
end
end
error = max(max(abs(told-t)));
told = t;
GS_iteration_SOR = GS_iteration_SOR + 1;
end
previous_dt = t;
error = 1;
end
time=toc;
plot=contourf(X,Y,t);
colormap(jet);
title(sprintf('Plot for SOR Jacobian Implicit method \n Tme=%d',time))
xlabel('X-axis')
ylabel('Y-axis')
GS_iteration_SOR
end
So from the above results, we can observe that the Number of Iterations taken and time required to solve Jacobian is higher than Gauss-Seidel solver which is, in turn, higher than Successive over-relaxation method. This is because of the fact Jacobi takes guess values throughout the process, so it takes much more iterations to get converged which proportionally increases the time required for the simulation to run. While the Gauss-Seidel method takes calculated values for (i-1) and (j-1) which is left and top values for calculation of the next set of temperature values. So the number of Iterations for this case is lower than the Jacobian method. On top of both of this solver, Successive over-relaxation method uses an alpha value which is based purely on element discretization sizes will converge much more rapidly than Gauss and Jacobi method.
In case of transient simulation, a number of iterations and time taken for Implicit is code is much lower than the time for explicit codes, this is because Implicit takes all values from previous time step while Explicit takes values from current time step.
Initial Temperature Matrix with boundary conditions and initial conditions for a number of divisions 10 in x- and y-axis:
Final Temperature Matrix:
This seems to be logical because at the lower bottom of the plot temperature is higher as it goes up it reduces. But this resultant temperature value is not accurate because of poor mesh. If we go for some fine mesh say 100 division in x- and y-axis the value will be better. The plot for this case is shown below:
If there are 250division in x- and the y-axis the plot will be like:
If there are 500 division in x- and the y-axis the plot will be like:
This might be weird. In general, as mesh size is reduced the results should be accurate but ironically it’s behaving another way. The main reason is we have to have an eye on tolerance limit as well. Because in all the cases mentioned above the tolerance is 1e-4. So the elements will converge much rapidly if mesh size is reduced. Thus there should be a proper balance for providing mesh size and tolerance value. So the plot for a number of divisions 250 in x- and y-axis with the tolerance of 1e-8 is shown below:
But this case takes 83576 iterations while 10 divisions take just 223 iterations to converge.
WHICH IS BEST TO SOLVE?
In my case, the fastest steady-state simulation is Successive over-relaxation method suing Gauss-Seidel approach.
For a number of divisions=500
For a number of divisions=10
Final Temperature Value:
For transient simulation too, fastest is Successive Overrelaxation Implicit approach.
Thus from the plots of Temperature contours and Temperature matrix we can see both the results are exactly the same. But if we look about computational time and Number of Iterations:
|
SOLVER |
TIME(s) |
NUMBER OF ITERATIONS |
|
Steady State |
0.0144 |
146 |
|
Transient Explicit SOR |
0.27553 |
162603 |
|
Transient Implicit SOR |
0.818265 |
276367 |
Steady-state solver solves the equation quicker as there is no play of temperature derivative term in this case. On the other hand transient case is having temperature derivative as well. So it takes some more time than steady-state. But the explicit method has lower time and number of iterations than an implicit method because it takes all values at previous time steps.
STABILITY ANALYSIS:
The unstable solution looks like:
Thus the unstable solution is mainly because of the CFL factor which is `(alpha*deltat)/(deltax^2). For the previous challenges, alpha value is 1e-3;deltat=0.1 and deltax=0.111, which gives CFL value as 8.1e-4. But for the two images above, alpha is 1e-1 and time interval is 1. Now the CFL number is equal to 8.1 which gives an unstable solution. If the number goes beyond 0.5 in SOR and 0.25 in Gauss&Jacobian solution becomes unstable.
The final temperature matrix will be like:
NaN values are infinity. This occurs because of the matrix is trying to converge with each iteration but still, the values are too higher so it goes to Infinity. It can be well explained by the following image:
This matrix is obtained after the first time step, in the 9th row we can see the temperature in the range of 6.97e3 which is practically impossible for this case. Since maximum temperature is 900K all temperature in the plate should be below this value. So this causes an unstable solution. This case going forward time step will produce temperature matrix like:
It goes in the scale of e20 and some values are in negative which is practically impossible. So this produces unstable contour at last as well.
VIDEOS:
UNITS: TEMPERATURE IN KELVIN
2D TRANSIENT HEAT CONDUCTION WITH 10 NODES IN X- AND Y-AXIS:
2D TRANSIENT HEAT CONDUCTION WITH 25 NODES IN X- AND Y-AXIS:
UNSTABLE SOLUTION FOR CFL=60.2
UNSTABLE SOLUTION FOR CFL=2.02
UNSTABLE SOLUTION FOR CFL=245