Solving the steady and unsteady 2D heat conduction problem
Aim - To solve 2D Transient Heat conduction Equation
Theory -
- For an unsteady state the 2-D heat conduction equation is represented as
∂T∂t=α⋅(∂2T/∂x2+∂2T/∂y2)
- For a steady state condition, the heat conduction equation is represented as
(∂2T/∂x2+∂2T/∂y2)=0
Because these equations are partial differential equations, numerical discretization techniques can be utilised to solve them numerically.
The problem's boundary conditions are as follows -
The equations can be solved in either implicit or explicit form for transient analysis.
The ∂T∂t">∂T/∂t is discretized using the forward time marching scheme as
Tn+1i,j−Tni,j / Δt
Δt
The subscript (i,j) denotes the nodes and the superscript represents the time step.
Since the right hand side of the equation is a second order partial derivative, the central differencing scheme is used for discretizing the terms as follows :
∂2T/∂x2=Ti-1,j-2⋅Ti,j+Ti+1,jΔx2">∂2T∂x2=Ti−1,j−2⋅Ti,j+Ti+1,jΔx2∂2�∂�2=��-1,�-2⋅��,�+��+1,�Δ�2
∂2T/∂y2=Ti,j-1-2⋅Ti,j+Ti,j+1Δy2">∂2T∂y2=Ti,j−1−2⋅Ti,j+Ti,j+1Δy2
Depending on which time step these second order derivated are computed, they can be termed as implicit and explicit.
The subscripts are replaced as shown in the diagram for clarity -
Explicit Form -
The equation has the explicit form when computed at the prior time step n
Tn+1p−TnpΔt=α⋅[TnL−2⋅TnP+TnRΔx2+TnT−2⋅TnP+TnBΔy2]
Rearranging,
Tn+1p=Tnp+α⋅(Δt)⋅[TnL−2⋅TnP+TnRΔx2+TnT−2⋅TnP+TnBΔy2]
Lets Consider
α⋅(Δt)Δx2=k1
and
α⋅(Δt)Δy2=k2">α⋅(Δt)Δy2=k2
α⋅(Δt)Δy2=k2">Substituting these values in the equation,
Tpn+1=Tpn+k1⋅(TLn-2⋅TPn+TRn)+k2⋅(TTn-2⋅TPn+TBn)">Tn+1p=Tnp+k1⋅(TnL−2⋅TnP+TnR)+k2⋅(TnT−2⋅TnP+TnB)
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">Tn+1p=Tnp⋅[1−2⋅k1−2⋅k2]+k1⋅(TnL+TnR)+k2⋅(TnT+TnB)
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">These equations can be used directly to compute the LHS, since all the values in the RHS is known
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">Implicit Form -
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">If the of the second order derivative terms are computed at the current time step (n+1), then the equation is of the implicit form.
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">Tn+1p−TnpΔt=α⋅[Tn+1L−2⋅Tn+1P+Tn+1RΔx2+Tn+1T−2⋅Tn+1P+Tn+1BΔy2]
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">Rearranging,
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">Tn+1p=Tnp+α⋅(Δt)⋅[Tn+1L−2⋅Tn+1P+Tn+1RΔx2+Tn+1T−2⋅Tn+1P+Tn+1BΔy2]
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">Lets Consider
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">α⋅(Δt)Δx2=k1
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">and
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">α⋅(Δt)Δy2=k2
Tpn+1=Tpn⋅[1-2⋅k1-2⋅k2]+k1⋅(TLn+TRn)+k2⋅(TTn+TBn)">Subtituting these values in the equation,
Tpn+1=Tpn+k1⋅(TLn+1-2⋅TPn+1+TRn+1)+k2⋅(TTn+1-2⋅TPn+1+TBn+1)">Tn+1p=Tnp+k1⋅(Tn+1L−2⋅Tn+1P+Tn+1R)+k2⋅(Tn+1T−2⋅Tn+1P+Tn+1B)
Tpn+1=Tpn+k1⋅(TLn+1+TRn+1)+k2⋅(TTn+1+TBn+1)1+2⋅k1+2⋅k2">Tn+1p=Tnp+k1⋅(Tn+1L+Tn+1R)+k2⋅(Tn+1T+Tn+1B)1+2⋅k1+2⋅k2
Tpn+1=Tpn+k1⋅(TLn+1+TRn+1)+k2⋅(TTn+1+TBn+1)1+2⋅k1+2⋅k2">
Because the RHS of the implicit form of the problem is still unknown, the system of equations is transformed to matrix form and then solved. Iterative solvers can also be used to solve the implicit form in different ways. The iterative solvers create guesses for the unknown variables and then compute the error value; the guess values are adjusted using a transformation matrix and a constant until the residual error falls below a set limit. The residual error, also known as the tolerance, is specified by the user prior to the start of the iterations.
Tpn+1=Tpn+k1⋅(TLn+1+TRn+1)+k2⋅(TTn+1+TBn+1)1+2⋅k1+2⋅k2">The iterative solvers' various methodologies are listed below -
Tpn+1=Tpn+k1⋅(TLn+1+TRn+1)+k2⋅(TTn+1+TBn+1)1+2⋅k1+2⋅k2">1. Jacobi Method
Tpn+1=Tpn+k1⋅(TLn+1+TRn+1)+k2⋅(TTn+1+TBn+1)1+2⋅k1+2⋅k2">The Jacobi iterative method is an iterative procedure that is used to find solutions to systems of linear equations. This method employs an approximation (guess) value, which is subsequently used to solve the complete problem. The process is iterated until it converges. The Jacobi technique maintains the values of Tni obtained in the nth iteration until the full (n+1)th iteration is calculated.
2. Gauss Siedel method
With the Gauss-Seidel method, we use the new values Tin+1">Tn+1i
as soon as they are known. For example,
once we have computed T1n+1">T1n+1 from the first equation, its value is then used in the second equation to obtain the new T2n+1">Tn+12
and so on. In simpler words the difference between the Gauss–Seidel and Jacobi methods is that the Jacobi method uses the values obtained from the previous step while the Gauss–Seidel method always applies the latest updated values during the iterative procedures.
3. Successive over relaxation (SOR) method
This iterative method uses a value ω">ω
(omega) to guess the value for next iteration. By employing the below mentioned formula the values are computed.
Tin+1=Tin⋅(1-ω)+ω⋅(TGS)">Tn+1i=Tni⋅(1−ω)+ω⋅(TGS)
In this TGS">TGS
is the value computed using the gauss-siedel iteration method formula. In this manner the gauss siedel method is extended to the relaxation methods.
- Transient Analysis using the explicit scheme
Iterative solvers are not employed because the RHS is totally understood. The temperature at each node is computed in this scenario due to numerical discretization. The MATLAB programme is shown below. The sole parameter that must be monitored is the CFL, which is represented in the programme by k1 and k2. For the analysis to be stable, the sum of both should be less than or equal to 0.5. The CFL value in this situation is 0.4536
- CODE
close all
clear all
clc
% input parameters
nx=10; % number of nodes in the x-direction
ny=10; % number of nodes in the y-direction
dt=0.002 ; % timestep
alpha = 1.4; % thermal diffusivity
nt = 500; % number of timesteps
T = 300*ones(nx,ny);
% initial conditions
for i = 1:nx
T(1,i) = 600;
T(ny,i)=900;
end
for j=1:ny
T(j,1) = 400;
T(j,nx) = 800;
end
T(1,1)= (T(1,2)+ T(2,1))/2;
T(1,nx) = (T(1,nx-1)+T(2,nx))/2;
T(ny,1) = (T(ny-1,1) + T(ny,2))/2;
T(ny,nx) = (T(ny,nx-1)+T(ny-1,nx))/2;
%creating a 2D mesh
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
Told = T;
%CFL values for stability
k1 = alpha*(dt/(dx^2));
k2 = alpha*(dt/(dy^2));
CFL = k1+k2;
% Time loop
for k = 1:nt
for i=2:nx-1
for j=2:ny-1
T(i,j) = (Told(i,j) + (k1*(Told(i+1,j)-(2*Told(i,j))+Told(i-1,j))) + (k2*(Told(i,j-1)-(2*Told(i,j)) +Told(i,j+1))));
end
end
Told = T;
% plotting the heatmap for each iteration.
colormap(jet);
A=contourf(x,y,T);
clabel(A);
set(gca,'YDir','reverse');
xlabel('X-Direction');
ylabel('Y-Direction');
grid on ;
title (' Transient Analysis of 2D-heat conduction equation using explicit scheme');
B = sprintf(' Number of timesteps = %d, timestep (dt) = %d',k,dt);
subtitle(B);
pause(0.0001);
endThe graphic is set up in a time loop to detect changes in the temperature distribution at each timestep. Even before the prescribed number of timesteps, the solution reaches a steady state. If the timestep size (dt) is changed, the time required to attain the steady state will likewise change. The solution reaches a stable state at roughly 100 to 150 iterations for the present timestep of 0.002. Aside from the timestep, node size and alpha (thermal diffusivity) also influence the time required to reach a steady state.
- OUTPUT -
- Transient Analysis using the implicit scheme
Because the RHS of the equation is unknown in this case, iterative solutions are used. The three iterative approaches mentioned above are employed in the MATLAB programme below.
- CODE
close all
clear all
clc
% input parameters
nx=10; % number of nodes in the x-direction
ny=10; % number of nodes in the y-direction
dt=0.0001 ; % timestep
alpha = 1.4; % thermal diffusivity
T = 300*ones(ny,nx); % initial guess values for temperature
% initial conditions
for i = 1:nx
T(1,i) = 600;
T(ny,i)=900;
end
for j=1:ny
T(j,1) = 400;
T(j,nx) = 800;
end
T(1,1)= (T(1,2)+ T(2,1))/2;
T(1,nx) = (T(1,nx-1)+T(2,nx))/2;
T(ny,1) = (T(ny-1,1) + T(ny,2))/2;
T(ny,nx) = (T(ny,nx-1)+T(ny-1,nx))/2;
%creating a 2D mesh
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
Told = T;
% selection of the types of solver
iteration = 3;
% tolerance for convergence of temperture
tol = 1e-4;
k1 = alpha*(dt/(dx^2));
k2 = alpha*(dt/(dy^2));
Tprev = T;
iter = 1 ;
%Joacobi method
if iteration ==1
for k=1:1500 % time loop
error = 9e9;
while(error > tol) % convergence loop
for i=2:ny-1
for j=2:nx-1
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
T(i,j) = (term1*Tprev(i,j)) + (term2*(Told(i+1,j)+Told(i-1,j))) + (term3*(Told(i,j-1)+Told(i,j+1)));
end
end
error = max(max(abs(Told-T))); % calculating maximum error value
Told = T; % updating the old temperature values
iter = iter +1;
end
Tprev = T;
end
end
%Gauss-siedel method
if iteration ==2
for k=1:1500 % time loop
error = 9e9;
while(error > tol) % convergence loop
for i=2:ny-1
for j=2:nx-1
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
T(i,j) = (term1*Tprev(i,j)) + (term2*(Told(i+1,j)+T(i-1,j))) + (term3*(T(i,j-1)+Told(i,j+1)));
end
end
error = max(max(abs(Told-T))); % calculating error values
Told = T; % updating the old temperature values
iter = iter +1;
end
Tprev = T;
end
end
%SOR method
if iteration ==3
omega = 1.013; % Value of omega
for k=1:1500 % time loop
error = 9e9;
while(error > tol) % convergence loop
for i=2:ny-1
for j=2:nx-1
term1 = (1/(1+(2*k1)+(2*k2)));
term2 = k1*term1;
term3 = k2*term1;
T(i,j) = ((1-omega)*Told(i,j)) + omega*((term1*Tprev(i,j)) + (term2*(Told(i+1,j)+T(i-1,j))) + (term3*(T(i,j-1)+Told(i,j+1))));
end
end
error = max(max(abs(Told-T))); % calculating error value
Told = T; % updating the old temperature values
iter = iter +1;
end
Tprev = T;
end
end
% plotting the temperature distribution
figure(1)
colormap(jet);
A=contourf(x,y,T);
clabel(A);
set(gca,'YDir','reverse');
xlabel('X-Direction');
ylabel('Y-Direction');
grid on ;
title(' Transient analysis of 2D heat equation using implicit scheme(Jacobi method)');
B = sprintf('number of timestep = %d , Number of iteration = %d',k,iter);
subtitle(B);
All three iterative solvers are employed in the code. The 'iteration' argument is used to determine the solver solver. The iterative solver iteration=1 employs the Jacobi approach. Iteration=2 is the gauss siedel method iterative solver, and iteration=3 is the successive over relaxation technique iterative solver.
The calculations and temperature distributions for various approaches are provided below:
Jacobi method
Formula -
Tn+1p=T(prev)np+k1⋅(T(old)n+1R+T(old)n+1L)+k2⋅(T(old)n+1T+T(old)n+1B)1+2⋅k1+2⋅k2
T(old) represents the previous iteration's guess values. These values are updated after each loop, not during it. The temperature values from the preceding timestep are represented by T(prev). It is constant at every given time step. Other elements to consider are the CFL considerations.
k1=α⋅ΔtΔx2">k1=α⋅ΔtΔx2 and k2=α⋅ΔtΔy2">k2=α⋅Δtk2=α⋅ΔtΔy2">Δy2
- OUTPUT PLOT
α⋅(Δt)Δy2=k2">�⋅(Δ�)Δ�2=�2
Gauss-siedel method
Formula
Tn+1p=T(prev)np+k1⋅(T(old)n+1R+Tn+1L)+k2⋅(Tn+1T+T(old)n+1B)1+2⋅k1+2⋅k2
In this scenario, the left and top point values are calculated during the iteration, and the guess value is changed within an iteration, but in Jacobi, the old guess values are used, and therefore the guess value is not updated inside an iteration.
- OUTPUT PLOT
Successive over relaxation method
Formula -
Tn+1p=(1−ω)⋅T(old)np+ω⋅T(prev)np+k1⋅(T(old)n+1R+T(old)n+1L)+k2⋅(T(old)n+1T+T(old)n+1B)1+2⋅k1+2⋅k2)
- OUTPUT PLOT
Comparison of no. of iterattions between different methods
This data is based on running the iterative solution for 1500 timesteps at 0.0001 second intervals. The number of iterations required for convergence is on the rise. Therefore, Jacobi method > Gauss-siedel method > SOR method
Steady state analysis of 2D heat conduction problem
For a steady state condition, the heat conduction equation is represented as,
(∂2T/∂x2+∂2T/∂y2)=0
- CODE
close all
clear all
clc
% input parameters
nx=10; % number of nodes in the x-direction
ny=10; % number of nodes in the y-direction
dt=1e-4 ; % timestep size
alpha = 1.4; % thermal diffusivity
T = 300*ones(nx,ny);
% initial conditions
for i = 1:nx
T(1,i) = 600;
T(ny,i)=900;
end
for j=1:ny
T(j,1) = 400;
T(j,nx) = 800;
end
T(1,1)= (T(1,2)+ T(2,1))/2;
T(1,nx) = (T(1,nx-1)+T(2,nx))/2;
T(ny,1) = (T(ny-1,1) + T(ny,2))/2;
T(ny,nx) = (T(ny,nx-1)+T(ny-1,nx))/2;
%creating a 2D mesh
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
Told = T;
% selection of the iterative solver method
iteration = 1;
% tolerance for convergence of temperture
tol = 1e-4;
error = 1e4;
k = ((2*(dy^2 + dx^2))/(dy^2*dx^2));
%Jacobi method
if iteration == 1
iter = 1 ;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term_1 = ((Told(i+1,j) + Told(i-1,j))/(k*dx^2));
term_2 = ((Told(i,j-1) + Told(i,j+1))/(k*dy^2));
T(i,j) = term_1 + term_2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter +1;
end
end
%Gauss-siedel method
if iteration == 2
iter = 1 ;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term_1 = ((Told(i+1,j) + T(i-1,j))/(k*dx^2));
term_2 = ((T(i,j-1) + Told(i,j+1))/(k*dy^2));
T(i,j) = term_1 + term_2;
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter +1;
end
end
%succesive over-relaxation method
omega=1.5;
if iteration == 3
iter = 1 ;
while(error > tol)
for i=2:nx-1
for j=2:ny-1
term_1 = ((Told(i+1,j) + T(i-1,j))/(k*dx^2));
term_2 = ((T(i,j-1) + Told(i,j+1))/(k*dy^2));
T(i,j) = ((Told(i,j)*(1-omega)) + omega*(term_1 + term_2));
end
end
error = max(max(abs(Told-T)));
Told = T;
iter = iter +1;
end
end
% plotting the temperature distribution
figure(1)
colormap(jet);
A=contourf(x,y,T);
clabel(A);
set(gca,'YDir','reverse');
xlabel('X-Direction');
ylabel('Y-Direction');
grid on ;
title('Steady State Analysis of 2D Heat conduction');
B = sprintf('Iterative method = %d , Number of iterations = %d , Error = %d',iteration,iter,error);
subtitle(B);To solve the PDE, this code employs three iterative solvers. The 'iteration' argument is used to determine the solver solver. The iterative solver iteration=1 employs the Jacobi approach. Iteration=2 is the gauss siedel method iterative solver, and iteration=3 is the successive over relaxation technique iterative solver.
The temporal loop is not identified since it is a stable state, as opposed to the unsteady analysis code. Because it is a steady-state equation, the formula changes.
When the three iterative solutions are applied, the formulas and output graphic are as displayed.
- Jacobi method
Formula
Where
- Gauss-Siedel method
- Successive over relaxation (SOR) method
Formula
OUTPUT PLOT -
Comparison of no. of iterattions between different methods
It may be concluded that the SOR approach requires the fewest iterations. In the Jacobi approach, the maximum number of iterations is observed. As a result, the trend in terms of the number of iterations required for convergence is, Jacobi method > Gauss-siedel method > SOR method.