Solving and Simulating the steady & unsteady 2-D Heat Conduction equation in MATLAB
Solving the Steady and Unsteady 2-D Heat Conduction Equation in MATLAB
Problem Specification:-
- Solving the steady-state 2-D heat conduction equation implicitly using point-iterative methods; such as Jacobi, Gauss-Seidel, and Successive Over-Relaxation (SOR).
- Solving the unsteady (transient)-state 2-D heat conduction equation using the explicit and implicit schemes (implementing Jacobi, Gauss-Seidel, and Successive Over-Relaxation methods).
Assumptions:-
- Assume that the domain is a unit square.
- Assume the number of grid points (nodes) along the x-direction is equal to the number of grid points along the y-direction.
- The boundary condition for steady and transient case: Left boundary temperature=400K, Right boundary temperature=800K, Top boundary temperature=600K, and Bottom boundary temperature=900K.
- Absolute error criteria (tolerance) =1e-4
- Relaxation factor, `ω`=1.2
- Thermal diffusivity, `α`=1.4
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1. Solving the Steady-State 2-D Heat Conduction equation implicitly using point-iterative methods; such as Jacobi, Gauss-Seidel and Successive Over-Relaxation (SOR)
THEORY:-
Steady-state 2-D heat conduction equation:`(∂^2 T)/(∂x^2)+(∂^2 T)/(∂y^2 )=0 or ∇^2 T=0`
To solve the above equation, we need to discretize the above equation. As we can see that the above equation is an elliptic partial differential equation where `x` and `y` play the role of spatial coordinate. We also know it as ‘Laplace Equation’ which in heat transfer describes the steady-state temperature distribution in a particular domain.
Since the above equation contains the second-order derivative of temperature and the temperature inevitably going to influence all the nearby points. So, to calculate the temperature at any point; we have to gather info from all the neighboring points as shown in Fig. (1). That means we have to use the Central Difference Scheme to discretize the above equation.
Fig. (1)
`(T_L-2T_P+T_R)/(∆x^2 )+(T_T-2T_P+T_B)/(∆y^2 )=0`
Taking `2T_P` terms one side and rest of the terms another side:
`(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 )` ………(I)
Let, `k=[(2(∆x^2+∆y^2 ))/(∆x^2 ∆y^2 )]`
So, `T_P.k=(T_L+T_R)/(∆x^2 )+(T_T+T_B)/(∆y^2 )``T_P=1/k [(T_L+T_R)/(∆x^2 )+(T_T+T_B)/(∆y^2 )]`
We’ll start with an initial value for temperature and then we would use iterative solvers like Jacobi, Gauss-Seidel, or SOR; to find the solution for a system of equation.
But, for this problem, since we are assuming the number of grid points in the x-direction is equal to the number of grid points in the y-direction, `nx=ny`; and the grid spacing `∆x=∆y=s`(let). Hence the equation (I) can be written as:
`T_P [(2(s^2+s^2 ))/(s^2.s^2 )]=(T_L+T_R)/s^2 +(T_T+T_B)/s^2`
`T_P [(4(s^2 ))/(s^2.s^2 )]=(T_L+T_R+T_T+T_B)/s^2`
`T_P=(T_L+T_R+T_T+T_B)/4` or, `T_(i,j)=(T_(i-1,j)+T_(i+1,j)+T_(i,j+1)+T_(i,j-1))/4` ………(II)
Equation (II) is the final equation that we’ll solve for this problem. We will use three different iterative solvers to solve the above equation; i.e. Jacobi, Gauss-Seidel, and Successive Over-Relaxation (SOR). Each of the solvers solves the equation iteratively until the iterative error falls below the specified tolerance value (which is 1e-4 for this problem). The major difference between these iterative solvers is how it updates the values for the successive iteration through the values from the previous iterations.
Note: The superscript `k` in these below equations represents the iteration number.
Jacobi Iterative Solver: `T_(i,j)^(k+1)=(T_(i-1,j)^k+T_(i+1,j)^k+T_(i,j-1)^k+T_(i,j+1)^k)/4`
In Jacobi solver, all the values on the right-hand side term are taking the information from the previous term as shown in the above equation.
Gauss-Seidel Iterative Solver: `T_(i,j)^(k+1)=(T_(i-1,j)^(k+1)+T_(i+1,j)^k+T_(i,j-1)^(k+1)+T_(i,j+1)^k)/4`
In Gauss-Seidel solver, the values (like `T_(i-1,j)^(k+1)` and `T_(i,j-1)^(k+1)`) that have already been updated in the current iteration are used and the rest of the values (like `T_(i+1,j)^k` and `T_(i,j+1)^k`) are taken from the previous iteration.
Successive Over-Relaxation (SOR) Iterative Solver: `T_(i,j)^(k+1)=(1-ω) T_(i,j)^k+ω[(T_(i-1,j)^(k+1)+T_(i+1,j)^k+T_(i,j-1)^(k+1)+T_(i,j+1)^k)/4]`
We can use the SOR solver on Jacobi or Gauss-Seidel method. This is an enhanced iterative technique that can enhance the performance of the Jacobi or Gauss-Seidel method.
Here in this problem, we’re applying SOR on the Gauss-Seidel method. That’s why we have included the Gauss-Seidel equation inside the parenthesis on the right-hand side of the above equation.
Here ‘`ω`’ is a scaling factor known as the Relaxation factor. If ‘`ω>1`’, then it makes the changes a bit more aggressive; that way our iterative solver is over tuned and this situation called ‘over-relaxation’. If ‘`ω<1`’, then it makes our solver less aggressive; that is what we called as ‘under-relaxation’. Here for this problem, we’re taking ‘`ω>1`’ (`i.e. ω=1.2`); so that our result can converge at a faster rate.
Before writing the code, it needs some clarification for the problem statement. In this problem, the domain is a square with a unit length. We divide the domain into an equal number of grid points (nodes) (i.e. `nx=ny`) and grid spacing (`∆x=∆y=s`) for both x & y-directions. The domain also has a specified boundary condition—as shown in Fig. (2)—to get a unique solution to the heat conduction problem. We define boundary conditions as:-
Left boundary: T=400K
Right boundary: T=800K
Top boundary: T=600K
Bottom boundary: T=900K
Fig. (2)
Note: The values of temperature at the corner nodes are calculated by taking the average of the values from the surrounding edges.
MATLAB code used to solve ‘Steady State 2-D Heat Conduction Equation’:-
% Solving the steady-state 2D Heat conduction PDE implicitly using different point iterative techniques
% (d^2T/dx^2)+(d^2T/dy^2)=0
%Final Discretized equation that we're going to solve
% T(i,j)=((T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1))/4
clear all
close all
clc
%variable inputs and defining domain
%grid points along x-direction
nx=10;
%grid points along y-direction
ny=nx;
%Length of domain
L_x=1; L_y=1; %length along x & y-direction respectively
%creating mesh
x=linspace(0,L_x,nx);
y=linspace(0,L_y,ny);
%relaxation factor,w=1.2
w=1.2;
%Initializing system temperature matrix
T=300*ones(nx,ny);
%Applying Boundary Conditions:
T(end,2:end-1)=600; %Top temperature
T(1,2:end-1)=900; %Bottom temperature
T(2:end-1,1)=400; %Left temperature
T(2:end-1,end)=800; %Right temperature
%For corner nodes; the temperature can be calculated as an average of surrounding nodes
T(end,1)=(600+400)/2;
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,end)=(600+800)/2;
[xx,yy]=meshgrid(x,y);
%creating copy of temperature, T
Told=T;
%operator for calling different solvers
%(for Jacobi; iterative_solver=1, for Gauss-Seidel; iterative_solver=2 and for SOR; iterative_solver=3)
iterative_solver=1;
error=9e9;
%tolerance value=1e-4
tol=1e-4;
iteration_counter=1;
%Jacobi Solver:-
tic;
if iterative_solver==1
%convergence loop
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=(Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1))/4;
end
end
error=max(max(abs(Told-T)));
Told=T; %updating value of temperature for each iteration
iteration_counter=iteration_counter+1;
end %end of convergence loop
toc;
time=toc;
%Plots for Jacobi Iterative Solver
figure(1)
[x_y,temp]=contourf(xx,yy,T);
colorbar
colormap('jet')
clabel (x_y,temp)
title_text=sprintf('For JACOBI: Number of Iterations=%d, Total Time=%f seconds',iteration_counter, time);
title(title_text)
xlabel('X')
ylabel('Y')
end %end of Jacobi iterative solver
%Gauss-Seidel Solver:-
tic;
if iterative_solver==2
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1))/4;
end
end
error=max(max(abs(Told-T)));
Told=T; %updating value of temperature for each iteration
iteration_counter=iteration_counter+1;
end %end of convergence loop
toc;
time=toc;
%Plots for Gauss-Seidel Iterative Solver
figure(1)
[x_y,temp]=contourf(xx,yy,T);
colorbar
colormap('jet')
clabel (x_y,temp)
title_text=sprintf('For Gauss-Seidel: Number of Iterations=%d, Total Time=%f seconds',iteration_counter, time);
title(title_text)
xlabel('X')
ylabel('Y')
end %end of Gauss-Seidel iterative solver
%Successive Over-relaxation (SOR) solver:-
tic;
if iterative_solver==3
%convergence loop
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=(1-w)*Told(i,j)+w*((T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1))/4);
end
end
error=max(max(abs(Told-T)));
Told=T; %updating value of temperature for each iteration
iteration_counter=iteration_counter+1;
end %end of convergence loop
toc;
time=toc;
%Plots for SOR Iterative Solver
figure(1)
[x_y,temp]=contourf(xx,yy,T);
colorbar
colormap('jet')
clabel (x_y,temp)
title_text=sprintf('For SOR: Number of Iterations=%d, Total Time=%f seconds',iteration_counter, time);
title(title_text)
xlabel('X')
ylabel('Y')
end %end of SOR iterative solver
RESULTS & DISCUSSIONS:-
The above plots show the contours of variation of temperature with the position. As it is, the steady-state solution; there is no change of temperature w.r.t time, but there is a significant variation of temperature w.r.t space (x & y-direction).
From the above plots, we can see that the Jacobi solver takes ‘208-iterations’ to get converge; whereas the Gauss-Seidel solver takes almost half of that number of iterations, which is ‘112- iterations’ to get converge. The Successive Over-Relaxation (SOR) solver takes only ’75-iterations’ to get converged. The total computational time taken by each solver is Jacobi=0.222876 seconds; Gauss-Seidel=0.003918 seconds, and SOR=0.005740 seconds.
We can deduce it from the above-discussed data that the Jacobi solver is the slowest one as it takes much more time and the number of iterations to get converged. This happens because Jacobi solver updates the current iteration with the values from the previous (old) iterations.
The total number of iterations for the Gauss-Seidel solver suggests that it is faster than Jacobi solver as it gets converged in a lesser number of time and iterations and this happens because the Gauss-Seidel solver updates the current iteration values with some values already calculated in the current iteration.
The SOR solver provides a greater edge over the Gauss-Seidel or Jacobi solver as it determines how aggressively the current iteration gets updated from the values in the current iteration, which makes it much faster than the Gauss-Seidel or Jacobi solver.
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
2. Solving the Unsteady (Transient)-State 2-D heat conduction equation using Explicit and Implicit schemes
THEORY:-
The 2-D heat conduction equation is a partial differential equation that governs the heat transfer through a medium by thermal conduction. The equation is defined as:
`(∂T)/(∂t)+α[(∂^2 T)/(∂x^2) +(∂^2 T)/(∂y^2) ]=0`
where `α` = thermal diffusivity of the material
In the above equation `(∂T)/(∂t)` the term governs the transient or unsteady behavior of the equation; that denotes the change of temperature w.r.t time.
To discretize the above equation; for the time derivative term we’ll use a forward difference scheme and for the space derivative terms we’ll use a central difference scheme. We can solve transient problems in two ways, either ‘explicitly’ or ‘implicitly’; depending on the time step at which we choose the spatial derivatives.
First, we’ll discuss the explicit scheme where all the spatial derivatives are taken from the previous time-step. Then after we’ll discuss the implicit scheme.
Discretizing the unsteady-state 2-D heat conduction equation using ‘explicit’ scheme:
Note: Superscript n in these below equations denotes the time-step at which we take the particular value of temperature. For nodal notations refer to Fig. (1).
`(T_P^(n+1)-T_P^n)/(∆t)=α[(T_L^n-2T_P^n +T_R^n)/(∆x^2 )+(T_B^n-2T_P^n +T_T^n)/(∆y^2 )]`
`T_P^(n+1)=T_P^n+α.∆t[(T_L^n-2T_P^n +T_R^n)/(∆x^2 )+(T_B^n-2T_P^n +T_T^n)/(∆y^2 )]` ………(III)
For this problem, we’re assuming `∆x=∆y=s` (let).
On putting the above assumption in equation (III); we get:
`T_P^(n+1)=T_P^n+(α.∆t)/ s^2 [T_L^n-2T_P^n +T_R^n+T_B^n-2T_P^n +T_T^n ]`
`T_P^(n+1)=T_P^n+(α.∆t)/ s^2 [T_L^n-4T_P^n +T_R^n+T_B^n +T_T^n ]`
Let, `(α.∆t)/(s^2) =k_1` and take ‘`T_P^n`’ common in the above equation:
`T_P^(n+1)=T_P^n (1-4k_1)+k_1 [T_L^n +T_R^n+T_B^n +T_T^n ]`
Or, `T_(i,j)^(n+1)=T_(i,j)^n (1-4k_1 )+k_1 [T_(i-1,j)^n +T_(i+1,j)^n+T_(i,j-1)^n +T_(i,j+1)^n ]` ………(IV)
The equation (IV) is the final discretized equation that we will solve for the explicit scheme.
For stability, the coefficient of the term `T_(i,j)^n` on the right-hand side of the equation can’t be negative (always greater than or equal to 0). This comes from the Courant-Friedrich-Lewy (CFL) condition for stability. The term `k_1=(α.∆t)/( s^2 )` is called the CFL number. From this condition, we can calculate the maximum time-step that can be chosen for solving the transient 2-D heat conduction equation explicitly.
`1-4k_1≥0 or 1-4[(α.∆t)/(s^2)]≥0`
`(α.∆t)/s^2 ≤1/4`
Or, `∆t≤[s^2/(4.α)]`
The above-discussed stability condition must be satisfied to get a converged solution. As for the explicit scheme, time-step: `∆t=[s^2/(4.α)]` is chosen.
Discretizing the unsteady-state 2-D heat conduction equation using ‘implicit’ scheme:
`(T_P^(n+1)-T_P^n)/(∆t)=α[(T_L^(n+1)-2T_P^(n+1) +T_R^(n+1))/(∆x^2 )+(T_B^(n+1)-2T_P^(n+1) +T_T^(n+1))/(∆y^2 )]`
`T_P^(n+1)=T_P^n+α.∆t[(T_L^(n+1)-2T_P^(n+1) +T_R^(n+1))/(∆x^2 )+(T_B^(n+1)-2T_P^(n+1) +T_T^(n+1))/(∆y^2 )]` ………(V)
For this problem, we’re assuming `∆x=∆y=s`(let).
On putting the above assumption in equation (V); we get:
`T_P^(n+1)=T_P^n+(α.∆t)/ s^2 [T_L^(n+1)-2T_P^(n+1) +T_R^(n+1)+T_B^(n+1)-2T_P^(n+1) +T_T^(n+1) ]`
`T_P^(n+1)=T_P^n+(α.∆t)/s^2 [T_L^(n+1)-4T_P^(n+1) +T_R^(n+1)+T_B^(n+1) +T_T^(n+1) ]`
Let, `(α.∆t)/(s^2) =k_1` and take ‘`T_P^n`’ common in the above equation:
`T_P^(n+1)+k_1 (4T_P^(n+1))=T_P^n+k_1 [T_L^(n+1) +T_R^(n+1)+T_B^(n+1) +T_T^(n+1) ]`
`(1+4k_1)T_P^(n+1)=T_P^n+k_1 [T_L^(n+1) +T_R^(n+1)+T_B^(n+1) +T_T^(n+1) ]`
`T_P^(n+1)=(T_P^n)/(1+4k_1 )+k_1/(1+4k_1 ) [T_L^(n+1) +T_R^(n+1)+T_B^(n+1) +T_T^(n+1) ]`
Or, `T_(i,j)^(n+1)=(T_(i,j)^n)/(1+4k_1 )+k_1/(1+4k_1 ) [T_(i-1,j)^(n+1) +T_(i+1,j)^(n+1)+T_(i,j-1)^(n+1) +T_(i,j+1)^(n+1) ]` ………(VI)
Since the coefficient of term, `T_P^n` or `T_(i,j)^n` on the right-hand side of equation (VI) can never be negative. Hence, the implicit scheme is unconditionally stable for any value of `k_1`.
In equation (VI), the right-hand side values inside the parenthesis are also unknown; to solve those values we will use iterative solvers (Jacobi, Gauss-Seidel, and Successive Over-Relaxation).
Applying Jacobi solver:
`T_(i,j)^(n+1,k+1)=(T_(i,j)^n)/(1+4k_1 )+k_1/(1+4k_1 ) [T_(i-1,j)^(n+1,k) +T_(i+1,j)^(n+1,k)+T_(i,j-1)^(n+1,k) +T_(i,j+1)^(n+1,k) ]`
Applying Gauss-Seidel solver:
`T_(i,j)^(n+1,k+1)=(T_(i,j)^n)/(1+4k_1 )+k_1/(1+4k_1 ) [T_(i-1,j)^(n+1,k+1) +T_(i+1,j)^(n+1,k)+T_(i,j-1)^(n+1,k+1) +T_(i,j+1)^(n+1,k) ]`
Applying Successive Over-Relaxation (SOR) solver:
`T_(i,j)^(n+1,k+1)=(1-ω) T_(i,j)^n+ω[T_(i,j)^n/(1+4k_1 )+k_1/(1+4k_1 ) (T_(i-1,j)^(n+1,k+1) +T_(i+1,j)^(n+1,k)+T_(i,j-1)^(n+1,k+1) +T_(i,j+1)^(n+1,k))]`
Note: The superscript ‘n’ denotes the time-step and superscript ‘`k`’ denotes the iteration within a particular time-step. The temperature term (`T_(i,j)^n`) missing the superscript ‘`k`’ shows that the values of temperature are taken from the previous time-step; which are constant throughout the current time-step and do not iterate.
Before writing the code, it needs some clarification for the problem statement. In this problem, the domain is a square with a unit length. We divide the domain into an equal number of grid points (nodes) (`i.e. nx=ny`) and grid spacing (`∆x=∆y=s`) for both x & y-directions. The domain also has a specified boundary condition—as shown in Fig. (2)—to get a unique solution to the heat conduction problem. We define boundary conditions as:-
Left boundary: T=400K
Right boundary: T=800K
Top boundary: T=600K
Bottom boundary: T=900K
MATLAB code used to solve ‘Unsteady (Transient) State 2-D Heat Conduction Equation’:-
% Solving the unsteady (transient) state 2D Heat conduction PDE implicitly using different point iterative techniques
% (dT/dt)+alpha[(d^2T/dx^2)+(d^2T/dy^2)]=0
clear all
close all
clc
%variable inputs and defining domain
%grid points along x-direction
nx=10;
%grid points along y-direction
ny=nx;
%Length of domain
L_x=1; L_y=1; %length along x & y-direction respectively
%creating mesh
x=linspace(0,L_x,nx);
y=linspace(0,L_y,ny);
%grid spacing, dx=dy=s(let)
s=x(2)-x(1);
%thermal diffusivity=1.4
alpha=1.4;
%relaxation factor,w=1.2
w=1.2;
%Initializing system temperature matrix
T=300*ones(nx,ny);
%Applying Boundary Conditions:
T(end,2:end-1)=600; %Top temperature
T(1,2:end-1)=900; %Bottom temperature
T(2:end-1,1)=400; %Left temperature
T(2:end-1,end)=800; %Right temperature
%For corner nodes; the temperature can be calculated as an average of surrounding nodes
T(end,1)=(600+400)/2;
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,end)=(600+800)/2;
[xx,yy]=meshgrid(x,y);
%creating copy of temperature, T
T_prev_dt=T;
Told=T;
%operator for calling different solvers (for Explicit scheme; iterative_solver=1, for Implicit scheme; Jacobi: iterative_solver=2, for Gauss-Seidel: iterative_solver=3, and for SOR: iterative solver=4)
iterative_solver=1;
%tolerance value=1e-4
tol=1e-4;
iteration_counter=1;
%EXPLICIT SCHEME:-
tic;
if iterative_solver==1
%time-step
dt=(s^2)/(4*alpha);
k1=(alpha*dt)/(s^2);
for nt=1:1400 %time loop
error=9e9;
%convergence loop
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T_parenthesis=(Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1));
T(i,j)=(1-(4*k1))*T_prev_dt(i,j)+(k1*T_parenthesis);
end
end
error=max(max(abs(Told-T)));
Told=T; %updating value of temperature for each iteration
T_prev_dt=T; %updating previous time-step temperature's matrix.
iteration_counter=iteration_counter+1;
end %end of convergence loop
end %end of time loop
toc;
time=toc;
%Plots for Explicit scheme
figure(1)
[x_y,temp]=contourf(xx,yy,T);
colorbar
colormap('jet')
clabel (x_y,temp)
title_text=sprintf('For EXPLICIT SCHEME: Number of Iterations=%d, Total Time=%f seconds',iteration_counter, time);
title(title_text)
xlabel('X')
ylabel('Y')
end %end of Explicit scheme
%Implicit Scheme
%Jacobi Solver:-
tic;
if iterative_solver==2
%time-step
dt=0.02;
k1=(alpha*dt)/(s^2);
for nt=1:1400 %time loop
error=9e9;
%convergence loop
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T_parenthesis=Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1);
T(i,j)=T_prev_dt(i,j)/(1+(4*k1))+((k1)/(1+(4*k1)))*T_parenthesis;
end
end
error=max(max(abs(Told-T)));
Told=T; %updating value of temperature for each iteration
T_prev_dt=T; %updating value of temperature matrix for previous time step
iteration_counter=iteration_counter+1;
end %end of convergence loop
end %end of time loop
toc;
time=toc;
%Plots for Jacobi Iterative Solver
figure(1)
[x_y,temp]=contourf(xx,yy,T);
colorbar
colormap('jet')
clabel (x_y,temp)
title_text=sprintf('For IMPLICIT SCHEME; JACOBI SOLVER: Number of Iterations=%d, Total Time=%f seconds',iteration_counter, time);
title(title_text)
xlabel('X')
ylabel('Y')
end %end of Jacobi iterative solver
%Gauss-Seidel solver:-
tic;
if iterative_solver==3
%time-step
dt=0.02;
k1=(alpha*dt)/(s^2);
for nt=1:1400 %time loop
error=9e9;
%convergence loop
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T_parenthesis=T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1);
T(i,j)=T_prev_dt(i,j)/(1+(4*k1))+((k1)/(1+(4*k1)))*T_parenthesis;
end
end
error=max(max(abs(Told-T)));
Told=T; %updating value of temperature for each iteration
T_prev_dt=T; %updating value of temperature matrix for previous time step
iteration_counter=iteration_counter+1;
end %end of convergence loop
end %end of time loop
toc;
time=toc;
%Plots for Gauss-Seidel Iterative Solver
figure(1)
[x_y,temp]=contourf(xx,yy,T);
colorbar
colormap('jet')
clabel (x_y,temp)
title_text=sprintf('For IMPLICIT SCHEME; GAUSS-SEIDEL SOLVER: Number of Iterations=%d, Total Time=%f seconds',iteration_counter, time);
title(title_text)
xlabel('X')
ylabel('Y')
end %end of Gauss-Seidel iterative solver
%Successive Over-Relaxation (SOR) solver:-
tic;
if iterative_solver==4
%time-step
dt=0.02;
k1=(alpha*dt)/(s^2);
for nt=1:1400 %time loop
error=9e9;
%convergence loop
while (error>tol)
for i=2:nx-1
for j=2:ny-1
T_parenthesis=T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1);
T_Gauss_Seidel=T_prev_dt(i,j)/(1+(4*k1))+((k1)/(1+(4*k1)))*T_parenthesis;
T(i,j)=(1-w)*Told(i,j)+w*T_Gauss_Seidel;
end
end
error=max(max(abs(Told-T)));
Told=T; %updating value of temperature for each iteration
T_prev_dt=T;
iteration_counter=iteration_counter+1;
end %end of convergence loop
end %end of time loop
toc;
time=toc;
%Plots for SOR Iterative Solver
figure(1)
[x_y,temp]=contourf(xx,yy,T);
colorbar
colormap('jet')
clabel (x_y,temp)
title_text=sprintf('For IMPLICIT SCHEME; SOR SOLVER: Number of Iterations=%d, Total Time=%f seconds',iteration_counter, time);
title(title_text)
xlabel('X')
ylabel('Y')
end %end of SOR iterative solver
RESULTS & DISUCSSIONS:-
From the above plots, we can see that the ‘explicit scheme’ takes ‘1607-iterations’ to converge. Whereas ‘implicit scheme’; Jacobi, Gauss-Seidel, and SOR solver take ‘1620-iterations’, ‘1534-iterations’, and ‘1494-iterations’ respectively.
The total computational time by each method is, Explicit Scheme=0.021542 seconds; Implicit Scheme: Jacobi solver=0.025442 seconds, Gauss-Seidel solver=0.024747 seconds, and SOR solver= 0.026009 seconds.
We can deduce it from the above-discussed data that the Jacobi solver is the slowest one as it takes much more time and the number of iteration; to get converged. Whereas, the SOR solver is the fastest amongst all, as it is an advanced method that determines how aggressively the current iteration gets updated from the values in the current iteration, which makes it faster than the other.
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
CONCLUSION:-
- We have solved unsteady (transient) 2-D heat conduction equations both explicitly and implicitly. When we compare both the schemes; we can observe that the explicit scheme took a paltry amount of time to get to the solution as compared to the implicit scheme. This happens because in an explicit scheme; each discretized equation at each time-step contains only one unknown and that can be solved in a straightforward manner; requiring fewer iterations and time. But, in an implicit scheme, we get unknowns through a simultaneous solution of discretized equation applied at all grid points (nodes) arrayed at a particular time-step. Because of this, we need to solve many simultaneous algebraic equations. The implicit scheme involves the manipulation of large matrices, and because it involves complex calculations; implicit scheme requires much more time than the explicit scheme.
- On comparing the results of the steady-state solution and unsteady-state solution, we can notice that for the steady-state, the time and number of iterations required to get the converged solution are much less than compared to the unsteady-state. This is because the governing equation for the steady-state 2-D heat conduction equation contains an only spatial variation of temperature with no thermal diffusivity; also as being steady-state it does not include time derivative term. These all make steady-state equations a simple partial differential equation to solve. Since we’re solving steady-state implicitly by point iterative solvers; we always get a stable solution. Whereas, unsteady-state 2-D heat conduction equation contains temperature varying w.r.t both time and spatial derivative terms. Hence, an unsteady-state solution takes a long time to get a converged solution as compared to a steady-state.
- We have used different point-iterative solvers, such as Jacobi, Gauss-Seidel, and Successive Over-Relaxation (SOR); to solve both the steady and unsteady state heat 2-D heat conduction equation. After comparing each of the iterative-solvers in terms of performance, we find that Jacobi iterative solver took a larger number of iterations and time to get converged as compare to the Gauss-Seidel and SOR iterative solvers. Whereas, the Gauss-Seidel iterative solver takes a lesser number of iterations and time as compared to the Jacobi iterative solver but larger than SOR iterative solver. The reason for all of this is; in the Jacobi iterative solver, we use old guess values for each iteration, and for the Gauss-Seidel iterative solver we use updated values (values that have been already calculated in the previous time-steps). The SOR iterative solver is the fastest amongst all the other iterative-solvers that we have discussed in this project because the SOR iterative solver is devised by applying extrapolation to Gauss-Seidel iterative solver. This extrapolation takes the form of a weighted average between previous iteration values and current iteration values for each component. But the fundamental idea behind using SOR with a full efficacy is to choose a value of relaxation factor; that will eventually speed up the rate of convergence; that iterates the solution at a faster rate. So, by choosing the value of relaxation factor greater than 1 (but close to 1), can speed up the rate of convergence.