Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
------------------------------------------------------------------- PART – A: STEADY STATE HEAT CONDUCTION EQUATION --------------------------------------------------------------------
1. Aim:
To solve the steady state heat conduction equation using iterative implicit techniques (Jacobi, Gauss Seidel and Successive Over-Relaxation). The Boundary conditions for the problem are as follows:
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
Absolute error criterion is 1e-4.
2. Theory and Governing Equations:
The general equation for heat conduction equation is:
`del/(delx)(k(delT)/(delx))+del/(dely)(k(delT)/(dely))+del/(delz)(k(delT)/(delz)) + S = rhoC_p(delT)/(delt)`
The equation for 2D heat conduction equation (without heat source) is:
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=(rhoC_p)/k(delT)/(delt)`
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=1/alpha(delT)/(delt)`
If the time derivative term becomes zero, we get steady state equation as:
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
After applying second order central difference to spatial derivative, we get:
`(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(dx^2)+(T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(dy^2)=0`
After rearranging the terms, we get:
`T_(i,j)=1/k((T_(i-1,j)+T_(i+1,j))/dx^2+(T_(i,j-1)+T_(i,j+1))/dy^2)` where `k=2((dx^2+dy^2)/(dx^2dy^2))`
The above equation can be solved by converging iterative methods which are as follows:
- For Jacobi method, for the RHS we use all temperature values from values from old convergence iteration. z denotes the iteration number from which temperature value is taken for calculation.
`T_(i,j)^(z+1)=1/k((T_(i-1,j)^z+T_(i+1,j)^z)/dx^2+(T_(i,j-1)^z+T_(i,j+1)^z)/dy^2)` - For Gauss-Seidel method, we use the new values for the nodes where we have already calculated Temperature in current convergence iteration and old ones for where calculation is still to be done. z denotes the iteration number from which temperature value is taken for calculation.
`T_(i,j)^(z+1)=1/k((T_(i-1,j)^(z+1)+T_(i+1,j)^z)/dx^2+(T_(i,j-1)^(z+1)+T_(i,j+1)^z)/dy^2)` - For SOR method, we use either Jacobi or Gauss-Seidel method, and scale up the convergence by a relaxation factor. If relaxation factor is less than 1, it is called under-relaxation. If it is greater than 1, it is called over-relaxation. `omega` is the relaxation factor in the equation.
`T_(i,j)^(z+1)=(1-omega)*T_(i,j)^z + omega*(Tgs)_(i,j)`
3. MATLAB Code:
clear all
close all
clc
global nx ny dx dy tol err k; %Global scope set to variables so they can be accessed inside functions
nx=15; %Number of grid points in x direction
ny=15; %Number of grid points in y direction
L=1; %Length of domain in both x and y directions
dx=L/(nx-1); %Grid Spacing in x direction
dy=L/(ny-1); %Grid Spacing in y direction
T=298*ones(ny,nx); %Initializing temperature matrix with room temperature at all nodes
T(1,:)=600; %Top Boundary Condition
T(ny,:)=900; %Bottom Boundary Condition
T(:,1)=400; %Left Boundary Condition
T(:,nx)=800; %Right Boundary Condition
T(1,1)=(600+400)/2; %Top Left Corner Condition
T(1,nx)=(600+800)/2; %Top Right Corner Condition
T(ny,1)=(900+400)/2; %Bottom Left Corner Condition
T(ny,nx)=(900+800)/2; %Bottom Right Corner Condition
T_initial=T; %Copy of initial temperature conditions
err=10^6; %Error is the max of absolute differences between values in successive iterations
tol=10^-4; %Tolerance is the maximum error that can be allowed between successive iterations
omega=1.64; %Relaxation factor for Successive Over-Relaxation Method
x=linspace(0,L,nx); %X-coordinate of nodes
y=linspace(0,L,ny); %Y-coordinate of nodes
k=2*(dx^2+dy^2)/(dx^2*dy^2);
%------------------------------------------------ JACOBI METHOD ------------------------------------------------
function [T,tJ,iter_J] = jacobi(T)
global nx ny dx dy tol err k;
iter_J=0; %Initializes iteration counter for Jacobi
tic; %Starts timer
while(err>tol) %CONVERGENCE LOOP
T_old=T; %Initializes temperature matrix corresponding to old convergence iteration with T
for i=2:nx-1 %NODAL LOOP
for j=2:ny-1
T(i,j)=((T_old(i-1,j)+T_old(i+1,j))/dx^2+(T_old(i,j-1)+T_old(i,j+1))/dy^2)/k;
endfor
endfor
err=max(max(abs(T_old-T)));
iter_J=iter_J+1;
endwhile
tJ=toc; %Stops timer and elapsed time for jacobian saved in tJ
endfunction
%--------------------------------------------- GAUSS SEIDEL METHOD ---------------------------------------------
function [T,tGS,iter_GS] = gauss(T)
global nx ny dx dy tol err k;
iter_GS=0; %Initializes iteration counter for Gauss Seidel
tic; %Starts timer
while(err>tol) %CONVERGENCE LOOP
T_old=T; %Initializes temperature matrix corresponding to old convergence iteration with T
for i=2:nx-1 %NODAL LOOP
for j=2:ny-1
T(i,j)=((T(i-1,j)+T_old(i+1,j))/dx^2+(T(i,j-1)+T_old(i,j+1))/dy^2)/k;
endfor
endfor
err=max(max(abs(T_old-T)));
iter_GS=iter_GS+1;
endwhile
tGS=toc; %Stops timer and elapsed time for Gauss Seidel saved in tGS
endfunction
%-------------------------------------- SUCCESSIVE OVER-RELAXATION METHOD --------------------------------------
function [T,tSOR,iter_SOR] = sor(T,omega)
global nx ny dx dy tol err k;
iter_SOR=0; %Initializes iteration counter for SOR
tic; %Starts timer
while(err>tol) %CONVERGENCE LOOP
T_old=T; %Initializes temperature matrix corresponding to old convergence iteration with T
for i=2:nx-1 %NODAL LOOP
for j=2:ny-1
T_GS=((T(i-1,j)+T_old(i+1,j))/dx^2+(T(i,j-1)+T_old(i,j+1))/dy^2)/k;
T(i,j)=(1-omega)*T_old(i,j)+omega*T_GS;
endfor
endfor
err=max(max(abs(T_old-T)));
iter_SOR=iter_SOR+1;
endwhile
tSOR=toc; %Stops timer and elapsed time for SOR saved in tSOR
endfunction
%------------------------------------- CALLING THE FUNCTIONS AND PLOTTING --------------------------------------
for i=1:3
err=10^6; %Global error value reset before calling each function
if i==1
[T,time,iter]=jacobi(T_initial);
elseif i==2
[T,time,iter]=gauss(T_initial);
elseif i==3
[T,time,iter]=sor(T_initial,omega);
endif
figure(i); %Opens a new figure window for each method
clabel(contourf(x,fliplr(y),T,9)); %Plots the contour and inserts 9 contour labels
colormap(jet); %To change color map type to jet
colorbar; %Shows the color legend as a bar
xlabel("X-direction"), ylabel("Y-direction");
title({"Temperature Distribution for steady state equation";sprintf("Iterations: %d Execution Time: %f",iter,time)})
endfor4. Explanation of MATLAB Code:
- Initially, all the constants and parameters like number of grid points (nx, ny), Length of domain (L), Initial temperature conditions (T), error (err), tolerance (tol) and relaxation factor for SOR method (omega) are defined. From these values, the grid spacing (dx, dy) and row matrices containing x and y-coordinate of nodes (x, y) are defined.
- Next, we define three functions, one each for obtaining steady state solution using Jacobi, Gauss Seidel and SOR method. They take the initial temperature conditions as parameters (T_initial) and return the temperature at steady state (T), execution time (tJ/tGS/tSOR) and number of iterations performed (iter_J/iter_GS/iter_SOR).
- In the function the iteration counter is initialized to 0 (iter) and the timer is started for calculating execution time (tic). A convergence loop is created using while statement that keeps iterating till the error value is no longer greater than tolerance.
- Inside the convergence loop, a nodal loop is created using two nested for loops to access the x and y coordinates of the nodes. Within the nested loop, the discretized steady state equation corresponding to each method is solved using temperature values at old iteration (T_old) to obtain new values. This is repeated till the solution is converged such that error is less than tolerance.
- The error is the maximum value among the absolute difference of current and previous iteration values. Iteration counter is increased by 1 (iter). After the while loop ends, the timer is stopped (toc).
- After all three functions have been defined, a ‘for’ loop running from i=1 to 3 is used to call each of the three functions using if condition particular to each i value. The error value is reset before calling each function.
- A new figure window is opened after calling each function and a contour is plotted of Temperature for each (x, y) value. The y-coordinates are flipped using fliplr command while plotting because if we don’t flip, the first value of y will be plotted at bottom and last at top and plot will be flipped upside down. Colormap command is used to select the type of color map in contour. Colorbar shows the contour values corresponding to a color as a bar. The axis labels and title have been mentioned accordingly.
5. Output:
1. Steady state equation solved using Jacobi method (Implicit):
2. Steady state equation solved using Gauss-Seidel method (Implicit):
3. Steady state equation solved using Successive Over-Relaxation method applied on Gauss-Seidel (Implicit):
6. Inference from results:
From the above contour plots, we see that all methods give the same Temperature Distribution and steady state has been achieved successfully through each one of them. However, we can see that the number of iterations and execution time is the lowest for SOR method, highest for Jacobi method and it is in between the former two for Gauss-Seidel method. These results were expected as Jacobi method uses all values from previous iterations whereas Gauss- Seidel uses updated values for the points where Temperature has already been computed. And hence convergence is achieved in lesser number of iterations. And in the case of SOR the basic equation is similar to Gauss-Seidel but a relaxation factor greater than 1 has been used to speed up or scale up the convergence and hence even lesser iterations are required to achieve it. The most suitable relaxation factor was used by manually varying the values between 1 and 2 and the one which took the least number of iterations to converge (1.64 in this case) was chosen for showing the results to display the computational difference between the methods.
Number of Iterations: Jacobi > Gauss-Seidel > SOR
Execution Time : Jacobi > Gauss-Seidel > SOR
----------------------------------------------------------------- PART – B: TRANSIENT STATE HEAT CONDUCTION EQUATION ------------------------------------------------------------------
1. Aim:
To solve the transient state heat conduction equation using explicit techniques and iterative implicit techniques (Jacobi, Gauss Seidel and Successive Over-Relaxation). The Boundary conditions for the problem are as follows:
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
Absolute error criterion is 1e-4.
2. Theory and Governing Equations:
The general equation for heat conduction equation is:
`del/(delx)(k(delT)/(delx))+del/(dely)(k(delT)/(dely))+del/(delz)(k(delT)/(delz)) + S = rhoC_p(delT)/(delt)`
The equation for 2D heat conduction equation (without heat source) is:
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=(rhoC_p)/k(delT)/(delt)`
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=1/alpha(delT)/(delt)`
After applying forward difference to time derivative and second order central difference to spatial derivative, we get:
`(T_(i-1,j)-2T_(i,j)+T_(i+1,j))/(dx^2)+(T_(i,j-1)-2T_(i,j)+T_(i,j+1))/(dy^2)=1/alpha((T_(i,j)^(n+1)-T_(i,j)^n)/dt)`
After rearranging the terms, we get:
`k1*(T_(i-1,j)-2T_(i,j)+T_(i+1,j))+k2*(T_(i,j-1)-2T_(i,j)+T_(i,j+1))=T_(i,j)^(n+1)-T_(i,j)^n`
where `k1=alphadt/dx^2` and`k2=alphadt/dx^2`
As we can see from above equation that in the RHS, one temperature is at (n+1)th time step and other is at nth time step. But the time step for the temperatures in LHS hasn't been indicated. So the above equation can be solved in following two ways depending on which temperatures are considered for calculation:
1. Explicit Method: If we take the temperatures in LHS at nth time step, then the values are already known and we can put the values and directly get the values at (n+1)th time step. This method of solving by taking all temperatures in LHS from previous time step is called EXPLICIT METHOD. The equation for explicit method is:
`k1*(T_(i-1,j)^n-2T_(i,j)^n+T_(i+1,j)^n)+k2*(T_(i,j-1)^n-2T_(i,j)^n+T_(i,j+1)^n)=T_(i,j)^(n+1)-T_(i,j)^n`
Rearranging,
`T_(i,j)^(n+1)=(1-2k1-2k2)*T_(i,j)^n+k1*(T_(i-1,j)^n+T_(i+1,j)^n)+k2*(T_(i,j-1)^n+T_(i,j+1)^n)`
2. Implicit Method: If we take the temperatures in LHS at (n+1)th time step, then the values are not known and we need to obtain the values by using different convergence iteration methods. This method of solving by using all temperatures in LHS from current time step is called IMPLICIT METHOD. The equation for implicit method is:
`k1*(T_(i-1,j)^(n+1)-2T_(i,j)^(n+1)+T_(i+1,j)^(n+1))+k2*(T_(i,j-1)^(n+1)-2T_(i,j)^(n+1)+ T_(i,j+1)^(n+1))=T_(i,j)^(n+1)-T_(i,j)^n`
Rearranging,
`T_(i,j)^(n+1)=(T_(i,j)^n+k1*(T_(i-1,j)^(n+1)+T_(i+1,j)^(n+1))+k2*(T_(i,j-1)^(n+1)+T_(i,j+1)^(n+1)))/(1+2k1+2k2)`
This implicit equation needs to be solved by converging iterative methods which are as follows:
- For Jacobi method, we use all temperature values from values from old convergence iteration. z denotes the iteration number from which temperature value is taken for calculation. For achieving convergence in the (z+1)th iteration at (n+1)th time step, equation is:
`(T_(i,j)^(z+1))^(n+1)=((T_(i,j))^n+k1*(T_(i-1,j)^z+T_(i+1,j)^z)^(n+1)+k2*(T_(i,j-1)^z+T_(i,j+1)^z)^(n+1))/(1+2k1+2k2)`
- For Gauss-Seidel method, we use the new values for the nodes where we have already calculated Temperature in current convergence iteration and old ones for where calculation is still to be done. z denotes the iteration number from which temperature value is taken for calculation. For achieving convergence in the (z+1)th iteration at (n+1)th time step, equation is:
`(T_(i,j)^(z+1))^(n+1)=((T_(i,j))^n+k1*(T_(i-1,j)^(z+1)+T_(i+1,j)^z)^(n+1)+k2*(T_(i,j-1)^(z+1)+T_(i,j+1)^z)^(n+1))/(1+2k1+2k2)`
- For SOR method, we use either Jacobi or Gauss-Seidel method, and scale up the convergence by a relaxation factor. If relaxation factor is less than 1, it is called under-relaxation. If it is greater than 1, it is called over-relaxation. `omega` is the relaxation factor in the equation. For achieving convergence in the (z+1)th iteration at (n+1)th time step, equation is:
`(T_(i,j)^(z+1))^(n+1)=(1-omega)*T_(i,j)^z + omega*(Tgs_(i,j))^(n+1)`
3. MATLAB Code:
clear all
close all
clc
global nx ny nt k1 k2 tol err; %Global scope set to variables so they can be accessed inside functions
nx=15; %Number of grid points in x direction
ny=15; %Number of grid points in y direction
L=1; %Length of domain in both x and y directions
dx=L/(nx-1); %Grid Spacing in x direction
dy=L/(ny-1); %Grid Spacing in y direction
dt=0.001; %Time step size
nt=500; %Number of time steps
alpha=1.25; %Thermal Diffusivity
T=298*ones(ny,nx); %Initializing temperature matrix with room temperature at all nodes
T(1,:)=600; %Top Boundary Condition
T(ny,:)=900; %Bottom Boundary Condition
T(:,1)=400; %Left Boundary Condition
T(:,nx)=800; %Right Boundary Condition
T(1,1)=(600+400)/2; %Top Left Corner Condition
T(1,nx)=(600+800)/2; %Top Right Corner Condition
T(ny,1)=(900+400)/2; %Bottom Left Corner Condition
T(ny,nx)=(900+800)/2; %Bottom Right Corner Condition
T_initial=T; %Copy of initial temperature conditions
err=10^6; %Error is the max of absolute differences between values in successive iterations
tol=10^-4; %Tolerance is the maximum error that can be allowed between successive iterations
omega=1.2; %Relaxation factor for Successive Over-Relaxation Method
k1=alpha*dt/(dx^2); %CFL number in x-direction
k2=alpha*dt/(dy^2); %CFL number in y-direction
x=linspace(0,L,nx); %X-coordinate of nodes
y=linspace(0,L,ny); %Y-coordinate of nodes
%------------------------------------------------ JACOBI METHOD ------------------------------------------------
function [T,tJ,iter_J] = jacobi(T)
global nx ny nt k1 k2 tol err;
iter_J=0; %Initializes iteration counter for Jacobi
tic; %Starts timer
for p=1:nt %TIME LOOP
T_prv=T; %Initializes temperature matrix corresponding to previous time step with T
err=10^6; %Initializes error value reset for next time step
while(err>tol) %CONVERGENCE LOOP
T_old=T; %Initializes temperature matrix corresponding to old convergence iteration with T
for i=2:nx-1 %NODAL LOOP
for j=2:ny-1
T(i,j)=(T_prv(i,j) + k1*(T_old(i-1,j)+T_old(i+1,j)) + k2*(T_old(i,j-1)+T_old(i,j+1))) / (1+2*k1+2*k2);
endfor
endfor
err=max(max(abs(T_old-T)));
iter_J=iter_J+1;
endwhile
endfor
tJ=toc; %Stops timer and elapsed time for Jacobi saved in tJ
endfunction
%--------------------------------------------- GAUSS SEIDEL METHOD ---------------------------------------------
function [T,tGS,iter_GS] = gauss(T)
global nx ny nt k1 k2 tol err;
iter_GS=0; %Initializes iteration counter for Gauss Seidel
tic; %Starts timer
for p=1:nt %TIME LOOP
T_prv=T; %Initializes temperature matrix corresponding to previous time step with T
err=10^6; %Initializes error value reset for next time step
while(err>tol) %CONVERGENCE LOOP
T_old=T; %Initializes temperature matrix corresponding to old convergence iteration with T
for i=2:nx-1 %NODAL LOOP
for j=2:ny-1
T(i,j)=(T_prv(i,j) + k1*(T(i-1,j)+T_old(i+1,j)) + k2*(T(i,j-1)+T_old(i,j+1))) / (1+2*k1+2*k2);
endfor
endfor
err=max(max(abs(T_old-T)));
iter_GS=iter_GS+1;
endwhile
endfor
tGS=toc; %Stops timer and elapsed time for Gauss Seidel saved in tGS
endfunction
%-------------------------------------- SUCCESSIVE OVER-RELAXATION METHOD --------------------------------------
function [T,tSOR,iter_SOR] = sor(T,omega)
global nx ny nt k1 k2 tol err;
iter_SOR=0; %Initializes iteration counter for SOR
tic; %Starts timer
for p=1:nt %TIME LOOP
T_prv=T; %Initializes temperature matrix corresponding to previous time step with T
err=10^6; %Initializes error value reset for next time step
while(err>tol) %CONVERGENCE LOOP
T_old=T; %Initializes temperature matrix corresponding to old convergence iteration with T
for i=2:nx-1 %NODAL LOOP
for j=2:ny-1
T_GS=(T_prv(i,j) + k1*(T(i-1,j)+T_old(i+1,j)) + k2*(T(i,j-1)+T_old(i,j+1))) / (1+2*k1+2*k2);
T(i,j)=(1-omega)*T_old(i,j) + omega*T_GS;
endfor
endfor
err=max(max(abs(T_old-T)));
iter_SOR=iter_SOR+1;
endwhile
endfor
tSOR=toc; %Stops timer and elapsed time for SOR saved in tSOR
endfunction
%----------------------------------------------- EXPLICIT METHOD -----------------------------------------------
function [T,tEX,iter_EX] = explicit(T)
global nx ny nt k1 k2 tol err;
iter_EX=0; %Initializes iteration counter for Explicit
tic; %Starts timer
for p=1:nt %TIME LOOP
T_prv=T; %Initializes temperature matrix corresponding to previous time step with T
for i=2:nx-1 %NODAL LOOP
for j=2:ny-1
T(i,j)=(1-2*k1-2*k2)*(T_prv(i,j)) + k1*(T_prv(i-1,j)+T_prv(i+1,j)) + k2*(T_prv(i,j-1)+T_prv(i,j+1));
endfor
endfor
iter_EX=iter_EX+1;
endfor
tEX=toc; %Stops timer and elapsed time for Explicit saved in tEX
endfunction
%------------------------------------- CALLING THE FUNCTIONS AND PLOTTING --------------------------------------
for i=1:4
if i==1
[T,time,iter]=jacobi(T_initial);
elseif i==2
[T,time,iter]=gauss(T_initial);
elseif i==3
[T,time,iter]=sor(T_initial,omega);
elseif i==4
[T,time,iter]=explicit(T_initial);
endif
figure(i); %Opens a new figure window for each method
clabel(contourf(x,fliplr(y),T,9)); %Plots the contour and inserts 9 contour labels
colormap(jet); %To change color map type to jet
colorbar; %Shows the color legend as a bar
xlabel("X-direction"), ylabel("Y-direction");
title({"Temperature Distribution for transient state equation";sprintf("Iterations: %d Execution Time: %f",iter,time)})
endfor4. Explanation of MATLAB Code:
- Initially, all the constants and parameters like number of grid points (nx, ny), Length of domain (L), time step size (dt), number of time steps(nt), thermal diffusivity (alpha), initial temperature conditions (T), error (err), tolerance (tol) and relaxation factor for SOR method (omega) are defined. From these values, the grid spacing (dx, dy), CFL number for x and y direction (k1, k2) and row matrices containing x and y-coordinate of nodes (x, y) are defined.
- Next, we define four functions, one each for obtaining transient state solution using Jacobi, Gauss Seidel, SOR method and Explicit method. They take the initial temperature conditions as parameters (T_initial) and return the temperature at steady state (T), execution time (tJ/tGS/tSOR/tEX) and number of iterations performed (iter_J/iter_GS/iter_SOR/iter_EX).
- In the function the iteration counter is initialized to 0 (iter) and the timer is started for calculating execution time (tic). A time loop is first created to control the time step variable and to help calculate Temperature at each intermediate time step.
- A convergence loop is created inside the time loop using ‘while’ statement that keeps iterating till the error value is no longer greater than tolerance. The error value is reset inside the time loop before starting the while loop. The convergence ‘while’ loop is not required in explicit method as known values are only used.
- Inside the convergence loop, a nodal loop is created using two nested for loops to access the x and y coordinates of the nodes. Within the nested loop, the discretized transient state equation corresponding to each method is solved to obtain new temperature values. For the first three methods (implicit), this is repeated till the solution is converged such that error is less than tolerance. In case of explicit method, there is no repetition as convergence is not required. The solution at all nodal points is calculated only once for each time step in explicit method.
- The error is the maximum value among the absolute difference of current and previous iteration values. Iteration counter is increased by 1 (iter) after the nodal loop. Where the time loop ends, the timer is stopped (toc).
- After all four functions have been defined, a ‘for’ loop running from i=1 to 4 is used to call each of the four functions using if condition particular to each i value.
- A new figure window is opened after calling each function and a contour is plotted of Temperature for each (x, y) value. The y-coordinates are flipped using fliplr command while plotting because if we don’t flip, the first value of y will be plotted at bottom and last at top and plot will be flipped upside down. Colormap command is used to select the type of color map in contour. Colorbar shows the contour values corresponding to a color as a bar. The axis labels and title have been mentioned accordingly.
5. Output:
1. Transient state equation solved using Jacobi method (Implicit):
2. Transient state equation solved using Gauss-Seidel method (Implicit):
3. Transient state equation solved using Successive Over-Relaxation method applied on Gauss-Seidel (Implicit):
4. Transient state equation solved using Explicit Method:
6. Inference from results:
From the above contour plots, we see that even though the equation solved is of transient state, eventually steady state is achieved through each of the methods. And all of them show the same temperature distribution and are hence validate each other. However, we can see that the number of iteration and execution time for Explicit method is less than Implicit methods. Among the implicit methods, number of iteration and execution time was least for SOR method, highest for Jacobi method and it is in between the former two for Gauss-Seidel method. These results were expected as in explicit method, the value at current step is calculated using previous time step which are already known. So, unlike implicit method, there is no need for iterative convergence. Hence it takes less time and lesser iterations. Among implicit methods, Jacobi method uses all values from previous iterations whereas Gauss-Seidel uses updated values for the points where Temperature has already been computed. And the computed values are closer to the required solution and hence convergence is achieved in lesser number of iterations. In the case of SOR the basic equation is similar to Gauss-Seidel but a relaxation factor greater than 1 has been used to speed up or scale up the convergence and hence even lesser iterations are required to achieve it. The most suitable relaxation factor was used by manually varying the values between 1 and 2 and the one which took the least number of iterations to converge (1.2 in this case) was chosen for showing the results to display the computational difference between the methods.
Number of Iterations: Jacobi > Gauss-Seidel > SOR > Explicit
Execution Time : Jacobi > Gauss-Seidel > SOR > Explicit
7. Conclusion:
After performing above challenge for both steady state and transient state, it can be concluded that among implicit solvers, if we use a suitable value for relaxation factor, SOR can be most efficient and Jacobi is least efficient. And if we compare explicit and implicit, although in this case explicit achieved results faster with lesser iterations, it may not always be possible. Implicit solvers work for any combination of values of diffusivity coefficient (alpha), grid spacing (dx) and time step size (dt). But this is not the case with explicit solver. Explicit solvers are stable only if`(alphadt/dx^2)+(alphadt/dy^2)+(alphadt/dz^2)<=0.5`. If this value goes above 0.5, it becomes unstable and the solution blows up. Because of these factors, the explicit method puts a limitation on the value of dt and dx. And hence it may not ideally be the best choice. Considering all these results, SOR method is the most efficient and versatile solving technique.