Solving the steady and unsteady 2D heat conduction problem
In this report, the emphasis has been laid about solving the 2-D heat conduction equation ( Steady State and Transient State) using several iterative methods. The methods implemented are :
1.) Jacobi
2.) Gauss-Seidel
3.) Successive Overrelaxation ( SOR )
But first of all, let's understand why do we need iterative methods to solve discretized differential equations. Consider a Linear system of discretized differential equations of the form
AU = F
where A is the coefficient matrix, U is the unknown vector, and F is the right-hand side function. One way to solve these types of linear equations is with Direct Methods which include :
1.) Cramers rule
2.) Gaussian Elimination
However, if the matrix is of billions of elements like in the engineering scale problems, it is computationally very expensive to solve these systems using the direct methods. One indirect way to solve these systems using the iterative techniques which uses an initial guess and then iterates it over and over again until the two consecutive iterative differ not more than the tolerance, i.e. we say the solution has converged up to the tolerance level if,
`norm(u^{k+1} - u^k) \leq \epsilon`
where `normu` is the standard defined norm
`u^k` is the solution at the `k^(th)` iterative step
`epsilon` is the specified tolerance limit
However, these Iterative methods themselves aren't that easy to set up. One very important thing to keep in mind while setting these matrices is that the matrix A must be diagonally dominant or should satisfy the Scarborough criteria i.e.
`|a_{ii}| \geq sum_(i \nej) |a_{ij}|`
meaning that the absolute value of the diagonal element in a particular row must be greater than the absolute value of all the other remaining elements. One more way to think about this is that if A is a positive definite symmetric matrix and if the following condition is satisfied :
`\rho(R)<1` , where `\rho(R)` is the spectral radius of R
then it will guarantee the convergence of the iterative method. Also, an important fact to note here is that as the spectral radius of A goes near to 0 it will accelerate the convergence.
Now, with a basic understanding of the mathematics behind iterative methods we are ready to develop a general form of the iterative methods :
Consider the following basic form of an iteration :
`u^(k+1) = Ru^k +g`
where `u^(k+1)` is the value at the `(k+1)^(th)` iteration and `u^k`is the value at the `k^(th)` iteration and the R is the iteration matrix according to the scheme. According to the different types of schemes the matrix R can have different versions.
JACOBI ITERATION
Consider the discretized governing equation in the matrix form AU=F.
Let A be decomposed as A = D - L - U where D = Diagonal Matrix, L = Strictly Lower Triangular Matrix, and U = Strictly Upper Triangular Matrix.
We can write the equation in the form :
`(D-L-U)u = F`
`Du = (L+U)u + F` , taking the inverse on both sides we get :
`u = D^(-1)(L+U)u + D^(-1)F`
which can be furthur simplified to :
`u = R_ju+g`
where `R_j` is the JACOBI ITERATION matrix. In iteration form this can be written as :
`u^(k+1) = R_ju^k + g`
This method will converge only if `rho(R_j) < 1`.
GAUSS - SEIDEL ITERATION
A similar analysis can be performed for the Gauss-Seidel Iteration.
We can write :
`(D-L-U)u = F`
`(D-L)u = Uu + F`
`u = (D-L)^(-1)Uu + (D-L)^(-1)F`
which can be simplified to :
`u = R_Gu+g`
where `R_G` is the GAUSS-SEIDEL ITERATION matrix. In Iteration form this can be written as :
`u^(k+1) = R_Gu^k + g`
Similarly this method will converge only if `rho(R_G) < 1`.
SUCCESSIVE OVER-RELAXATION
A more elegant modification can be made to basic iterative techniques. The new iteration can be taken as a weighted average of the new and the old approximation according to the relaxation factor `alpha`.
In equation form :
`u^(k+1) = (alpha)(R_Gu^k + g) + (1-alpha)u^k`
simplifying this a little bit more we get :
`u^(k+1) = [(1-alpha)I + alphaR_G]u^k + alphag`
where the term `[(1-alpha)I + alphaR_G]` can be labelled as a new matrix `R_(alpha)`, so that the iteration can be written as :
`u^(k+1) = R_(alpha)u^k + alphag`
Now with a basic knowledge of the iteration techniques, we are now in a position to apply these to the 2-D Conduction problem :
`frac{partial^2T}{partialx^2} + frac{partial^2T}{partialy^2} = 0`
Consider the following 5-point stencil as shown below :
Discretizing the equation using the finite difference we get :
`frac{T_{i+1,j}-2T_{i,j}+T_{i-1,j}}{Deltax^2}+frac{T_{i,j+1}-2T_{i,j}+T_{i,j-1}}{Deltay^2} = 0`
Simplifying the equation as :
`T_(i,j) leftarrow frac( T_{i+1,j}+ T_{i-1,j}+ T_{i,j+1}+ T_{i,j-1} )(4)`
According to the several Iterative Techniques discussed the following schemes were implemented in MATLAB to solve the 2-D Heat Conduction equation. The domain selected for the problem is a 2-D unit Square so that `Deltax = Deltay`. The number of subintervals used for the discretization are taken as 10. The code for the same has been attached below :
clear all
close all
clc
%
% In this program I shall be demonstrating on how to solve the 2-D steady
% state equation with various iteration techniques
%
% The equation is of the form :
%
% Txx + Tyy = 0
%
% With the boundary conditions specified as :
%
% T_top = 600K
% T_bottom = 900K
% T_left = 400K
% T_right = 800K
%
% The domain of this problem is a 2-D unit square with equal number of grid
% points in the x and the y direction
%
% Defining the number of sub-intervals
n = input('Input the number of sub-intervals you want : ');
fprintf('\n');
% Defining the interval spacing
h = 1/n;
% Defining the initial guess of the temperature on the domain
T_initial = 300*ones(n);
% Defining the boundary conditions on the Initial Guess
T_initial(1,2:end-1) = 600;
T_initial(end,2:end-1) = 900;
T_initial(2:end-1,1) = 400;
T_initial(2:end-1,end) = 800;
% Defining the Initial Guess as the average of its nearest neighbours
T_initial(1,1) = 500;
T_initial(1,end) = 700;
T_initial(end,1) = 650;
T_initial(end,end) = 850;
% Creating a copy of T_initial
T_new = T_initial;
% inputing the type of solver to be employed
fprintf(' Types of iterative solvers available :');
fprintf('\n');
fprintf(' 1 - Jacobi Iterative Solver');
fprintf('\n');
fprintf(' 2 - Gauss Seidel Iterative Solver');
fprintf('\n');
fprintf(' 3 - Successive Relaxation solver');
fprintf('\n');
iterative_solver = input('Enter the type of solver you want to employ : ');
% Jacobi Method
if iterative_solver == 1
% Defining the tolerance value
tol = 1e-4;
% Defining the initial Error
error = 9e9;
% Defining the iteration count Jacobi_iteration count
Jacobi_count = 0;
% Main LOOP
while(error>tol)
% rows
for ii = 2:n-1
% columns
for jj = 2:n-1
% storing the values in the new variable while keeping the
% old ones intact per relaxation sweep
T_new(ii,jj) = T_initial(ii+1,jj)/4 + T_initial(ii-1,jj)/4 + ...
T_initial(ii,jj+1)/4 + T_initial(ii,jj-1)/4;
end
end
% Defining the error
error = max(max(abs(T_new-T_initial)));
% Updating the values of the t_initial for the next iteration
T_initial = T_new;
% Updating the iteration count
Jacobi_count = Jacobi_count + 1;
% Plotting the surface to see how the plot looks like
contourf(T_initial,'ShowText','on');
colormap(jet);
colorbar;
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title([ 'After ' num2str(Jacobi_count) ' Jacobi Iterations' ]);
set(gca,'XTick',[],'YTick',[]);
% pausing the loop to see the relaxation happening
pause(0.01)
end
% printing a message to screen as to how many JACOBI iterations it
% took for the solution to converge
fprintf([ ' It took ' num2str(Jacobi_count) ' Jacobi Iterations for the solution to converge']);
fprintf('\n');
% Employ the Gauss-Seidel Iteration
elseif iterative_solver == 2
% Defining the tolerance value
tol = 1e-4;
% Defining the initial Error
error = 9e9;
% Defining the iteration count Jacobi_iteration count
Gauss_seidel_count = 0;
% Main LOOP
while(error>tol)
% rows
for ii = 2:n-1
% columns
for jj = 2:n-1
% storing the values in the T after each iteration
T_new(ii,jj) = T_new(ii+1,jj)/4 + T_new(ii-1,jj)/4 + ...
T_new(ii,jj+1)/4 + T_new(ii,jj-1)/4;
end
end
% Defining the error
error = max(max(abs(T_new-T_initial)));
% Updating the values of the T_initial for the next iteration to
% check for the error convergence
T_initial = T_new;
% Updating the iteration count
Gauss_seidel_count = Gauss_seidel_count + 1;
% Plotting the surface to see how the plot looks like
contourf(T_initial,'ShowText','on');
colormap(jet);
colorbar;
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title([ 'After ' num2str(Gauss_seidel_count) ' Gauss-Seidel Iterations' ]);
set(gca,'XTick',[],'YTick',[]);
% pausing the loop to see the relaxation happening
pause(0.01)
end
% printing a message to screen as to how many JACOBI iterations it
% took for the solution to converge
fprintf([ ' It took ' num2str(Gauss_seidel_count) ' Gauss Seidel Iterations for the solution to converge']);
fprintf('\n');
% for successice Over-Relaxation
elseif iterative_solver == 3
% Inputing the relaxation parameter
alpha = input('Please Input the relaxation parameter : ');
% this is the SOR Iterative Solver
% Define the error :
error = 9e9;
% Define the tolerance
tol = 1e-5;
% Initialising the SOR count to keep track of the iterations
sor = 0;
% Intialising T_new_1
T_new_1 = T_initial;
% Main LOOP
while(error>tol)
% rows
for ii = 2:n-1
% columns
for jj = 2:n-1
% storing the values in the T after each iteration
T_new_1(ii,jj) = T_new(ii+1,jj)/4 + T_new(ii-1,jj)/4 + ...
T_new(ii,jj+1)/4 + T_new(ii,jj-1)/4;
% Adding this using the relaxation parameter
T_new = (alpha)*T_new_1 + (1-alpha)*T_new;
end
end
% Finding the error to check for the convergence
error = max(max(abs(T_new-T_initial)));
% Updating tthe T_initial for the next iteration step
T_initial = T_new;
% Updating the iteration count
sor = sor + 1;
% Plotting the surface to see how the plot looks like
contourf(T_initial,'ShowText','on');
colormap(jet);
colorbar;
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title([ 'After ' num2str(sor) ' Successive Over-relaxation Iterations using relaxation parameter ' ...
num2str(alpha)]);
set(gca,'XTick',[],'YTick',[]);
% pausing the loop to see the relaxation happening
pause(0.01)
end
% printing a message to screen as to how many JACOBI iterations it
% took for the solution to converge
fprintf([ ' It took ' num2str(sor) ' Successive Over-relaxation Iterations for the solution to converge with relaxation parameter ' ...
num2str(alpha)]);
fprintf('\n');
end
JACOBI ITERATION
The iteration uses the fact that all the values in a matrix/vector need to be updated before they can be used again for the next iteration. In mathematical form this can be written as :
`T_(i,j)^(k+1) leftarrow frac( T_{i+1,j}^k+ T_{i-1,j}^k+ T_{i,j+1}^k+ T_{i,j-1}^k )(4)`
As it is evident to calculate the `T_(i,j)` all the values are used from the previous time-step. The following results were obtained in MATLAB :
As it is clear it took about 207 JACOBI ITERATIONS to converge to the solution, i.e. to reduce the maximum norm of the error below the tolerance.
GAUSS-SEIDEL METHOD
A simple modification can be made to the JACOBI METHOD, instead of using the grid values from the previous step iteration they can be used as soon as they are calculated. In mathematical form :
`T_(i,j)^(k+1) leftarrow frac( T_{i+1,j}^k+ T_{i-1,j}^(k+1)+ T_{i,j+1}^k+ T_{i,j-1}^(k+1) )(4)`
As it is seen the values at the points `T_(i-1,j)` and `T_{i,j-1)` are used at `(k+1)^(th)` iteration step to calculate the value of `T_(i,j)` at the `(k+1)^(th)` iteration step i.e. the values are updated as soon as they are calculated and hence we could expect a faster convergence as the values are updated quickly and hence we could reduce the steps to reach to the convergence. This hypothesis is in fact in agreement to the results obtained in MATLAB :
As it is clear from the graph obtained, it took only 111 GAUSS-SEIDEL ITERATION steps to converge the solution below the tolerance, which are merely about half as many required for the JACOBI METHOD.
SUCCESSIVE OVER-RELAXATION
A simple modification can be made to the simple iterative methods i.e. a weighted average can be taken between the new and the old iteration values as can be stored in the new variable. In mathematical form :
`T_(i,j)^(k+1) leftarrow alpha[frac( T_{i+1,j}^k+ T_{i-1,j}^(k+1)+ T_{i,j+1}^k+ T_{i,j-1}^(k+1) )(4)] + (1-alpha)T_(i,j)^(k)`
If put in simple words this method can be thought of as pushing the new approximate towards the exact_solution. Let's assume that the JACOBI/GAUSS-SEIDEL Iteration pushes the old approximation towards the original solution by 1 unit, then the relaxation method pushes the solution towards the exact solution by roughly `alpha` times. If `alpha` is selected to be too large it can cause a phenomenon called Over-Shooting i.e. instead of pushing the approximation towards the solution it can shoot it away from the solution.
In this method, the best selection of the relaxation parameter `alpha` is a must as without a proper selection of `alpha` the method will diverge. It can be shown from rigorous mathematical calculations that if `0\leqalpha\leq2` the iteration matrix has a spectral radius less than one. However, still, a careful selection for `alpha` is a must to accelerate the convergence process. The following results were obtained in MATLAB for the employment of the SOR technique.
The following outputs were obtained in the command window :
As is clear it took merely 104 iterations for the solution to converge using the relaxation parameter as 1.4 Now we shall see some more results with different relaxation parameters:
As it is seen from the above outputs it took about :
- 97 iterations with `alpha = 1.5` to converge the solution
- 86 iterations with `alpha = 1.7` to converge the solution
- 97 iterations with `alpha = 1.9` to converge the solution
- 148 iterations with `alpha = 2` to converge the solution
From the data obtained above, it can be concluded that the optimum acceleration parameter,`alpha` is between 1.6 and 1.8.
Consider one more case but this time with the relaxation parameter as 4 which is much higher than the prescribed range.
It may seem that the solution has converged after 46 iterations using relaxation parameter as 2.3. But if we look at the plot :
We can see the solution has blown up and MATLAB isn't even able to plot it properly.
As we can see the iterations in the SOR method are comparatively lower than that of JACOBI but if the relaxation parameter is not selected properly it may take even more time to converge than the GAUSS-SEIDEL iteration.
CONCLUSION: As it is seen from the numerical computation performed the 2-D steady-state heat conduction equation we see that the gauss-seidel method gave us the best convergence result i.e. the least number of iteration steps taken followed by SOR followed by the JACOBI-METHOD. The reason can be thought of as like this: In JACOBI ITERATION the new values are calculated from the previous values first completely before utilizing them for the next iteration. This causes an increase in the number of computations required. In Gauss-Seidel the new values are utilized as soon as they are calculated which leads to accelerated convergence as one is using the improved values to improve the new approximation in one iteration and hence it leads to computational advantage as well as less. memory storage requirement. In the SOR technique, we see that we got a higher iteration number than the GAUSS-SEIDEL for some of the relaxation parameters. This is because in pushing the approximation towards the solution it may sometimes push it a little forward and hence the approximation will try to come back to the exact solution and it will keep on oscillating around the solution till the point the frequencies in the error gets damped. It can be though of much like a potential mountain where a particle left from above will settle down at equilibrium ( i.e. the valley ) only if it has the right momentum, if it has too excessive momentum it will bounce off and will never return to the valley. If it has too low a momentum it will take a longer time to come to the valley. If it has slightly more momentum than normal than the particle will oscillate about the valley before stopping at the equilibrium. Much of what described above happens in the iterative techniques. Therefore a careful selection of the relaxation parameter `alpha` is a must to avoid over-shooting as it will take an unnecessary number of iterations to converge.
TRANSIENT STATE ANALYSIS
Now we are ready to analyze the transient systems. The only difference between the steady-state and transient state analysis is that in the study state the solution will not change with time. However, in the transient systems the solution changes with time and after a very long time it will reach steady-state. The transient form of the 2-D heat conduction equation without the presence of the source term can be simplified to :
`frac(partialT)(partialt) = alpha(frac(partial^2T)(partialx^2)+frac(partial^2T)(partialy^2))`
Then in the discretized form, with `Deltax` and `Deltay` the grid spacing along the x and the y-directions respectively and `Deltat` as the time step, this equation can be written as:
`frac(T_(current) - T_(old))(Deltat) = alpha(frac(T_L - 2T_P + T_R)(Deltax^2) +frac(T_T - 2T_P + T_B)(Deltay^2))`
With the grid point (i,j) and the time step as n this equation can be further simplified to :
`frac(T_(i,j)^(n+1) - T_(i,j)^n)(Deltat) = alpha(frac(T_(i-1,j) - 2T_(i,j) + T_(i+1,j))(Deltax^2) +frac(T_(i,j+1) - 2T_(i,j) + T_(i,j-1))(Deltay^2))`
In the Explicit method, all the values at the right-hand-side are from the previous time step and hence are known. So we can isolate the new time step term `T_(i,j)^(n+1)` to the left and take everything else to the right and we get an equation of the form :
`T_(i,j)^(n+1) = T_(i,j)^n+ alphaDeltat( frac{T_(i+1,j)^n - 2T_(i,j)^n + T_(i-1,j)^n }{Deltax^2}+frac{T_(i,j+1)^n - 2T_(i,j)^n + T_(i,j-1)^n}{Deltay^2})`
Simplifying it a little bit, we get :
`T_(i,j)^(n+1) = T_(i,j)^n( 1-2 frac(alphaDeltat)(Deltax^2) -2 frac(alphaDeltat)(Deltay^2)) + alphaDeltat(frac{T_(i+1,j)^n + T_(i-1,j)^n }{Deltax^2})+ alphaDeltat(frac{T_(i,j+1)^n + T_(i,j-1)^n}{Deltay^2})`
Now with the values known from the previous time step, we can easily substitute them in this equation and calculate the new time step value at each grid point. One more important thing to keep in mind is the selection of `alpha`,`Deltat`,`Deltax` and `Deltay`. If they are not selected properly the solution may never reach steady-state, Instead, it will keep growing and will blow up. For convenience the convergence criteria for the explicit method can be checked as :
`frac(alphaDeltat)(Deltax^2) + frac(alphaDeltat)(Deltax^2) \leqfrac(1)(2)`
The above condition is called the CFL stability criteria and the solution from the explicit form will converge only if this condition is satisfied.
The following code was implemented in MATLAB to solve the equation using the Explicit method and the following results were obtained:
clear;
close all;
clc;
%
% In this program I will try to solve the 2-D transient Heat Conduction
% Problem using Explicit Method :
%
% The equation is of the form :
%
% alpha(Txx + Tyy) = Ttt
%
% With the boundary conditions specified as :
%
% T_top = 600K
% T_bottom = 900K
% T_left = 400K
% T_right = 800K
%
% The domain of this problem is a 2-D unit square with equal number of grid
% points in the x and the y direction
%
% Defining the number of sub-intervals
nx = input('Input the number of sub-intervals you want along the y-direction : ');
fprintf('\n');
ny = input('Input the number of sub-intervals you want along the x-direction : ');
fprintf('\n');
% Defining the interval spacing in x-direction
h = 1/nx;
% Defining the interval spacing in y-direction
h_1 = 1/ny;
% Defining the time step spacing
dt = input(' Input the time step value you want to solve the problem : ');
% Defining the initial guess of the temperature on the domain to start
% iteration in the time loop
T_previous_dt = 0*ones(nx,ny);
% Defining the boundary conditions on the Initial Guess
T_previous_dt(1,2:end-1) = 600;
T_previous_dt(end,2:end-1) = 900;
T_previous_dt(2:end-1,1) = 400;
T_previous_dt(2:end-1,end) = 800;
% Defining the Initial Guess as the average of its nearest neighbours
T_previous_dt(1,1) = 500;
T_previous_dt(1,end) = 700;
T_previous_dt(end,1) = 650;
T_previous_dt(end,end) = 850;
% Creating a copy of T_initial
T_new_dt = T_previous_dt;
% Defining the alpha
alpha = input(' Enter the thermal diffusivity(alpha) you would like to use : ');
fprintf('\n');
% As the spacing in the x and the y-direction is the same
CFL_constant = alpha*dt/h^2 + alpha*dt/h_1^2;
% Defining the constants
K1 = alpha*dt/h^2;
K2 = alpha*dt/h_1^2;
% As it is an Explicit method we don't need any iterative methods to solve
% it as the temperature values are already known from the previous time
% step
%
% doing the sanity check to ensure the convergence of the method using the
% CFL criteria :
if CFL_constant > 0.5
error([ 'Please enter the time step and the grid size again as the solution won''t converge' ...
'\n The CFL stability criteria is not satisfied !!!! \n']);
% If the criteria is satisfied begin with the main loop
else
% Defining the Initial difference
difference = 9e9;
% tolerance
tol = 1e-4;
% Time loop count
nt = 0;
% the main Time loop
%
% The steady state will be reached when the difference between the old
% and new values is negligible which will be satisfied by this while
% loop
%
% setup movie
mov = VideoWriter('Transient State 2-D conduction using Explicit method.avi');
mov.FrameRate = 15;
open(mov);
while(difference>tol)
% the space loop
for ii = 2:nx-1
% columns
for jj = 2:ny-1
% computing the value of T at the new time step
T_new_dt(ii,jj) = T_previous_dt(ii,jj) + K1 * ( T_previous_dt(ii+1,jj) + ...
T_previous_dt(ii-1,jj) -2*T_previous_dt(ii,jj) ) + ...
K2*(T_previous_dt(ii,jj+1) + T_previous_dt(ii,jj-1) -2*T_previous_dt(ii,jj));
end
end
% updating the time loop count
nt = nt+1;
% difference between previous and new iteration time step
difference = max(max(abs(T_previous_dt - T_new_dt)));
% Updating the previous time step value for the next iteration
T_previous_dt = T_new_dt;
% Plotting some figures
% Plotting the surface to see how the plot looks like
figure(1);clf;
[C,q] = contourf(T_new_dt,'ShowText','on');
colormap(jet);
clabel(C,q,'Fontsize',10);
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title([ 'After ' num2str(nt) ' Time loop Iterations' ]);
set(gca,'XTick',[],'YTick',[]);
% update movie frame
writeVideo(mov,getframe(gcf));
% pausing the loop to see the relaxation happening
pause(0.01);
end
end
% write out movie and close the handle
close(mov);
% Making the final Plot
figure(2);clf;
% Plotting the final steady state graph
[C,q] = contourf(T_new_dt,'ShowText','on');
colormap(jet);
clabel(C,q,'Fontsize',10);
xlabel('X-Axis');
ylabel('Y-Axis');
colorbar;
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title([ 'After ' num2str(nt) ' Time loop Iterations' ]);
set(gca,'XTick',[],'YTick',[]);
The time-step value was selected as 0.0001 sec along with an equal number of sub-intervals in the x and the y-direction. The following graph and animation were obtained:
As is seen it takes nearly 2500 iterations in the main time loop so that the two consecutive iterates differ not more than `10^(-3)`. Ideally, it would take an infinite amount of time for the solution to converge.
The following animation was obtained regarding the same :
IMPLICIT METHOD
A simple yet elegant modification can be made to the explicit method. Instead of the taking the values on the right-hand side of the equation at `n^(th)` we can take the values at the `(n+1)^(th)` time step which would give us better results as we are aware that as soon the information is transmitted it will be used in the calculation.
`frac(T_(i,j)^(n+1) - T_(i,j)^n)(Deltat) = alpha(frac(T_(i-1,j) - 2T_(i,j) + T_(i+1,j))(Deltax^2) +frac(T_(i,j+1) - 2T_(i,j) + T_(i,j-1))(Deltay^2))`
So, instead of calculation from the previous time step, why not use the current time step. This results in the equation of the form :
`T_(i,j)^(n+1) = T_(i,j)^n+ alphaDeltat( frac{T_(i+1,j)^(n+1) - 2T_(i,j)^(n+1) + T_(i-1,j)^(n+1) }{Deltax^2}+frac{T_(i,j+1)^(n+1) - 2T_(i,j)^(n+1) + T_(i,j-1)^(n+1) }{Deltay^2})`
Simplifying it a little bit we get :
`T_(i,j)^(n+1)( 1 + frac(alphaDeltat)(Deltax^2) +frac(alphaDeltat)(Deltay^2)) = T_(i,j)^n+ alphaDeltat( frac{T_(i+1,j)^(n+1) + T_(i-1,j)^(n+1) }{Deltax^2}+frac{T_(i,j+1)^(n+1) + T_(i,j-1)^(n+1) }{Deltay^2})`
Several Iterative Techniques can be applied to solve these equations :
1.) Jacobi
2.) Gauss-Seidel
3.) Successive-Over-relaxation
The following MATLAB code was implemented using these several Iterative solvers to solve the 2-D transient heat conduction equation.
clear;
close all;
clc;
%
% In this program I will try to solve the 2-D Transient Heat Conduction
% Problem using Explicit Method :
%
% The equation is of the form :
%
% alpha(Txx + Tyy) = Ttt
%
% With the boundary conditions specified as :
%
% T_top = 600K
% T_bottom = 900K
% T_left = 400K
% T_right = 800K
%
% The domain of this problem is a 2-D unit square with equal number of grid
% points in the x and the y-direction
%
% Defining the number of sub-intervals
nx = input('Input the number of sub-intervals you want along the y-direction : ');
fprintf('\n');
ny = input('Input the number of sub-intervals you want along the x-direction : ');
fprintf('\n');
% Defining the interval spacing in x-direction
h = 1/nx;
% Defining the interval spacing in the y-direction
h_1 = 1/ny;
% Defining the time step spacing
dt = input(' Input the time step value you want to solve the problem : ');
% Defining the initial guess of the temperature on the domain to start
% iteration in the time loop
T_previous_dt = 300*ones(nx,ny);
% Defining the boundary conditions on the Initial Guess
T_previous_dt(1,2:end-1) = 600;
T_previous_dt(end,2:end-1) = 900;
T_previous_dt(2:end-1,1) = 400;
T_previous_dt(2:end-1,end) = 800;
% Defining the Initial Guess as the average of its nearest neighbours
T_previous_dt(1,1) = 500;
T_previous_dt(1,end) = 700;
T_previous_dt(end,1) = 650;
T_previous_dt(end,end) = 850;
% Creating a copy of T_initial
T_new_dt = T_previous_dt;
% Creating a new copy of T for the use of the spacial iterations
T_new = T_new_dt;
% Defining the alpha
alpha = input(' Enter the thermal diffusivity(alpha) you would like to use : ');
fprintf('\n');
% Defining the constants
K1 = alpha*dt/h^2;
K2 = alpha*dt/h_1^2;
% As it is an Implicit method we need iterative methods to solve
% as the temperature values are unknown from the previous time
% step
%
% In this three types of Iterative Solvers have been implemented :
%
% 1.) Jacobi
% 2.) Gauss - Seidel
% 1.) Successive - Over - Relaxation
%
% inputing the type of solver to be employed
fprintf(' Types of iterative solvers available :');
fprintf('\n');
fprintf(' 1 - Jacobi Iterative Solver');
fprintf('\n');
fprintf(' 2 - Gauss Seidel Iterative Solver');
fprintf('\n');
fprintf(' 3 - Successive Relaxation solver');
fprintf('\n');
iterative_solver = input('Enter the type of solver you want to employ : ');
%%
% Jacobi Method
if iterative_solver == 1
% tolerance
tol = 1e-4;
% Defining the state error to allow to reach the steady-state
state_error = 9e9;
% Time counter loop index
ct = 0;
% the main Time loop
%
% The steady-state will be reached when the difference between the old
% and new values are negligible which will be satisfied by this while
% loop
%
% setup movie
mov = VideoWriter('Transient State 2-D conduction using Implicit method with the help of JACOBI ITERATIONS.avi');
mov.FrameRate = 10;
open(mov);
% Main Time Loop
while(state_error>tol)
% Initialising the T for the space iteration matrix as the average
% as to allow for the faster convergence using iterative techniques
T_initial = 0.5*(T_new_dt + T_previous_dt);
% Defining the error in the JACOBI ITERATION
error = 9e9;
% the main loop for the convergence of the iterative solver
while(error>tol)
% rows
for ii = 2:nx-1
% columns
for jj = 2:ny-1
% storing the values in the new variable while keeping the
% old ones intact per relaxation sweep
T_new(ii,jj) = T_previous_dt(ii,jj)/(1+2*K1+2*K2) + K1*( T_initial(ii+1,jj) + ...
T_initial(ii-1,jj) )/(1+2*K1+2*K2) + ...
K2*( T_initial(ii,jj+1) + T_initial(ii,jj-1) )/(1+2*K1+2*K2);
end
end
% Defining the error to converge for the JACOBI LOOP
error = max(max(abs(T_new - T_initial)));
% Updating the T_old value
T_initial = T_new;
end
% Finding the state error to allow for convergence in steady state
state_error = max(max(abs(T_previous_dt - T_new)));
% Updating the time loop counter
ct = ct + 1;
% Updating the previous time step value once the solution has
% converged in a particular time step
T_previous_dt = T_new;
% Plotting the surfaces in the multiples of 20's only
if rem(ct,20) == 0
% Plotting the surface to see how the plot looks like
contourf(T_initial,'ShowText','on');
colormap(jet);
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title({[ 'After ' num2str(ct*dt) ' seconds'];...
['Using the Jacobi Iteration']});
set(gca,'XTick',[],'YTick',[]);
% update movie frame
writeVideo(mov,getframe(gcf));
% pausing the loop to see the relaxation happening
pause(0.001)
end
end
% write out the movie and close the handle
close(mov);
% Plot the final graph after about 1400 iterations
% Plotting the surface to see how the plot looks like
figure(2);clf;
contourf(T_initial,'ShowText','on');
colormap(jet);
colorbar;
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title({[ 'After ' num2str(ct*dt) ' seconds'];['Using the Jacobi Iteration']});
set(gca,'XTick',[],'YTick',[]);
% Employ the Gauss-Seidel Iteration
elseif iterative_solver == 2
% tolerance
tol = 1e-4;
% Defining the state error to allow to reach the steady-state
state_error = 9e9;
% Time counter loop index
ct = 0;
% the main Time loop
%
% The steady-state will be reached when the difference between the old
% and new values are negligible which will be satisfied by this while
% loop
%
% setup movie
mov = VideoWriter('Transient State 2-D conduction using Implicit method with the help of GAUSS-SEIDEL ITERATIONS.avi');
mov.FrameRate = 10;
open(mov);
% Main Time Loop
while(state_error>tol)
% Initialising the T for the space iteration matrix
T_initial = 0.5*(T_new_dt + T_previous_dt);
% Defining the error in the GAUSS-SEIDEL ITERATION
error = 9e9;
% the main loop for the convergence of the iterative solver
while(error>tol)
% rows
for ii = 2:nx-1
% columns
for jj = 2:ny-1
% storing the values in the new variable while keeping the
% old ones intact per relaxation sweep
T_new(ii,jj) = T_previous_dt(ii,jj)/(1+2*K1+2*K2) + ...
K1*( T_new(ii+1,jj) + T_new(ii-1,jj) ) /(1 + 2*K1 + 2*K2 ) + ...
K2*( T_new(ii,jj+1) + T_new(ii,jj-1) ) /(1 + 2*K1 + 2*K2);
end
end
% Defining the error to converge for the JACOBI LOOP
error = max(max(abs(T_new - T_initial)));
% Updating the T_old value
T_initial = T_new;
end
% Updating the state error to allow for convergence of steady state
state_error = max(max(abs(T_previous_dt - T_new)));
% Updating the previous time step value once the solution has
% converged in a particular time step
T_previous_dt = T_new;
% Updating the loop counter
ct = ct + 1;
% Plotting the surfaces in the multiples of 30's only
if rem(ct,20) == 0
% Plotting the surface to see how the plot looks like
contourf(T_initial,'ShowText','on');
colormap(jet);
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title({[ 'After ' num2str(ct*dt) ' Seconds'];...
'Iterations using the Gauss-Seidel ITERATION' });
set(gca,'XTick',[],'YTick',[]);
% update movie frame
writeVideo(mov,getframe(gcf));
% pausing the loop to see the relaxation happening
pause(0.001)
end
end
% write out the movie and close the handle
close(mov);
% Plot the final graph after about 1400 iterations
% Plotting the surface to see how the plot looks like
figure(2);clf;
contourf(T_initial,'ShowText','on');
colormap(jet);
colorbar;
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title({[ 'After ' num2str(ct*dt) ' Seconds'];...
'Iterations using the Gauss-Seidel ITERATION' });
set(gca,'XTick',[],'YTick',[]);
% for successice Over-Relaxation
elseif iterative_solver == 3
% tolerance
tol = 1e-4;
% Defining the state error to allow to reach the steady-state
state_error = 9e9;
% Time counter loop index
ct = 0;
% the main Time loop
%
% The steady-state will be reached when the difference between the old
% and new values are negligible which will be satisfied by this while
% loop
%
% Entering the relaxation parameter
w = input('Enter the relaxation parameter omega : ');
fprintf('\n');
% setup movie
mov = VideoWriter('Transient State 2-D conduction using Implicit method with the help of SOR Technique.avi');
mov.FrameRate = 10;
open(mov);
% Main Time Loop
while(state_error>tol)
% Initialising the T for the space iteration matrix
T_initial = 0.5*(T_new_dt + T_previous_dt);
% Defining the error in the SOR ITERATION
error = 9e9;
% the main loop for the convergence of the iterative solver
while(error>tol)
% rows
for ii = 2:nx-1
% columns
for jj = 2:ny-1
% storing the values in the new variable while keeping the
% old ones intact per relaxation sweep
T_new(ii,jj) = T_previous_dt(ii,jj)/(1+2*K1+2*K2) + K1*( T_initial(ii+1,jj) + ...
T_initial(ii-1,jj) )/(1+2*K1+2*K2) + ...
K2*(T_initial(ii,jj+1) + T_initial(ii,jj-1) )/(1+2*K1+2*K2);
% Actual T_new will be the average of the two
T_new(ii,jj) = w*T_new(ii,jj) + (1-w)*T_initial(ii,jj);
end
end
% Defining the error to converge for the SOR LOOP
error = max(max(abs(T_new - T_initial)));
% Updating the T_old value
T_initial = T_new;
end
% Updating the state error
state_error = max(max(abs(T_previous_dt - T_new)));
% Updating the previous time step value once the solution has
% converged in a particular time step
T_previous_dt = T_new;
% Updating the loop counter
ct = ct + 1;
% Plotting the surfaces in the multiples of 20's only
if rem(ct,20) == 0
% Plotting the surface to see how the plot looks like
contourf(T_initial,'ShowText','on');
colormap(jet);
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title({[ 'After ' num2str(ct*dt) ' Seconds'];...
'Iterations using the SOR ITERATION' ; ...
['With relaxation parameter \alpha = ' num2str(alpha)]});
set(gca,'XTick',[],'YTick',[]);
% update movie frame
writeVideo(mov,getframe(gcf));
% pausing the loop to see the relaxation happening
pause(0.001);
end
end
% write out the movie and close the handle
close(mov);
% Plot the final graph after about 1400 iterations
% Plotting the surface to see how the plot looks like
figure(2);clf;
contourf(T_initial,'ShowText','on');
colormap(jet);
colorbar;
% making the plot look a bit nicer
set(gca,'Fontsize',20);
set(gca,'YDir','reverse');
title({[ 'After ' num2str(ct*dt) ' Seconds'];...
'Iterations using the SOR ITERATION' ; ...
['With relaxation parameter \alpha = ' num2str(alpha)]});
set(gca,'XTick',[],'YTick',[]);
else
% Printing an error message
error('Please select from the choices above only !!!!');
endThe following results were obtained for each iterative solver:
JACOBI SOLVER :
After about 0.29 seconds in the main time loop, the following graph and animation of the converged result were obtained.
The following animation was obtained :
GAUSS-SEIDEL SOLVER:
After about 0.29 seconds in the main time loop, the following graph and animation of the converged result were obtained.
The following animation was obtained :
SUCCESSIVE-OVER-RELAXATION (SOR) SOLVER
After about 0.29 seconds in the main time loop, the following graph and animation of the converged result were obtained for relaxation parameter `alpha = 1.4`:
The following animation was obtained:
CONCLUSION: As it can be seen that in all the three methods we are reaching to the same solution in about 0.3 seconds. It reaches the steady-state so fast because of the thermal diffusivity is quite high as compared to ones encountered in everyday life. Well, the main difference between the application of these methods is the way it calculates the Temperature values at each step. Generally, Gauss-Seidel would lead to the faster calculation of the values per spatial loop than the Jacobi and the successive-over-relaxation methods. However, a careful selection of the relaxation parameter alpha in the SOR technique can lead to very fast convergence. The reasons for these have already been discussed in the previous section on the steady-state analysis of the 2-D Heat Conduction equation.