Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Problem Statement: To solve the 2D Heat Conduction Eqn. for steady state and unsteady state by explicit and impicit method with given boundary conditions. Using implicit method, solve using different iterative solvers i.e. Jacobi, Gauss-Seidal, SOR, etc. for both steady and transient state heat conduction eqn. The accuracy of the solution obtained in both methods(implicit & explicit) should be in accordance with the tolerance value mentioned.
Boundary Conditions: 1. Top Boundary=600 K
2. Bottom Boundary=900K
3. Left Boundary=400 K
4. Right Boundary=800 K
5. Tolerance Value=>error=1e-4
Introduction: Explicit Method:
When a direct computation is made for an unknown quantity in terms of known quantity, this computation is called explicit. This method requires stability criteria which when satisfied leads to a stable solution. Because of the stability criteria, small time steps have to be taken to obtain stable solutions. With decrease in the time step, the no. of computations increases and hence the computational time required for solution- converging is less and solution accuracy is more.
Implicit Method: When a computation is made for a dependent variable defined by a set of equations and a matrix is involved to find the unknown variable, the computation is called implicit method. This method is unconditionally unstable as a result and no stability criteria is needed to obtain a stable solution and the computational time for solution-converging is high and therefore accuracy of solution is less compared to explicit method because of the usage of simultaneous equations and guess values as the initial values to figure out the approximated values. If the guess values are closer to the desired values, the computational time is less and since there's no necessity for stability criteria, larger time steps can be used.
In order to solve for dependent variables, 3 iterative solvers are used:
1. Jacobi iterative solver
2. Gauss-Seidal iterative solver
3. Sucessive Over Relaxation iterative solver.
2D Heat Conduction Eqn.: `(delT)/(delt)=alpha.((del^2T)/(delx^2)+(del^2T)/(dely^2))`
The time derivative term `(delT)/(delx)`is hyperbolic in nature and can be discretized using Forward Differencing Scheme(FDS) and the spatial derivatives are elliptical in nature and can be discretized using Central Differencing Scheme(CDS).
Discretization of time derivative using FDS is done because the time derivative is varying in time direction which can be solved by marching method.
`((delT)/(delt))^ni,j=((Ti,j)^(n+1)-(Ti,j)^n)/(trianglet)+O(trianglet)`
Discretization of spatial derivatives using CDS since the spatial derivatives are elliptical in nature and used information from neighbouring points to obtain the solution. Below given is the discretization for 'x' and 'y' spatial derivatives:
`((del^2T)/(delx^2))^ni,j=((Ti+1,j)^n-2(Ti,j)^n+(Ti-1,j)^n)/(trianglex^2)+O(trianglex^2)`
`((del^2T)/(dely^2))^ni,j=((Ti,j+1)^n-2(Ti,j)^n+(Ti,j-1)^n)/(triangley^2)+O(triangley^2)`
From the PDE Heat Conduction Eqn.:
`((Ti+1,j)^n-2(Ti,j)^n+(Ti-1,j)^n)/(trianglex^2)+((Ti,j+1)^n-2(Ti,j)^n+(Ti,j-1)^n)/(triangley^2)=1/alpha(((Ti,j)^(n+1)-(Ti,j)^n)/(trianglet))`
where 'n' is current time step and 'n+1' is next time step; 'i' and 'i+1' are the neighbouring points in 'x' direction and 'j' and 'j+1' are the neighbouring points in 'y' direction. `trianglex`is grid spacing in 'x' direction and `triangley`is grid spacing in 'y' direction.
Iterative Solvers:
1. Jacobi Iterative Solver: In this method, the temperature at the 'n+1'th iteration is found using old values of termperature defined at 'n'th iterations. Using the approximation obtained iterations procedure, it continues until the desired accuracy of the solution is obtained.
2. Gauss-Seidal Iterative Solver: In this type, temperature at 'n+1' iterations uses values of 'n+1' iterations that have already been computed. By this way, old approximated values are overwritten with new approximate values for the iterations performed. Unlike Jacobi, this solver requires one storage vector which can be overwritten with each computation. This helps the solution converge faster that Jacobi.
3. Sucessive Over Relaxation Iterative Solver: In this type of iterative solver, relaxation factor is used to increase the convergence rate. The formula is given below:
`T^(n+1)=T^n.(1-omega)+omega.T(GS)`or`Tjacobi`
where `omega`is the relaxation factor and ranges from 0 to 2. When it is between 0 and 1, it is under relaxed state. This condition is used when the temperature gradient is high and helps in faster convergence because the time step doesn't overshoot. When it is between 1 and 2, it is in over-relaxed state. Here, temp. gradient is low and larger time steps help in faster convergence i.e. obtain steady state.
`T^n`is the value of temp. for previous iterations. `T(GS) or Tjacobi`is the value of temp. for 'n+1'th iteration using jacobi or gs iterative solver.
Steady State 2D Conduction Eqn:
Steady state is where the fluid properties do not change wrt. to time. From the 2D Heat Conduction Eqn. , since the system is steady-state, time derivative is zero and diffusivity is constant.
`(delT)/(delt)=0`
Hence, the eqn. for steady-state is: `(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
Discretization for steady-state: `(T_L-2.T_P+T_R)/(trianglex^2)+(T _T-2.T_P+T_B)/(triangley^2)=0`
`T_P=T_i,j`
`T_L=T_(i-1,j)`
`T_R=T_(i+1,j)`
`T _T=T_(i,j+1)`
`T_B=T_(i,j-1)`Solving for T(i,j), we get:
`T_P=K_1(T_L+T_R)+K_2(T_T+T_B)`
`K_1=(trianglex^2)/(2(trianglex^2+triangley^2)`;
`K_2=(trianglex^2)/(2(trianglex^2+triangley^2)`
Solving the 2D Steady Heat Conduction Eqn. implicitly:
Formula for Jacobi Iterative Solver:
`T_P^(n+1)=K_1(T_L+T_R)^n+K_2(T_T+T_B)^n`
Formula for Gauss-Seidal Iterative Solver:
`T_P^(n+1)=K_1(T_L+T_R)^(n+1)+K_2(T_T+T_B)^(n+1)`
Formula for Successive Over Relaxation Iterative Solver:
`T_P^(n+1)=T_P^n.(1-omega)+omega.(K_1(T_L+T_R)^(n+1)+K_2(T_T+T_B)^(n+1))`
Transient State: Unsteady state is where the fluid properties change wrt. to time. The formula for 2D Transient Heat Conduction Eqn. is given below:
`(delT)/(delt)=alpha.((del^2T)/(delx^2)+(del^2T)/(dely^2))`
Discretization for unsteady-state:
`(T_P^(n+1)-T_P^n)/(trianglet)-alpha.((T_L-2.T_P+T_R)/(trianglex^2)+(T_T-2.T_P+T_B)/(triangley^2))=0`
`T_P=T_i,j`
`T_L=T_(i-1,j)`
`T_R=T_(i+1,j)`
`T_T=T_(i,j+1)`
`T_B=T_(i,j-1)`
Solving the 2D Transient Heat Conduction Eqn. explicitly:
`T_P^(n+1)=(1-2.K_1-2.K_2).T_P^n+K_1(T_L+T_R)^n+K_2(T_T+T_B)^n`
`K_1=alpha.(trianglet)/(trianglex^2)`
`K_2=alpha.(trianglet)/(triangley^2)`
Stability Criteria: `K_1+K_2<=0.5`is the criteria for 2D heat conduction eqn. when solved explicitly.
Solving the 2D Transient Heat Conduction Eqn. implicitly:
`T_P^(n+1)=1/K.[T_P^n+K_1(T_L+T_R)^(n+1)+K_2(T_T+T_B)^(n+1)]`
`K=1+2.K_1+2.K_2`
`K_1=alpha.(trianglet)/(trianglex^2)`
`K_2=alpha.(trianglet)/(triangley^2)` here, 'n+1' is the next time step and 'n' the preceding time step.
Formula for Jacobi Iterative Solver:
`T_P^(n+1)=[T_P^n+K_1(T_L^n+T_R^n)+K_2(T_T^n+T_B^n)]/(1+2.K_1+2.K_2)`
Formula for Gauss-Seidal Iterative Solver:
`T_P^(n+1)=[T_P^n+K_1(T_L^(n+1)+T_R^n)+K_2(T_T^n+T_B^(n+1))]/(1+2.K_1+2.K_2)`
Formula for Successive Over Relaxation Iterative Solver:
`T_P^(n+1)=T_P^n.(1-omega)+omega.[(T_P^n+K_1(T_L^(n+1)+T_R^n)+K_2(T_T^n+T_B^(n+1)))/(1+2.K_1+2.K_2)]`
Boundary Conditions:
Boundary Conditions are defined at the boundaries of independent variables of the function. In this problem, mesh size of 10x10 is considered for which the boundary conditions are defined at the extreme rows and columns of the mesh.
| 500 | 600 | 600 | 600 | 600 | 600 | 600 | 600 | 600 | 700 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 400 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 300 | 800 |
| 650 | 900 | 900 | 900 | 900 | 900 | 900 | 900 | 900 | 850 |
For the mesh size, 10x10 shown above, the corner boxes are average of the two boundaries intersecting at that point. The values other than the boundaries are taken as guess values which in this case is 300K.
Main Program:
% code to study the behaviour of steady state and transient state analysis
% using explicit and implicit method
clear all
close all
clc
% defining spatial domain
nx=10; ny=10;
% specifying length of domain
L=1;
% defining space domain
x=linspace(0,L,nx);
y=linspace(0,L,ny);
% grid size
dx=x(2)-x(1);
dy=y(2)-y(1);
% defining temperature
t=300*ones(nx,ny);
% [t]=300*meshgrid(x,y);
t(1:end,1)=400; % left boundary
t(end,1:end)=900; % bottom boundary
t(1,1:end)=600; % top boundary
t(1:end,end)=800; % right boundary
% edge temperature
t(1,1)=(t(1,2)+t(2,1))/2;
t(end,1)=(t(end,2)+t((end-1),1))/2;
t(1,end)=(t(1,(end-1))+t(2,end))/2;
t(end,end)=(t((end-1),(end))+t(end,(end-1)))/2;
% tolerance
tol=1e-4;
% making a copy of initial temperature matrix
told=t;
t_pre=t;
% coefficient obtained from derivation for steady state equation
k=2*(dx^2+dy^2)/(dx^2*dy^2);
% diffusivity coefficient in the 2D heat conduction eqn
alpha=1.4;
% calculation of time step
dt=(0.5)/((alpha/dx^2)+(alpha/dy^2));
nt = round(1.4/dt);
% coefficient multipied for x diffuse term
k1=dt*alpha/dx^2;
% % coefficient multipied for y diffuse term
k2=dt*alpha/dy^2;
% relaxation factor for successive over relaxation for steady state
% heat conduction
a=1.5;
% relaxation factor for successive over relaxation for transient heat
% conduction
a1=1.15;
% function to calculate the temp. of steady state using jacobi iterative
% solver
[error_j,iter1]=jacobi(nx,ny,t,told,k,dx,dy,tol);
% function to calculate the temp. of steady state using gauss seidal
% iterative solver
[error_gs,iter2]=gs(nx,ny,t,told,k,dx,dy,tol);
% function to calculate the temp. of steady state using successive over
% relaxation iterative solver
[error_sor_gs,iter3]=sor_gs(nx,ny,t,told,k,dx,dy,tol,a);
% function to calculate the temp. of transient heat conduction explicitly
[t,nt]=transient_heat_conduction(nx,ny,k1,k2,told,t,nt);
% function to calculate the temp. of unsteady state using jacobi iterative
% solver
[error_jac,iter_jac]=transient_heat_conduction_jacobi(nx,ny,k1,k2,told,t,tol,t_pre);
% function to calculate the temp. of unsteady state using gauss seidal iterative
% solver
[error_gs_imp,iter_gs_imp]=transient_heat_conduction_gs(nx,ny,k1,k2,told,t,tol,t_pre);
% function to calculate the temp. of unsteady state using SOR iterative
% solver
[error_sor_gs_imp,iter_sor_gs_imp]=transient_heat_conduction_sor_gs(nx,ny,k1,k2,told,t,tol,t_pre,a1);Explanation of Main Program:
1. Specify the no. of grid points in the spatial domain 'x' and 'y' direction. Here, we're using 10 grid points in 'x' and 'y' direction.
2. Specify the length of domain in 'x' and 'y' space domain.
3. Define the spatial domain for the length and grid points specified.
4. Calculate the grid spacing in 'x' and 'y' space domain.
5. Define temperature with guess values 300K for no. of grid points specified.
6. Apply the boundary conditions as mentioned in the above problem. The accuracy of the solution is measured with computed approximated solution and compared with the tolerance value so that the converged solution is of desired accuracy.
7. A copy of the initial temperature matrix is made in told and t_pre which is used while calculating the temperature values of 'n+1' time steps explicitly and implicitly.
8. Calculate the value of K,K1 and K2 and here, the value of `alpha`is taken as 1.4.
9. Functions are called to calculate the temperature values at 'n+1' time step for both steady and unsteady state using different iterative solvers and furthermore, explicit method is used to calculate temperature for unsteady state heat conduction equation.
Explanation of Functions Seen in Main Program:
1. In the functions, dummy error is assigned to error variable so that 'while loop' condition would satisfy and iterations are calculated.
2. The error is calculated and updated for each iteration. The error calculated is coded only to take the absolute values since negative values would not satisfy the 'while loop' condition and the iterations will stop and results obtained will not converge.
3. No. of iterations taken for the solution to converge is stored in the 'iter' variable. This is done by incrementing the variable by 1 after the end of each iteration. As the iteration proceeds to calculate temperature at 'n+1' time steps, temperature at 'n'th time step is required and the values must be updated for each iteration. If the temp. values are not updated for each iteration, the solution will not converge and the 'while loop' will run infinite no. of times.
4. Contourf plot is used to check if the solution converges or not and told stores the value of temp. of the preceding time step and t stores the temp. values of the 'n+1' time step. This logic is used for all the functions to calculate the temp. of steady state heat conduction eqn. for the desired accuracy.
Function to Calculate the Temperature of Steady State using Jacobi Iterative Solver:
% function to calculate the temp. of steady state using jacobi iterative
% solver
function[error_j,iter1]=jacobi(nx,ny,t,told,k,dx,dy,tol)
% defining iteration variable to store and display the number of
% iterations computed for the values to converge.
iter1=0;
% dummy error value
error_j=9e9;
% convergence loop
while(error_j>tol)
% iteration loop
for i=2:(nx-1)
for j=2:(ny-1)
% % formula to calculate the temp. at each grid point
%
t(i,j)=(1/k)*(((told((i-1),j)+told(i+1,j))/dx^2)+((told(i,(j+1))+told(i,(j-1)))/dy^2));
end
end
% error for each iteration
% error1=abs(told-t);
error_j=max(abs(told-t),[],'all');
% updating the temp. values for the next iteration
told=t;
% incrementing iter variable by 1 to display the no. of
% iterations taken for solution to converge
iter1=iter1+1;
end
% plot representing the converged results
figure(1)
[C,h]=contourf(t);
% labelling the plot here C represents the temp. values and
% h rotates the temp. value as it moves
clabel(C,h);
axis ij;
title('Number of iterations for Jacobi method for steady state heat conduction',iter1)
end
Function to Calculate the Temperature of Steady State using Gauss Seidal Iterative Solver:
% function to calculate the temp. of steady state using gauss seidal iterative
% solver
function [error_gs,iter2]=gs(nx,ny,t,told,k,dx,dy,tol)
% defining iteration variable to store and display the number of
% iterations computed for the values to converge.
iter2=0;
% dummy error value
error_gs=9e9;
% convergence loop
while(error_gs>tol)
% iteration loop
for i=2:(nx-1)
for j=2:(ny-1)
% formula to calculate the temp. at each grid point
t(i,j)=(1/k)*(((t((i-1),j)+told(i+1,j))/dx^2)+((told(i,(j+1))+t(i,(j-1)))/dy^2));
end
end
% error for each iteration
error1=abs(told-t);
error_gs=max(error1,[],'all');
% updating the temp. values for the next iteration
told=t;
% incrementing iter variable by 1 to display the no. of
% iterations taken for the solution to converge
iter2=iter2+1;
end
% plot representing the converged results
figure(2)
[C,h]=contourf(t);
% labelling the plot here C represents the temp. values and
% h rotates the temp. value as it moves
clabel(C,h);
axis ij;
title('Number of iterations for Gauss Seidal method for steady state heat conduction',iter2)
end
Function to Calculate the Temperature of Steady State using SOR Iterative Solver:
% function to calculate the temp. of steady state using SOR iterative
% solver
function [error_sor_gs,iter3]=sor_gs(nx,ny,t,told,k,dx,dy,tol,a)
% defining iteration variable to store and display the number of
% iterations computed for the values to converge.
iter3=0;
% dummy error value
error_sor_gs=9e9;
% convergence loop
while(error_sor_gs>tol)
% iteration loop
for i=2:(nx-1)
for j=2:(ny-1)
% formula to calculate the temp. at each grid point
t(i,j)=(1/k)*(((t((i-1),j)+told(i+1,j))/dx^2)+((told(i,(j+1))+t(i,(j-1)))/dy^2));
% successive over relaxation formula to calculate
% the temp. at each grid point
t(i,j)=(told(i,j))*(1-a)+a*(t(i,j));
end
end
% error for each iteration
error1=abs(told-t);
error_sor_gs=max(error1,[],'all');
% updating the temp. values for the next iteration
told=t;
% incrementing iter variable by 1 to display the no. of
% iterations taken for the solution to converge
iter3=iter3+1;
end
% plot representing the converged results
figure(3)
[C,h]=contourf(t);
% labelling the plot here C represents the temp. values and
% h rotates the temp. value as it moves
clabel(C,h);
axis ij;
title('Number of iterations for SOR method for steady state heat conduction',iter3)
endFunction to Calculate the Temperature of Unsteady State using Explicit Method:
% function to calculate thetemp. of unsteady state using explicit method
function [t,nt]=transient_heat_conduction(nx,ny,k1,k2,told,t,nt)
% defining iteration variable to store and display the number of
% iterations computed for the values to converge.
for k=1:nt
% iteration loop
for i=2:(nx-1)
for j=2:(ny-1)
H=told(i-1,j)-2*told(i,j)+told(i+1,j); % this line represents temp. values in 'x' direction
V=told(i,j-1)-2*told(i,j)+told(i,j+1); % this line represents temp. values in 'y' direction
% formula to calculate the temp. at each grid point
t(i,j)=told(i,j)+k1*H+k2*V;
end
end
% updating the temp. values for the next iteration
told=t;
end
% plot representing the converged results
figure(4)
[C,h]=contourf(t);
% labelling the plot here C represents the temp. values and
% h rotates the temp. value as it moves
clabel(C,h);
axis ij;
title('Number of iterations using transient heat conduction explicitly',nt);
end
Function to Calculate the Temperature of Unsteady State using Jacobi Iterative Solver:
% function to calculate the temp. of unsteady state using jacobi iterative
% solver
function [error_jac,iter_jac]=transient_heat_conduction_jacobi(nx,ny,k1,k2,told,t,tol,t_pre)
% defining iteration variable to store and display the number of
% iterations computed for the values to converge.
iter_jac=0;
% dummy error value
error_jac=9e9;
error_f=9e9;
% while loop to stop iterations when temp. values converge while
% computing for different time steps
while(error_f>tol)
% convergence loop each time step
while(error_jac>tol)
% iteration loop
for i=2:(nx-1)
for j=2:(ny-1)
% coefficient obatained while deriving for transient
% heat conduction eqn.
term1=(1+(2*k1)+(2*k2))^-1;
% k1 and k2 represent the cfl value for 'x' and 'y'
% domain respectively
H=told(i-1,j)+told(i+1,j); % this line represents temp. values in 'x' direction
V=told(i,j-1)+told(i,j+1); % this line represents temp. values in 'y' direction
% formula to calculate the temp. at each grid point
t(i,j)=(t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
end
end
% error for each iteration
error1=abs(told-t);
error_jac=max(max(error1,[],'all'));
% updating the temp. values for the next iteration
told=t;
% incrementing iter variable by 1 to display the no. of
% iterations taken for the solution to converge
iter_jac=iter_jac+1;
end
% error defined to obtain desired accuracy for convergence
% rate
error_f=max(max(abs(t_pre-t)));
t_pre=t;
% resetting error value for the convergence loop to compute
% for next time step
error_jac=9e9;
end
% plot representing the converged results
figure(5)
[C,h]=contourf(t);
% labelling the plot here C represents the temp. values and
% h rotates the temp. value as it moves
clabel(C,h);
axis ij;
title('Number of iterations for implicit jacobi',iter_jac);
end
Function to Calculate the Temperature of Unsteady State using Gauss Seidal Iterative Solver:
% function to calculate the temp. of unsteady state using gauss seidal iterative
% solver
function [error_gs_imp,iter_gs_imp]=transient_heat_conduction_gs(nx,ny,k1,k2,told,t,tol,t_pre)
% defining iteration variable to store and display the number of
% iterations computed for the values to converge.
iter_gs_imp=0;
% dummy error value
error_gs_imp=9e9;
error_f=9e9;
% while loop to stop iterations when temp. values converge while
% computing for different time steps
while(error_f>tol)
% convergence loop each time step
while(error_gs_imp>tol)
% iteration loop
for i=2:(nx-1)
for j=2:(ny-1)
% coefficient obatained while deriving for transient
% heat conduction eqn.
term1=(1+(2*k1)+(2*k2))^-1;
% k1 and k2 represent the cfl value for 'x' and 'y'
% domain respectively
H=t(i-1,j)+told(i+1,j); % this line represents temp. values in 'x' direction
V=t(i,j-1)+told(i,j+1); % this line represents temp. values in 'y' direction
% formula to calculate the temp. at each grid point
t(i,j)=(t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
end
end
% error for each iteration
error1=abs(told-t);
error_gs_imp=max(error1,[],'all');
% updating the temp. values for the next iteration
told=t;
% incrementing iter variable by 1 to display the no. of
% iterations taken for the solution to converge
iter_gs_imp=iter_gs_imp+1;
end
% error defined to obtain desired accuracy for convergence
% rate
error_f=max(max(abs(t_pre-t)));
t_pre=t;
% resetting error value for the convergence loop to compute
% for next time step
error_gs_imp=9e9;
end
% plot representing the converged results
figure(6)
[C,h]=contourf(t);
% labelling the plot here C represents the temp. values and
% h rotates the temp. value as it moves
clabel(C,h);
axis ij;
title('Number of iterations for implicit gauss seidal',iter_gs_imp);
end
Function to Calculate the Temperature of Unsteady State using SOR Iterative Solver:
% function to calculate the temp. of unsteady state using SOR iterative
% solver
function [error_sor_gs_imp,iter_sor_gs_imp]=transient_heat_conduction_sor_gs(nx,ny,k1,k2,told,t,tol,t_pre,a1)
% defining iteration variable to store and display the number of
% iterations computed for the values to converge.
iter_sor_gs_imp=0;
% dummy error value
error_sor_gs_imp=9e9;
error_f=9e9;
% while loop to stop iterations when temp. values converge while
% computing for different time steps
while(error_f>tol)
% convergence loop each time step
while(error_sor_gs_imp>tol)
% iteration loop
for i=2:(nx-1)
for j=2:(ny-1)
% coefficient obatained while deriving for transient
% heat conduction eqn.
term1=(1+(2*k1)+(2*k2))^-1;
% k1 and k2 represent the cfl value for 'x' and 'y'
% domain respectively
H=t(i-1,j)+told(i+1,j); % this line represents temp. values in 'x' direction
V=t(i,j-1)+told(i,j+1); % this line represents temp. values in 'y' direction
% formula to calculate the temp. at each grid point
% t(i,j)=(t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V);
t(i,j)=(1-a1)*told(i,j)+a1*((t_pre(i,j)*term1)+(term1*k1*H)+(term1*k2*V));
end
end
% error for each iteration
error1=abs(told-t);
error_sor_gs_imp=max(error1,[],'all');
% updating the temp. values for the next iteration
told=t;
% incrementing iter variable by 1 to display the no. of
% iterations taken for the solution to converge
iter_sor_gs_imp=iter_sor_gs_imp+1;
end
% error defined to obtain desired accuracy for convergence
% rate
error_f=max(max(abs(t_pre-t)));
t_pre=t;
% resetting error value for the convergence loop to compute
% for next time step
error_sor_gs_imp=9e9;
end
% plot representing the converged results
figure(7)
[C,h]=contourf(t);
% labelling the plot here C represents the temp. values and
% h rotates the temp. value as it moves
clabel(C,h);
axis ij;
title('Number of iterations for SOR',iter_sor_gs_imp);
end
Output Graphs:
Inference from Plots:
A. Steady State Heat Conduction Eqn.:
Iterative Solver No. of Iterations
Jacobi 207
Gauss Seidal 111
SOR 28
B. Unsteady State Heat Conduction Eqn.:
Iterative Solver No. of Iterations
Jacobi 1982
Gauss Seidal 1458
SOR 1148
Conclusion: Implicit method can be used to solve problems when converged results are needed with less computation times as they are unconditionally stable. With increases in time step sizes, the convergence rate increases which leads to steady state with less no. of iterations. Explicit method results in time-accurate solutions as it is governed by stability criteria and if its not satisfied, it leads to an unstable solution and requires more no. of iterations to acheive steady state condition. But, in the code provided above, time step size remains constant which is calculated wrt. to stability criteria. Hence, for this time step size, the explicit scheme leads to faster convergence rates as time step size is small whereas for implicit scheme, the no. of iterations increases which results in a slower convergence rate.