Transient and steady heat conduction in a square plate - review of iterative methods
In this project we will compute the temperature distribution over an aluminium square plate both in transient and steady-state cases by solving the heat conduction equation using 3 iterative methods (Jacobi, Gauss-Seidel and Successive Over-Relaxation). Then, we will compare the performance of these methods.
In 2D, the heat conduction equation for the transient case when no heat is applied is:
`(del T)/(del t) = alpha *((del^2 T)/(del x^2) + (del^2 T)/(del y^2))`
with alpha being a possitive real number called "thermal diffusivity" with units of `m^2//s`. For Aluminium, alpha is approximately equal to `1*10^-4 m^2//s`
For the steady-state case, the time-derivative disappears and the equation is simply:
`(del^2 T)/(del x^2) + (del^2 T)/(del y^2) = 0`
To solve this equations numerically, we will use finite difference schemes to discretize the derivatives.
In the transient case, to ensure stability for any time-step size, we will need to use implicit schemes, for which we have to solve a linear system `A*x = b` in every time-step. When the mesh is very big, solving the system will be time consuming. Iterative methods have been developed in the past to accelerate the resolution of large linear systems. We will review some of these methods.
1. TRANSIENT HEAT CONDUCTION
We use central scheme for the spatial second derivatives, and first order backward scheme for the time derivative, resulting in an implicit solver.
If we discretize the square plate in equally spaced Nx points in the x direction and Ny points in the y direction resulting in a grid of `N_x*N_y` points, then the coordinates x,y of each point are given by:
`x_i = i*Deltax`
`y_j = j*Deltay`
with `Deltax = L/(N_x-1)` and `Deltay = L/(N_y-1)` and i,j being the integer indexes that uniquely localize every point.
If we define the constants:
`k_1 = (alpha*Delta t)/(Delta x^2)`
`k_2 = (alpha*Delta t)/(Delta y^2)`
Then, the temperature at every point of the square plate of spatial indexes `(i,j)` can be calcuated for the next time step (n+1) from the previos time step (n) by solving the linear equation:
`T_(i,j)^(n+1) = T_(i,j)^(n) + k_1*(T_(i-1,j)^(n+1) - 2T_(i,j)^(n+1) +T_(i+1,j)^(n+1) )+k_2*(T_(i,j-1)^(n+1) - 2T_(i,j)^(n+1) +T_(i,j+1)^(n+1) )`
If we write this equations as a system `A*x = b`, the matrix A would have the terms:
`(1+2k_1+2k_2)` in the diagonal
`k_1` for the left and right points of each point in the grid.
`k_2` for the top and bottom points of each point in the grid.
And the independent term b, will have the terms `T_p^(n)` for each point p.
To solve this `N_x N_y xx N_xN_y` linear system, there are direct methods such as Gauss reduction that solve the system with exact precision, and there are also iterative methods which consist on guessing the solution and applying operations that will make the guess closer and closer to the exact solution, until we get the desired precision. This methods are much more efficient than direct methods when the system is very large.
Some of the most famous methods are Jacobi, Gauss-Seidel and SOR:
Jacobi method
For the system `sum_(j)a_(ij) x_i = b_i`, the (k+1) guess of the Jacobi method would be:
`x_i^((k+1)) = 1/a_(ii) (b_i - sum_(j!=i)a_(ij)x_j^((k)))`
Gauss-Seidel method
Similarly, the GS method would be:
`x_i^((k+1)) = 1/a_(ii) (b_i - sum_(j<i)a_(ij)x_j^((k+1)) - sum_(j>i)a_(ij)x_j^((k)))`
Successive Over-Relaxation method
This method adds a parameter that needs to be tuned, called "relaxation factor" (here noted omega) that is a real number between 0 and 2:
`x_i^((k+1)) = (1-omega)x_i^((k)) + omega/a_(ii) (b_i - sum_(j<i)a_(ij)x_j^((k+1))-sum_(j>i)a_(ij)x_j^((k)))`
Notice that the SOR method with omega = 1 is equivalent to the GS method.
Let's apply these methods to the following problem:
Square plate of length Lx = Ly = 0.1 m
Boundary conditions: applied temperatures:
- Left: 400 K
- Right: 800 K
- Top: 600 K
- Bottom: 900 K
Initial conditions: homogeneous temperature of 300 K
Material: aluminium (alpha = 1e-4 m^2/s)
Duration of the simulation: tfin = 15 s
We will initially discretize the domain in a grid of Nx = Ny = 50 points for each dimension so 2500 points total. The time discretization will be Nt = 500 points
This code will set up the problem for all the methods we are going to use:
%% Transient 2D heat conduction equation
% @Angel Gonzalez 2019
%% Parameters
Lx = 0.1; % x length [m]
Ly = 0.1; % y length [m]
alpha = 1e-4; %Heat diffusivity [m^2/s] approximate value for Aluminium
tfin = 15; %Final time [s]
T0 = 300; %Initial temperature of domain [K]
Tleft = 400; %Temperature at left boundary [K]
Ttop = 600; %Temperature at top boundary [K]
Tright = 800; %Temperature at right boundary [K]
Tbottom = 900; %Temperature at bottom boundary [K]
Nx = 50;
Ny = 50;
Nt = 500;
errconv = 1e-4; %Error criteria for convergence loop
%Derived parameters
dx = Lx/(Nx-1);
dy = Ly/(Ny-1);
dt = tfin/(Nt-1);
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
%% Initial conditions and Boundary Conditions
T = nan(Nx, Ny, Nt);
%I.C.
T(:,:,1) = T0*ones(Nx,Ny);
%B.C.
T(1,:,:) = Tleft;
T(end,:,:) = Tright;
T(:,1,:) = Tbottom;
T(:,end,:) = Ttop;Then, for each time iteration we will apply the iterative method until the error (the maximum difference between each iteration) is less than 1e-4 Kelvin.
For the Jacobi method, the loop will be:
%% Jacobi
iterJac=1;
tic;
for i = 2:Nt
T_old = T(:,:,i-1);
T_new = T_old;
err = inf;
while err > errconv
T_new(2:Nx-1,2:Ny-1) = 1/(1+2*k1+2*k2)*( ...
T(2:Nx-1,2:Ny-1, i-1) ...
+ k1*( T_old(1:Nx-2, 2:Ny-1) + T_old(3:Nx, 2:Ny-1) ) ...
+ k2*( T_old(2:Nx-1, 1:Ny-2) + T_old(2:Nx-1, 3:Ny) ) );
err = max(max(abs(T_new - T_old)));
T_old = T_new;
iterJac=iterJac+1;
end
T(:,:,i) = T_new;
end
tocJac = toc;
Tjac = T;The Gauss-Seidel loop would be:
%% Gauss-Seidel
iterGS=1;
tic;
for i = 2:Nt
T_old = T(:,:,i-1);
T_new = T_old;
err = inf;
while err > errconv
for j = 2:Nx-1
for k = 2:Ny-1
T_new(j,k) = 1/(1+2*k1+2*k2)*( ...
T(j,k, i-1) ...
+ k1*(T_new(j-1,k) + T_new(j+1,k)) ...
+ k2*(T_new(j,k-1) + T_new(j,k+1)) );
end
end
err = max(max(abs(T_new - T_old)));
T_old = T_new;
iterGS=iterGS+1;
end
T(:,:,i) = T_new;
end
tocGS = toc;
Tgs = T;And finally the SOR loop:
%% SOR
om = 1.5;
iterSOR = 1;
tic;
for i = 2:Nt
T_old = T(:,:,i-1);
T_new = T_old;
err = inf;
while err > errconv
for j = 2:Nx-1
for k = 2:Ny-1
T_new(j,k) = (1-om)*T_old(j,k) ...
+ om/(1+2*k1+2*k2)*( ...
T(j,k, i-1) ...
+ k1*(T_new(j-1,k) + T_new(j+1,k)) ...
+ k2*(T_new(j,k-1) + T_new(j,k+1)) );
end
end
err = max(max(abs(T_new - T_old)));
T_old = T_new;
iterSOR = iterSOR+1;
end
T(:,:,i) = T_new;
end
tocSOR = toc;
Tsor = T;
All of the simulations give the same results (within a 1e-4 K error). The number of iterations is a significant indicator of speed. Jacobi is the slowest with 13653 iterations, then GS with 8493 and the best is SOR with 5905. This is reflected in machine time of execution, which in my machine is 0.3822s for Jacobi, 0.3065s for GS and 0.2857s for SOR.
The animation of SOR results can be seen in the following links:
1) for the Temperature contour plot:
2) for the Heat flux quiver plot:
The code to generate this animations is:
%% Refactor the Temperature matrices
for i = 1:Nt
Tjac(:,:,i) = flipud(rot90(Tjac(:,:,i)));
Tgs(:,:,i) = flipud(rot90(Tgs(:,:,i)));
Tsor(:,:,i) = flipud(rot90(Tsor(:,:,i)));
end
%% Animation of Temperature and Heat
% 1. Temperature contour
[xx, yy ] = meshgrid(linspace(0,Lx,Nx), linspace(0,Ly,Ny));
figure('Position', [0 0 1024 768]);
for i = 1:Nt
[x1,y1] = contourf(xx,yy, Tsor(:,:,i), 12);
clabel(x1,y1);
M1(i) = getframe();
end
% 2. Heat quiver
xxquiv = xx;
yyquiv = yy;
xxquiv(1,:) = [];
xxquiv(end,:) = [];
xxquiv(:,1) = [];
xxquiv(:,end) = [];
yyquiv(1,:) = [];
yyquiv(end,:) = [];
yyquiv(:,1) = [];
yyquiv(:,end) = [];
figure('Position', [0 0 1024 768]);
for i = 1:Nt
[qX, qY] = gradient(Tsor(:,:,i));
qX = -qX;
qY = -qY;
quiver(xxquiv,yyquiv, qX(2:end-1,2:end-1), qY(2:end-1,2:end-1), 2, 'Color', [1,0,0],'LineWidth', 1.5);
hold on
contour(xx, yy, Tsor(:,:,i), 12);
hold off
pause(1e-3);
M2(i) = getframe();
end
2. STEADY HEAT CONDUCTION
In the steady-state case, the solution does not depend on time and is equivalent to the transient problem when time tends to infinity.
The same iterative methods can be applied, with only one convergence loop and no need for a time loop. The formulas are slightly different as the matrix to solve A x = b has changed.
The common setup is very similar than the previous one:
%% Steady-state 2D heat conduction equation
% @Angel Gonzalez 2019
%% Parameters
Lx = 0.1; % x length [m]
Ly = 0.1; % y length [m]
Tleft = 400; %Temperature at left boundary [K]
Ttop = 600; %Temperature at top boundary [K]
Tright = 800; %Temperature at right boundary [K]
Tbottom = 900; %Temperature at bottom boundary [K]
T0 = 400; %Initializing temperature [K]
Nx = 50;
Ny = 50;
errconv = 1e-4; %Error criteria for convergence loop
%Derived parameters
dx = Lx/(Nx-1);
dy = Ly/(Ny-1);
a1 = 1/(dx^2*2*(1/dx^2+1/dy^2));
a2 = 1/(dy^2*2*(1/dx^2+1/dy^2));
%% Boundary Conditions
T = nan(Nx, Ny);
%B.C.
T(1,:) = Tleft;
T(end,:) = Tright;
T(:,1) = Tbottom;
T(:,end) = Ttop;The Jacobi loop:
%% Jacobi
T_old = T;
T_old(2:Nx-1, 2:Ny-1) = T0*ones(Nx-2, Ny-2);
T_new = T_old;
err = inf;
iterJac=1;
tic;
while err > errconv
T_new(2:Nx-1,2:Ny-1) = ...
+ a1*( T_old(1:Nx-2, 2:Ny-1) + T_old(3:Nx, 2:Ny-1) ) ...
+ a2*( T_old(2:Nx-1, 1:Ny-2) + T_old(2:Nx-1, 3:Ny) );
err = max(max(abs(T_new - T_old)));
T_old = T_new;
iterJac=iterJac+1;
end
tocJac = toc;
Tjac = T_new;The GS loop:
%% Gauss-Seidel
T_old = T;
T_old(2:Nx-1, 2:Ny-1) = T0*ones(Nx-2, Ny-2);
T_new = T_old;
err = inf;
iterGS =1;
tic;
while err > errconv
for i = 2:Nx-1
for j = 2:Ny-1
T_new(i,j) = ...
+ a1*( T_new(i-1,j) + T_new(i+1,j) ) ...
+ a2*( T_new(i,j-1) + T_new(i,j+1) );
end
end
err = max(max(abs(T_new - T_old)));
T_old = T_new;
iterGS = iterGS+1;
end
tocGS = toc;
Tgs = T_new;The SOR loop:
%% SOR
om = 1.5;
T_old = T;
T_old(2:Nx-1, 2:Ny-1) = T0*ones(Nx-2, Ny-2);
T_new = T_old;
err = inf;
iterSOR = 1;
tic;
while err > errconv
for i = 2:Nx-1
for j = 2:Ny-1
T_new(i,j) = (1-om)*T_old(i,j) ...
+ om*(...
a1*( T_new(i-1,j) + T_new(i+1,j) ) ...
+ a2*( T_new(i,j-1) + T_new(i,j+1) ) );
end
end
err = max(max(abs(T_new - T_old)));
T_old = T_new;
iterSOR = iterSOR +1;
end
tocSOR = toc;
Tsor = T_new;The results are equivalent (within a 1e-4 K error) for all iterative methods. The plot generated for SOR is:
3. PERFORMANCE COMPARISON OF ITERATIVE METHODS
The number of iterations required to converge to the solution for a desired precision is an important measure of the method's performance and computational cost.
For the SOR method, there is an extra degree of freedom to choose, the "omega" value, which has to be between 0 and 2. Let's repeat the calculations in steady state for multiple values of omega to see how the performance varies with it:
%% SOR omega list
omlist = linspace(0.1, 1.99, 190);
tocSORlist = nan(size(omlist));
iterSORlist = nan(size(omlist));
for n = 1:length(omlist)
om = omlist(n);
T_old = T;
T_old(2:Nx-1, 2:Ny-1) = T0*ones(Nx-2, Ny-2);
T_new = T_old;
iter = 1;
tic;
err = inf;
while err > errconv
for i = 2:Nx-1
for j = 2:Ny-1
T_new(i,j) = (1-om)*T_old(i,j) ...
+ om*(...
a1*( T_new(i-1,j) + T_new(i+1,j) ) ...
+ a2*( T_new(i,j-1) + T_new(i,j+1) ) );
end
end
err = max(max(abs(T_new - T_old)));
T_old = T_new;
iter = iter +1;
end
tocSORlist(n) = toc;
iterSORlist(n) = iter;
end
figure()
semilogy(omlist, iterSORlist);
hold on
plot([omlist(1), omlist(end)], [iterJac, iterJac]);
plot([omlist(1), omlist(end)], [iterGS, iterGS]);
title('Number of iterations, (Nx=Ny=50)')
xlabel('omega')
ylabel('iterations');
legend({'SOR', 'Jacobi', 'Gauss-Seidel'})The resulting plot is:
We can see that GS is better than Jacobi and that for the range (1, 1.99] the SOR method has a lower number of iterations than both Jacobi and GS, resulting in better performance, but is worse for the rest of values. The optimal omega value is 1.84
However, an even more accurate measure of performance is the number of floating-point operations of a simulation (number of FLOPS). A FLOP is an addition, substraction, multiplication or division of two floating point numbers.
Observing the code, for each iteration, inside the for loops, in Jacobi we have 10 flops, for GS we have also 10 flops, but for SOR we have 13 flops due to the extra operations with the omega parameter. These operations are repeated for all the internal points so 28*28 = 2304 times due to the for loops. Additionally, the calculation of the error involves 2500 more flops.
Hence, the number of total flops for all methods can be computed:
%% FLOPS calculation
flopsJacperIter = 10*(Nx-2)*(Ny-2) + Nx*Ny;
flopsGSperIter = 10*(Nx-2)*(Ny-2) + Nx*Ny;
flopsSORperIter = 13*(Nx-2)*(Ny-2) + Nx*Ny;
flopsJac = flopsJacperIter*iterJac;
flopsGS = flopsGSperIter*iterGS;
flopsSOR = flopsSORperIter*iterSOR;
flopsSORlist = flopsSORperIter*iterSORlist;
%% Plot
figure()
semilogy(omlist, flopsSORlist);
hold on
plot([omlist(1), omlist(end)], [flopsJac, flopsJac]);
plot([omlist(1), omlist(end)], [flopsGS, flopsGS]);
title('Number of floating-point operations (FLOPS), (Nx=Ny=50)')
xlabel('omega')
ylabel('flops');
legend({'SOR', 'Jacobi', 'Gauss-Seidel'})
The resulting plot is:
The results are very similar to the iterations plot, but due to the increased number of operations per cycle in the SOR loop, it is now only better than GS for the range (1.17, 1.99).
CONCLUSION
Iterative methods solve A x = b more efficiently for large systems, which is convenient for CFD problems. The GS method tends to be more performant than the Jacobi method. The SOR method is normally more performant than both, but only for a specific range of omega values that is a-priori unknown and always between (0,2). Guessing the right omega value is an important skill.