Solving the steady and unsteady 2D heat conduction equation
Aim: To solve a heat conduction equation of a 2D plate surrounded by constant temperature by the Finite Difference Method. Here we will consider both the Steady-state and Transient state(Unsteady) case and compare which iterative method gives better results with fewer number iteration. Also, we will understand the effect of CFL numbers in the case of transient cases.
Assumptions:
1. Absolute error criteria is 1e-4
2. Consider an equal number of grid points or nodes along the X-axis and Y-axis
3. No internal heat generation
4. Boundary conditions are as follows
- Top Boundary = 600 k
- Bottom Boundary = 900 k
- Left Boundary = 400 k
- Right Boundary = 800k
Governing Equation:`(del^2T)/(delx^2)+(del^2T)/(dely^2)+(del^2T)/(delz^2)+q/k=(rhoc_p)/k (delT)/(delt)`
Here.
q= Rate of Heat Generation
k= Thermal Conductivity
`rho`= Density
`c_p`=Specific heat capacity at constant pressure
`rhoc_p`=volumetric heat capacity
Considering 2D case and no heat generation constant pressure
`(del^2T)/(delx^2)+ (del^2T)/(dely^2)=1/alpha (delT)/(delt)`
`(delT)/(dt)=alpha[(del^2T)/(delx^2)+(del^2T)/(dely^2)]`
Numerical Discretization:
1. Steady-State Case:
In the case of steady-state, the time derivative term of temperature becomes zero. So the governing equation becomes :
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
Now Discretizing each of the terma by finite difference method we get,
`(del^2T)/(delx^2)=(T_(i,j+1)^n+2T_(i,j)^n+T_(i,j-1)^n)/((Deltax)^2)`
`(del^2T)/(dely^2)=(T_(i+1,j)^n+2T_(i,j)^n+T_(i-1,j)^n)/((Deltay)^2)`
After rearranging the terms we get,
`T_(i,j)=(T_(i+1,j) + T_(i-1,j)+ T_(i,j+1)+T_(i,j-1))/4`
This simultaneous set of eqaution can be solved using the followiing iterative method:
- Jacobi:`T_P=T_(old_(T))+T__(old_(B))+T_(old_(L))+T_(old_(R))`
- Gauss -Seidal:`T_P=T_T+T_(old_(B)+T_L+T_(old_(R)))`
- Successive Order Reduction:`T_(p_(SOR))=(1-w)T_(p_(old))+wT_(p_(GS))`
Here the subscript old refers to the value in the previous iteration and the subscripts T, B R and L refders to Top, Bottom,Right and Left respectively.
2. Transient State Case:
In the case of the transient case the time derivative term of temperature is considered. So the governing equation becomes:
`alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=(delT)/(delt)`
Here,
`alpha=k/(rhoc_p)=1.4m^2/s`(considered arbitarily)
Now Discretizing each of the terms by finite difference method we get,
`(delT)/(delt)=(T_l^(n+1)-T_l^n)/(Deltat)`
Now we get two cases, Explicit and Implicit based on the discretization of the remaining two terms.
Considering the previous time step will give rise to explicit terms and considering the same time step will give rise to the implicit case.
Explicit Case:
`(del^2T)/(delx^2)=(T_(i,j+1)^n +2T_(i,j)^n+T_(i,j-1)^n)/(Deltax)^2`
`(del^2T)/(dely^2)=(T_(i+1,j)^n+2T_(i,j)^n+T_(i-1,j)^n)/(Deltay)^2`
After calculation consider `Delta x=Delta y` and k1=k2
`T_p =T_(p_(old))[1-4k]+k(T_(T_(old))+T_(B_(old))+T_(L_(old))+T_(R_(old)))`
Implicit Case:
`(del^2T)/(delx^2)= (T_(i,j+1)^(n+1)+2T_(i,j)^(n+1)+T_(i,j-1)^(n+1))/(Deltax^2)`
`(del^2t )/(dely^2)=(T_(i+1,j)^(n+1)+2T_(i,j)^(n+1)+T_(i-1,j)^(n+1))/(Deltay)^2`
Now just like the staeadt state equations this set of implicit case can be solved by the following iterative methods.
Here `T_(P_(prev))`= values of `T_p`in the previous time step
- jacobi
- gauss- seidal
- Successive Over Relaxation.
User-Defined Functions:
1. linear_eqn_solver_jacobi
function linear_eqn_solver_jacobi(T,Told,Tprev,nx,ny,nt,x,y,k1,k2,tol,type)
%%Transient Implicit (Jacobi)
if type == 11
tic;
iteration = 1;
for k= 1:nt %time loop
error=9e9;
while(error > tol) % Convergence loop
for i= 2:nx-1 % Space loop
for j= 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
V = Told(i-1,j) + Told(i+1,j);
H = Told(i,j-1) + Told(i,j+1);
T(i,j) = (Tprev(i,j)*term1) + (k1*term1*H) + (k2*term1*V);
end
end % end of space loop
error = max(max(abs(Told-T)));
Told =T;
iteration = iteration +1;
end % end of convergence loop
Tprev =T;
end % end of time loop
time = toc;
f1=figure(1)
contourf(x,y,T,'--','ShowText', 'on')
colorbar
colormap('jet')
text=sprintf("2D Transient Heat Conduction-Implicite-Jacobi No. of iteration %f Time of simulation %f Sec",iteration,time);
title(text)
set(gca,'YDIR','reverse')
%Steady State (Jacobi)
elseif type== 21
iteration =1;
tic;
error=9e9;
while (error> tol) %Convergence loop
for i= 2:nx-1 %Space Loop
for j=2:ny-1
term1 = Told(i,j-1) +Told(i,j+1);
term2 = Told(i-1,j) + Told(i+1,j);
T(i,j) = (1/4)*(term1+term2);
end
end % End of space loop
error=max(max(abs(Told-T)));
Told=T;
iteration =iteration+1;
end % End of convergence loop
time=toc;
f1=figure(1)
contourf(x,y,T,'--','ShowText','on')
colorbar
colormap('jet')
text = sprintf("2D Steady State Heat Conduction - Jacobi No. of iteration %f Time of Simulation %f Sec",iteration,time);
title(text)
set(gca,'YDIR','reverse')
end
end
2.linear_eqn_solver_Gs
function linear_eqn_solver_Gs(T,Told,Tprev,nx,ny,nt,x,y,k1,k2,tol,type)
%Transient Implicit(Gauss seidal)
if type == 12
iteration = 1;
tic;
for k=1:nt %Time loop
error = 9e9;
while(error > tol) % convergence loop
for i=2:nx-1 %Space Loop
for j=2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
V = T(i-1,j) + Told(i+1,j);
H = T(i,j-1)+ Told(i,j+1);
T(i,j) = (Tprev(i,j)*term1)+(k1*term1*H)+(k2*term1*V);
end
end %end of space loop
error =max(max(abs(Told-T)));
Told=T;
iteration = iteration+1;
end % End of convergence loop
Tprev= T;
end% Enf of time loop
time= toc;
f2 = figure(2)
contourf(x,y,T,'--','ShowText','on')
colorbar
colormap('jet')
text =sprintf("2D Transient Heat Conduction - Implicit - GS NO. of iteration %f Time of simulation %f Sec", iteration, time);
title(text)
set(gca,'YDIR','reverse')
%%Steady State Gauss Seidal
elseif type ==22
iteration =1;
error= 9e9;
tic;
while(error > tol ) %convergence loop
for i=2:nx-1 % Space loop
for j=2:ny-1
term1 = T(i,j-1) + Told(i,j+1);
term2 = T(i-1,j) + Told(i+1,j);
T(i,j) = (1/4)*(term1 +term2);
end
end %end of space loop
error = max(max(abs(Told-T)));
Told =T;
iteration =iteration+1;
end % End of convergence loop
time = toc;
f2= figure(2)
contourf(x,y,T,'--','ShowText','on')
colormap('jet')
text = sprintf("2D ASteady State Heat Conduction - GS No. of iterations %f Time of Simulation %f Sec", iteration,time);
title(text)
set(gca,'YDIR','reverse')
end
end
3.linear_eqn_solver_SOR
function [iteration] = linear_eqn_solver_SOR(T,Told,Tprev,x,y,nx,ny,nt,k1,k2,tol,type)
%%% Transient Implicit (SOR)
w=1.1;
if type == 13
iteration =1;
tic;
for k =1:nt %time loop
error=9e9;
while(error > tol) %convergence loop
for i=2:nx-1 %space loop
for j= 2:ny-1
term1 = (1+(2*k1)+(2*k2))^-1;
V = T(i-1,j) + Told(i+1,j); % Ttop(updated) + Tbottom
H = T(i,j-1) + Told(i,j+1); %Tleft(updated) + Tright
T_GS(i,j) = (Tprev(i,j)*term1) + (k1*term1*H) + (k2*term1*V);
T(i,j) = ((1-w)*Told(i,j)) + (w*T_GS(i,j));
end
end % end of space loop
error = max(max(abs(Told-T)));
Told=T;
iteration = iteration+1;
end %end of convergence loop
Tprev= T;
end % end of time loop
time = toc;
f3 = figure(3)
contourf(x,y,T,'--','ShowText','on')
colorbar
colormap('jet')
text =sprintf("2D Transient Heat Conduction - Implicit - SOR NO. of iteration %f Time of simulation %f Sec", iteration, time);
title(text)
set(gca,'YDIR','reverse')
%%% Steady State SOR
elseif type == 23
iteration = 1;
error=9e9;
tic;
while(error > tol) %convergence loop
for i=2:nx-1 % Space loop
for j=2:ny-1
term1 = T(i,j-1) + Told(i,j+1);
term2 = T(i-1,j) + Told(i+1,j);
T_GS(i,j) = (1/4)*(term1 + term2);
T(i,j) = ((1-w)*Told(i,j)) + (w*T_GS(i,j));
end
end %end of space loop
error= max(max(abs(Told-T)));
Told=T;
iteration= iteration+1;
end % end of convergence loop
time=toc;
f3 = figure(3)
contourf(x,y,T,'--','ShowText','on')
colorbar
colormap('jet')
text =sprintf("2D Steady State Heat Conduction -SOR NO. of iteration %f Time of simulation %f Sec", iteration, time);
title(text)
set(gca,'YDIR','reverse')
end
end
4. linear_eqn_solver_exlpicit
function linear_eqn_solver_explicit(T,Told,nx,ny,nt,x,y,k1,k2)
%%Transient Explicit
iteration= 1;
tic;
for k=1:nt %time loop
for i=2:nx-1 %time loop
for j=2:ny-1
term1=(1+(2*k1)+(2*k2))^-1;
H = Told(i+1,j) + Told(i-1,j);
V = Told(i,j+1) + Told(i,j-1);
T(i,j) = (Told(i,j)*term1) + (k1*H) + (k2*V);
end
end % end if space loop
error = max(max(abs(Told-T)));
Told =T;
iteration=iteration+1;
end %end of time loop
time =toc;
f4=figure(4)
contourf(x,y, T,'--','ShowText','on')
colorbar
colormap('jet')
text=sprintf("2D Transient Heat Conduction- Explicit- No. of iterations %f Time of Simulation %f Sec",iteration,time);
title(text)
set(gca,'YDIR','reverse')
end
MAIN CODE:
clear all
close all
clc
% Inputs
L=1;
nx=11;
x=linspace(0,L,nx);
dx=x(2)-x(1);
W=1;
ny=nx;
y=linspace(0,W,ny);
dy=y(2)-y(1);
prompt = 'type='; %transient case =1 and Steady State =2
type= input(prompt);
if type ==1
prompt = 'CFL=';
CFL = input(prompt); % Fails for Explicit if CFL >0.5
% CFL CFLx+CFLy = alpha*dt/dx^2 + alpha*dt/dy^2=2*alpha*dt/dx^2
else
CFL =0;
end
alpha =1.4;
% Finding dt & nt
dt = CFL*dx^2/(2*alpha);
time=.03*60; % mins
nt = time/dt;
% Error Limit Inputs
error=9e9;
tol=1e-4;
%Boundary Conditions
T=300*ones(nx,ny);
T(:,1) = 400; % left
T(:,end) = 800;% right
T(1,:) = 600; %top
T(end,:) = 900;%bottom
T(1,1) = (400+600)/2;
T(end,end) = (800+900)/2;
T(end,1) = (400+900)/2;
T(1,end) = (600+800)/2;
Told=T;
Tprev=T;
k1=(alpha*dt/(dx^2));
k2=(alpha*dt/(dy^2));
% Calling the Solvers
if type ==1 % Transient case
linear_eqn_solver_jacobi(T,Told,Tprev,nx,ny,nt,x,y,k1,k2,tol,11);
linear_eqn_solver_Gs(T,Told,Tprev,nx,ny,nt,x,y,k1,k2,tol,12);
linear_eqn_solver_SOR(T,Told,Tprev,x,y,nx,ny,nt,k1,k2,tol,13);
linear_eqn_solver_explicit(T,Told,nx,ny,nt,x,y,k1,k2); % Transient Explicit
elseif type ==2 % Steady State Case
linear_eqn_solver_jacobi(T,Told,Tprev,nx,ny,nt,x,y,k1,k2,tol,21);
linear_eqn_solver_Gs(T,Told,Tprev,nx,ny,nt,x,y,k1,k2,tol,22);
linear_eqn_solver_SOR(T,Told,Tprev,x,y,nx,ny,nt,k1,k2,tol,23);
endAlgorithm:
- Firstly we have defined the size of the 2D rectangle onto which we want to calculate the temperature contour. Then divided the geometry into equal grid points or nodes of equal length and width. After that, we calculated the node size, dx & dy which are equal in magnitude that is dx=dy.
- Then we accept the type of simulation we want to run that is whether a steady-state or transient state. Pressing"1" signifies transient state for which we take in another user input for CFL number and pressing "2" means steady state.
- From CFL number time step dt is automatically calculated. Here the total simulation time is set to 0.03 minutes and from dt and total time, we calculate the total number of time steps,nt.
- In Step 4 we define the boundary conditions which are as follow: Top boundary = 600 k Bottom Boundary = 900 k Left Boundary = 400 k Right Boundary = 800 k
- Then as per the type of simulation, the calculations are being done by the 4 user defined functions namely , linear_eqn_solver_jacobi(),linear_eqn_solver_Gs(),linear_eqn_solver_SOR(),linear_eqn_solver_explicit() and print the graphs along with the iteration number and simulation time.
- Among this functions for the first 3 linear_eqn_solver_jacobi(),linear_eqn_solver_Gs(),linear_eqn_solver_SOR() solves the Transient state and Steady state is solved by next 3 mention under the type 2.
The main difference between these solvers is that the discretized equations are solved by different iterative methods like Jacobi, Gauss-seidel, or Successive overrelaxation methods.
OUTPUT:
1.Transient-State Analysis:
EXPLICIT
2.STEADY-STATE.
Discussion of Results:
1. In the case of both implicit Transient and Steady-state SOR converges faster than others. SOR>GS>Jacobi. This can will be understood by the iteration number and simulation time.
2. Between Implicit and Explicit methods of Transient, there is a limitation of explicit methods on the CFL number. The solutions blow up for CFL greater than 0.5
3. Between the fastest method of the steady-state and transient that is SOR, steady-state converges faster than transient.