2d Heat Conduction 1

Solving the Steady and Unsteady 2D Heat Conduction Problem

-By Mandeep P Shetty

In the current project, I have written a Matlab code to solve the steady-state and transient heat conduction equation in a 1m*1m square plate domain. I have considered the domain to contain air with a thermal diffusivity value of 1m^2/s. Dirichlet boundary conditions are specified at the four edges of the plate as:

1. Left edge : 400K
2. Top edge : 600K
3. Right edge : 800K
4. Bottom edge : 900K

The aim of this project is to simulate the steady state temperature distribution within the plate by solving both the steady state and transient heat conduction equations using finite difference method and to understand the nuances in solving the steady state (elliptic) and transient (parabolic) heat conduction equations. The square domain is discretized with a structured mesh having 16 nodes along the x and the y axis. The time domain is discretized with time step, dt = 1e-3. The error tolerance value is taken to be 1e-4 and here I am running the simulation (transient) for 5000 time steps.

The 3D anisotropic heat conduction governing equation for a solid can be given as :

`ρC_P ∂T/∂t= ∂/∂x (k_x ∂T/∂x)+ ∂/∂y (k_y ∂T/∂y)+ ∂/∂z(k_z ∂T/∂z)`

Where T(x,y,z,t) is the temperature field in the solid, is the density of the material, is the specific heat, k(x,y,z) is the conductivity in different direction (anisotropic conductivity), t is the time, and (x,y,z) are the space dimension.

Simplifying this equation by making certain assumptions like 2D and isotropic conductivity, we get:

`∂T/∂t+ α((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))=0`

Which is the 2D isotropic heat conduction equation where `α= k/(ρC_p )` is the thermal diffusivity of the material. This equation is parabolic in nature.

The steady-state heat conduction equation can be given as:

`((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))=0`

Which is elliptic in nature. Elliptic equations take information from all direction and so they don’t have a specific marching direction.

Solving these PDEs by discretizing using FDM

Discretizing the steady-state equation by using the central differencing technique for second-degree space derivative

`((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))=0`

(Note: We can use second order central differencing scheme to discretize spatial derivatives to take care of the elliptic nature of the equation)

`(T_L-〖2T〗_P+ T_R)/(∆x^2 )+ (T_T-〖2T〗_P+ T_B)/(∆y^2 )=0`

`(2T_P)/(∆x^2 )+ (2T_P)/(∆y^2 )= (T_L+ T_R)/(∆x^2 )+ (T_T+ T_B)/(∆y^2 )`

`T_p (2(∆x^2+ ∆y^2 )/(∆x^2* ∆y^2 ))= (T_L+ T_R)/(∆x^2 )+ (T_T+ T_B)/(∆y^2 )`

`T_P= 1/w ((T_L+ T_R)/(∆x^2 )+ (T_T+ T_B)/(∆y^2 ))`

Where `w= (2(∆x^2+ ∆y^2 )/(∆x^2* ∆y^2 ))`

Applying this equation at ever nodal points in the domain results in a system of linear equations and we can employ one of the iterative methods to solve these equations.

The transient heat conduction equation is parabolic in nature. The parabolic equation has a marching direction, (which is the time dimension here). This equation can be numerically solved using the FDM by implementing implicit or explicit methods.

To employ the explicit method we utilizing the Forward-(in)Time-Central-(in)-Space (FTCS) scheme, which results in the discretized transient heat conduction equation as:

`∂T/∂t+ α((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))=0`

`(T_P^(n+1)- T_P^n )/∆t+ α〖((T_R-2T_P+T_L)/(∆x^2 )+(T_T-2T_P+T_B)/(∆y^2 ))〗^n=0`

`T_P^(n+1)= T_P^n+ α*∆t*〖((T_R-2T_P+ T_L)/(∆x^2 )+ (T_T-2T_P+T_B)/(∆y^2 ))〗^n`

The method is explicit because it describes the unknown dependent variable at time step (n+1) in terms of know variables at time step (n). Therefore, the dependent variable can be solved directly. The explicit method is conditionally stable. Von Neumann stability analysis reveals the stability criterion to be CFL < 0.25, (across one dimension) or CFL < 0.5 across both the dimension.

Utilizing the implicit Backward-(in)Time-Central-(in)Space (BTCS) scheme, which is first order in time and second order in space. The discretized transient heat conduction equation can be derived as:

`∂T/∂t+ α((∂^2 T)/(∂x^2 )+ (∂^2 T)/(∂y^2 ))=0`

`T_P^(n+1)= T_P^n+ α*∆t*〖((T_R-2T_P+ T_L)/(∆x^2 )+ (T_T-2T_P+T_B)/(∆y^2 ))〗^(n+1)`

Let, `k_1= (α*∆t)/(∆x^2 )` and `k_(2 )= (α*∆t)/(∆y^2 )`

`T_P^(n+1) (1+2k_1+2k_2 )- k_1 T_L^(n+1)-k_1 T_R^(n+1)-k_2 T_T^(n+1)-k_2 T_B^(n+1)=T_P^n`

Let,

`a= 1⁄((1+2k_1+2k_2))`

`b=k_1⁄((1+2k_1+2k_2))`

`c=k_2⁄((1+2k_1+2k_2))`

`∴ T_P^(n+1)=a*T_P^n+b*(T_L^(n+1)+T_R^(n+1) )+c*(T_T^(n+1)+T_B^(n+1) )`

The method is implicit because it describes the unknown dependent variable at time step (n+1) in terms of know and unknow variables at time step (n) and (n+1) respectively. Therefore, the implicit scheme results in a set of linear equations which we solve here using iterative methods like

1. Jacobi method
2. Gauss-Seidal (GS) method
3. Successive Over Relaxation (SOR) method, with a relaxation factor, .

Implicit methods are unconditionally stable and so we can use larger time steps.

Results

Solving the discretized steady state heat equation:

1. Jacobi method

2. GS method

3. SOR method

Solving the discretized transient heat equation (till steady state) by the implicit scheme:

1. Jacobi method

2. GS method

3. SOR method

Solving the discretized transient heat equation (till steady state) by the explicit scheme:

Solution verification using Ansys APDL:

Conclusions

Discretizing the steady state equation results in a system of linear equations which can be solved by employing an iterative technique with appropriate boundary conditions. The transient equations can be solved by utilizing either the implicit or explicit scheme. The implicit scheme results in a set of linear equations which can be solved utilizing an iterative scheme and the explicit scheme results in equations which can be solved directly at every time step for every nodal value.

Among the 3 iterative methods used here to solve the set of linear equations got from discretizing the steady-state heat equation, SOR method, with the right value for, results in the least number of iterations, while employing the Jacobi method results in the most number of iterations.

The parabolic transient heat conduction equation can be solved using either the implicit or the explicit time discretization schemes. Employing the explicit discretization scheme gives the unknown temperature at a node point in terms of know values of its neighbors, which can be solved directly without using any iterative methods. The explicit scheme is conditionally stable, which means we have to use smaller time steps. The condition for stability can be given in terms of the CFL criterion. This method is useful when we need a transient solution but takes longer time to reach the steady state solution.

The implicit discretization scheme results in a set of linear equations which are solved using each of the 3 iterative methods mentioned earlier. Again, the SOR method (with the proper value for ), results in the least number of iterations and the Jacobi method results in the highest number of iterations. The implicit solution is unconditionally stable, which means we can utilize large time steps. This results in the solution converging to the steady-state solution at a faster rate but the quality of the transient part of the solution is not dependable.

Matlab Code

Main Code :

clear all
close all
clc
%Solver for 2D heat conduction equation in a square plate.
%The equation we are solving for is the heat conduction equaiton:
%dT/dt = alpha(d^2T/dx^2+d^2T/dy^2), where the material is considered to
%have isotropic conductivity and the heat source terms are absent.

%----------------%%Inputs%%---------------%
%Dimensions of the 2D square plate
L = 1; %Dimension of the unit square plates. Length = width
W = L;

% Material conditions (taken for stainless steel 304)
% roh = 7900; %Density of the material, kg/m^3
% Cp = 502.416; %Specific heat of the material, J/kg-K
% k = 14.9; %Conductivity of the material (Taken to be isotoropic),
global alpha;
alpha = 1.1; %k/(roh*Cp) %Thermal diffusivity

%Number of grid points
global nx;
global ny;
global dt;

nx = 16;
ny = nx; %Number of node points in x and y direction are equal.
dt = 1e-3; %Time steps

%Initial condition
Tini = 298; %Initial condition.(Consider the inital plate temparature to be 298K)

%Boundary Conditions
Tl = 400; %Left plate edge temperature
Tt = 600; %Top plate edge temperature
Tr = 800; %Right plate edge temperature
Tb = 900; %Bottom plate edge temperature

%Grid sizing
global dx
global dy

dx = (L/(nx-1)); %Space step size in x direction
dy = (W/(ny-1)); %Space step size in y direction

global error_lim
error_lim = 9e9; %Initializing the error term.

global tol;
tol = 1e-4; %Tolarance ie if the error is less than tol we can take it as converged result.

%Mesh generation

T = Tini*ones(nx,ny); %Mesh creation and initilization. For the Transient explicit method the initialization values act as the true iniial ocndition
%For Transient explicit method and the steady state ellliptic equation
%method the initialization values act as the inital Guess values for the
%itteratioin to begin.

%Need to consider x y coordinates of the grid points not row and columns of
%a matrix

T(:,1) = Tl; %Applying boundary conditions, Dirchlit BC
T(1,:) = Tb; %Applying boundary conditions, Dirchlit BC
T(:,end) = Tr; %Applying boundary conditions, Dirchlit BC
T(end,:) = Tt; %Applying boundary conditions, Dirchlit BC

% T(1,1) = (Tb + Tl)/2; %Corner value, left-top. The corner corrctions dosent matter.
% T(1,end) = (400 + 600)/2; %Corner value, right-top. The corner corrctions dosent matter.
% T(end,end) = (Tt + Tr)/2; %Corner value, right-bottom. The corner corrctions dosent matter.
% T(end, 1) = (Tb + Tr)/2; %Corner value, left-bottom. The corner corrections dosent matter.

TS = 5000; %Number of time steps.

flag = menu('Steady state or Transient solver?','Transient','Steady State');
if flag == 1
list = {'Jacobi','Gauss-Seidel','SOR','Explicit'};
num_met = listdlg('PromptString','Choose the solver:','ListString',list,'SelectionMode','single','Name','Solvers');
else
list = {'Jacobi','Gauss-Seidel','SOR'};
num_met = listdlg('PromptString','Choose the solver:','ListString',list,'SelectionMode','single','Name','Solvers');
end

switch num_met
case 1
[Tsol, iter,time_steps, error, header] = jacobi(T,flag, TS);
case 2
[Tsol, iter,time_steps, error, header] = gausss(T,flag, TS);
case 3
[Tsol, iter,time_steps, error, header] = sor(T,flag, TS);
case 4
[Tsol, time_steps, header] = explicit(T, flag, TS);
end

cfl = (alpha*dt)/(dx^2)+ (alpha*dt)/(dy^2); %CFL number value for 2D problem

%Plotting
figure(2)
[M,c] = contourf(Tsol);
colorbar;
title(header)
clabel(M,c);
colormap(jet);

Function definitions:

jacobi():

function [T, jacobi_iter,time_steps, error, header] = jacobi(T, flag, TS)
global tol dx dy nx ny dt alpha error_lim;
error = error_lim;
time_steps = 0;
jacobi_iter = 0;
Told = T; %Told is the value from last itteration.

%Steady state
if flag == 2

a = ( (dx^2*dy^2)/(2*(dx^2 + dy^2)) );
while(error > tol)
jacobi_iter = jacobi_iter + 1;
for i = 2:(nx - 1)
for j = 2:(ny - 1)
T(i,j) = a *( ((Told(i,j-1)+Told(i,j+1))/(dx^2)) + ((Told(i+1,j)+Told(i-1,j))/(dy^2)) ); %Here T is a matrix, not grid point.Matrix are indexed
%by (row_number, column_number) and grid points have x
%and y coordinates. So here think in terms of rows and
%columns and not x and y coordinates ie (i,j) is row
%and column value not x and y coordinates

end
end
error = max(max(abs(T - Told)));
Told = T;
end
end

%Transient Implicit method
if flag == 1

k1 = (alpha*dt)/dx^2;
k2 = (alpha*dt)/dy^2;
term1 = (1 + 2*k1 + 2*k2)^-1;
term2 = k1*term1;
term3 = k2*term1;

%do = 1;
%while(do) %Time loop has to continue till the steady state is reached.Check condition to enter the time loop. Using a do-while structure so that atleast once the loop is entered before the condition is checked.
for time_step = 1:TS
Tprev = T; %Tprev is the value from previous time step.
time_steps = time_steps + 1;
error = error_lim; %This is required as the error variable gets changed during each time step and without this assignment for the next
%time step will take the changed value.

while(error > tol)
jacobi_iter = jacobi_iter + 1; %To get the total number of itteration (cummulative itterations is itterations @ each time step.)
Told = T;
for i = 2:(ny - 1)
for j = 2:(nx - 1)
T(i,j) = term1*Tprev(i,j) + term2*(Told(i,j-1)+Told(i,j+1)) + term3*(Told(i+1,j)+Told(i-1,j));
end
end
error = max(max(abs(T - Told)));

end

% if(max(max(abs(Tprev-T))) < 1e-3)
% do = 0;
% end
end
end
header = sprintf('Jacobi method; %d iterations; %d time steps',jacobi_iter, time_steps);

end

gauss():

function [T, gauss_iter, time_steps,error, header] = gausss(T, flag, TS)
global tol dx dy nx ny dt alpha error_lim;
error = error_lim;
time_steps = 0;
gauss_iter = 0;

%Steady state
if flag == 2
Told = T;
a=( (dx^2*dy^2)/(2*(dx^2 + dy^2)) );
while(error > tol)
gauss_iter = gauss_iter + 1;
for i = 2:(nx - 1)
for j = 2:(ny - 1)
T(i,j) = a *( ((T(i,j-1)+T(i,j+1))/(dx^2)) + ((T(i+1,j)+T(i-1,j))/(dy^2)) ); %Here T is a matrix, not grid point.Matrix are indexed
%by (row_number, column_number) and grid points have x
%and y coordinates. So here think in terms of rows and
%columns and not x and y coordinates ie (i,j) is row
%and column value not x and y coordinates
end
end
error = max(max(abs(T - Told)));
Told = T;
end
end

%Transient Implicit method
if flag == 1

k1 = (alpha*dt)/dx^2;
k2 = (alpha*dt)/dy^2;
term1 = (1 + 2*k1 + 2*k2)^-1;
term2 = k1*term1;
term3 = k2*term1;
Told = T; %Told is the value from last itteration.

%do = 1;
%while(do) %Time loop has to continue till the steady state is reached.Check condition to enter the time loop. Using a do-while structure so that atleast once the loop is entered before the condition is checked.
for time_step = 1:TS
Tprev = T; %Tprev is the value from previous time step.
time_steps = time_steps + 1;
error = error_lim; %This is required as the error variable gets changed during each time step and without this assignment for the next
%time step will take the changed value.
while(error > tol)
gauss_iter = gauss_iter + 1; %To get the total number of itteration (cummulative itterations run @ each time step.)
for i = 2:(ny - 1)
for j = 2:(nx - 1)
T(i,j) = term1*Tprev(i,j) + term2*(T(i,j-1)+T(i,j+1)) + term3*(T(i+1,j)+T(i-1,j));
end
end
error = max(max(abs(T - Told)));
Told = T;
end
% if(max(max(abs(Tprev-T))) < 1e-3)
% do = 0;
% end
end
end
header = sprintf('GS method; %d iterations; %d time steps',gauss_iter, time_steps);

end

sor():

function [T, sor_iter, time_steps,error, header] = sor(T, flag, TS) %This SOR method is used to enhance the GS method.
global tol dx dy nx ny dt alpha error_lim;
error = error_lim;
time_steps = 0;
sor_iter = 0;
omega = 1.2;
global test; %Testing

%Steady state
if flag == 2
Told = T;
a=( (dx^2*dy^2)/(2*(dx^2 + dy^2)) );
while(error > tol)
sor_iter = sor_iter + 1;
for i = 2:(nx - 1)
for j = 2:(ny - 1)
T(i,j) = a *( ((T(i,j-1)+T(i,j+1))/(dx^2)) + ((T(i+1,j)+T(i-1,j))/(dy^2)) ); %Here T is a matrix, not grid point.Matrix are indexed
%by (row_number, column_number) and grid points have x
%and y coordinates. So here think in terms of rows and
%columns and not x and y coordinates ie (i,j) is row
%and column value not x and y coordinates
T(i,j) = Told(i,j)*(1-omega) + omega*(T(i,j));
end
end
error = max(max(abs(T - Told)));
Told = T;
end
end

%Transient Implicit method
if flag == 1

k1 = (alpha*dt)/dx^2;
k2 = (alpha*dt)/dy^2;
term1 = 1/(1 + 2*k1 + 2*k2);
term2 = k1*term1;
term3 = k2*term1;
Told = T; %Told is the value from last itteration.

%do = 1;
%while(do) %Time loop has to continue till the steady state is reached.Check condition to enter the time loop. Using a do-while structure so that atleast once the loop is entered before the condition is checked.
for time_step = 1:TS;
Tprev = T; %Tprev is the value from previous time step.
time_steps = time_steps + 1;
error = error_lim; %This is required as the error variable gets changed during each time step and without this assignment for the next
%time step will take the changed value.
test_iter = 0;
while(error > tol)
sor_iter = sor_iter + 1; %To get the total number of itteration (cummulative itterations run @ each time step.)
test_iter = test_iter + 1;
for i = 2:(ny - 1)
for j = 2:(nx - 1)
test(time_steps,test_iter) = 1; %Testing
T(i,j) = term1*Tprev(i,j) + term2*(T(i,j-1)+T(i,j+1)) + term3*(T(i+1,j)+T(i-1,j));
T(i,j) = Told(i,j)*(1-omega) + omega*(T(i,j));
end
end
error = max(max(abs(T - Told)));
Told = T;
end
% if(max(max(abs(Tprev-T))) < 1e-3)
% do = 0;
% end
end
end
header = sprintf('SOR method; %d iterations; %d time steps',sor_iter, time_steps);

end

explicit():

function [T, time_steps, header] = explicit(T, flag, TS)
global dx dy nx ny dt alpha;
%t = 300*10^3;

%Transient explicit method
if flag == 1

time_steps = 0;


do = 1;
%while(do) %Time loop has to continue till the steady state is reached.Check condition to enter the time loop. Using a do-while structure so that atleast once the loop is entered before the condition is checked.
for time_step = 1:TS;
time_steps = time_steps + 1;
Told = T;

for i = 2:(ny - 1)
for j = 2:(nx - 1)

T(i,j) = (Told(i,j)) + ((alpha*dt)*(((Told(i,j+1)-(2*Told(i,j))+Told(i,j-1))/dx^2) + ((Told(i+1,j)-(2*Told(i,j))+Told(i-1,j))/dy^2)));
end
end

if(max(max(abs(Told-T))) < 1e-3)
do = 0;
end
end
end
header = sprintf('Explicit method; %d time steps',time_steps);
end